sm79n1-3903.html

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

Piera Carpi 1, Alberto Santojanni 1, Fortunata Donato 1, Sabrina Colella 1, Vanja Čikeš Keč 2,
Barbara Zorica 2, Iole Leonori 1, Andrea De Felice 1, Vjekoslav Tičina 2, Tomaž Modic 3,
Polona Pengal 3, Enrico Arneri 4

1 CNR-ISMAR, Largo Fiera della Pesca, 60125, Ancona, Italy. E-mail: piera.carpi@an.ismar.cnr.it
2 Institute of Oceanography and Fisheries, Šetalište I. Meštrovića 63, 21000 Split, Croatia.
3 Fisheries Research Institute of Slovenia, Zupanciceva 9, 1000 Ljubljana, Slovenia.
4 FAO-AdriaMed, Viale delle Terme di Caracalla, Rome, Italy.

Summary: 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: Laurec-Sheperd 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.

Keywords: Engraulis encrasicolus, Adriatic Sea, age-structured stock assessment, ICA, natural mortality, stock collapse.

Evaluación conjunta de la población de la anchoa del Mar Adriático septentrional y central: comparación entre dos métodos de dinámica poblacional

Palabras clave: Engraulis encrasicolus, mar Adriático, modelos analíticos de evaluación de stock, ICA, mortalidad natural, colapso de la población.

Citation/Como citar este artículo: Carpi P., Santojanni A., Donato F., Colella S., Čikeš Keč V., Zorica B., Leonori I., De Felice A., Tičina V., Modic T., Pengal P., Arneri E. 2015. A joint stock assessment for the anchovy stock of the northern and central Adriatic Sea: comparison of two catch-at-age models. Sci. Mar. 79(1): 57-70. doi: http://dx.doi.org/10.3989/scimar.03903.29A

Editor: W. Norbis.

Received: June 5, 2013. Accepted: November 14, 2014. Published: January 8, 2015.

Contents

IntroductionTop

Fishery

Anchovy (Engraulis encrasicolus L.) is one of the most important commercial species in the Adriatic Sea (, ), accounting—together with sardine—for a high percentage of the overall fish production of the surrounding countries, mainly Italy and Croatia (). 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 (). 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 (); 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 (, , ).

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 (, , ). This separation is now less evident, partly due to the changed market situation with the breakdown of the former Yugoslavia (). 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 (, ).

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

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 () and acoustic surveys (, ). The Croatian coast, up to the midline, has been covered since 2003 by the Croatian national pelagic monitoring programme PELMON (). In order to ensure comparability of acoustic estimates, an intercalibration between the two research vessels conducting these surveys was performed ().

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 (), with peaks in May-June and August-September (, , ). 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 (). 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 (). 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 (, ).

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 and , 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) (). Because the Adriatic Sea is an important habitat for anchovies (), their key role as a prey element has been documented through a process-oriented model (Ecopath with Ecosim) applied to the Adriatic ecosystem by , ).

State of the art and aim of the study

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

Up to 2009, Laurec-Shepherd VPA and the De Lury depletion model were mainly used to assess the Adriatic stock of anchovy (, , ): 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 (, ).

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

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 (). 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 (, ), so its performance is not always satisfactory.

In this paper, integrated catch at age (ICA) () 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 (, ). Moreover, in ICA, tuning indices may be either age-structured or based on non–age-structured measures of spawning stock biomass (SSB) (), which can be useful in the case of small pelagic data.

ICA has been successfully applied to other small pelagic short-lived species (, ) 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.

Materials and methodsTop

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

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

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.

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.

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 (, , ); 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 (, ,).

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

 $M=[ B t2−t1 ]ln(t2−t1)+A,t>0$ (1)

and (ii) the empirical equation proposed by (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, t0=-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 assessment was carried out taking into account a conventional birth date on 1 June (split-year), as previously done by . 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 (, ). 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 (E=0.4).

Virtual population analysis

The analyses were performed using the MAAF-VPA software package developed by . Several model configurations were tested, varying the way in which the terminal F is estimated. The model allows two types of Laurec-Shepherd tuning (, ), 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 = Fbar(n)·k. For example, if the last true age (oldest) is set equal to 3, n equal to 2 (so Fbar 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 (Fx) 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 (). 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 (, , ). A separable period assumes constant selection-at-age, requiring estimation of a fishing mortality made by a year effect (Fy) and an age effect (Sa), such that F at age is a multiplicative combination of these (Fy,a=Fy·Sa). 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 () and, finally, such statistical models give better estimates of the sizes of cohorts still being fished, for which VPA can be unreliable (). ICA was performed using the software developed by .

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

RESULTSTop

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 Linf value (Table 3).

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

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

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

Table 4. – Log catchability, standard error and slope of the regression.

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 (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 Fa,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.

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.

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

Table 6. – ICA distribution statistics for acoustic survey (ages 0 to 4+).

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

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

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

DiscussionTop

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 : 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 (). 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 () and from the egg-production method used by .

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; 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 (), ICA still offers more reliable estimates of the most recent cohorts, weighting options for surveys and catch data according to belief about validity (). 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 ().

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 (, , ). 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 ); several authors have pointed out that changes in the environmental conditions and exceptional events in that period played a determining role (, , ).

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 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 (, ). 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 (). , 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 2000 and 2002. 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.

ACKNOWLEDGEMENTSTop

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.

ReferencesTop

Abella A., Caddy J.F., Serena F. 1998. Estimation of the parameters of the Caddy reciprocal M-at-age model for the construction of natural mortality vectors. In: Lleonart J. (ed.). Dynamique des populations marines. Zaragoza: CIHEAM. p. 191-200. (Cahiers Options Méditerranéennes; n. 35). Deuxième Réunion du Groupe de Travail DYNPOP, 1996/10/02-05, Genova (Italy).

Alheit J., Mollmann C., Dutz J., et al. 2005. Synchronous ecological regime shifts in the central Baltic and the North Sea in the late 1980s. ICES J. Mar. Sci. 62: 1205-1215.
http://dx.doi.org/10.1016/j.icesjms.2005.04.024

Alheit J., Pohlmann T., Casini M., et al. 2012. Climate variability drives anchovies and sardines into the North and Baltic Seas. Prog. Oceanogr. 96(1): 128-139.
http://dx.doi.org/10.1016/j.pocean.2011.11.015

Antonakakis K., Giannoulaki M., Machias A., et al. 2011. Assessment of the sardine (Sardina pilchardus Walbaum, 1792) fishery in the eastern Mediterranean basin (North Aegean Sea). Medit. Mar. Sci. 12: 333-357.
http://dx.doi.org/10.12681/mms.36

Azzali M., De Felice A., Luna M., et al. 2002. The state of the Adriatic Sea centered on the small pelagic fish populations. PSZNI Mar. Ecol. 23: 78-91.
http://dx.doi.org/10.1111/j.1439-0485.2002.tb00009.x

Beverton R.J., Holt S.J. 1956. A review of methods for estimating mortality rates in exploited fish populations, with special reference to sources of bias in catch sampling. Rap. Proc.-verb. Réun. Cons. Inter. Exp. Mer 140: 67-83.

Checkley D., Alheit J., Oozeki Y., et al. (eds) 2009. Climate change and small pelagic fish. Cambridge Univ. Press. 372 pp.

Cingolani N., Santojanni A. 2003. Manual of the Recorder. AdriaMed Training Course on Data Collection and Biological Sampling System on Small Pelagics. AdriaMed Occasional Papers. No.6 (Rev.1). GCP/RER/010/ITA/OP-06-Rev.1, Termoli, 2003: 53 pp.

Cingolani N., Giannetti G., Arneri E. 1996. Anchovy fisheries in the Adriatic Sea. Sci. Mar. 60: 269-277.

Cingolani N., Santojanni A., Arneri E., et al. 2008. Programma nazionale per la raccolta dei dati alieutici -ModuloH- Campionamento biologico delle catture: composizione per età e per lunghezza (risorse demersali e piccoli pelagici). Relazione finale CNR-ISMAR.

Coll M., Santojanni A., Palomera I., et al. 2007. An ecological model of the Northern and Central Adriatic Sea: analysis of ecosystem structure and fishing impacts. J. Mar. Syst. 67: 119-154.
http://dx.doi.org/10.1016/j.jmarsys.2006.10.002

Coll M., Santojanni A., Palomera I., et al. 2009. Food-web changes in the Adriatic Sea over the last three decades. Mar. Ecol. Prog. Ser. 381: 17-37.
http://dx.doi.org/10.3354/meps07944

Darby C.D., Flatman S. 1994. Virtual Population Analysis: Version 3.1 (Windows/Dos) user guide. Info. Tech. Ser. No1. Lowestoft.

Daskalov G.M., Mamedov E.V. 2007. Integrated fisheries assessment and possible causes for the collapse of anchovy kilka in the Caspian Sea. ICES J. Mar. Sci. 64: 503-511.
http://dx.doi.org/10.1093/icesjms/fsl047

De Oliveira J.A.A., Darby C.D., Roel B.A. 2010. A linked separable–ADAPT VPA assessment model for western horse mackerel (Trachurus trachurus), accounting for realized fecundity as a function of fish weight. ICES J. Mar. Sci. 67: 916-930.
http://dx.doi.org/10.1093/icesjms/fsp290

FAO 2012. Fishery Statistical Collections: Global Production Statistics 1950-2010 (online query).

Fréon P., Cury P., Shannon L., et al. 2005. Sustainable exploitation of small pelagic fish stocks challenged by environmental and ecosystem changes: a review. Bull. Mar. Sci. 76: 385-462.

GFCM 2008. Report of the ninth meeting of the working group on small pelagics. Technical report General Fisheries Commission for the Mediterranean, Sub-Committee on Stock Assessment (SCSA).

GFCM 2013a. GFCM e-glossary: http://www.gfcm.org/gfcm/topic/17107/en

Giannoulaki M., Iglesias M., Tugores M. P., et al. 2013. Characterizing the potential habitat of European anchovy Engraulis encrasicolus in the Mediterranean Sea, at different life stages. Fish Oceanogr. 22(2): 69-89.
http://dx.doi.org/10.1111/fog.12005

Gislason H., Daan N., Rice J. C., et al. 2009. Size, growth, temperature and the natural mortality of marine fish. Fish Fish. 11: 149-158.
http://dx.doi.org/10.1111/j.1467-2979.2009.00350.x

Grbec B., Dulčić J., Morović M. 2002. Long-term changes in landings of small pelagic fish in the eastern Adriatic - Possible influence of climate oscillations over the Northern Hemisphere. Clim. Res. 20: 241-252.
http://dx.doi.org/10.3354/cr020241

Hutching J.A. 2001. Influence of population decline, fishing and spawner variability on the recovery of marine fishes. J. Fish. Biol. 59: 306-322.
http://dx.doi.org/10.1111/j.1095-8649.2001.tb01392.x

Hutchings J.A., Reynolds J.D. 2004. Marine fish population collapses: consequences for recovery and extinction risk. BioScience. 54: 297-309.
http://dx.doi.org/10.1641/0006-3568(2004)054[0297:MFPCCF]2.0.CO;2

IREPA 2009. Osservatorio economico sulle strutture produttive della pesca marittima in Italia 2009. Istituto di ricerche economiche per la pesca e l’acquacoltura (IREPA). 180 pp.

Kimura D.K., Chikuni, S. 1987. Mixtures of empirical distributions: an iterative application of the Age-Length Key. Biometrics. 43: 23-35.
http://dx.doi.org/10.2307/2531945

Klanjšček J., Legović T. 2007. Is anchovy (Engraulis encrasicolus, L.) overfished in the Adriatic Sea? Ecol. Model. 201: 312-316.
http://dx.doi.org/10.1016/j.ecolmodel.2006.09.020

Laurec A., Shepherd J.G. 1983. On the analysis of catch and effort data. J. Cons. Int. Explor. Mer. 41: 81-84.
http://dx.doi.org/10.1093/icesjms/41.1.81

Leonori I., De Felice A., Campanella F., et al. 2011. Assessment of small pelagic fish biomass in the western Adriatic Sea by means of acoustic methodology. In: Brugnoli E., Cavarretta G., Mazzola S., et al. (eds), Marine research at CNR Vol. DTA/06 Fishery and Sea Resources. Rome: pp. 2019-2029.

Leonori I., Tičina V., Felice A.D., et al. 2012. Comparison of two research vessels’ properties in the acoustic surveys of small pelagic fish. Acta Adriat. 53: 389-398.

Magoulas A., Castilho R., Caetano S., et al. 2006. Mitochondrial DNA reveals a mosaic pattern of phylogeographical structure in Atlantic and Mediterranean populations of anchovy (Engraulis encrasicolus). Mol. Phylogenet. Evol. 39: 734-746.
http://dx.doi.org/10.1016/j.ympev.2006.01.016

Mannini P., Massa F. 2000. Support paper prepared for the first AdriaMed Coordination Committee. FAO-MiPAF Scientific Cooperation to Support Responsible Fisheries in the Adriatic Sea. GCP/RER/010/ITA/TD-01: 64 pp.

Marčeta B. 2001. Status of Slovene research and fishery on small pelagics. ADRIAMED Tech. Doc. 3: 24-29.

Morello E.B., Arneri E. 2009. Anchovy and sardine in the Adriatic Sea - An ecological review. Oceanogr. Mar. Biol. 47: 209-256.

Needle C.L. 2000. The Ins and Outs of ICA. Technical Report 04/00 Marine Laboratory, Aberdeen. 80 pp.

Needle C.L. 2003. Course on Fish Stock Assessment Techniques. Workshop Course on Fish Stock Assessment Techniques. ICES Headquarters, Copenhagen. 182 pp.

Patterson K. 1992. Fisheries for small pelagic species: an empirical approach to management targets. Rev. Fish Biol. Fisher. 2: 321-338.
http://dx.doi.org/10.1007/BF00043521

Patterson K.R., Melvin G.D. 1996. Integrated Catch at Age Analysis - Version 1.2. Technical Report 58 The Scottish Office - Agriculture, Environment and Fisheries Department.

Pauly D. 1980. On the interrelationships between natural mortality, growth parameters, and mean environmental temperature in 175 fish stocks. J. Conseil. 39: 175-192.
http://dx.doi.org/10.1093/icesjms/39.2.175

Payne M.R., Hatfield E.M.C., Dickey-Collas M., et al. 2009. Recruitment in a changing environment: the 2000s North Sea herring recruitment failure. ICES J. Mar. Sci. 66: 272-277.
http://dx.doi.org/10.1093/icesjms/fsn211

Pope J., Shepherd J.G. 1985. A comparison of the performance of various methods for tuning VPA’s using effort data. J. Cons. Int. Explor. Mer. 42: 129-151.
http://dx.doi.org/10.1093/icesjms/42.2.129

Quinn T.J., Deriso R.B. 1999. Quantive Fish Dynamics. Oxford Univ. Press, 542 pp.

Rampa R., Arneri E., Belardinelli A., et al. 2005. Length at first maturity of the Adriatic anchovy (Engraulis encrasicolus L.). Document presented at the General Fisheries Commission for the Mediterranean (GFCM), Scientific Advisory Committee (SAC), Sub Committee on Stock Assessment (SCSA). Rome.

Regner S. 1985. Ecology of planktonic stages of the anchovy, Engraulis encrasicolus (Linnaeus, 1758), in the central Adriatic. Acta Adriat. 26: 5-113.

Regner S. 1996. Effects of environmental changes on early stages and reproduction of anchovy in the Adriatic Sea. Sci. Mar. 60: 167-177.

Ruiz J., Gonzalez-Quiros R., Prieto L., et al. 2009. A Bayesian model for anchovy (Engraulis encrasicolus): the combined forcing of man and environment. Fish Oceanogr. 18: 62-76.
http://dx.doi.org/10.1111/j.1365-2419.2008.00497.x

Santojanni A. 2009. Comments on “Is anchovy (Engraulis encrasicolus L.) overfished in the Adriatic Sea?” by Klanjšček and Legović [Ecol. Model. 201 (2007) 312–316]. Ecol. Model. 220: 430-433.
http://dx.doi.org/10.1016/j.ecolmodel.2008.05.003

Santojanni A., Arneri E., Barry C., et al. 2003. Trends of anchovy Engraulis encrasicolus L.) biomass in the northern and central Adriatic Sea. Sci. Mar. 67: 327-340.

Santojanni A., Cingolani N., Arneri E., et al. 2005. Stock assessment of sardine (Sardina pilchardus Walb.) in the Adriatic Sea, with an estimate of discards. Sci. Mar. 69: 603-617.

Santojanni A., Arneri E., Bernardini V., et al. 2006a. Effects of environmental variables on recruitment of anchovy in the Adriatic Sea. Clim. Res. 31: 181-193.
http://dx.doi.org/10.3354/cr031181

Santojanni A., Cingolani N., Arneri E., et al. 2006b. Use of an exploitation rate threshold in the management of anchovy and sardine stocks in the Adriatic Sea. Biol. Mar. Med. 13: 98-111.

Santojanni A., Leonori I., Carpi P., et al. 2011. Report of the working group on stock assessment of small pelagic species. Technical report Scientific Advisory Committee, Sub-Committee on Stock Assessment, Chania, Greece.

SARDONE 2010. Improving assessment and management of small pelagic species in the Mediterranean. Deliverable D18: Final report of WP5, Improved Stock assessments for sardine and anchovy in the Mediterranean Sea (CNR).

Shepherd J. G. 1999. Extended survivors analysis: an improved method for the analysis of catch-at-age data and abundance indices. ICES J. Mar. Sci. 56: 584-591.
http://dx.doi.org/10.1006/jmsc.1999.0498

Sinovčić G., Zorica B. 2006. Reproductive cycle and minimal length at sexual maturity of Engraulis encrasicolus (L.) in the Zrmaja River estuary (Adriatic Sea, Croatia). Estuar. Coast Shelf Sci. 69: 439-448.
http://dx.doi.org/10.1016/j.ecss.2006.04.003

Tičina V. 2003. Pelagic resources of the Adriatic Sea. Croatian International Relations Review 9(32): 33-35.

Tičina V., Katavić I., Dadić V., et al. 2006. Acoustic estimates of small pelagic fish stocks in the eastern part of Adriatic Sea. Proceedings of GFCM-SAC WG on Identification of Reference Points for the management of fishery resources, Roma, 20-21 April 2004. Biol. Mar. Med. 13: 124-136.

Tudela S. 1999. Morphological variability in a Mediterranean, genetically homogeneous population of the European anchovy, Engraulis encrasicolus. Fish. Res. 42: 229-243.
http://dx.doi.org/10.1016/S0165-7836(99)00052-1

Watanabe Y., Zenitani H., Kimura R. 1995. Population decline of the Japanese sardine Sardinops melanostictus owing to recruitment failures. Can. J. Fish. Aquat. Sci. 52: 1609-1616.
http://dx.doi.org/10.1139/f95-154

Zorica B., Vilibić I., Čikeš Keč V., et al. 2013. Environmental conditions conducive to anchovy (Engraulis encrasicolus) spawning in the Adriatic Sea. Fish. Oceanogr. 22: 32-40.
http://dx.doi.org/10.1111/fog.12002