A joint stock assessment for the anchovy stock of the northern and central Adriatic Sea : comparison of two catch-at-age models

Anchovy (Engraulis encrasicolus, L.) is one of the most important commercial species of the northern and central Adriatic Sea, as well as one of the most productive fisheries in the whole Mediterranean. In the Adriatic Sea the stock of anchovy is shared between Italy, Croatia and Slovenia. A joint stock assessment was carried out using catch data from all the fleets for the time interval 1975-2009. Analyses were performed using estimates of natural mortality at age obtained by means of two different methods and two population dynamics methods based on the analysis of catch-at-age data: LaurecSheperd virtual population analysis (VPA) and integrated catch-at-age (ICA), both tuned to acoustic estimates of abundance. Gislason’s estimates for natural mortality appeared to be more realistic and were thus preferred for short-lived species. The general trend of biomass and fishing mortality is similar for the two models, highlighting the major collapse of the stock in 1987. Nevertheless, ICA has enough flexibility to combine all the data available without adding too much complexity in comparison with a VPA approach and seems to perform better in terms of the spawning stock biomass/recruitment relationship and diagnostics (i.e. the retrospective pattern). For the stock status, the exploitation rate from ICA is higher than the suggested threshold of 0.4 proposed by Patterson for small pelagic species.


Fishery
Anchovy (Engraulis encrasicolus L.) is one of the most important commercial species in the Adriatic Sea (Cingolani et al. 1996, IREPA 2009), accounting-together with sardine-for a high percentage of the overall fish production of the surrounding countries, mainly Italy and Croatia (FAO 2012).Anchovy is distributed in the whole Adriatic Sea, but the area investigated in the present work covers the northern-central part of the basin, also known as Geographical Sub-Area 17, where the majority of the catches are taken (Fig. 1).The geographical sub-areas (GSAs) are geographically defined zones in the Mediterranean, Black Sea and connecting waters that were established to facilitate data collection and to monitor and assess fisheries resources in a geo-referenced manner (GFCM 2013a).GSA 17 is enclosed between Italy, Slovenia and Croatia and its southern limit is Vieste to the west and the Croatian-Montenegro border to the east.
In Italy anchovy are fished all year round by pelagic trawlers (known in Italy as volante) and in the warm season (late spring, summer and early autumn) by purse seiners (known in Italy as lampara).The pelagic trawlers usually observe a closure period in August.In Croatia anchovy is fished all year round during the daytime with pelagic trawl-boats in the shallow areas and during the night with purse seine with artificial lights (Tičina 2003); although a small number of pelagic trawlers were active in the past, in the most recent years their number has fallen almost to zero.Before 2006 the Croatian fishery was open all year round; after then, a closure period of one month starting from December 15 was introduced.In Slovenia the last two commercial midwater pelagic trawlers stopped working in 2012, and in 2013 only a small fleet of purse seiners was actively operating (Marčeta 2001, Tičina 2003, Santojanni et al. 2011).
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 (Grbec et al. 2002, Cingolani and Santojanni 2003, Santojanni et al. 2005).This separation is now less evident, partly due to the changed market situation with the breakdown of the former Yugoslavia (Marčeta 2001).In particular, in the past 30 years two main events have driven the anchovy markets in the area and had a strong influence on the landings.The first took place in Italy from 1975 to 1985, when a regulation from the European Economic Community (EEC) charged the Azienda di stato per gli Interventi nel Mercato Agricolo (Gazzetta Ufficiale delle Comunità Europee, N. L20/1; A.I.M.A.) to buy the unsold catches from the Italian fishermen for a price that was sometimes higher than the market price.In 1982, for example, the A.I.M.A. withdrew a total of 45000 tons of anchovy and sardine for a total price of about €4.5 millions.This is the main reason that brought the Italian landings of anchovy in the Adriatic to reach the threshold of 54000 tons in 1980 (Santojanni et al. 2005 and references therein).The second event occurred in the 1990s, due to the war and subsequent partition of Yugoslavia that brought the reduction of the stock to a 70% decrease in catches in the eastern Adriatic area.The recorded landings reached their historical minimum (298 t in 1996), while the Italian share of pelagic landings in the total Adriatic catches increased considerably (Mannini andMassa 2000, Marčeta 2001).
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 (Istituto di Scienze Marine (CNR-ISMAR) of Ancona database).In the last ten years, with some fluctuations, the catches have increased again to higher levels on both the western and eastern side of the Adriatic coasts.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 (Santojanni et al. 2003) and acoustic surveys (Azzali et al. 2002, Leonori et al. 2011).The Croatian coast, up to the midline, has been covered since 2003 by the Croatian national pelagic monitoring programme PELMON (Tičina et al. 2006).In order to ensure comparability of acoustic estimates, an intercalibration between the two research vessels conducting these surveys was performed (Leonori et al. 2012).

Biology
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 (Regner 1996), with peaks in May-June and August-September (Sinovčić and Zorica 2006, Morello and Arneri 2009, Zorica et al. 2013).Recent surveys performed in the Adriatic Sea within the Sardone Project found larvae as early as February and as late as November, indicating an extension of the spawning season in the last few years (SARDONE 2010).The main spawning activities are located all along the Italian coast line, from the Gulf of Trieste to the Gargano peninsula (outside the lower boundary of GSA 17); the Gulf of Trieste and the river Po delta are the areas where the largest number of eggs occur (Regner 1996).The migration patterns of anchovy can be summarized as an inshore spawning migration towards the Italian coast at the beginning of the spring and an offshore trophic migration during the colder months (Tičina 2003, Morello andArneri 2009).
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 Tudela (1999) and Magoulas et al. (2006), we considered the Adriatic anchovy as a unique stock.
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) (Checkley et al. 2009).Because the Adriatic Sea is an important habitat for anchovies (Giannoulaki et al. 2013), their key role as a prey element has been documented through a process-oriented model (Ecopath with Ecosim) applied to the Adriatic ecosystem by Coll et al. (2007Coll et al. ( , 2009)).

State of the art and aim of the study
The most detailed and recent published stock assessment on this species dates back to 2003and comprises data from 1975to 1996. However, Santojanni et al. (2006a) related the biomass obtained by two more recent stock assessment models (virtual population analysis [VPA] and the De Lury method) to some environmental factors acting in the Adriatic Sea.Fur-thermore, updated stock assessments are presented and discussed during the annual meeting held by the Scientific Advisory Committee of the General Fisheries Commission for the Mediterranean (GFCM 2013b).
Up to 2009, Laurec-Shepherd VPA and the De Lury depletion model were mainly used to assess the Adriatic stock of anchovy (Cingolani et al. 1996, Santojanni et al. 2003, 2006a): both these models were implemented using landings from the entire northern and central Adriatic but biological data from the Italian commercial fleets; the tuning index was represented by CPUE data from one Italian port (Porto Garibaldi), which, until the late 1990s, was one of the most important in terms of number of boats and landings; the natural mortality assumed in the models was not varying by age and was set to 0.6 (GFCM 2008, SARDONE 2010).
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 ad hoc for the small pelagic fisheries (Cingolani and Santojanni 2003).
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 (GFCM 2009).Furthermore, acoustic data collected during acoustic surveys performed by the Croatian Institute of Oceanography and Fisheries (IOF) in Split become available and were joined to the acoustic data from the western acoustic survey.
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 (Shepherd 1999, Needle 2003), so its performance is not always satisfactory.
In this paper, integrated catch at age (ICA) (Patterson and Melvin 1996) was applied to the data: this method uses separable VPA with weighted tuning indices, takes into account errors in the catch-at-age data, and estimates selectivity at age based on assumptions for the fully recruited age and for the selectivity at the last true age (Patterson andMelvin 1996, Needle 2000).Moreover, in ICA, tuning indices may be either age-structured or based on non-age-structured measures of spawning stock biomass (SSB) (Needle 2000), which can be useful in the case of small pelagic data.
ICA has been successfully applied to other small pelagic short-lived species (Daskalov andMamedov 2007, Antonakakis et al. 2011) and has been considered a good stock assessment tool because it is relatively robust to noisy data and integrates various sources of information without adding too much complexity.
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.

Data
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 (Fig. 2).
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 (Table 1).
Since otolith reading has not always been performed, for some years there are no ALKs available.An iterative algorithm from Kimura and Chikuni (1987) (iterative age-length key [IALK]) was aimed at filling those gaps: this method reconstructs the age structure from catch at length data, using an available age-length key as a template.This calculation did not perform well with the data because the derived ALKs did not resemble well the observed ones.Therefore, in cases in which no ALKs were available, it was preferred to use the ALK from the closest year.The catch-at-age matrix built using the available ALKs shows a good level of internal consistency (i.e.well defined cohorts), the only flaw being in the years from 1985 to 1992 (Fig. 3).
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 (Table 2).
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 Table 1.-Summary of the biological data available by year and by area, collected through commercial samples (length-frequency distribution (LFD) for years 1982-1983 was reconstructed through linear interpolation between previous and following years).ALK = age-length key.
Year 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 Total Adriatic Landings (t) Year 1987Year 1988Year 1989Year 1990Year 1991Year 1992Year 1993Year 1994Year 1995Year 1996Year 1997Year 1998 Total 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.

Methodology
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 (Beverton and Holt 1956, Pauly 1980, Quinn and Deriso 1999); for example, asymptotic length has been proved to be negatively correlated with natural mortality, so small and short-lived species tend to have higher mortality rates.For the Adriatic anchovy a value of 0.6, equal for all the age classes, was previously assumed (Santojanni et al. 2003(Santojanni et al. , 2006a,b),b).
Since the spawning takes place mostly in springsummer (Regner 1985), the assessment was carried out taking into account a conventional birth date on 1 June (split-year), as previously done by Santojanni et al. (2003).Consequently, all data were shifted by six months in order to have each year compounded by the time interval ranging from 1 June to 31 May of the following year.The smallest mature specimen in the Adriatic is found to be around 6.7 cm for males and 7.1 cm for females, most of the specimens being mature at the end of the first year of life (Rampa et al. 2005, Sinovčić andZorica 2006).However, according to maturity data collected from biological samples, it is not rare to find age 1 animals with immature gonads.For this reason, the proportion of mature specimens was set equal to 0.5 for age 0, 0.75 for age 1 and equal to 1 for the older ages.
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 Patterson (1992) (E=0.4).

Virtual population analysis
The analyses were performed using the MAAF-VPA software package developed by Darby and Flatman (1994).Several model configurations were tested, varying the way in which the terminal F is estimated.The model allows two types of Laurec-Shepherd tuning (Laurec andShepherd 1983, Pope andShepherd 1985), which differ in the approach for the estimation of F in the last true age (Fs) (i.e. the last age before the plus group).In one case (option 1), the oldest age Fs is derived as a proportion of the average fishing mortality on n younger ages in the same year: Fs = F 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)(1977)(1978)(1979)(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 (Fig 4B).
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 (Fig. 4A, B, C).
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 (Darby and Flatman 1994).The settings chosen for the model are the default ones, with the standard error of the mean F logarithm set equal to 0.5 and the mean applied over 5 years.
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 analysis
Integrated catch at age (ICA) is a methodology that bases its algorithm on a separability assumption (Patterson and Melvin 1996, Needle 2003, De Oliveira et al. 2010).A separable period assumes constant selection-at-age, requiring estimation of a fishing mortality made by a year effect (F 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 (Patterson and Melvin 1996) and, finally, such statistical models give better estimates of the sizes of cohorts still being fished, for which VPA can be unreliable (Needle 2003).ICA was performed using the software developed by Patterson and Melvin (1996).
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).

Natural mortality
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 (Table 3).
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).

VPA results
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 (Table 4).Both the mid-year biomass and the mid-year SSB (Fig. 5) increase constantly after the collapse in 1987, reaching, at the end of the time series, values comparable to the highest recorded before the collapse.
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 (Fig. 6).
The selectivity pattern (Table 5) of the associated fishery for each age class was calculated as the ratio between F 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.
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 (Fig. 7).
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) (Fig. 8).
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 Patterson (1992) as a reference point based on fishing mortality.

ICA results
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 (Table 6).The year and age marginal total residuals for the separable analysis show no clear pattern and are reasonably close to 0, indicating a good fit to the separable constraint (Fig. 9).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 (Fig. 10).
The trend of F by age is represented in Figure 11: the main difference from the VPA model concerns the F from ages 2 to 4, which is overall higher in ICA than in VPA.In particular, the average F from age 2 to 4 estimated from ICA is higher in the time span from 1999 to 2003 than in the year of the collapse; by contrast, for VPA the values in 1986-1987 are the highest of the time series.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 (Fig. 3).
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 Table 7 (it should be noted that the selectivity for age 2 is equal to 1 since this is the reference age, and that we assumed the selectivity for the oldest ages to be 1.6 the selectivity of the reference 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 (Fig. 12).The retrospective results improve for the biomass and are still reasonable for F values (Fig. 13).
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.

DISCUSSION
The present work aims to update and improve the stock assessment of Engraulis encrasicolus in the Adriatic Sea.The last published assessment for this area was performed in 2003 by Santojanni et al. (2003): the author used VPA and the De Lury model, tuned with CPUE data from one representative port on the Italian Adriatic coast and an M value of 0.6 constant for all the age classes.In the present work we tested the performance of VPA vs ICA for the Adriatic stock and we explore how the use of an M-vector affects the results.Furthermore, a fishery-independent tuning index from acoustic data covering the whole area was used in the models.
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 (Abella et al. 1998).On the other hand, Gislason's equation takes into account the asymptotic length and the growth rate K; it also strongly relates M to the body length, raising it by a factor equal to -1.6: this allows the population structure to be taken into consideration.The comparison between these estimates and Gislason's estimates, as well as the exploration of the results given by the VPA and ICA runs, suggests the conclusion that ProdBiom is more suitable for long-lived species in which contribution to the net biomass at sea is more related to the growth rate than for the quick and massive recruitment of new specimens such as short-lived small pelagics.
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 (Azzali et al. 2002) and from the egg-production method used by Regner (1996).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 (Fig. 2), probably because the biomass was in any case available to the fishing fleet.The observed pattern may be due to a change in the age-length frequency distribution: in fact, in the cohorts before 1996 a greater percentage of the oldest ages (4+) with the same length distribution was observed, while after 1996 few of age 4 or older were recorded: only after 2002 did this number slowly start to increase (Fig. 3).This distribution contributed to the decrease in biomass observed in the model.The fixed terminal F method used in the VPA allowed this effect to be softened and increased the biomass estimated for this period, bringing it to a level reasonably higher than that estimated for the years of collapse.On the other hand, ICA needed no additional assumptions (e.g.fixed terminal F) to maintain the biomass in the period 1999-2003 above the values estimated for the collapse period.
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; Darby and Flatman (1994) suggest avoiding fleets with only a few years of data.This pattern will be monitored in the future and some improvements will be adopted if it does not change.
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 (Patterson and Melvin 1996), ICA still offers more reliable estimates of the most recent cohorts, weighting options for surveys and catch data according to belief about validity (Needle 2003).The disadvantage of this flexibility is the greater computational demand of the ICA algorithm, which causes a lower robustness with a higher demand in terms of skills and judgment by the analyst (Shepherd 1999).
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 (Watanabe et al. 1995, Fréon et al. 2005, Payne et al. 2009).In the case of the Adriatic stock, it is strongly believed that the major collapse of the 1987 was mainly due to a failure in the recruitment, whose causes cannot be attributed only to overfishing (as hypothesized by Klanjšček and Legović (2007)); several authors have pointed out that changes in the environmental conditions and exceptional events in that period played a determining role (Regner 1996, Santojanni et al. 2006a, Santojanni 2009).
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 (Alheit et al. 2005(Alheit et al. , 2012)).The fast recovery observed in the biomass trend after the collapse in 1987 is not surprising: for example, several studies on fish population collapses and extinction risk report a high resilience of Clupeidae in comparison with other marine species (Hutching 2001, Hutchings andReynolds 2004).Fluctuations of small pelagics have been observed all over the world and can be driven by changes in environmental conditions, changes in the prey-predator composition, decadal fluctuations dominated by the alternation between small pelagic stocks (e.g.anchovy and sardine), or most likely a combination of all these factors (Checkley et al. 2009).Ruiz et al. (2009), in the Gulf of Cadiz, report a similar case, with a sudden collapse of anchovy in 1994, probably due to a recruitment failure, and a quite fast recovery in the following years.
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 (Fig. 14).
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 2000and 2002. Regner (1996) ) hypothesized the same regarding the collapse in 1987.In his paper he writes: "it seems that during periods of massive blooms fish were forced to reproduce in relatively "clean" oligotrophic waters of the middle and southern Adriatic instead of in its main spawning areas".A more exhaustive study is, however, necessary to address these questions and it is not among the purposes of this paper to answer.
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.
From 2000 to 2009, Italian catches fluctuated between 20400 and 35400 t, Croatian catches were around 8000 t, and Slovenian catches varied between 58 t in 2003 and 438 t in 2005 (FAO 2012).Overall, in GSA 17 average landings of anchovy were about 29000 tfrom 1975 to 2009 and about 35000 tons in the last 10 years (CNR-ISMAR Database; official state statistics), and together with sardine (Sardina pilchardus) it accounted for approximately 41% of total Adriatic catches (average 1992-2008) (FAO 2012).

Fig. 1 .
Fig. 1. -Major fishing ports for small pelagic fish along the Adriatic side of the Italian peninsula and along the Slovenian and Croatian coasts.In dark grey: Geographical Sub-Area 17 (GSA 17), according to the subdivision established from the General Fishery Commission for the Mediterranean (GFCM).

Fig. 3 .
Fig. 3. -Log of the numbers at age in the cohorts from commercial landings from 1976 to 2009.The data presented have already been converted to split-year (see text); the year of each plot refers to the cohort, e.g. the 1972 cohort is represented only from specimens that in 1976 are of age 4+, while, for example, the 2000 cohort is represented from specimens that are of age 0 in 2000, of age 1 in 2001, of age 2 in 2002, and so on.

Fig. 4 .
Fig. 4. -A, trend in landings from 1976 to 2009 for the entire Adriatic fleet (dashed line) and the Porto Garibaldi fleet (full line); B, trend in effort (expressed as numbers of fishing days at sea) for the entire Adriatic fleet (dashed line) and the Porto Garibaldi fleet (full line); C, F vector on the oldest age estimated from Equation 1 for the entire time series.

Fig. 5 .
Fig. 5. -Stock biomass (full line) and spawning stock biomass (dashed line) obtained by Laurec-Shepherd VPA using F vector at the oldest age and shrinking the terminal F to the mean.
Fig. 7. -Stock recruitment relationship for the time period 1976-2009 obtained by the model with fixed terminal F and shrinkage to the mean.Spawners in year x are related to the recruitment in year x+1 (SSB, spawning stock biomass, R, recruitment).

Fig. 8 .
Fig. 8. -Retrospective analysis carried out for the years from 2006 to 2009 for the VPA run.Trend in biomass (left axis) and trend in F (right axis).

Fig. 9 .
Fig. 9. -ICA model catch residuals, displayed by age (to the left) and by year (to the right).

Fig. 10 .
Fig. 10.-Stock biomass (full line) and spawning stock biomass (dashed line) estimated by the means of ICA using 10 years of separability.

Fig. 12 .
Fig. 12. -Stock recruitment relationship for the time period 1976-2009 obtained by the ICA model.Spawners in year x are related to the recruitment in year x+1 (SSB, spawning stock biomass, r, recruitment).Beverton-Holt curve is also plotted (estimated a=58397871, estimated b=294720).

Fig. 13
Fig. 13.-Retrospective analysis carried out for the years from 2006 to 2009 for the ICA run.Trend in biomass (left axis) and trend in F (right axis).

Fig. 14 .
Fig. 14. -Time series from 1976 to 2009 of mean length in the catch obtained from the catch-weighted length-frequency distribution: linear regressions for the whole time series (1976-2009) and for the period 1996-2005 are drawn.Equations and n (number of years) for both regressions are also shown.No data were collected from 1982 to 1984.

Table 2
. -Total biomass estimates and length-frequency distribution (LFD) of the anchovy stock derived from acoustic surveys.The areas covered from acoustic surveys are the following: western, North Adriatic, from the Gulf of Trieste up to Giulianova, available since 1976; western, mid-Adriatic, from Giulianova up to Vieste, available since 1987; eastern Adriatic, from cape Savudrija to cape Oštro, available since 2003.

Table 3 .
-M-vector obtained by the means of the ProdBiom method and Gislason's equation.

Table 5 .
-Selectivity at age for the time period 1998-2009 obtained by the VPA model with fixed terminal F and shrinkage to the mean.
Fig. 6. -Fishing mortality vector by year for ages 0, 1 and 2-4, respectively, obtained by the VPA model with fixed terminal F and shrinkage to the mean.

Table 7 .
-ICA selectivity estimates with standard deviation and coefficient of variation (CV).