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 and Massa 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. 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).

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).

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 and Arneri 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. (2007, 2009).

The most detailed and recent published stock assessment on this species dates back to 2003 and comprises data from 1975 to 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. Furthermore, 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 and Melvin 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 and Mamedov 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.

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 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 (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, 2006a,b).

In the present work, two options for estimating a natural mortality vector by age were analysed; the methods explored were (i) the ProdBiom software (Abella et al. 1998) (Eq. 1):

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

and (ii) the empirical equation proposed by Gislason (2009) (Eq. 2):

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₀=-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 (Cingolani et al. 2008).

Since the spawning takes place mostly in spring-summer (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ć and Zorica 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).

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 and Shepherd 1983, Pope and Shepherd 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-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 (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).

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).

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.

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.