Assessing Northwest Pacific Fishery Stocks Using Two New Methods: The Monte Carlo Catch-MSY (CMSY) and the Bayesian Schaefer Model (BSM)

The CMSY and Bayesian Schaefer model (BSM) methods were applied to assess data-limited fishery stocks in the Japan Sea and surrounding areas of the Northwest Pacific. Ten stocks including 4 fish species and 5 cephalopod species were assessed; the CMSY method was used in 3 stocks with catch data only, and the BSM method in 7 stocks with both catch time series and catch per unit effort (CPUE) data available. The two methods estimated the maximum intrinsic rate of population increase (r) and carrying capacity of each stock, which allowed the computation of maximum sustainable yield (MSY), and exploited biomass relative to the biomass at maximum sustainable yield (B/BMSY). All 10 stocks were overfished, if to a different extent, and one, the spear squid (Heterololigo bleekeri) has collapsed. The reference points estimated here may be used as indicator for fishery management in this ecoregion.


INTRODUCTION
Global marine fishery catches fluctuated from 75 to 85 million tonnes since the late 1980s (FAO, 2018). From 1950 to 2012, both nominal fishing effort and capacity have increased from 1950 to 2012, especially in Asia and developing countries (Bell et al., 2016). Meanwhile, fishing fleets increased from 1.7 to 3.7 million vessels between 1950 to 2015, and effective catch per unit of effort (CPUE) has decreased substantially since 1950 (Rousseau et al., 2019).
Major stock collapses due to overfishing began to occur in the 1970s, and this became worse in the 1980s and 1990s (Pauly, 2008). Due to excessive fishing effort, the percentage of fish stocks being overfished within a decade after a fishery was fully developed increased from 26% in the 1950s to 35% in the 1980s. Nearly 50% of 900 important exploited species were overfished, collapsed or abandoned due to overfishing in 1999 (Froese and Kesner-Reyes, 2002), and the number of sustainably exploited stocks decreased from 90% in 197490% in to 67% in 201590% in (FAO, 2018. Large predatory fishes are strongly depleted, with their current biomass at only 10% of their pre-industrial levels (Myers and Worm, 2003).
While global fish fisheries declined, the populations of commercially exploited cephalopods have increased due to reduced predation by fish, thus allowing an increased fishing pressure to generate higher catches (Caddy and Rodhouse, 1998;Arkhipkin et al., 2015). Cephalopod fisheries increased from 1 million metric tons (mt) in the 1970s to over 4.3 mt in 2007 and fluctuating thereafter (Jereb et al., 2010;Arkhipkin et al., 2015;FAO, 2018). In the Japan Sea (which is surrounded by the Korean Peninsula, the Japanese islands and the Russian coast, forming a semi-enclosed marginal sea), small pelagic fish and cephalopods jointly contributed 52% of the total catch (Zhang et al., 2004). The growing interest in cephalopod fisheries is mainly due to the high value of their catch, and, for fisheries scientists, to their response to ecosystem changes (Rodhouse, 2001).
The lack of scientific knowledge about the development, fluctuation and status of the exploited stocks has often prevented the implementation of fisheries management plans (Ludwig et al., 1993). To help overcome this state of affairs, the CMSY and Bayesian state-space Schaefer surplus production model (BSM) methods were developed to allow the assessment of data-limited fishery stocks (Froese et al., 2017). These methods have been applied to 397 stocks in 14 European ecoregions by Froese et al. (2018).
The CMSY method relies mainly on catch time series and some ancillary information ("priors"), while the BSM method relies on catch time series and (relative) abundance data, such as catch/effort (CPUE) data. The CMSY and BSM methods generate estimates of the intrinsic growth rate of a population (r) along with an estimate of its carrying capacity (k); from these, time series of biomass (B) and fishing mortality (F) can be computed, including the biomass (B MSY ) from which maximum sustainable yield (MSY) can be extracted given F MSY . These results can help understand the status of exploited fishery stocks and be used to design fishery management plans.
The Northwest Pacific is a very productive fishing area, from which 22% of the global catch was taken in 2017; however, only a small fraction of the stocks therein are exploited optimally (FAO, 2018(FAO, , 2019.
In this contribution, we assessed 10 data-limited stocks in the Japan Sea and surrounding areas of the Northwest Pacific using the CMSY and BSM methods, to provide reference points for the fishery management in this ecoregion. We also compared the results obtained by the CMSY and BSM in 7 stocks with both catch and CPUE time series.

MATERIALS AND METHODS
The CMSY method was applied to assess the status of fishery stocks in the Northwest Pacific (Froese et al., 2017). In addition, the Bayesian state-space Schaefer surplus production model (BSM) Millar and Meyer, 1999) that is part of CMSY R-code was applied to account for variability in both population dynamics (process error) and measurement and sampling (observation error) (Thorson et al., 2014;Froese et al., 2017). The probability distributions of the parameters were sampled by the JAGS software with the Markov chain Monte Carlo method (Plummer, 2003). Three sampling chains were included in basic parameter settings, with a chain length of 60,000 steps each and a burn-in phase of 30,000 steps; every 10th value was used to reduce autocorrelation (Froese et al., 2017). All estimated posterior parameters were assumed to be approximately log-normally distributed; the median was used as central value with 95% confidence intervals approximated by 2.5th and 97.5th percentiles to find values at which test statistics attain less than 0.05 significance (Gelman et al., 1995;McAllister et al., 2001;Owen, 2013;Froese et al., 2017). All data files and R-code of the method were available in Supplement of Froese et al. (2017).
Catch and abundance data used in this contribution were derived from the published literature on 10 stocks, including 4 fish species and 5 cephalopod species ( Table 1). All 10 stocks were located in the Japan Sea, or in neighboring areas of the Northwest Pacific Ocean (Figure 1). Catch time series were the only input for 1 fish stock and 2 cephalopod stocks, which were assessed with the CMSY method ( Table 1). The other 7 stocks, including 3 fish and 4 cephalopod stocks with catch and CPUE time series were assessed with the BSM method ( Table 1).
The CMSY and BSM methods estimate parameters, including MSY, B MSY and F MSY based on the most probable r-k pairs filtered by a Monte Carlo test (Froese et al., 2017). The viable r-k pairs found by the CMSY and BSM methods generate triangularshaped clouds in the plot's log space, in which the most probable r-k pair (and its approximate 95% confidence limits) is found in the tip of the triangle (Figure 2).
A Bayesian state-space implementation of the BSM is then developed to verify the r and k according to catch and abundance data (i.e., biomass and CPUE) in Eq. 1: where B t is the biomass in year t, B t+1 is the exploited biomass in year t + 1, r is the intrinsic rate of population increase, k is the carrying capacity (i.e., the mean unexploited stock size), and C t is the catch data in year t. When a stock is severely depleted with its biomass below 0.25k, Eq. 1 is modified to account for depensation or reduced recruitment in Eq. 2: (2) where the term 4rB t /k expresses the assumption the intrinsic population growth rate declines linearly with biomass below half the biomass associated with MSY.
Resilience information in FishBase 1 and SeaLifeBase 2 were translated into the prior of r-ranges ( Table 2).
The maximum catch (C) was divided by the upper and lower bounds of r; these values were then used to inform the lower and upper bounds of k: Eq. 3 accounts for the stock with lower prior biomass at the end of the given time series, while as for the high biomass, it was modified: k low = 2 max (C) /r high ; k high = 12 max (C) /r low (4) All available r values were assigned to bins in log-space with equal width, of which the most probable was derived from the 75th  Demersal 1975-20061975-2006Tian et al., 2011 Todarodes pacificus (Japanese flying squid) Japan Sea Bentho-pelagic 19791995Fang and Chen, 2018Northwest Pacific 19791990Fang and Chen, 2018 Uroteuthis edulis (Swordtip squid) Japan Sea Demersal 1975-2006-Tian et al., 2011 FIGURE 1 | The areas of the 10 fishery stocks (see Table 1  percentile of the mid-values of occupied bins. If the r value was higher than the 50th percentile of mid-values of occupied bins, the most probable k value is derived from a linear regression: The range and density of the viable r values fitted inversely in Eq. 6 and a uniform distribution between 0.001 irf and 0.02 irf was used to describe the standard deviation of r in log-space: where irf is an inverse range factor to determine the prior density of r, r high and r low are defined as above. The k values are assumed to have a log-normal distribution, the standard deviation of which was assumed to be a quarter of the distance between the central value and the lower bound of the k range (McAllister et al., 2001).
The abundance estimation is attainable for data-limited stocks by a catchability coefficient q in Eq. 7, by which the Schaefer model used to transform the CPUE into biomass: where CPUE t and B t are mean catch per unit effort and biomass in year t, respectively, and q is the catchability coefficient. The dynamic of abundance as CPUE is expressed by Eq. 8: where the variables and parameters are defined as in Eqs. 1 and 7, and the prior q is derived from Eq. 9: where Y is the equilibrium yield for B, and other parameters are defined as in Eq. 1. The lower and higher priors for q are derived from Eqs 10 and 11 for stocks with recent high biomass: FIGURE 2 | Viable r-k pairs of spear squid (Heterololigo bleekeri) (Japan Sea) (A) and Japanese sandfish (Arctoscopus japonicus) (Japan Sea) (B) obtained from the CMSY (gray) and the BSM (black) methods. The most reliable r-k pair and its approximate 95% confidence limits are indicated by a black cross for the CMSY method (A,B); for the BSM method, the corresponding cross is gray (B).
where q low and q high are the lower and upper prior catchability coefficient for stocks with high recent biomass, respectively, r pgm is the geometric mean of r, r high is the upper prior range for r, CPUE mean is the mean CPUE over the last 5 or 10 years, and C mean is mean catch over the same period. For stocks with recent low biomass, the multipliers were changed to 0.5 for q low and to 1.0 for q high . Mean catch and CPUE were applied to species with medium and high resilience over the past 5 years and to species with low or very low resilience over the past 10 years. B start /k and B end /k were the priors of relative biomass at the start and the end of each time series, and their ranges were estimated depending on the assumed depletion level (Froese et al., 2017;FAO, 2019). B start /k of the Japanese amberjack with available fishery data starting from the 1950s was set as 0.9 to 1.0 (Tian et al., 2012), and the other 9 stocks were set as 0.5-0.9 due to the increased depletion in recent decades ( Table 3; Pauly et al., 2002;Tian et al., 2011;Watson et al., 2013Watson et al., , 2014Arkhipkin et al., 2015;FAO, 2019). The B end /k of spear squid was set as 0.01 to 0.2 due to high depletion (Tian, 2009;Arkhipkin et al., 2015), with the other 9 stocks set as 0.01 to 0.4, according to the increasing fishing pressure and efforts in recent years (Tian et al., 2011;Watson et al., 2013Watson et al., , 2014Arkhipkin et al., 2015;FAO, 2019).
The relative total biomass in the last year (B/B MSY ) was used to assess the exploitation level and the stock status estimated by the CMSY and BSM methods ( Table 4; Palomares et al., 2018).

RESULTS
The B/B msy values in the last year for all 10 stocks assessed were less than 1, indicating the occurrence of overfishing ( Table 5). All stocks were depleted, but the trends were different (Figures 3, 4).
All 4 fish stocks were subject to ongoing overfishing, as shown in Figures 3A-D. The Japanese sandfish and blackthroat seaperch had similar trends in estimated biomass, with a substantial decrease since 1975 (Figures 3A,B). The relative biomass (B/B MSY ) were 0.54 and 0.51 for them at the end of the time series, respectively ( Table 5). The Okhotsk atka mackerel has been subject to ongoing overfishing since the early 1980s with B/B MSY reaching 0.76 at the end of time series (Figure 3C and Table 5). The Japanese amberjack was in a healthy condition and had a biomass above the one that can produce MSY before the 2000s, but was overfished after that ( Figure 3D and Table 5).
Of the 6 cephalopod stocks, 3 were overfished, 2 were endangered by reduced recruitment (0.2 < B/B MSY < 0.5) and 1 has collapsed (B/B MSY < 0.2; Table 5). The spear squid was severely depleted and subject to unsustainable exploitation (B/B MSY = 0.12; Table 5). The catches and exploitation levels were different for golden cuttlefish Sepia esculenta stock and swordtip squid stock from the Japan Sea; however, the estimated biomasses for both stocks have declined sharply to below B MSY since the 1980s (Figures 4B,E). The value of B/B MSY of neon flying squid was below 0.5 in 2010, which lasted for several years (Figure 4A). Two Japanese flying squid stocks have been overfished since the late 1980s and late 1990s, respectively (Figures 4C,D).
The CMSY results could be reproduced using the BSM method (Table 6). Overall, the 5th-95th percentile ranges were narrower for r, k, MSY, and B/B MSY from BSM than that from CMSY. Thus, only the estimates r of okhotsk atka mackerel (P. azonus) and MSY of golden cuttlefish (S. esculenta) had wider ranges in BSM. The 95% confidence limits of the CMSY estimates for r, k and MSY of all 7 stocks included the most probable BSM estimate.

DISCUSSION
Many stock assessment models require plenty of data making their implementation generally limited to valuable species or very abundant species (Harley et al., 2011(Harley et al., , 2014; however, less attention is paid to other species (Pauly, 2006;Costello et al., 2012;Costello and Ovando, 2019). However, science-based catch or effort limits are vital information for effective fishery management (Melnychuk et al., 2017). MSY is often used as a reference point for fishery stock assessment, and the status of the fishery is commonly reported in terms of B/B MSY (Schaefer, 1954;Costello and Ovando, 2019). Understanding of the nature of stock dynamics is of great significance for the prevention of stock depletion and for rebuilding depleted stocks. In this study, two recent data-limited stock assessment methods were applied to estimate the exploitation status of 10 fishery stocks in the Northwest Pacific. A key advantage of using these two methods was that they work in data limited situations (catch time series for CMSY, catch time series CPUE for BSM), yet produce results that can support policy and management decisions at national and regional levels. Also note that before applying BSM method, we standardized the CPUE data by accounting for "technological creep" (Palomares and Pauly, 2019), i.e., the gradual increase of effective over nominal effort.

Demersal Fishes
Fisheries catches in the Japan Sea reached its peak in the 1980s, followed by an abrupt decrease due to the collapse of the Pacific sardine Sardinops sagax, and the decline of demersal fish and invertebrate fisheries resulting from overfishing (Tian et al., 2006). The variations in demersal fish communities were strongly associated with the oceanographic conditions; the Tsushima Warm Current experienced changes from a colder to a warmer regime in the late 1980s, which had a strong effect on the ecosystems of the region and the fish communities therein (Tian et al., 2006(Tian et al., , 2011. Japanese sandfish, a cold-water species, was one of the most important commercial resources in Japan (Watanabe et al., 2011: Yoon et al., 2018. In the northern and western Sea of Japan, the Japanese sandfish stock decreased sharply after the mid-1970s and was severely depleted during the 1980s (Watanabe et al., 2011). This confirms the results of this study, which showed that Japanese sandfish stock decreased rapidly from 1974 on. Blackthroat seaperch, a high-value species exhibited a similar dynamic, and also ended up being overfished in the 1980s. The catches of Japanese amberjack showed an increasing trend from the 1950s, mainly due to an increase of effort by the purse seine fishery, which resulted in the this stock being overfished in the Japan Sea (Tian et al., 2012).

Cephalopods
As finfish landings have declined, invertebrate fisheries, especially for cephalopods, have grown and become more important in global fisheries (Anderson et al., 2011: Rodhouse et al., 2014. Therefore, more attention is being paid to invertebrate fisheries, which also need careful management to avoid the same fate as many finfish fisheries (Anderson et al., 2011;Doubleday et al., 2016). It was hypothesized that the growth of cephalopod catches might be related to the release of predation and competition pressure as a result of fish stock depletions (Caddy and Rodhouse, 1998). Cephalopod populations increased from the early 1950s to mid-2010 but declined after that (Doubleday et al., 2016;FAO, 2018). Cephalopods, particularly squids, are generally shortlived. Their abundance and distribution are greatly influenced by oceanographic conditions (Sakurai et al., 2000;Chen et al., 2006). For instance, the regime shift of the Tsushima Warm Current in the late 1980s had a strong effect on the Japan Sea, which may have affected the abundance and distribution of squids (Tian et al., 2006Tian, 2007). Since the late 1980s, the effort of Japanese fleets targeting squid has decreased continuously (Arkhipkin et al., 2015). The swordtip squid stock was not overfished in the southwest of the Japan Sea in the 1980s, but a recent stock assessment indicated that it remained at a low level and the total allowable catch (TAC) should be reduced (Yoda and Fukuwaka, 2013;Arkhipkin et al., 2015).
The catch of spear squid from the southwestern Sea of Japan was generated mainly by pair trawlers, and it decreased from a peak of 13,700 t in 1977 to 16 t in 2003, and remained low level thereafter (Tian, 2009;Arkhipkin et al., 2015). The abundance of southern stocks decreased abruptly around 1990 both in the Japan Sea and the wider North Pacific (Arkhipkin et al., 2015), which was closely associated with the rising water temperature (Tian et al., 2013). The spear squid fishery collapsed in the late 2000s in the Japan Sea (Tian, 2007(Tian, , 2009, and also declined in the wider North Pacific (Bower and Ichii, 2005).
Neon flying squid is widely distributed from subtropical to cold temperate waters, but is only commercially fished in the Pacific Ocean (Arkhipkin et al., 2015). The fishing mortality of the winter/spring cohort exerted by the jigging fishery appeared to be sustainable, but the biomass decreased from 2001 to 2005 (Chen et al., 2008). The western winter-spring cohort in the Northwest Pacific was over-exploited in 2010 (Ding et al., 2019). Our results showed that the Neon flying squid stock was overfished since the 2000s, in agreement with the literature.
The autumn cohort and winter cohort of the Japanese flying squid, of great commercial importance in Japan, are commonly assessed separately. The former was mainly fished from June to December in the Japan Sea, while the latter is fished from July to December in the Pacific and from January to March in the Sea of Japan (Arkhipkin et al., 2015). This stock was managed since 1998, based on relatively high TACs issued by the Japanese Fisheries Agency (Arkhipkin et al., 2015). The fishery for Japanese flying squid in the Sea of Japan experienced a rapid decline in the late 1970s. However, it has partly recovered since the mid-1980s due to the implementation of management measures, such as the implementation of Total Allowable Catch since 1998 (Fang and Chen, 2018).

Comparison of Results Between CMSY and BSM
Of 7 stocks, the most probable BSM estimates of r, k, and MSY were within the 95% confidence limits of the CMSY estimates, suggesting a good agreement between two methods. Froese et al. (2017) applied BSM and CMSY to 28 data-limited stocks with catch and CPUE available and indicated that estimates were not significantly different for 25 of them. As expected, the 5th-95th percentile ranges of most CMSY estimates were wider than those provided by the BSM method. Estimates from available catch and CPUE data sets can thus be corroborated, and combining them can lead to narrower confidence intervals, indicating the use of CPUE estimate is an important step in reducing the uncertainty of r, k and MSY as estimated from a catch time series.

CONCLUSION
This contribution assessed 10 commercial fishery stocks in the Northwest Pacific, especially the Japan Sea, using the CMSY and BSM methods. The results indicated that all 10 stocks were overfished at different levels according to the values of B/B MSY , with one squid stock (the spear squid) having collapsed. More studies should be conducted to understand the fisheries in the Northwest Pacific region, which are essential for future fishery management.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article can be obtained in the Supplementary Material.

AUTHOR CONTRIBUTIONS
QR wrote the first draft. ML revised the manuscript. QR and ML performed the data analyses. Both authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by the National Programme on Global Change and Air-Sea Interaction (GASI-02-PAC-YD spr/sum/aut). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.