Assessments of 16 Exploited Fish Stocks in Chinese Waters Using the CMSY and BSM Methods

Sixteen marine fish species (populations) exploited by Chinese fisheries were assessed, using published time series of catch and the CMSY and BSM methods. Given the catch times series as inputs, some ancillary information and reasonable constraints, carrying capacity, maximum sustainable yield, and likely time series of biomass and exploitation rate were estimated. The results show that one (7%) of the assessed species was severely depleted, four species (27%) were fully/overfished, six (40%) were outside of safe biological limits, one species (7%) was recovering and three species (20%) were in a healthy state at the end year of their assessment. However, one species, Pacific sardine (Sardinops sagax), could not be assessed using CMSY, as the exceedingly large fluctuations of its biomass were mainly environmentally driven. These results correspond with previous knowledge on the status of fish populations along the coast of China, where overfishing is rampant. Based on these assessments, some of the benefits that would result from a reduction of the excessive fishing effort are outlined.


INTRODUCTION
China is the country with the world's largest marine fisheries catch (FAO, 2018). It is widely agreed that China's domestic fisheries resources are overexploited (Shen and Heino, 2014). However, it seems that overfishing has gradually changed the structure and function of marine ecosystems of China's coastal seas (Zhai and Pauly, 2020) and that the state of its domestic resources is the main reason for its current emphasis on distant-water fishing (Mallory, 2013).
China's fisheries management system has been gradually improving since the 1980s (Huang and He, 2019). Its most powerful regulations are the "double control" system and summer fishing moratoria (Shen and Heino, 2014). The former are regulations of both the total number of marine engine-powered fishing vessels and their total engine power; the latter have been implemented since 1995 and extended from 3 to 4 months and more in many areas (Ministry of Agriculture and Rural Affairs of the People's Republic of China [MARA], 2018). Also, some other regulations and programs are being implemented including vessel buyback, alternative employment opportunities for fishers (Song, 2007)  China has carried out a large number of fishery resource surveys and stock assessment work in recent years, but still faces difficulties due to a lack of historical data. Thus, there are still deficiencies regarding China's fish stock assessments: the range of species covered by assessments is narrow, evaluation methods are limited and the results are not usually expressed as B/B MSY , i.e., the ratio of stock biomass to the biomass that can produce the maximum sustainable yield (MSY) or other management-relevant indicators. Until recently, very few stock assessments were conducted that allowed for estimating the potential for stock rebuilding (Villasante et al., 2013), although this dire situation is now being overcome (see Zhai and Pauly, 2019;Ji et al., 2019;Liang et al., 2020a,b;Wang et al., 2020;Zhang et al., 2020).
Here, the newly developed CMSY method (Froese et al., 2017), catch time series and some Supplementary Information are used to infer likely biomass trajectories for 16 species of exploited fish along Chinese coastlines. This method was shown to be adequate for the assessment of hundreds of fish stocks in and around Europe (Froese et al., 2017(Froese et al., , 2018Demirel et al., 2020). Moreover, a number of contributions show that the method can be used for stock assessment in China (Ji et al., 2019;Liang et al., 2020b;Wang et al., 2020).
Also, a Bayesian state-space implementation of the Schaefer model (BSM; Froese et al., 2017) is applied here with the same catch time series and priors, and additional catch/effort (C/f, or CPUE) and/or biomass data, to assess the biomass and exploitation rate, and thus to assess the stock status and extent of overfishing in Chinese and adjacent waters. This contribution is a further example of the study of data-poor fisheries and our results provide information toward a science-based management of China's fisheries.

Data Sources
The main features of the catch time series and additional information available for 16 fish populations investigated here are given in Table 1, while the data sources are shown in Supplementary Table 1. All species investigated here are important commercial fish in China since the 1950s that have more than 10 years of catch data.

The CMSY and BSM Methods
The principle of the CMSY method is that, given catch time series and using a Monte Carlo approach, multiple biomass trajectories of the biomass of the stock in question are traced, and the parameters retained that generated the biomass trajectory (or trajectories) that is (are) compatible with the time series of catches, and a number of constraints (Froese et al., 2017). Here, "compatible" means that, among other things, the stock does not crash, i.e., its biomass does not drop to zero.
When, as is the case here, relative abundance data, such as CPUE or spawning stock biomass (SSB), are also available, a Bayesian state-space implementation of the Schaefer production model (BSM; Millar and Meyer, 1999) is applied to estimate the intrinsic rate of population growth (r) and unexploited stock size (or carrying capacity; k) of each stock (Froese et al., 2017).
The CMSY and BSM methods are based on the logic of the surplus-production model of Schaefer (1954Schaefer ( , 1957. Thus, they assume that from 1 year (t) to the next (t + 1), the biomass (B t ) follows the equation: where r is the intrinsic rate of population growth, k the carrying capacity, and C t is the catch in year t. However, Eq. (1) is slightly modified when the biomass falls below 0.25k, to allow for depensation or reduced recruitment when stock size is severely depleted (Froese et al., 2017): (2) where the term 4rB t /k ensures a linear decline of recruitment below half of the biomass capable of generating MSY.

Prior Information
The R code implementing the CMSY method incorporates a routine which estimates wide (uniform) priors for k (Froese et al., 2017), whose output were here accepted: where k low and k high are the lower and upper bounds of the prior range of k, max(C) is the maximum catch in the time series and r low and r high are the lower and upper bound of the range of r-values that the CMSY method explores. This is expressed by where the variables and parameters are defined as in Eq. (3). As stated in Froese et al. (2017), when running the BSM method, the estimated standard deviation of r in log-space is described by a uniform distribution between 0.001 irf and 0.02 irf by irf = 3/ r high − r low (5) where irf is an inverse range factor to determine r range, and r high and r low is provided in Table 2.
The k estimation by BSM also assumes that k has a log-normal distribution, while the mean of k provides the reasonable central value. The standard deviation of k is assumed to be a quarter of the distance between the central value and the lower bound of the k-range (McAllister et al., 2001).
Additionally, the BSM method allows the estimation of the catchability coefficient q, which relates biomass to CPUE, when the latter is available. For this, priors are defined as q low = 0.25r pgm CPUE mean /C mean ; q high = 0.5r high /CPUE mean (6) where q low and q high are the lower prior and higher prior for the catchability coefficient respectively; r pgm is the geometric mean of the prior range for r; CPUE mean is the mean CPUE over the last 5 or 10 years; C mean is mean catch over the same period.
For stocks with low recent prior biomass, the ranges of multipliers are 0.25-0.5 for q low and 0.5-1 for q high .  1989CPUE (1979-1997 Coilia mystus (Osbeck's grenadier anchovy) 2 Yangtze River and Estuary, and ECS 1960-2007CPUE (1995 Ilisha elongata ( The intrinsic rate of population growth (r) is largely determined by the size and age at first maturity of individuals of the species in question, and FishBase 1 (Froese and Pauly, 2019), based on these and other traits, provides ranges of likely values ("priors") for all fish species (Froese et al., 2000). Also, a resilience classification based on the stock fecundity was proposed by Musick (1999). The available resilience ranges and those selected for the 16 stocks studied here are given in Table 2.
Additionally, two constraints are required for the CMSY method, i.e., the faction of carrying capacity that the biomass was exhibiting just before the first annual catch in the available time series is subtracted (B start /k; see Eq. 1) and the faction of carrying capacity that the biomass reached at the end of the catch time series, again expressed as a fraction of k (B end /k). The first of these constraints is a function of the fishing pressure to which the population (or stock) in question was exposed to prior to the period for which catch data are available, and which would have reduced the biomass below carrying capacity.  Table 3 presents suggested ranges of B start /k and B end /k for the stocks in Table 1. Note that for some species that were already exploited in the 1950s, we selected a very low depletion (B start /k = 0.8-1) for the start of the series because of the low technology that was deployed at the time (Zhan et al., 1986;Yang, 1988;Zhang and Hua, 1990;Yu and Zheng, 2000;Liu, 2005;Xu and Liu, 2007;Wang et al., 2011;Xu et al., 2014;Yan et al., 2014). For some other species, we selected a low range (B start /k = 0.4-0.8), in view of the fact that from the 1970s to the 1990s, a massive increase of fishing effort occurred, and the technological sophistication of that effort also increased (Zhu, 1992;Qiu, 1995;Lu et al., 1998;Bo et al., 2005).
There are two exceptions here, Larimichthys polyactis and Sardinops sagax, were set at the strong (0.01-0.4) and very strong (0.01-0.2) depletion, respectively (Table 3). For L. polyactis, the catch of the stock in the final year increased to 2∼3 times that of the start year (Lin et al., 2004;Yan et al., 2014). Thus, we have a good reason to assume that L. polyactis was already strongly overfished in the mid-1950s.
The B start /k for S. sagax was set at a very strong depletion (0.01-0.2; see Table 3) to force its biomass at the start of the assessment period to be very low. This biomass subsequently increased due to environmental factors (Wei and Li, 1986;Lluch-Belda et al., 1989), a topic to which we will return.
The second constraint, the expected biomass at the end relative to carrying capacity (B end /k), should be roughly reflective of the fishing pressure which generated the available catch time series. Table 3 presents suggested ranges of B end /k, i.e., the fraction of the population initially present that is left in the water, set based on the default rules of CMSY method (Froese et al., 2017).
The last constraint, i.e., the ratio of biomass reached during an intermediate year to carrying capacity (B int /k), can be used if such knowledge is available. We set "default" B int /k for all stocks in this study except for Engraulis japonicus (Japanese anchovy; int. year is 2009, prior range is 0.2-0.6; Li et al., 2015). The CMSY method identifies the r and k values generating viable biomass trajectories (if any), then output the geometric means and confidence intervals of these estimates. As well, various parameters that can be derived from r and k, notably MSY = 0.25r·k, F MSY = 0.5r and B MSY = 0.5k (Schaefer, 1954;Ricker, 1975) are also estimated.
The BSM method can run with additional constraints, such as SSB or catch/effort (C/f, or CPUE) (Froese et al., 2017), which represent (relative) abundance. In this case, the viable biomass trajectories will take account of the (relative) abundance data, even if they pertain to a shorter period than the catch data ( Table 1). In the process, an estimate of the catchability coefficient q = (C/f )/B is also produced as an average for the period with CPUE data (see Eq. 6).
Four of the sixteen stocks in the study, i.e., E. japonicus, Coilia mystus, L. polyactis and Muraenesox cinereus with more than one SSB or CPUE time series, we used a routine "Bcrumb, " to interpolate and average SSB or CPUE for them. This routine was developed as a component of JARA (Just Another Redlist Assessment) described in Winker and Sherley (2019). Moreover, considering the technological improvement of the fishery, the CPUE is corrected by 2% increases every year according to Palomares and Pauly (2019), given gear and other technological improvements of the industrial sector.
The ratio B/B MSY in the final year, which is often used to express stock status, is estimated by both CMSY as well as BSM. For management purposes, the more precise results of BSM may be preferred over the results of a CMSY assessment. Table 4 presents suggestion of fish stock status based on B/B MSY and

RESULTS
Altogether 16 fish populations were analyzed by CMSY method, with 13 also assessed using the BSM methods, as CPUE/biomass data were available for them. The results of estimated r-k for Scomberomorus niphonius and E. japonicus are shown in Figure 1 as examples, while all other results are shown in Table 5 and Supplementary Figure 1(except for S. sagax).
In Figure 1A, featuring Japanese Spanish mackerel (S. niphonius) and Figure 1B, featuring Japanese anchovy (E. japonicus), the viable r-k pairs (gray) were obtained by the CMSY method, and the solid crosses identify the most probable r-k pairs, along with their 95% confidence intervals. In Figure 1B, the BSM method used (relative) abundance data (i.e., SSB/CPUE), and the viable r-k pairs it estimated are in black. The dotted cross identifies the most probable r-k pair, and its 95% confidence intervals. The overlaps between the two clouds of dots (and hence the closeness of the two crosses) implies that the results of the two methods are similar, and thus more credible than if they did not overlap.
The biomass trends resulting from the CMSY and BSM analyses are shown in the different panels of Figure 2. In each panel, the key item is the biomass trajectory (solid black line) and its 95% confidence intervals (dotted line). The horizontal dashed line represents B MSY , while the vertical solid lines show the priors on biomass range at the start and end of time series. Finally, the open dots indicate corrected CPUE trends of the 13 stocks for which such data was available.
As shown in Figure 2, the biomass of all species generally declined, although a few species, notably L. polyactis, substantially recovered. The mean B/B MSY for 15 stocks is 0.59 ± 0.10 (SE) and mean B end /k is 0.29 ± 0.05 (SE).
The results demonstrate that the exploitation rate (i.e., F) for 9 species of the 15 are higher than F MSY . The average value for 15 species is 1.39 ± 0.26 (SE). The species Clupea pallasi have the highest exploitation rate of 4.33 (Table 5 and Supplementary Figure 2).
According to the definition of fish stock status (see Table 4), the average status of China's coastal fish stock assessed here is fully/overfished. Three of the assessed stocks were healthy, one Frontiers in Marine Science | www.frontiersin.org FIGURE 1 | Examples of the identification of viable and best trajectories (as defined by values of r and k), for the populations of (A) Japanese Spanish mackerel (S. niphonius) and (B) Japanese anchovy (E. japonicus). In (A), the viable r-k pairs (gray) were obtained by the CMSY method, and the solid crosses identify the most probable r-k pairs, along with their 95% confidence intervals; in (B), the viable r-k pairs (black) were obtained by the BSM method, the dotted cross identifies the most probable r-k pair, and its 95% confidence intervals. was recovering, six were outside of safe biological limits, four were fully/overfished, one was severely depleted (see Figure 3).

DISCUSSION
Time series of catches such as analyzed here are scarce for Chinese waters; however, the few that are available provided the impression of massive overfishing, consistent with other analyses (Liang and Pauly, 2017b;Zhai and Pauly, 2019;Liang et al., 2020a,b). All of the thirteen stocks with CPUE or SSB included in this contribution were well fitted by the CMSY and BSM methods, i.e., the most likely estimated r-k pair are included within the 95% confidence intervals of each other (see Supplementary Figure 1). Similarly, Ji et al. (2019) found a good match between the estimates of population parameters of Trichiurus lepturus in the Yellow and Bohai Seas obtained by using CMSY and those from the Schaefer and Fox models. Thus, the CMSY method can for assessing data-limited fisheries, although the related BSM methods, which additional information (e.g., CPUE data) produces narrower confidence intervals . The intrinsic rate of population growth (r) is key in defining the response of populations to challenges such as fisheries (Cheung et al., 2005;Anderson et al., 2008). FishBase (see text  Table 1, as estimated by the CMSY method (solid black lines) and their uncertainty (gray dashed lines), related abundance series (gray dots; SSB and corrected CPUE by 2% every year), when available, and vertical lines represented priors (from Tables 3, 4). Except for S. sagax (see text), all of these trajectories are due to excessively high catches. footnote 1; Froese and Pauly, 2019) provides estimated ranges of r for fish populations based on life-history traits of fish populations which are very useful, as they compensate for the lack of estimates from local data.
The uncertainties of catch, SSB or CPUE data, which may result, e.g., from discards not having been considered, or unreported catch issues, will lead to errors in the assessment. In these cases, the local experience of fisheries experts should be accounted for when screening the data sets and results.
The catches of most of the species investigated here exhibited an explosive growth in the 1990s and/or the 2000s, then stagnated and declined. This feature basically coincides with China's fishery development history. Fisheries of China developed gradually from the 1950s to the 1970s, then rapidly in 1980s and 1990s, with its total catch peaking in the 2000s (Anon, 1963(Anon, -2019Shen and Heino, 2014 ; Supplementary Figure 3). Even though China's national catch statistics were over-reported in the 1990s, resulting in distortions in the catch and CPUE data (Watson and Pauly, 2001), the rapid growth of China's  (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017). Red area: stocks that are being overfished or are outside of safe biological limits; yellow area: recovering stocks; green area: stocks subject to sustainable fishing pressure and of a healthy stock biomass that can produce high yields close to MSY. overall fishing effort in undeniable, and it had a strong effect on the structure of its coastal ecosystems (Zhai and Pauly, 2019).
Some important fish species, such as Decapterus muroadsi, Branchiostegus japonicus, and Larimichthys crocea, were already strongly exploited in the 1950s. Thus, the stock of large yellow croaker (L. crocea), a prized species in China, was quickly depleted when China developed its industrial fisheries; it finally collapsed in 1974, once Chinese and South Korean fleet reached its overwintering grounds (Xu and Liu, 2007), rendering it commercially extinct (Zhao et al., 2002). Although protection measures for the spawning grounds of large yellow croaker were put in place by China in the early 1980s, the stock was severely depleted, and does not seem able to recover (Liu and Sadovy de Mitcheson, 2008).
Another important commercial species of China, yellow croaker (L. polyactis), showed first signs of overexploitation in the 1960s, but they were masked by the offshore expansion of the fisheries (Liang and Pauly, 2017a) and large catches of other species (Fei, 1980). The catch of yellow croaker initially decreased, but later appeared to increase to more than two times than that in the 1950s. This stock recovery was the result of China's and South Korea's protective measures, notably a 6-year total ban on fishing in its main spawning ground in the southern Yellow Sea in the 1980s (Lee and Midani, 2014;Yan et al., 2014). However, the recovery of yellow croaker came at a cost, as the adult individuals became much smaller, and reached maturity at smaller sizes (Dieckmann et al., 2005;Yan et al., 2014), as can be expected on theoretical grounds (Pauly, 1984(Pauly, , 2019. Japanese anchovy (E. japonicus), which had one of three populations found not to be overfished, is a very common and important small pelagic fish, especially as food for other commercial species. As a low-value species (Yu and Zheng, 2000), it was largely a bycatch species until the mid-1980s, when it became the target species of reduction (i.e., fishmeal) fisheries (Zhu et al., 1990;Tang et al., 2002) and thus began to show up in China fisheries statistics and fisheries research. While this research showed that Japanese anchovy was overfished around 2006 (Li et al., 2015), its stock seems to have recovered since, possibly because the top predators in China's coastal seas have been fished out, which reduced predation pressure on small pelagic fish.
Japanese flounder (Paralichthys olivaceus), another stock that is currently in good shape, was depleted in the early 1990s. Then, management measures, which include restocking (via the release of juveniles) were put in place and strengthened, which allowed the population to recover (Zhou, 2011).
Osbeck's grenadier anchovy (C. mystus), is the only stock found to be recovering. This anchovy is an amphidromous species that is very popular with consumers along the Yangtze River; in fact, it was one of the most important resource species of the Yangtze Estuary before the 1960s (Wang and Ni, 1984;Zheng, 2012). It became even more valuable in the last 30 years, as the populations of associated species, such as Reeves shad (Tenualosa reevesii) and the noodlefish (Salanx prognathus) strongly declined (Ni, 1999). The environment of C. mystus is very vulnerable to human activities and terrigenous pollution. Thus, we laud the decision to prohibit the harvesting of Osbeck's grenadier anchovy in the Yangtze River and its estuary, as published on February 1, 2019 (Anon, 2019), and the 10 years' fishing ban for the entire Yangtze River Basin that began in January 2020. These decisions, if vigorously implemented, may help the stock to recover, in spite of the pollution and degradation of its habitat (Yang et al., 2012). Note, finally that the stock of C. mystus in the Min River Estuary was found to be in good shape , thus suggesting a strong influence of local conditions. Pacific sardine (S. sagax), also known as "South American pilchard" (see text footnote 1; Froese and Pauly, 2019), ranges across the entire Pacific and beyond, and has supported many huge fisheries, notably off Chile and Peru, California (Steinbeck, 1945;Cisneros-Mata et al., 1995;Yáñez et al., 2001), and Japan (Watanabe et al., 1995). The stocks of Pacific sardine fluctuate strongly in response to climatic events, often in similar ways across very distant regions (Lluch-Belda et al., 1989, 1992Deyle et al., 2013), which led to high catches in China in the 1970s and 1980s (Wang, 1985). As the 2-parameter CMSY model cannot readily accommodate massive changes in biomass due to environmental fluctuations, Pacific sardine treated as if it was overfished in 1950-1975. This allowed its biomass to increase until the 1980s (as occurred in reality), then to decline due to strong fishing pressure and the fading of the environmental conditions that had generated the massive increase. However, we did not report on it "MSY" or other statistics, which would be unrealistic.

CONCLUSION
In conclusion, most stocks we studied have biomasses that are much lower than B MSY , i.e., the biomass associated with MSY. Thus, China loses millions of tonnes every year of potential catch to overfishing, and huge sums in the form of subsidies to fisheries that exploit overfished coastal stocks (Mallory, 2016). Although mariculture and the release of juveniles are maintaining several stocks, this is not a sustainable proposition; also, these measures lead to genetic diversity losses (Wang et al., 2012). Aquaculture should not be expected to maintain wild fish populations. Prudent fisheries management, on the other hand, can do this, and counter fisheries resources degradation. China has now taken measures to control its excessive fishing effort and we hope that the present study will serve as supporting scientific evidence toward fishing effort reduction.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
LZ was responsible for data collection, formal analysis, and writing original draft. CL was responsible for data curation, reviewing, and editing the manuscript. DP was responsible for conceptualization, methodology, and supervision. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
LZ would like to thank the China Scholarship Council (CSC). CL thanks the Qingdao Postdoctoral Applied Research Project. DP acknowledges support by the Sea Around Us, funded by a number of philanthropic foundations. The authors also thank Rainer Froese for his comments on an early draft of this contribution and Evelyn Liu for the figures.