Environmental influences on commercial oceanic ommastrephid squids : a stock assessment perspective

1 College of Marine Sciences, Shanghai Ocean University, Shanghai 201306, China. (JW) E-mail: jtwang@shou.edu.cn. ORCID-iD: http://orcid.org/0000-0003-3332-1209 (XC) (Corresponding author) E-mail: xjchen@shou.edu.cn. ORCID-iD: http://orcid.org/0000-0003-2529-0852 2 School of Marine Sciences, University of Maine, Orono, ME 04469, USA. (KT) E-mail: kisei.tanaka@maine.edu. ORCID-iD: http://orcid.org/0000-0002-1901-6972 (JC) E-mail: jie.cao@maine.edu. ORCID-iD: http://orcid.org/0000-0002-9291-9182 (YC) E-mail: ychen@maine.edu. ORCID-iD: http://orcid.org/0000-0003-0655-3864 3 Collaborative Innovation Centre for Distant-water Fisheries, Shanghai 201306, China. 4 National Engineering Research Centre for Oceanic Fisheries, Shanghai Ocean University, Shanghai 201306, China. 5 Key Laboratory of Sustainable Exploitation of Oceanic Fisheries Resources, Ministry of Education, Shanghai Ocean University, Shanghai 01306, China.


INTRODUCTION
Cephalopods are considered to be one of the three major marine fishery resources that are not fully exploited (Chen 1996).Global annual catches of cephalopods have risensteadily over the last fewdecades, with the highest catch being 4.3 million tonnes in 2007 and 4.77 million tonnesin 2014.The annual catches of cephalopods accounted for 3.04-4.75% of the global capture fishery yields (animals) from 1991 to 2012 (FAO 2016).Commercially targeted cephalopods include species of the families Sepiidae, Loliginidae, Ommastrephidae and Octopodidae, with Ommastrephidaehaving the largest proportion (Chen et al. 2012).The Ommastrephidae squid species such as Todarodes pacificus, Ommastrephes bartramii, Nototodarus sloani, Illex argentinus and Dosidicus gigas are commercially important species.O. bartramii, I. argentinus and D. gigas are traditionally targeted by the Chinese squid-jigging fleets, and their annual production reached 400000 to 700000 tonnes in 2013 and 2014.
The neon flying squid O. bartramii is widely distributed in the northwest Pacific Ocean (Murata 1990).The Japanese squid-jigging fleet has commercially exploited this species since 1974, and it was later exploited by South Korea, China and Chinese Taipei (Wang and Chen 2005).The O. bartramii population is composed of four stocks: the central stock of the autumn cohort, the eastern stock of the autumn cohort, the western stock of the winter-spring cohort and the central-eastern stock of the winter-spring cohort (Yatsu 1992, Yatsu et al. 1997).Of these stocks, the western winter-spring cohort has become an important target in the regions of 150-165°E (Chen et al. 2008b).Since 2000, this stock has mainly been captured by China and Chinese Taipei.The western winter-spring cohort migrates from subtropical waters to the subarctic boundary during the first half of the summer and then moves northward into the sub-Arctic domain in August-November.This species matures gradually in autumn and it seems thatit starts the spawning migration in October and November (Murata and Hayase 1993).
The spawning and feeding grounds for the western winter-spring cohort of O. bartramii are located at 20-30°N, 130-170°E and38-46°N, 150-165°E (Bower 1997, Bower andIchii 2005).The suitable sea surface temperature (SST) for the spawners ranges from 21 to 25°C during the spawning season (January-April) (Bower 1997, Bower andIchii 2005) and favourable SSTs for the feeding ground vary between August and November: 15 to 19°C in August, 14 to 18°C in September, 10 to 13°C in October and 12 to 15°C in November (Chen and Tian 2005).
Illex argentinus, the Argentinian shortfin squid, is a major component of commercial squid fisheries in the southwest Atlantic Ocean.I. argentinus is widely distributed across the Patagonian Shelf and adjacent oceanic waters between 22 and 54°S (Wang andChen 2005, Chen et al 2008b).The population structure of I. argentinus is complex and based on length at maturity, area and timing of spawning, and the distribution of early juveniles and adult life stages, divided into four stocks: the South Patagonia Stock (SPS), the Bonaerensis-Northpatagonia Stock, the Summer Spawning Stock and the Southern Brazil Stock (Hatanaka 1986).The Falkland Islands fishery consists almost exclusively of winter spawners and fishery managers assume that it target a single stock (SPS) south of 45°S (Brunetti and Ivanovic 1992).China (Lu and Chen 2013), Chinese Taipei (Chen and Chiu 2009) and the Falkland Islandshave exploited the SPS cohort of I. argentinus since 2000 (Waluda et al. 1999).
Previous studies from this region suggest that the inferred hatching area (the region where winterhatched I. argentinus paralarvae are most commonly found) of SPS is found between 32 and 39°S and 49 and 61°W from June to July (a peak hatching season).Waluda et al. (2001a) found that suitable SST for hatching ranges from 16 to 18°C.The species feeds in the waters from 42 to 54°S and 54 to 62°W from January to June, and the preferred SST for feeding is 7 to 14°C (Lu et al. 2013a).
Dosidicus gigas, the jumbo flying squid, supports a major fishery in the eastern Pacific Ocean and its range of distributed extends from 37-40°N to 45-47°S, between California and southern Chile, extending westward from the coast to a maximum of 125-140°W at the equator (Nigmatullin et al. 2001).D. gigas fisheries occur off the coast of Peru (Taipe et al. 2001), and Central America (Ichii et al. 2002) and in the Gulf of California, Mexico (Nevárez-Martínez et al. 2000).D. gigas's population structure is not well defined, but based on individual mantle lengths it may include three stocks: large, medium and small sizes (Argüelles et al. 2001).Small-scale fisheries, particularly off Peru, Chile and Mexico, traditionally exploited these stocks (Nevárez-Martínez et al. 2000, Rocha andVega 2003).However, since the early 1990s there has been an increase in industrial vessels, largely from East Asia (Japan, Korea and China), fishing off Peru and Central America, with a resulting increase in squid catches from the eastern Pacific Ocean (Ichii et al. 2002).The spawning ground of D. gigas is defined in the region from 5 to 14°N and 85 to 95°W and the favourable SST for spawning is from 24 to 28°C in September.The feeding ground of D. gigas is found from 20°N to 20°S and from 70 to 110°W, based on the distribution of fishing fleets, and the favourable SST is from 17 to 22°C in July (Ichii et al. 2002, Waluda andRodhouse 2006).
Recent studies suggest that variations in the oceanographic environmental variables, such as SST, sea surface height (SSH), and chlorophyll a (Chl a), play an important role in regulating the distribution and abundance of squid populations (Waluda et al. 2001b, 2006, Chen et al. 2007).Here we focus on SST influencing abundances of ommastrephid squid.For example, Chen et al. (2007) found that the SST rose by 1°C in a La Niña year compared with the normal level, and a higher SST resulted in unfavourable recruitment conditions for O. bartramii.In an El Niño year, however, the SST was almost the same as in a normal year; thereby leading to favourable conditions for recruitment.SST data analyses for I. argentinus showed that variations in the availability of optimum sea temperature for larval development during the spawning season explained about 55% of the recruitment variability (Waluda et al. 2001a).Moreover, D. gigas abundance was strongly influenced by the mesoscale variability linked to the ENSO, with low levels of upwelling during the very strong El Niño leading to low catches of D. gigas off Peru (Waluda et al. 2006).
The majority of ommastrephid squid, such as O. bartramii, I. argentinus and D. gigas, are short-lived species with about aone-year life span, and the biomass of these species in the subsequent year is largely composed of new recruits.As the annual recruitment variability of ommastrephid squid is largely driven by the environment, there is a growing need to incorporate environmental variability in the assessment of ommastrephid stocks.However, to date few studies of ommastrephid stock assessment have included environmental factors.
In this study, we modified a surplus production model by incorporating known environmental factors of spawning and feeding grounds, and then used this model to evaluate the environmental influence on the populations of O. bartramii, I. argentinus and D. gigas.This study shows the importance of including environmental variables in the assessment of squid populationdynamics.

Fishery data
Daily catch (tonnes, t), effort (days fished, d), fishing dates and fishing locations (longitude and latitude) were obtained from the Chinese commercial jigging fleets operating in the northwest Pacific Ocean, the southeast Pacific Ocean and the southwest Atlantic Ocean.Table 1 lists the spatial and temporal ranges ofommastrephid squid-jigging fishery.The spatial scale used for aggregating the fishery data is 0.5° latitude by 0.5° longitude.
There is virtually no by-catch in the squid-jigging fishery (Wang and Chen 2005).All the Chinese jigging vessels in a given fishing area areequipped with the same type of jiggers, and had identical fishing power and operational behaviour.Meanwhile, Chinese jigging vessels operated the same number of jigging machines; the fishing effort extracted from different vessels is consistent (Wang and Chen 2005).Therefore, we consider catch per unit effort (CPUE, t/d) of the squid-jigging vessels as a reliable indicator of stock abundance on the fishing grounds (Chen et al. 2008a).The formula for calculating monthly nominal CPUE in one fishing unit of 0.5×0.5°was where CPUE y,m,i is the monthly CPUE at i fishing unit in month m and year y; C y,m,i is the monthly catch (t) at i fishing unit in month m and year y; and E y,m,I is the number of fishing days at i fishing unit in month m and year y.
Annual O. bartramii catches by China accounted for approximately 80% of the total catches of this species (Wang and Chen 2005).Therefore, we extrapolated Chinese catch (divided by 80%) to estimate the total annual catch of O. bartramii from 2003 to 2013 in the northwest Pacific.The annual catch of I. argentinus was obtained from mainland China and Chinese Taipei, which landed almost all the catch from the SPS cohort in the Falkland Island Interim Conservation Zone, and from the Falkland Islands, which landed almost all the catch in the Exclusive Economic Zone.The Food and Agriculture Organization (FAO) provided annual catches of D. gigas (Fig. 1).

Environmental data
The favourable SSTs on the spawning and feeding grounds for three species are shown in Table 2. Monthly global SST, SSH, and Chl a concentration data from 2000 to 2013 were downloaded from the Live Access Server of the National Oceanic and Atmospheric Administration OceanWatch (http://Oceanwatch.pifsc.noaa.gov/las/servlets/data).The spatial resolution of SST, SSH and Chl a concentration data were 0.1×0.1°,0.25×0.25°and 0.05×0.05°,respectively.All the environmental data were then converted to a 0.5×0.5°grid for each month in order to correspond to the spatial grid of CPUE.

Standardizing yearly CPUE by the generalized additive model
Many environmental factors such as SST, SSH and Chl a determine the distribution of squid populations (Tian et al. 2009a), and those factors are used to forecast the fishing grounds for squids.To remove possible effects of environmental factors in fishing processes, we standardized nominal CPUE.In a similar approach, a generalized additive model (GAM)was used to standardize yearly CPUE of ommastrephid squids by Wang et al. (2015) and discussed in Tian et al. (2009b)  where s is a spline smoother function; and constant c is fixed 10% of mean CPUE (Campbell 2004); Var ε=σ 2 and E(ε)=0.

Development of environmentally dependent surplus production models
Bakun and Csirke (1998) suggested that environmental factors tended to be more important than factors such as predators and diseasein influencing the dynamics of squid recruitment.Therefore, we assumed that suitable SSTs in the presumed spawning grounds during the spawning seasons (referred to as Ps) and presumed feeding grounds during the feeding seasons (referred to as Pf) have effects on parameters K and r in surplus production models, respectively.
The following four surplus production models were proposed: (3) where B t is the biomass in t year; K is the carrying capacity; r is the intrinsic rate of stock growth; u t is independent and identically distributed IID N(0, τ 2 ); and Ps t and Pf t , are annual environmental indices for K and r, which are the averages of monthly Ps and Pf in year t, calculated as the number of fishing spatial grids with optimal SST divided by the total number of the fishing grids on the spawning and feeding grounds.Equation ( 3) is the original Schaefer surplus production model (referred to as SP).Equation ( 4) is the surplus production model with added Ps effect (referred to as Ps-EDSP).Equation ( 5) is the surplus production model with added Pf effect (referred to as Pf-EDSP).Finally, Equation ( 6) is the surplus production model with the consideration of Ps and Pf effects (referred to as Ps-Pf-EDSP).The standardized CPUE, I t , is assumed to be proportional to B t .This assumption can be written as where q is the catchability coefficient and v t is independent and identically distributed IID N(0, τ 2 ).

Parameter estimation
To fit the models to observed data for estimating model parameters, we used the maximum likelihood (ML) method.A likelihood function was used to calculate the goodness of fit between the observed data and the predicted values.We assumed that the observational errors followed the log-normal distribution, and the likelihood function was written as

∏ (
) To estimate uncertainties for model parameters, we developed an algorithm that incorporates a bootstrap approach into the ML.The following steps describe this bootstrap-based ML method: -estimate the model parameters using the ML method; -calculate the predicted CPUEs using the MLestimated parameters; -calculate the residuals between observed and predicted log CPUEs; -randomly sample the residuals with replacement to add to the predicted logarithm CPUEs to yield pseudo-observed CPUEs; -apply the ML method to the pseudo-observed CPUEs to estimate bootstrapped estimates; -repeat steps 4 to 5 1500 times to simulate 1500 sets of pseudo-CPUE data and to subsequently estimate the corresponding 1500 sets of bootstrapped parameters; and -calculate the mean value and 95% confidence intervals for each parameter using the 1500 bootstrapped estimates.
The meanvalues and 95% confidence intervals derived from the 1500 sets of bootstrapped parameters were used as parameter estimates and their associated uncertainties.The choice of 1500 bootstrap runs was made because stable estimates of uncertainly could be achieved when it reached 1500 runs.The best model for each species was selected based on the Akaike information criterion (AIC) of models and the CVs of the estimated parameters.

Biological reference points
Biological reference pointswere estimated, including maximum sustainable yield (MSY), fishing mortality at MSY level (F MSY ), and biomass at MSY level (B MSY ).F 0.1 and B MSY were used as the target reference points, and F MSY and 25% of B MSY were used as limit reference points (i.e.F lim and B lim ) in our proposed management system (Wang et al. 2015).

GAM-standardized CPUEs for ommastrephid squids
All the covariates (year, month, longitude, latitude, SST, SSH and Chl a) were statistically significant in a preliminary analysis, with no interaction terms being considered.The comparison of GAM-standardized CPUEs and nominal CPUEs showed temporal fluctuation in the species abundance (Fig. 2).For O. bartramii, the major difference was observed in 2007, with the highest nominal CPUE and lowest GAM-standardized CPUEs.For I. argentinus, the trend of GAM-standardized CPUEs was more stable than the trend of nominal CPUEs.The lowest and highest GAM-standardized CPUEs were observed in 2006 and 2009, respectively.A similar trend between nominal CPUEs and GAMstandardized CPUEs was observed for D. gigas.The highest CPUE was observed in 2004 and the lowest CPUE was observed in 2007.
For D. gigas, the smallest MSE, CV and AIC were 0.82, 11% and 33.17, respectively, for the Pf-EDSP model and the Pf-EDSP was identified as the best model for D. gigas (Table 6).Means r, K and q estimated with the bootstrapped ML method were 1.72, 270×10 4 t and 0.027×10 -4 , respectively.The estimated MSY ranged from 100×10 4 to 340×10 4 , with a mean value of 119×10 4 t (Fig. 5).The management reference points for ommastrephid squids using the original surplus production model and its optimal EDSP model were plotted in Figures 6, 7 and 8.For O. bartramii, the peak annual catch, (14×10 4 t in 2007) was near MSY, and all other annual catches were less than MSY.The fishing mortality rates from 2003 to 2013 were less than F 0.1 (=0.22) except in 2007 and 2008.The biomasses from 2003 to 2013 were higher than B MSY (Fig. 6).
For I. argentinus, the peak annual catch (54×10 4 t in 2007) was more than MSY, and other annual catches were all less than MSY.The fishing mortality rates from 2000 to 2010 were less than F 0.1 (=0.17) except in 2007 for the Ps-EDSP model.The biomasses during these years was higher than B MSY (Fig. 7).
For D. gigas, the peak annual catch (92×10 4 t in 2012) was less than MSY.The fishing mortality rates from 2003 to 2012 were less than F 0.1 (=0.77), except in 2012 for the Pf-EDSP model.The biomasses during the modelled period were higher than B MSY (Fig. 8).Overall, the assessment results for each species based on the original and its optimal surplus production models show that the fisheries of ommastrephid squid (O.bartramii, I. argentinus and D. gigas) are currently in good condition as their biomasses are all above their corresponding B MSY values.

DISCUSSION
Many studies indicated that variability in marine ecosystems, especially SST, causes inter-annual variability in the ommastrephid squid stocks (Waluda et al. 1999, 2006, Chen et al. 2007), but the inter-annual variability cannot be controlled and predicted effectively.Thus, assessment of ommastrephid squid stocks by incorporating environmental factors is difficult to implement.To our knowledge, there is only one study about assessing the western winter-spring cohort of neon flying squid with Bayesian inferences using SST on spawning ground and feeding ground based on surplus production model, and the model with spawning ground SST was the optimal model (Wang et al. 2015).
The longevity of ommastrephid squids is roughly one year, so the stock of a given year is comprised almost entirely of new recruits (Rosenberg et al. 1990, Pierce andGuerra 1994).It is therefore essential to allow the suitable spawning habitat (for spawners to spawn and the larvae to survive) to ensure a high probability of successful recruitment in the following year.Furthermore, ommastrephid squids are fast-growing   species (Boyle and Rodhouse 2005) and a suitable feeding habitat is vital for their growth in stock biomass.Thus, how to measure the habitat of the spawning and feeding grounds when assessing ommastrephid stocks is very important.This study incorporates the environmental indices (Ps and Pf) to reflect inter-annual variability in environmental conditions of spawning grounds and feeding grounds.On spawning grounds, climate changes and some large-scale or mesoscale oceanographic processes may drive the inter-annual variability of Ps, which can influence the recruitment of ommastrephid squids.For example, ENSO events and the strength of the North Equatorial Current and the Kuroshio Current change Ps to influence O. bartramii recruitment in the northwest Pacific (Chen et al. 2007).The relative strength of the Brazil current and the Falkland current may influence the Ps to affect the recruitment of I. argentinus in the southwest Atlantic (Waluda et al. 1999).ENSO events also strongly influence D. gigas recruitment, as weak upwelling during strong El Niño years lead to very low catches in the southeast Pacific (Waluda and Rodhouse 2006).On feeding grounds, ommastrephid squids are mainly distributed in waters with a suitable SST (Chen et al. 2007), which are mainly located in transition areas, including meandering eddies, frontal zones and so on, yielding highly productive habitats that serve as feeding grounds for ommastrephid squids (Zainuddin et al. 2006).The best surplus production models for the three-ommastrephid stocks (the surplus production model with Ps for O. bartramii and I. argentinus, the surplus production model with Pf for D. gigas) are all models incorporating environmental factors.This is consistent with our hypothesis that considering environmental variability can improve the assessment of ommastrephid squids.
For each modelled ommastrephid stock, the original surplus production model and corresponding best model showed differences in the exploitation status.The O. bartramii fisheries arenot overfished, which was consistent with previous results (Chen et al. 2008a, Wang et al. 2015).Parameters r, K and q based on SP and Ps-EDSP showed little difference, but K in SP model was about 10% less than that in Ps-EDSP model (Table 4).The status for each year based on Ps-EDSP model was closer to the limit (F/F MSY ) than any other based on the SP model (Fig. 6).This finding was similar to the results obtained by different methods modelled by Bayesian inferences (Wang et al. 2015).
For the I. argentinus fishery, status for almost the whole year in the target zone was in a good state.The MSY was about 44×10 4 t and the results were consistent with Lu et al. (2013b).However, the status of all years based on the SP model was much tight and no differences were found between them relative to the status of all years based on the Ps-EDSP model (Fig. 7).Parameter r in the SP model was about four times that in the Ps-EDSP model, and K in the SP model was about 12% less than that in Ps-EDSP model (Table 5).
For the D. gigas fishery, the development status in the SP and the Pf-EDSP model was similar (Fig. 8), but parameters r and K were about 19 and 26% less than any other one in the Pf-EDSP model, respectively (Table 6).
The improvements in the results through incorporation of environmental factors benefit the assessment and management of ommastrephid squid stocks.For example, one benefit is another method of assessing ommastrephid squid stocks.Because assessment of fished cephalopod stocks is a relatively undeveloped field, only a few methods are currently used, such as fishing and acoustic surveys, egg and larval production estimates, mark-recapture, the Leslie-Delury method and the surplus production model (Boyle andRodhouse 2005, Wang et al. 2015); even fewer of those methods incorporate environmental variables.The other benefit of incorporating environmental factors into management procedures could be the potential increase in fishing effort when future recruitment is expected to be stable or increasing, and the reduction in fishing effort when future recruitment is expected to be poor (Basson 1999).With the effective forecasting of environmental conditions, the surplus production models developed will allow precautionary approaches for the management of ommastrephid squid fisheries.
The dynamic life history of ommastrephid squids poses challenges for the assessment of stocks.This study shows that considering environmental factors can improve their assessment.Because of the importance of thermal habitat in regulating the dynamics of ommastrephid squids, we suggest that environmental factors such as SST should be considered in the assessment of ommastrephid squid stocks.

Fig. 1 .
Fig. 1. -Annual catch for O. bartramii in the northwest Pacific, I. argentinus in Southwest Atlantic and D. gigas in the southeast Pacific.

Fig. 2 .
Fig. 2. -Nominal and GAM-standardized CPUEs for O. bartramii in the northwest Pacific in 2003-2013, I. argentinus in southwest Atlantic in 2000-2010 and D. gigas in the southeast Pacific in 2003-2012 based on Chinese commercial jigging fleets.

Fig. 3 .
Fig. 3. -Distribution of the parameters estimated with the bootstrapped ML in the Ps-EDSP model for O. bartramii from 2003 to 2013.The red lines are the mean value for each parameter and the dotted red lines are the 95% confidence intervals.

Fig. 4 .
Fig. 4. -Distribution of the parameters estimated with the bootstrapped ML in the Ps-EDSP model for the I. argentnius from 2003 to 2013.The red lines are the mean value for each parameter and the dotted red lines are the 95% confidence intervals.

Fig. 5 .
Fig. 5. -Distribution of the parameters estimated with the bootstrapped ML in the Pf-EDSP model for the D. gigas from 2003 to 2013.The red lines are the mean value for each parameter and the dotted red lines are the 95% confidence intervals.

Fig. 6 .
Fig. 6. -Development of the O. bartramii fishery from 2003 to 2013.A, based on the PS model; B, based on the Ps-EDSP model.Triangle represents the biomass status of the first year (2003) of study; square represents the biomass status of the last year (2013).

Fig. 7 .
Fig. 7. -Development of the I. argentinus fishery from 2000 to 2010.A, based on the PS model; B, based on the Ps-EDSP model.Triangle represents the biomass status of the first year (2000) of study; square represents the biomass status of the last year (2010).

Fig. 8 .
Fig. 8. -Development of the D. gigas fishery from 2003 to 2012.A, based on the PS model; B, based on the Pf-EDSP model.Triangle represents the biomass status of the first year (2003) of study; square represents the biomass status of the last year (2012).

Table 1
. -The spatial and temporal range of Chinese squid-jigging fishery for O. bartramii in the Northwest Pacific, I. argentinus in the Southwest Atlantic and D. gigas in the southeast Pacific.Stock operating time ranges of fishing grounds O. bartramii from June to November in 2003-2013 35-45°N, 140-165°E I. argentines from January to June in 2000-2010 30-45°S, 40-65°W D. gigas from January to December in 2003-2012 8-18°S, 70-95°W

Table 2 .
-The favourable SSTs on the spawning and feeding grounds for O. bartramii in the northwest Pacific, I. argentinus in the southwest Atlantic and D. gigas in the southeast Pacific.

Table 3 .
-Time series indices of Ps and Pf for O. bartramii, I. argentinus and D. gigas.

Table 4 .
-Summary of the estimates of parameters with the bootstrapped ML method from 1500 runs of bootstrap simulation for CPUE and catch data observed from 2003 to 2013 for O. bartramii.For each bootstrap run, mean square error (MSE) was calculated.

Table 5 .
-Summary of the estimates of parameters with the bootstrapped ML method from 1500 runs of bootstrap simulation for CPUE and catch data observed from 2000 to 2010 for I. argentinus.For each bootstrap run, mean square error (MSE) was calculated.

Table 6 .
-Summary of the estimates of parameters with the bootstrapped ML method from 1500 runs of bootstrap simulation for CPUE and catch data observed from 2003 to 2012 for D. gigas.For each bootstrap run, mean square error (MSE) was calculated.