Anchovy (

La anchoa (

Anchovy (

In Italy anchovy are fished all year round by pelagic trawlers (known in Italy as

Historically, the Italian catches of anchovy represented the largest part of the total landings, while the eastern fleet used to target mainly sardine, due to fish market preferences (

In the late 1980s the small pelagic fisheries suffered a major crisis and the Italian total catches dropped to the historic minimum of 5875 t (

Since 1975, the CNR-ISMAR of Ancona (Italy) has been conducting research on the biology and stock assessment of anchovy (together with sardine) in the northern and central Adriatic by means of stock assessment models (

European anchovies are small-bodied, short-lived pelagic species. They are serial batch spawners: the spawning takes place in the warmer months, between April and October (

Many studies have been carried out regarding the presence of a unique stock or the presence of different sub-populations living in the Adriatic Sea. This has several implications for the management, i.e. differences in the growth features between subpopulations could imply the need for ad hoc strategies in the management. In the present paper, following the suggestions and the findings in

The interest in this species, its stock and its distribution, is due not only to the importance it has for the fishery and the economy of the countries involved, but also to its ecological role: small pelagics occupy a key position in the trophic chain, responding to changes from levels below and above. With their fluctuations they affect their prey (plankton) and their predators (humans among them) (

The most detailed and recent published stock assessment on this species dates back to 2003 and comprises data from 1975 to 1996. However,

Up to 2009, Laurec-Shepherd VPA and the De Lury depletion model were mainly used to assess the Adriatic stock of anchovy (

In 2001, the Population Dynamics Unit and the Marine Acoustics unit of Ancona began to collaborate with the countries on the eastern Adriatic side (in particular Croatia and Slovenia): with the support of the Italian Ministry for Agriculture Food and Forestry Policies (MIPAAF) and the FAO-AdriaMed regional project, all the data available among the parties involved were joined in a common database developed

Catch-at-age data from all the countries were pooled together and the first stock assessment based on this new database was presented in occasion of the 2009 GFCM meeting (

Some drawbacks of the VPA model are that it does not take into account errors in abundance indices and catch data, it ignores correlation between values of q for different ages and it is not based on a formal statistical model (

In this paper, integrated catch at age (ICA) (

ICA has been successfully applied to other small pelagic short-lived species (

This paper addresses several goals. First, it meets the need for an updated stock assessment of small pelagics in the Adriatic Sea, including in the analysis the new data obtained from the collaboration of the countries involved, namely Italy, Croatia and Slovenia. Second, it presents a comparison between the simpler VPA and the more elaborate ICA. Finally, it improves former estimates of natural mortality and tests the effects on the assessment results.

The present assessment was carried out using data collected in both the western (Italy) and eastern side (Croatia and Slovenia) of GSA17. In particular, monthly and/or yearly landings of anchovy are available for all the main ports of Italy, Croatia and Slovenia starting from 1975: from that, annual total landings for the three countries were derived (

Biological data from commercial samples and abundance and biomass data from acoustic surveys have been collected, at least from the western side, from 1975-1976 to present. However, some years are missing in the time series and the data are not always homogeneously available throughout the area.

Biological data from commercial samples include (i) year-specific length frequency distribution of the landings from 1975 onwards for the western side, and from the 1998 for the eastern side; (ii) age-length keys (ALKs) by means of otolith readings available for both the western and eastern Adriatic since 1984 and 2001, respectively, with some gaps in the years; and (iii) catch-at-age given by the sum of independent samplings between the eastern and western side only since 1998: biological data from the western side have been assumed to be valid for the eastern side before 1998 (

Year | 1975 | 1976 | 1977 | 1978 | 1979 | 1980 | 1981 | 1982 | 1983 | 1984 | 1985 | 1986 |
---|---|---|---|---|---|---|---|---|---|---|---|---|

Total Adriatic Landings (t) | + | + | + | + | + | + | + | + | + | + | + | + |

Western LFD | + | + | + | + | + | + | + | + | + | + | + | + |

Western ALK | + | + | ||||||||||

Eastern LFD | ||||||||||||

Eastern ALK | ||||||||||||

Year | 1987 | 1988 | 1989 | 1990 | 1991 | 1992 | 1993 | 1994 | 1995 | 1996 | 1997 | 1998 |

Total Adriatic Landings (t) | + | + | + | + | + | + | + | + | + | + | + | + |

Western LFD | + | + | + | + | + | + | + | + | + | + | + | + |

Western ALK | + | + | + | + | ||||||||

Eastern LFD | + | |||||||||||

Eastern ALK | ||||||||||||

Year | 1999 | 2000 | 2001 | 2002 | 2003 | 2004 | 2005 | 2006 | 2007 | 2008 | 2009 | |

Total Adriatic Landings (t) | + | + | + | + | + | + | + | + | + | + | + | |

Western LFD | + | + | + | + | + | + | + | + | + | + | + | |

Western ALK | + | + | + | |||||||||

Eastern LFD | + | + | + | + | + | + | + | + | + | + | + | |

Eastern ALK | + | + | + | + | + | + | + | + |

Since otolith reading has not always been performed, for some years there are no ALKs available. An iterative algorithm from

Acoustic surveys in three different areas were carried out by the Marine Acoustics unit of the CNR-ISMAR of Ancona and by the Institute of Oceanography and Fisheries of Split (Croatia). Length-frequency distributions of the samples collected during the acoustic surveys are available for 1980 and from 1982 from the northwest Adriatic survey, and from 1987 from the mid-Adriatic survey. For the eastern Adriatic survey only data on the total stock biomass are available (

Total biomass | LFD | ||
---|---|---|---|

Western | North Adriatic | From 1976 to present | 1980 - From 1982 to present |

Mid-Adriatic | From 1987 to present | From 1987 to present | |

Eastern | From 2003 to present | From 2009 to present | |

Gaps (western Adriatic only) | North Adriatic | 1979 - 1984 - 1986 - 2002 - 2003 | |

Mid-Adriatic | 1990 - 1991 - 1996 - 2000 - 2003 - 2004 |

Biomass at length from the western acoustic surveys was turned into numbers at length and then split into age classes by means of ALKs reconstructed from the Italian commercial samples. The total biomass from the eastern acoustic survey was distributed into length classes and then split into age classes by means of length-frequency distributions and ALKs coming from the Croatian commercial samples. Hence, tuning indexes of acoustic abundance at age were obtained under the assumption of equivalent length-frequencies and length-at-age distribution between acoustic surveys and commercial catches recorded in the same area and in the same year.

The effect of natural mortality on both the models was tested. Natural mortality is closely related to the growth parameters of a given stock and changes as a function of age, being highest in the first years of life (

In the present work, two options for estimating a natural mortality vector by age were analysed; the methods explored were (i) the ProdBiom software (

$$M=[\frac{B}{t2-t1}]ln(t2-t1)+A,t>0$$ (1)

and (ii) the empirical equation proposed by

ln M = 0.66 – 1.69 ln L + 1.45 ln L_{∞} + 0.9 ln K (2)

The growth parameters used in the above equations (L_{∞}=16.1, k=0.4, t_{0}=-2.04, a=0.003, b=3.263) have been estimated through monthly sampling conducted along the western side of GSA17 and are therefore most likely representative of the majority of the population (

Since the spawning takes place mostly in spring-summer (

The oldest specimens found in the catches are 6 years old, but the bulk of animals are between age 0 and age 3, so age 4, 5 and 6 were grouped into a 4+ age-class. In the ICA model it was necessary to have an age 5 plus group due to an internal constraint of the Patterson’s software, which requires a minimum of 6 age classes. This inconsistency between the two models has been solved by giving a lower weight (equal to 0.05) to age class 5+.

Numbers at age estimated from western and eastern acoustic surveys were pooled together to obtain a unique index, aggregated in four age groups, from 2004 to 2009.

The retrospective pattern was tested for both VPA and ICA to investigate the results of sequential estimates of fishery variables (i.e. biomass and fishing mortality) as additional year of data are added. The temporal window tested goes from the third year of the acoustic survey (2007) to the last year in the data (2009). The stock recruitment relationship was also investigated, plotting the estimated SSB in the year x against the recruitment in year x+1.

Finally, the exploitation rate (E) was computed for each model as the ratio of fishing mortality (F) to total mortality (Z) (given by fishing mortality plus natural mortality (M), E=F/(F+M) or E=F/Z). The estimation was averaged between ages 1 to 3, which are the ages most represented in the catches. This exploitation rate has been related to the reference point for F proposed by

The analyses were performed using the MAAF-VPA software package developed by _{bar}(n)·k. For example, if the last true age (oldest) is set equal to 3, n equal to 2 (so F_{bar} is calculated as the mean of 2 younger ages) and k equal to 1.6, the estimated selectivity curve will be higher in the last age groups, which means that the fleet is targeting specifically older age classes. Various ratios of these values were tested. In the other case (option 2), a vector of annual F (F_{x}) at the last true age is required. This was calculated using the following proportion:

The average F for the time interval 1976-1980 (F_{(76-80)}) was derived from a previous assessment carried out following the standard procedure presented above (option 1). This specific time window (1976-1980) was chosen for two reasons: first, since it is before any sign of collapse, it can be considered as a sort of steady state that characterized the abundance of the stock well; second, in a VPA model estimates for the earlier years are more reliable than those from the recent one.

E(76-80) is an index of fishing effort (expressed in fishing days, not standardized) from 1976 to 1980: this effort estimate was obtained by averaging each year effort of the fleet from 1976 to 1980 in Porto Garibaldi (i.e. the most representative port for the small pelagic fish species back in those days: CNR-ISMAR Database) with the entire fleet fishing effort in the same time period (

Ex is the annual effort, expressed as number of fishing days per year, from 1976 to 2009, estimated by averaging the entire fleet fishing days with the effort from the Porto Garibaldi fleet: from 2000 onwards, to this estimation a 30% increasing was added in order to take into account the improvements in technology of fishing vessels, otherwise not taken into account by the raw fishing days.

This laborious estimation of F aims to overcome the deficiency of both the effort indices: Porto Garibaldi effort is in fact too low in the most recent years due to a loss of importance of this harbour, while the overall effort shows some strong fluctuations in the time series (

The effect of shrinkage was tested: in the Laure-Shepherd VPA the shrinkage allows the arithmetic mean of the fishing mortality values estimated for the five years preceding the final year to be incorporated with an appropriate weight. This procedure reduces the predicted variance at the cost of a bias towards the mean (

In this paper the results of the model with fixed terminal F vector and shrinkage will be presented, since it was the one that performed best in terms of predicted q, standard error, retrospective pattern and overall results. The performances of this will then be compared with the ICA model.

Integrated catch at age (ICA) is a methodology that bases its algorithm on a separability assumption (_{y}) and an age effect (S_{a}), such that F at age is a multiplicative combination of these (F_{y,a}=F_{y}·S_{a)}. The tuning indices in the separable model can be weighted manually or using an inverse-variance reweighting. Deviations from constant selectivity over the separable model period can be modelled in two ways, namely gradual and abrupt change. Moreover, catches at age are assumed to be measured with an error that has an independent lognormal distribution (

Both age 1 and age 2 were tested as a reference age for the separable constraint, i.e. the first age that is to be considered fully recruited. Age 1 values, although reasonable, cause some instabilities in the algorithm, with run time errors or weird/bad results in terms of biomass and coefficients of variation, which is why age 2 was used in the final model. Various assumptions for last age selectivity and for the separability period were investigated: The first was explored within a range between 0.6 and 1.6 years, and the second was explored within a range between 5 years and 15 years. The last age selectivity in the final model was set to 1.6, and a separability assumption of 10 years was chosen as the most consistent with the data available (it is reasonable to assume that in the last 10 years the fishing effort has been more or less stable).

The values estimated from ProdBiom were very low: even when experimenting ith different growth curves, we found that the model is really sensitive to changes in the input values, ranging from really low values, unrealistic for such short-lived species, to really high values (almost 3 for age 0). On the other hand, Gislason’s equation looks more stable, being influenced in the age 0 mortality mainly by the L_{inf} value (

ProdBiom | Gislason | |
---|---|---|

Age 0 | 0.68 | 0.94 |

Age 1 | 0.43 | 0.70 |

Age 2 | 0.34 | 0.59 |

Age 3 | 0.30 | 0.53 |

Age 4+ | 0.28 | 0.50 |

The results using M-vector from ProdBiom had several problems in the diagnostics of both models (evident trends in the catchability residuals, issues with the objective function minimization, very large residuals), and for this reason are not discussed here. On the other hand, the results obtained using Gislason’s values were satisfactory in terms of low observation errors (small differences between observed and predicted values) and low process errors (low standard error of the log catchability).

The shrinkage and the fixed terminal F appeared to be good instruments for stabilizing the estimation of both the biomass and the F in the final model, which otherwise fluctuated widely. The standard error of the mean Log catchability indicates a good fit for ages 0 and 1 (standard error below 0.5), and a slightly less good fit for age 2 (

Age 0 | Age 1 | Age 2 | ||
---|---|---|---|---|

VPA Run | Pred Log q | –7.14 | –5.98 | –5.27 |

S.E. (log q) | 0.413 | 0.384 | 0.789 | |

Slope | 0.144 | 0.099 | 0.388 |

Both the mid-year biomass and the mid-year SSB (

The F vector, in particular for ages 0 and 1, reaches its maximum in the years just before the collapse, with a value for age 1 equal to 0.97 in the year 1986. F_{(2-4)} increases again to higher values between 2000 and 2004, when a slight decrease in biomass is observed. From the 2005 onward, the F remains stable around values of 0.03, 0.2 and 0.4 for age 0, age 1 and age 2, respectively (

The selectivity pattern (_{a,y} and F_{(0-4),y}: the time period chosen for the calculus is the same as that used in the ICA for the separability assumption, i.e. from 1998 to 2009.

Age 0 | Age 1 | Age 2 | Age 3 | Age 4 | |
---|---|---|---|---|---|

Selectivity | 0.17 | 0.9 | 1.55 | 1.19 | 1.19 |

The SSB/recruitment relationship as obtained by Laurec-Shepherd tuning showed two clouds of points, one at the lower part of the curve including the years of the collapse from 1982 to 2002, and the other one on the upper part of the plot including the oldest years and the more recent ones. However, no clear relationship (i.e. Beverton-Holt nor Ricker) is recognizable (

The retrospective analysis shows some retrospective pattern in the biomass estimate, with a slightly systematic underestimation (graph on top), but it shows almost no pattern for F (graph on the bottom) (

Mean exploitation rate (F/Z ratio) was calculated for ages 1 to 3 and ranges between 0.24 and 0.58, reaching its maximum in 1986. The last year exploitation rate is 0.37, which is considered a reasonable value for small pelagics, below that suggested by

The statistic for the acoustic index by age show a good fit. In particular, the residual variance by age remains low and the assumption of log-normally distributed error is justified since skewness and Kurtosis absolute values are lower than 2 and the chi-square test value is highly significant for all ages (

Age | 0 | 1 | 2 | 3 | 4 |
---|---|---|---|---|---|

Variance | 0.0264 | 0.0223 | 0.0992 | 0.052 | 0.4486 |

Skewness | –0.8247 | –0.6043 | –0.105 | –0.0295 | –0.3809 |

Kurtosis | –0.24 | –0.3706 | –0.6577 | –0.7019 | –0.3403 |

Partial chi-square | 0.0064 | 0.0053 | 0.0245 | 0.0141 | 0.1459 |

Significance in fit | 0 | 0 | 0.0001 | 0 | 0.0025 |

Weight in the analysis | 0.2 | 0.2 | 0.2 | 0.2 | 0.2 |

What really made the difference between the various runs was the selectivity assumption: although the value of 1.6 for the last age may be questionable, overall it gave the best model fit, with no trend in the residuals and satisfactory diagnostic values. Parameter estimates concerning the catches and F also indicated a good model fit, as the coefficient of variation did not exceed 20% in most cases. The coefficient of variation for the catchability was around 30% for all ages.

The trend in biomass did not differ so much between the separability periods tested: the collapse in 1987 is strongly pronounced and after that the biomass tends to increase with a peak in 2004-2005 (

The trend of F by age is represented in

The analysis of the cohorts in the catch-at-age matrix justifies the trend in F and biomass observed: in the years between 1985 and 1992, with few exceptions, a relatively high number of small specimens and a concurrent low proportion of old age specimens were caught, a pattern that most likely contributed to the collapse observed in 1987. A second, strong reduction of old ages in the catch is observed from 1996 to 2001 and corresponds to a second decline in biomass with high F values for ages 2 to 4 (

The selection pattern of the associated fishery for each age-class with the coefficient of variation and the confidence intervals calculated for the separability period is shown in

Selection (S) by Age | ||||||
---|---|---|---|---|---|---|

Age | Estimate | CV | Lower 95% CL | Upper 95% CL | –s.e. | +s.e. |

0 | 0.0645 | 20 | 0.043 | 0.0967 | 0.0525 | 0.0793 |

1 | 0.4021 | 15 | 0.2956 | 0.5471 | 0.3437 | 0.4705 |

2 | 1 | Fixed: reference age | ||||

3 | 1.5434 | 12 | 1.2104 | 1.9679 | 1.3634 | 1.7471 |

4 | 1.6 | Fixed: last true age |

The shape of the relationship as obtained by the ICA model has a more discernible trend than the previous VPA model and the Beverton-Holt stock recruitment function was estimated (

The retrospective results improve for the biomass and are still reasonable for F values (

Mean exploitation rate (F/Z ratio) was also calculated for ages 1 to 3: it ranges between 0.27 in 1991 and 0.69 in 2000, the high second value being due to the high F estimated by ICA for age 3 in this year. The last year exploitation rate is 0.54, higher than the VPA estimate.

The present work aims to update and improve the stock assessment of

The final run for each model included Gislason’s estimate of natural mortality. Furthermore, an estimate of the F vector for the last age was used, together with the shrinkage to the mean F; in ICA, a separable period of 10 years and a selectivity at the last age equal to 1.6 were preferred.

There is no perfect method for estimating the natural mortality vector. We tested two models to calculate it. ProdBiom was originally developed in order to estimate the mortality at age of demersal species. In this method the growth parameters are used in order to calculate the overall biomass production and the overall losses in biomass due to a mean natural mortality, in an equilibrium situation without the influence of the fishery. The fraction between these two terms in an equilibrium situation should be equal to 1; the program fixes this ratio to 1 and achieves it by changing iteratively the two parameters, namely an asymptotic M and the parameter describing the concavity of the curve, used to calculate the final M. Hence, the life history is not taken strictly in consideration, and this method was built on hake (

Overall, two models give similar results in terms of trends, with the absolute biomass estimated by the VPA higher than the one estimated by the ICA. The collapse of 1987 is strong in both models, with a decrease in the stock size from the peak in 1978 ranging from 90% in VPA to 80% in ICA. This collapse is a focal point of the assessment: it was not only seen in landings data (and in the major crisis that hit fishermen in that period) but also in abundance estimates derived from acoustic surveys (

After that the trend in biomass recovers until 1997. In the following years the estimated biomass decreases again to a value slightly higher than the minimum of 1987. This reduction is consistent between the two models, even if it was not perceived by the fishery, as the landings were not affected (

The trend in F is similar for VPA and ICA, apart from the average F between age 2 and 4: in ICA the F shows higher values from 1999 to 2003 than in the years of the collapse. This contrast is due to the differences in the F calculation between the two methods: in the VPA run the terminal F is fixed, in ICA instead F is defined in terms of age and year components.

The stock recruitment relationship is better defined in ICA: in fact, the points corresponding to the earlier years of the time series are less dispersed and allows a Beverton-Holt curve to be drawn even if no plateau is achieved. On the other hand, in VPA two groups of points are discernible, corresponding to the years of higher biomass and years of lower biomass.

The retrospective results improved slightly in the ICA estimations, even if the pattern did not disappear completely: this may be due to the short time series of the survey index;

While VPA is strongly influenced by the choice of the F pattern, the strongest assumption in ICA regards the reference age together with the selectivity value: a careful examination of the diagnostic, though, will help in the final decision.

Overall, ICA seems to perform well with the data: the results are consistent with VPA, but ICA offers more options and several advantages. One of the limits of the Laurec-Shepherd VPA is the fact that only disaggregated tuning indices by age can be used. ICA, on the other hand, is more flexible, allowing, for example, the use of aggregated indices by age such as CPUE data, SSB index and larval indexes together with different tuning models. Though the Patterson software has some constraints due to intrinsic deficiency of the program, i.e. constraints on the length of the maximum separability period and on the value of the reference age (

One of the crucial points of the whole modelling regards the attempt to explain the collapse of 1987 and the reduction in biomass observed about ten years later. The consequences that a recruitment failure has on the following stock biomass have been documented for several small pelagic stocks all around the world (

Other anomalies of several fish stocks were recorded in the late 1980s around the world and some authors correlate them with climatic oscillations on decadal and multidecadal timescales, such as the North Atlantic Oscillation Index, which conspicuously changed from mainly negative to mainly positive values in that window of time (

The estimated reduction in biomass starting in 1999 must clearly be attributed to a diminishing of the oldest ages in the catch-at-age matrix. However, it is not clear what caused the disappearance of the oldest ages in those years. The weighted mean length during the whole time period has a negative trend (slope=b=–0.011) whose slope, however, is not significant (n=31, p-value of the slope=0.2798), implying that the mean length in the catches over the entire time series is somewhat stable. A strong decline is observed from 1993 until 2005 (slope=b=–0.14, n=10, p-value of the slope <0.05), corresponding to the second decline in biomass detected from the models: this decline in the mean length surely influenced the age structure and consequently played a determining role in the biomass estimation. After 2005 the length in the catches increases again to average values (

One hypothesis to explain the change in age composition of the catches could be the shift in distribution of the bigger animals: older specimens may have moved further from the main fishing ground. This is just a speculation, but is based on the fact that in 2001, 2003 and 2000 the highest peak in sea surface temperature since 1976 was recorded, respectively, in winter, spring-summer and autumn (AVHRR Pathfinder V5 dataset, NOAA NODC). Furthermore, the Po River flow rate had two extraordinary outflow rates in autumn 2000 and 2002.

In conclusion, ICA was found to be a good tool for evluating small pelagic stocks in the Adriatic Sea: this approach allows sufficient flexibility to combine all the data available without adding too much complexity in comparison with a VPA approach. Due to the importance of anchovy in the Adriatic Sea, further development will be attempted in order to include the influence of environment into the stock assessment and to understand better the dynamics behind the observed fluctuations.

We thank the Direzione Generale della Pesca Marittima del Ministero Italiano per le Politiche Agricole, Alimentari e Forestali (MIPAAF), which partially funded this research. We are also very grateful to the AdriaMed project, which helped in the cooperation part, Mr. N. Cingolani from CNR-ISMAR Ancona for the priceless contribution to the management of the data, Mr. A. Belardinelli, Dr. M. Panfili, Dr. S. Angelini from CNR-ISMAR Ancona who were involved in fisheries data collection, and the acoustic team, in particular Dr. F. Campanella. Finally, we would like to thank the scientific editor and one anonymous reviewer for all the suggestions and the help to improve the manuscript.