Stock assessment of the western winter-spring cohort of Ommastrephes bartramii in the Northwest Pacific Ocean using a Bayesian hierarchical DeLury model based on daily natural mortality during 2005-2015 ; Evaluación del stock de la cohorte occidental de invierno-primavera de Ommastrephes bartramii 

1 College of Marine Sciences, Shanghai Ocean University, Shanghai 201306, China. 2 Key Laboratory of Sustainable Exploitation of Oceanic Fisheries Resources, Ministry of Education, Shanghai Ocean University, Shanghai 201306, China. (QD) E-mail: dingqi@ysfri.ac.cn. ORCID iD: https://orcid.org/0000-0001-7181-7059 (XC) (Corresponding author) E-mail: xjchen@shou.edu.cn. ORCID iD: https://orcid.org/0000-0001-6311-6842 3 School of Aquatic and Fisheries Sciences, University of Washington, Box 355020, Seattle, WA, 98195, USA. (JC) E-mail: souketu@gmail.com. ORCID iD: https://orcid.org/0000-0002-9291-9182 4 National Engineering Research Centre for Oceanic Fisheries, Shanghai Ocean University, Shanghai 201306, China.


INTRODUCTION
The neon flying squid (Ommastrephes bartramii), one of the most abundant and economically important oceanic squid in the family Ommastrephidae, is widely distributed in subtropical and temperate waters of the global ocean (Roper et al. 1984, Arkhipkin et al. 2015. Only the North Pacific population is fished commercially. The North Pacific population is widely distributed between 20°N and 50°N, and comprises two cohorts (the autumn-spawning cohort and the winter-spring cohort) and four stocks: (1) the central stock of the autumn cohort, (2) the eastern stock of the autumn cohort, (3) the western stock of the winter-spring cohort, and (4) the central-eastern stock of the winter-spring cohort (Bower and Ichii 2005). The western winter-spring stock is the main fishing target of the Chinese squid-jigging vessels. This stock annually migrates between spawning grounds (subtropical waters) and feeding grounds (the sub-Arctic domain) Nakamura 1998, Ichii et al. 2006), and the latter grounds are where the fishery takes place (Fig.  1). Chinese squid-jigging fleets mainly operate in the waters of 150-165°E, where over 80% of the total catch occurs from August to November (Wang and Chen 2005, Chen et al. 2008b, Cao et al. 2015. Three fishery companies from the Chinese Mainland (Zhoushan marine fisheries corporation, Ningbo marine fisheries corporation, and Yantai marine fisheries corporation) conducted the first survey on O. bartramii in the Northwest Pacific Ocean in 1993 using squid-jigging vessels, and a large-scale jigging fishery by Chinese vessels started in 1994. Annual catch of O. bartramii in China gradually increased from 23000 t in 1994 to the maximum value 132000 t in 1999, and ranged from 80000 to 124000 t from 2000 to 2008. However, annual catches substantially decreased in the period 2009 to 2015, ranging from 31000 to 55000 t (Fig. 2). The number of vessels gradually increased from 98 to the maximum of 446 in 2000. Afterwards, it showed an overall downward tendency and decreased to 116 in 2015 (Fig. 2).
The neon flying squid matures at 7 to 10 months and has an estimated one-year lifespan (Yatsu et al. 1997). This unique life history limits the use of traditional age-or length-structured models to evaluate the effect of intensive commercial jigging fishery on this stock. DeLury depletion models (De Lury 1947) are often used for the stock assessment and management of short-lived species, and have been successfully applied in the stock assessment and management of squid and crustacean species (Young et al. 2004, Roa-Ureta and Arkhipkin 2007, Roa-Ureta 2012. However, estimates from depletion models can be misleading when sampling error variation masks trends in abundance (McAllister et al. 2004). Hierarchical modelling offers an approach that allows for variability in parameter estimates, while minimizing the potential for overspecifying the model. Previous studies indicated that a Bayesian hierarchical DeLury depletion model can produce more reliable estimates than similar non-hierarchical methods, and can better quantify the variability in parameters by combining all sets of time-series data hierarchically in the same model (McAllister et al. 2004, Cao et al. 2015.
There are currently no management measures for O. bartramii fisheries on the high seas of the North Pacific other than the driftnet moratorium (Arkhipkin et al. 2015). Because of the annual lifecycle of squid species and their high variability in abundance from year to year, fishery management for squid should implement real-time management, i.e. set target proportional escapement based on the catch and effort data during the period of depletion. To quantify the effect of commercial jigging fishery on the stock of the western winter-spring cohort, Chen et al. (2008a) conducted a stock assessment using a modified depletion model. Assessment results indicated that proportional escapement was close to the management target of 40% from 2000 to 2005, and fishing mortality of the squid-jigging fishery was considered to be sustainable. However, annual catches have decreased significantly, and the management target has not been reached since then. This paper used Bayesian hierarchical depletion models to estimate the abundance and daily fishing mortality of neon flying squid based on Chinese squidjigging fishery data during the period 2005-2015, which is compared with previous studies with 10-day natural mortality (Chen et al. 2008a, Cao et al. 2015. These results are combined with previous estimates to provide a 16-year time series of stock size and fishing pressure. The study aimed to evaluate the escapement rate for the western winter-spring cohort of neon flying squid in recent years and to provide an update of the exploitation status of this stock.

Catch data
The western winter-spring cohort of neon flying squid is mainly distributed in the area of 150-165°E in the Northwest Pacific Ocean (Fig. 1). The fishing season starts in June and ends in late November or early December, peaking from August to October. Therefore, the squid fishery can be treated as a depletion experiment from its peak to early December (Chen et al. 2008a). The annual catch from China comprises more than 80% of the total catch in this area. The squid fishery operates at night and there is no bycatch. Data on daily catch (kg), fishing effort (days fished) and fishing locations were obtained from the Chinese mainland commercial jigging fleets between August and December from 2005 to 2015.

Body weight
Body weight is used to convert catches (kg) to individuals for catch per unit effort (CPUE). Fifteen individuals were sampled on each day during the fishing season from 2010 to 2015, and the weight of each individual was measured. Not all days were included in each month due to data limitation. Comparisons using a Kruskal-Wallis test indicated no significant difference among years in mean body weight for each 10-day period. Therefore, we pooled the weight data across years to estimate the daily growth (Fig. 3).

Abundance index
In accordance with Cao et al. (2011), standardization of CPUE was conducted to eliminate the effects of spatial, temporal and environmental factors that bias CPUE as an index of abundance. A generalized linear model including five variables was considered to standardize the CPUE: where Date i is the date i, lon is longitude, lat is latitude, SST is sea surface temperature and SS H is sea surface height. The environmental data were obtained from the Ocean-Watch Data Library (http://oceanwatch.pifsc. noaa.gov).
We identified the best-fitting model using a stepwise approach based on the Akaike information criterion. The best fitted model was used in the standardization. The standardized CPUE of each year shows a prominent decreasing trend (Fig. 4).

Depletion period
We treated the squid fishery targeting the western winter-spring cohort (Fig. 1) as a depletion experiment.
Since the CPUE generally peaks in August, we defined the depletion period as between August and November each year. By doing so, we assume that the majority of squids migrate to the feeding ground in August, and from then the stock is considered as a closed stock. We used daily CPUE to track the depletion in this study, as opposed to the weekly or 10-day CPUE used in the previous studies.

Natural mortality
Natural mortality rate (M) of O. bartramii is largely unknown, despite its great impact on the stock assessment (Chen et al. 2008a). In order to address the uncertainty of M, a sensitivity analysis was conducted in this study. We applied a series of M values ranging from 0.003 to 0.01 (per day) with a 0.001 increment and compared the outcomes of the assessments. Note that M is assumed to be constant over time.

Depletion model
The DeLury model formulation of Rosenberg et al. (1990) was used to estimate the initial population size, N 0,y , and the catchability coefficient, q y , during each fishing season beginning on 1 August. The predicted CPUE on a given day d (beginning on 1 August) for year y was calculated as where q y is the catchability coefficient, N 0,y is the initial abundance of year y, C y,k is the catch in numbers taken in day k, and M is the daily rate of natural mortality, which is assumed to be constant over the fishing period. The likelihood function of the observed CPUE given the prediction from a model is assumed to be lognormal (Rosenberg et al. 1990): where CPUE y represents a vector of daily CPUE data for year y, ndk(y) is the number of fishing days in year y, CPUE y,d is the observed catch rate on day d, τ L is the precision that is assumed to be constant across years (τ L =1/σ 2 and σ was the standard deviation), and CPUE ! y,d is the predicted catch rate on day d for year y. Bayesian estimation requires specification of prior information for random variable parameters. We assumed that N 0,y and q y were uncorrelated, and there was a population of N 0,y and q y over years and that a value in any one year came from a particular parameter probability density function (pdf). Following Cao et al. (2015), we also assumed that the priors for N 0,y and q y follow lognormal distributions with a mean of μ N 0 and μ q and a precision of τ N 0 and τ q , respectively: μ q , μ N 0 , τ q , and τ N 0 were called hyper-parameters, which describe mean values and uncertainty for parameters for the entire population, respectively. Following the approach of McAllister et al. (2004), all prior and hyper-prior probability density functions (pdfs) applied in the hierarchical model were assumed to be normally distributed with a very large prior coefficient of variation of 100 in order to ensure that posterior pdfs were driven by the data. Considering that sampling error variation may masks trends in abundance, the depletion model was implemented only on the set of days on which a decline in catch rates was observed in a given fishing season, which is from the peak value of CPUE to the end of fishery (McAllister et al. 2004).

Calculation of the proportion not harvested
The escapement ratio is measured as a proportion of squid that would have been alive at the end of the fishing season. The proportional escapement was calculated following Brunetti et al. (2000). Stock size at time i+1 in year y was modelled as where N iy was the number of squid at time i, M was the instantaneous natural mortality, and C iy was the catch number if the catch occurred at the middle of each interval i to i+1. In general, the escapement ratio P iy is calculated as where N 0y is the number of squid at the start of the fishing season, which was 1 August in this study.

Sensitivity analyses of body weight
Given that body weight data were limited during the period 2005-2009, sensitivity analyses were conducted to evaluate the effect of between-year variation of body weight in the study period. We explored the sensitivity of the depletion model to the body weight values in 2005-2009 that were lower than 20%, lower than 10%, higher than 10%, and higher than 20% of the corresponding values calculated on the basis of the derived daily growth function in this study.  (2007) individuals. The central tendencies for q y also vary considerably, from 11.68×10 -6 (2006) to 117.13×10 -6 (2015). The posterior coefficients of variation (CV) in estimates of N 0,y ranged from 0.010 to 0.21 and the posterior CVs in estimates of q y ranged from 0.048 to 0.24. The ranges of posterior CVs for the marginal distribution for N 0,y and q y indicate that the data for each year provide different information. Murata and Shimazu (1982) estimated that M for the winter-spring cohort may be around 0.007 per day. Based on this value, the posterior mean values of N 0,y and q y ranged from 109.77 to 414.41 million individuals and 17.98×10 -6 to 69.71×10 -6 , respectively (Figs 5, 6). All the posterior distributions converged on the basis of the convergence diagnostic of Geweke (1992). Hierarchical modelling quantifies how a parameter varies across different members of a population. In this case, the population is the set of years in which DeLury assessments have been applied. In our study, the units are all sets of years (i.e. 2005-2015), the parameters  were assumed to come from lognormal distributions, and the distribution of estimated parameters was used to describe the variation of parameters across units of the population. The posterior estimated mean values for N 0,y and q y were 199.11 million individuals and 36.97×10 -6 under M=0.007, respectively.

Interannual trends in initial and final population estimates
Initial abundance and abundance at the end of fishing season in each year were derived from depletion models. For consistency, initial abundance was the stock size on 1 August , and the abundance at the end of fishing season was the stock size at the end of November. When the the day-to-day progression of estimated population size was calculated over the period 2005-2015, the initial abundance varied considerably across the years (Fig. 7).  (Fig. 7). Figure 8 shows the relationships between initial abundance in the current year and final stock size in the previous year based on M=0.007. The years 2001-2008 had high initial abundance and high population size at the end of previous fishing season, whereas these values have been at low levels since 2009.

Estimation of proportional escapement
The proportional escapement P iy for each fishing season over the period 2005-2015 are shown in Figure  9. As natural mortality rate M increased, the P iy also increased. The wide range of M (0.003 and 0.01) used in this study generated P iy estimates differing between 8.94% (2008) and 70.56% (2015), with an average of 38.58%, which was close to the management target of 40%. Based on M of 0.007, P iy ranged from 14.04% (in  For different M given in our analysis, the mean value of P iy throughout the study period changed from 13.74% in 2008 to 47.39% in 2014 (Fig. 9). This stock suffered from overexploitation for 2008 and 2010, and P iy of these two years was lower than 20% and 35%, respectively, for all the M values used in this study (Fig.  9). The mean values of P iy met the management target of 40% in the years 2011-2015.

Sensitivity analysis of body weight
Sensitivity analyses of body weight in [2005][2006][2007][2008][2009] were implemented in the Bayesian hierarchical depletion model for daily body weight values ranging from lower than 20% to higher than 20% of the values calculated on the basis of the derived daily growth function. The results suggested that only N 0,y was sensitive to daily body weight, while q y was quite robust to changes in body weight within the ranges tested (Table 1). Specifically, using body weight values in 2005-2009 that were lower than 20%, lower than 10%, higher than 10% and higher than 20% of the values in 2010-2015 that derived from the daily growth function, the estimated values of N 0,y for the years 2005-2009 were shown to have generally decreased with the increase in body weight (Table 1). With body weight values lower than 20%, lower than 10%, higher than 10% and higher than 20% of the derived value, the estimated N 0,y were higher than 25%, higher than 11%, lower than 9% and lower than 17% of the corresponding values estimated with body weight equal to the mean values in 2010-2015, respectively. In addition, with daily M values ranging from 0.003 to 0.01, the above percentages were almost constant.
It should be noted that proportional escapement was also quite robust to changes in body weight values within the ranges tested. Using M of 0.007, the fluctuation of proportional escapement was lower than 1% under the extreme body weight values used in this study (Table 1). Similarly, with daily M values ranging from 0.003 to 0.01, the above percentage were almost constant.

DISCUSSION
Given the records of catch and effort from the Chinese Squid-Jigging Technology Working Group, this study conducted a Bayesian hierarchical DeLury depletion model in which data from the years 2005-2015 were combined hierarchically in the same stock assessment model to estimate the population abundance and annual catchability coefficients of the western winterspring cohort of neon flying squid. The results indicated that O. bartramii may have been overexploited since 2008, resulting in low levels of stock size and annual catch since 2009. We further combine our results with previous estimates to provide a 15-year time series of stock size and fishing pressure. The Chinese squid-jigging fishery operated in similar waters on the traditional fishing ground and the fishing technique was almost the same during the period 2001-2015. However, annual catch from 2009 to 2015 was significantly lower than the catch from 2001 to 2008. The above catch trend was in accordance with the trend of initial abundance, i.e. initial abundances in the years 2001 to 2008 were obviously higher than those in the years 2009-2015 (Fig. 8).
Because of the annual lifecycle, managers have very little information on the potential size of the exploitable stock until shortly before recruitment, so the short lifespan of squids requires a different management approach to that taken for most finfish fisheries. Caddy (1993) proposed that the management of squids should be based on effort limitation, with the possibility of short-term adjustment of effort, and with the objective of allowing a maximum proportion (40%) of the catchable biomass to be removed each year. The scheme outlined by Caddy (1993) still remains valid today and has the advantage that it enables a fishery to be managed in real time, which is especially important in volatile stocks of short-lived species such as squids. The proportional escapement was much lower than 40% in the year 2008 (with an average of 13.74%), Table 1. -Estimated mean values of initial abundance, catchability coefficient, and proportional escapement of different body weight using a Bayesian hierarchical model with M=0.007. "Basic" refers to body weight values in 2005-2009 calculated on the basis of the derived daily growth function. "Lower than 20%", "Lower than 10%", "Higher than 10%", and "Higher than 20%" refer to body weight values in 2005-2009 lower than 20%, lower than 10%, higher than 10%, and higher than 20% of the corresponding values calculated on the basis of the derived daily growth function, respectively. Interannual variation in initial population and catchability coefficients was significant across the years. The significant between-year variations in initial abundance resulted in O. bartramii failing to meet the management target in 2010. Annual catch and number of fishing vessels were extremely similar between the years 2010 and 2013. However, initial abundance in 2013 was 1.3 times greater than that in 2010 under M=0.007, which led to the exploitation level meeting the management target of 40% annual escapement in 2013 but being overexploited in 2010. As they are short-lived ecological opportunists, the abundance and distribution of O. bartramii were extremely vulnerable to the multi-scale environmental conditions, especially when anomalous environmental conditions occurred (Chen et al. 2007, Cao et al. 2009, Yu et al. 2016. For example, Yu et al. (2015) found that strong linkage existed between anomalous environmental events and variability in the O. bartramii habitats, and the El Niño-Southern Oscillation (ENSO) event tended to result in unfavourable habitats and reductions in squid recruitment, whereas the La Niña event yielded more favourable habitats and yielded high squid catches. In this study, the lowest initial abundance occurred in the years 2009 and 2015, which coincided with the occurrence of an El Niño event. Heavy fishing pressure together with unfavourable environmental conditions may have resulted in the sharply decreased initial abundance of O. bartramii in 2009. Thus, variations in the above two parameters may be related to environmental conditions. In addition, many other characteristics, such as biological factors and fishing processes, may affect catchability in the squid-jigging fishery (Chen et al. 2008a). Without more data, it is hard to identify the main factors leading to the annual variation in catchability.
The quality of results obtained from a stock assessment depends on the assumptions made. One primary assumption from this analysis was the closed population which ignored the immigration and emigration during the modelled time period. However, climate changes such as the ENSO could result in interannual variation in the timing of peak recruitment migration to the fishing grounds (Chen et al. 2007, 2012, Yu et al. 2015, and the lagged migration may mask the decline in abundance. Given that the depletion period is relatively short, i.e. only around four or five months, this study used a per day time interval, which can provide more accurate assessments. The DeLury estimators are sensitive to the value of natural mortality, but M cannot always be determined with certainty (Chen et al. 2008a (Hendrickson et al. 1996), and Murata and Shimazu (1982) estimated that M for the winter-spring cohort of neon flying squid may be around 0.007 per day. To improve our understanding of impacts of M on the assessment, we conducted an extensive sensitivity analysis of the Bayesian hierarchical depletion model to per-day M values ranging from 0.003 to 0.01 in this study, and depletion model estimators were demonstrated to be sensitive to the value of natural mortality, M. With increasing M, estimated initial abundance increased, while estimated catchability decreased.
Bayesian hierarchical modelling approaches are suitable for modelling the variability of parameters underlying different sets of historical data by sharing information among populations probabilistically, and are also useful for years when sampling error variation masks trends in abundance, and the catch rate data are not well behaved (McAllister et al. 2004). The reliability of stock assessment is obviously subject to data quality. The data used in this study were fishery-dependent data from the Chinese mainland squid-jigging fleet, which accounted for >80% of the total catch from the study area. Considering that a small percentage of catch data were not available in this analysis, we standardized the CPUE before conducting stock assessment, which may be more reliable than using nominal CPUE. A Bayesian hierarchical model can better quantify the variability in parameters and provide more reliable estimates by combining all sets of time-series data hierarchically in the same stock assessment model. Furthermore, the posterior distribution derived from a Bayesian hierarchical model using all the historical data combined can serve as a prior distribution for the current year's assessment (Adkison and Su 2001, Mc-Allister et al. 2004, Gelman et al. 2014. The Bayesian hierarchical depletion model is helpful near the beginning of the fishing season when only a few days of catch rate data are available and estimates from the conventional estimation method for the DeLury model can be highly imprecise. However, the development of prior distributions for future stock assessment may be problematic if the range of hierarchical estimates from the historical years do not reflect the range in future years (McAllister et al. 2004). The estimators of the hierarchical model could be improved as more data become available in the future.
Considering the short lifespan and high annual variability in initial abundance, fishery management based on effort limitation is much flexible for conserving short-lived stocks such as O. bartramii at a sustainable level (Arkhipkin et al. 2015). This study showed that O. bartramii may suffer from overexploitation, and the sharply decreased initial abundance of O. bartramii from 2009 to 2015 might be related to the heavy fish-ing pressure together with the unfavourable environmental conditions. Although current fishing mortality of the squid-jigging fishery is generally considered to be sustainable, the catches and stock abundance still fluctuate at low levels in recent years, leading to a possible shift in hypothesis. The environmental conditions change and may not be as favourable as they were prior to 2008, so decreasing fishing effort may not lead to increasing abundance under current environmental conditions. This should be an interesting hypothesis for future studies. Significant between-year variations in initial abundance resulted in O. bartramii suffering from a certain degree of overexploitation in 2010. For short-lived species such as O. bartramii, recruitment success dominates annual variations in stock abundance and is extremely sensitive to multiscale environmental events such as the ENSO and La Niña (Sakurai et al. 2000, Yu et al. 2015, which may result in large annual variability in abundance. Improved knowledge of the influences of environmental conditions on stock dynamics would open the way to pre-recruitment prediction, at least in broad terms, of likely stock size in the coming season, thus improving the sustainable management of squid fishery resources.