Stock Assessment of Small Yellow Croaker (Larimichthys polyactis) Off the Coast of China Using Per-Recruit Analysis Based on Bayesian Inference

This paper presents a framework for quantifying uncertainty in per-recruit analysis for small yellow croaker (Larimichthys polyactis) fisheries in China, in which credible estimates of life history parameters from Bayesian inference were used to generate the distribution for a quantity of interest. Small yellow croakers were divided into five spatial groups. The status of each group was examined using a yield-per-recruit (YPR) model and a spawning stock biomass-per-recruit (SSBPR) model. The optimal length at first capture (Lcopt) was proposed to recover the biomass. The maximum observed age in the current stocks (3 years) and the maximum recorded age (≥20 years) were adopted in per-recruit analysis. Our results suggest that the framework can quantify uncertainty well in the output of per-recruit analysis for small yellow croaker. It is suited to other fish species. The SSBPR at FMSY (SSBPRMSY) is a better benchmark than the spawning potential ratio (SPR) at FMSY because SSBPRMSY had a unimodal distribution. The SSBPR analysis can lead to a more conservative Lcopt than the YPR analysis. The key factor influencing the assessment conclusions may be the growth parameters rather than the natural mortality rate for a stock with a younger maximum age. Overfishing likely occurred for all groups and recruitment overfishing may not occur if the maximum age is maintained at 3 years. Increasing lengths at first capture to the recommended values can help this population recover. However, Fcur is too high for small yellow croakers to attain the maximum recorded age. Both reducing fishing mortality rate and increasing length at first capture are needed to attain the maximum recorded age.


INTRODUCTION
Per-recruit models such as the yield-per-recruit (YPR) model and the spawning stock biomass-perrecruit (SSBPR) model are commonly used in fish stock assessments. These models can produce YPR and SSBPR at a fishing mortality rate (F), and biological reference points. Per-recruit models are defined mainly by fish life history parameters, such as the parameters on growth, length-weight relationship (LWR), and natural mortality. Thus, the output of a per-recruit model is primarily dependent on a combination of different types of life history parameters (Restrepo and Fox, 1988). The inherent uncertainty surrounding life history parameters can lead to uncertainty in per-recruit analysis (Chang et al., 2009). Incorrect estimates of biological reference points may be derived if uncertainty in parameters is ignored. Consequently, conclusions about the status of fish stock will be inaccurate (Helser et al., 2001;Chen and Wilson, 2002;Lin et al., 2015). Quantifying uncertainty in assessment results can enhance decision making for fishery resource management and help determine the best harvest strategies (Gavaris, 1993;Doll et al., 2017).
Random values from the distributions of input parameters must be generated to quantify uncertainty in the outcome of per-recruit analysis. Conventionally, these random values are produced by introducing random errors into each subset model with a set of particular estimates. A coefficient of variation (CV) or standard deviation (SD) is often specified for a subset model according to the available knowledge of parameter uncertainty (e.g., Chen and Wilson, 2002;Grabowski and Chen, 2004;Jiao et al., 2005;Chang et al., 2009). The estimated SDs from the observed data can also be used to introduce random errors (Lin et al., 2015). For correlated parameters, random errors can be introduced using the estimated standard errors and the correlation (Hart, 2013), and the covariance matrix constructed by the estimated SDs and the correlation (Lin et al., 2010). For instance, Hart (2013) took the standard error of 0.004 for K and 0.4 for L ∞ , and the correlation of −0.6 to simulate the negatively correlated parameters K and L ∞ in the von Bertalanffy growth model (von Bertalanffy, 1938). Recent advances in Bayesian estimation provide a new method to obtain random values for life history parameters. Bayesian analysis can incorporate multiple sources of information to reduce uncertainty in life history parameter estimates and obtain more realistic estimates of uncertainty (Pulkkinen et al., 2011;Froese et al., 2014;Romakkaniemi, 2015;Zhu et al., 2020). Possible parameter estimates derived from Bayesian inference are given as samples drawn from its posterior distribution using the Markov chain Monte Carlo (MCMC) technique. MCMC samples are a number of estimates for an independent parameter or pairs of estimates for the correlated parameters. MCMC samples can generate distributions for the output of a per-recruit model. Doll et al. (2017) demonstrated that uncertainty in YPR estimates could be measured using the posterior distributions of growth parameters, LWR parameters, and natural mortality rate (M) when a combination of exploitation rate and minimum length limit was given.
Small yellow croaker (Larimichthys polyactis) is a warmtemperate demersal fish species that is distributed throughout the northwest Pacific Ocean (Ma et al., 2017). It is a commercially valuable species in China and inhabits the Bohai Sea, Yellow Sea, and East China Sea (Figure 1). Biological characteristics of the stocks have significantly changed; average body size is smaller and individuals become mature at an earlier age (Jin, 1996;Li et al., 2011;Lin et al., 2011). The maximum age decreased from 21-23 years in the 1960s (Mao et al., 1987;Guo et al., 2006) to 5-7 years in the early and mid-1980s (Mao et al., 1987;Liu et al., 1999;Shui, 2003;Lin et al., 2004;Guo et al., 2006), to 2-4 years in the 1990s and mid-2000s (Liu et al., 1999;Shui, 2003;Lin et al., 2004;Guo et al., 2006;Yan et al., 2006;Chen et al., 2010;Shan et al., 2017). Changes in life history traits of small yellow croaker have been attributed to increased fishing effort (Shan et al., 2013). The population has been at high risk of being overfished in recent years (Ma et al., 2020). Fishing effort needs to be reduced to recover the population to biomass level that supports maximum sustainable yield (MSY) (Ma et al., 2020). However, it is difficult to control fishing effort in China (Shen and Heino, 2014). The fishing effort of small yellow croaker increased each year from 1968 to 2015, with the exception of 1991 (Ma et al., 2020). Measures are urgently needed to reduce fishing pressure on this population. In addition, it is unclear if recruitment overfishing has occurred.
Trawl selectivity experiments demonstrated that enlarging mesh size could reduce fishing pressure on small yellow croaker, especially for juveniles (Chen et al., 2018). Thus, increasing length at first capture (L c ) may be an efficient measure to recover the population of small yellow croaker. Researchers have recommended the optimal age (t copt ) or length (L copt ) at first capture for small yellow croaker in different regions based on the results from the YPR model (Jin, 1996;Liu et al., 1999;Lin et al., 2006;Zhang et al., 2010a;Gao et al., 2019). However, uncertainty in input parameters was not considered. The derived L copt varied among these studies. In addition, Zhou et al. (2020) suggested the fishing mortality rate to achieve MSY (F MSY ) or the spawning potential ratio (SPR) at F MSY (SPR MSY ) as the benchmark because the surplus production model takes into account population dynamics from one generation to the next. Thus, there may be a gap between L copt derived from the YPR analysis and the SSBPR analysis.
In a previous study, the probability distributions of life history parameters were estimated using the Bayesian hierarchical approach for five groups of small yellow croakers (Zhu et al., 2020). Based on this work, a framework was presented that quantified uncertainty in per-recruit analysis for small yellow croaker and examined the status of each group of small yellow croakers using the YPR model and the SSBPR model by incorporating uncertainty in parameters. L copt was estimated when SSBPR at F MSY (SSBPR MSY ) was regarded as the management goal. The objective of this study was: (1) to establish a procedure for quantifying uncertainty in per-recruit analysis for small yellow croaker based on Bayesian inference and (2) to further check the status of small yellow croaker and to recommend management measures by considering both the yield and the spawning stock biomass.

Methods and Data Sources
Previous studies have recommended to consider spatial heterogeneity in the studies for small yellow croaker due to the existence of geographic subpopulations for this species (Ying et al., 2011;Ma et al., 2020). However, division of the subpopulation for small yellow croakers is a controversial issue (Liu et al., 2013). Lin et al. (1965) divided small yellow croakers into three subpopulations based on body morphometrics in which the boundaries were 32.0 and 34.0 • N (Figure 1). Based on this division, Liu (1990) proposed a fourth subpopulation in the central Yellow Sea (between 34.0 and 37.5 • N) according to sex gland, which was supported by Wang et al. (2016) based on analysis of stable isotopic composition of otoliths. However, Xu and Chen (2010) argued that only two subpopulations exist according to their migration trajectory, and they defined a boundary of 34.0 • N. Zhang et al. (2014) proposed another division of four subpopulations according to otolith features and the boundaries were about 35.5, 30.0, and 27.0 • N.
In response to the controversy, small yellow croakers were divided into five groups according to the intensive investigation regions (Zhu et al., 2020), in which Group 3 was included in Group 4 (Figure 1). Some researchers studied small yellow croaker in the region between 32.0 and 34.0 • N (e.g., Zhang et al., 2010a;Liu et al., 2013) according to the subpopulation divisions proposed by Lin et al. (1965), whereas some investigated the population in the area between 31.0 and 34.0 • N (e.g., Lin et al., 2004;Yan et al., 2006). Thus, Groups 3 and 4 were regarded as two independent groups since it is unclear whether they belong to the same subpopulation (Zhu et al., 2020). Zhu et al. (2020)  used the Bayesian hierarchical models to estimate life history parameters and their uncertainty for each group based on 25 growth curves, 41 LWRs, and 16 natural mortality rates collected from the literature (Supplementary Table 1). The resulting estimates of each type of life history parameter were 30,0000 sets of MCMC samples for each of the five groups. Growth parameters, LWR parameters, and natural mortality rate were summarized and are shown in Supplementary Tables 2-4, in which highest density interval (HDI) that gave the span of the most credible values (Kruschke, 2014) was used to summarize the uncertainty of a parameter.
Based on this work, we developed a procedure for quantifying uncertainty in per-recruit analysis for small yellow croakers (Figure 2). It consisted of two steps: estimating the posterior probability distributions for life history parameters and conducting per-recruit analysis. In the first step, already developed and undertaken by Zhu et al. (2020), the posterior distributions of life history parameters are first obtained using Bayesian approaches by combining the priors and the existing estimates, which were used to conduct per-recruit analysis when newly observed data were not collected. Otherwise, the resulting posteriors served as the priors to update the parameters by combining with new data. Then, the latest knowledge about life history parameters was used to carry out per-recruit analysis. For the second, novel step, the credible MCMC samples of each parameter generated from the first step were used to conduct per-recruit analysis and to calculate the quantities of interest by combining with fishery parameters, such as F and L c .
We used the parameter values that had been obtained from the first step to undertake the second step. New data for estimating the growth parameters, LWR parameters, and natural mortality 1 | Number of the total parameter sets in per-recruit analysis and number of the parameter sets failing to define F max and F 20% .

Group
Total F max F 20% t max = 3 t max = 5 t max = 7 t max = 20-23 t max = 3 t max = 5 t max = 7 t max = 20-23 rate were not collected because almost all the recent estimates have been incorporated in Bayesian analysis in our previous work (Zhu et al., 2020). The MCMC samples inside the 95% HDI were taken as credible estimates to avoid the impact of extremely high or low values on per-recruit analysis. For the correlated parameters, a pair of estimates were taken together if the estimate of each parameter was located in its 95% HDI. The number of credible MCMC samples was different for different types of parameters. The same number of MCMC samples was taken according to the minimum number of credible MCMC samples among three types of parameters for a group. The number of parameter sets was enough to obtain the possible distributions for the output of a per-recruit model ( Table 1).

Per-Recruit Model
Considering implementation of summer-fishing-closed-season policy in China, a discrete YPR model and a discrete SSBPR model were constructed for small yellow croaker fisheries based on the work of Govender et al. (2006). YPR (expressed in mass g) was expressed as (1) where Y is the attained yield, R is the number of the recruit, and t max is the maximum age in the fishery in months a w and b w are LWR parameters (W = a w L b w ). L t is the predicted length at age t, which was calculated using von Bertalanffy growth model (von Bertalanffy, 1938). j is used to produce the cumulative sum of fishing mortality rate and natural mortality rate at each stage, and this term exists when t ≥ 2. Here, R was assumed to be the number at 1 month old and set to 1. S L t is the gear selectivity at length L t and was assumed to be knife-edged.
A t indicated whether a particular month corresponding to age t was open to fishing or not (Govender et al., 2006). It was 1 when the month was open to fishing. Otherwise, it was 0.
Spawning stock biomass-per-recruit (expressed in mass g) was expressed as where G t was the fraction of mature fish at age t and was also assumed to be knife-edged (Govender et al., 2006).
where t m was the age at 50% maturity. Spawning potential ratio is the ratio of SSBPR at a fishing mortality rate to the maximum SSBPR under unfished conditions (Goodyear, 1993). It was expressed as Spawning potential ratio has a maximum value of unity and declines toward zero with an increase in fishing mortality rate (Goodyear, 1993). Similarly, SPR MSY was expressed as (Zhou et al., 2020)

Biological Reference Points
Biological reference points derived from per-recruit models are expressed as fishing mortality rates, which are metrics of stock status (Gabriel and Mace, 1999). F max and F 0.1 are reference points derived from the YPR model. F max is the fishing mortality rate maximizing YPR and F 0.1 is where the slope of the YPR curve is 10% of the maximum slope (Grabowski and Chen, 2004). The reference point from the SSBPR model is expressed as F x% representing the fishing mortality rate that reduces SSBPR to x% of the maximum SSBPR in an unexploited state (Gabriel and Mace, 1999). F x% is exactly the value of SPR at a fishing mortality rate according to Eq. 5. Mace and Sissenwine (1993) advocated F 20% as a recruitment overfishing threshold for species with high and medium resilience and F 30% as a recruitment overfishing threshold for species with low resilience. Generally, small yellow croaker grows slowly, matures at an older age (2-3 years), and has low resilience (Ye, 1991;Jin and Deng, 2000). However, its resilience may increase when it matures earlier. At present, FishBase 1 gives the resilience of small yellow croaker as medium. Thus, we took F 20% as a recruitment overfishing threshold for small yellow croaker. We calculated F max , F 0.1 , and F 20% to examine stock status and check the impact of uncertainty in life history parameters on per-recruit analysis.

Impact of Uncertainty in Parameters on Biological Reference Points
To examine the impact of a type of history life parameter on uncertainty in per-recruit analysis, reference points were calculated under four scenarios. We used credible MCMC samples of all types of life history parameters to calculate reference points in Scenario 1. Credible MCMC samples of the growth parameters, LWR parameters, and natural mortality rate were adopted in Scenarios 2, 3, and 4, respectively, and point estimates were taken for other types of parameters. The variation in biological reference points in Scenarios 2, 3, and 4 was caused by the uncertainty of a type of life history parameter.

Effect of Maximum Age
The maximum age may have not exceeded 3 years for each group of small yellow croakers since the mid-2000s Guo et al., 2006;Chen et al., 2010;Yan et al., 2014;Hu et al., 2015;Shan et al., 2017). It was 2 years in some years Guo et al., 2006;Yan et al., 2014;Hu et al., 2015). In this study, a maximum age of 3 years was used in per-recruit analysis for all groups. The maximum age influences YPR and SSBPR for a given F according to Eqs 1, 3, and may further affect the derived reference points. We set t max = 3, 5, and 7 years to test the effect of maximum age. At present, it is unclear if small yellow croaker would reach an age of 21-23 years in the absence of fishing. The optimal maximum age (t optmax ) was 21 years for Groups 1 and 2, 23 years for Groups 3 and 4, and 20 years for Group 5 (Gao et al., 2019). We calculated reference points and derived L copt when t max = t optmax to further check the status of small yellow croaker and to explore management measures.

Optimal Length at First Capture
In past studies, the YPR contour was often used to determine t copt , which was transformed to L copt . In the YPR isopleth, age at first capture (t c ) was the vertical axis and F the horizontal axis. One line representing maximum YPR for a given F and the other line representing maximum YPR for a given t c were plotted in the contour (Caddy, 1984). At the current fishing mortality rate, any t c located between these two lines can be chosen as t copt . Thus, t copt has a range and one value is usually chosen subjectively.
In this study, we chose SPR MSY and SSBPR MSY as the target reference points to which to aim to recover the biomass of small yellow croakers and calculated their distributions. As shown in Figure 2, L copt was found to achieve these distributions.

Fishing Mortality Rate and Other Parameters
Fishing mortality rate has an increasing trend for small yellow croaker fisheries which is in accordance with increasing fishing effort (Supplementary Table 5). Considering that fishing effort has consistently increased and fishing mortality rate also has a level of uncertainty, current fishing mortality rates (F cur ) were assigned the mean of the latest two estimates for Group 1 (1.73 year −1 ) and Group 5 (2.73 year −1 ). F cur took the latest estimate for Group 2 (1.56 year −1 ), Group 3 (2.38 year −1 ), and Group 4 (2.38 year −1 ) because other values were estimated four or more years ago for those groups. Groups 3, 4, and 5 had much higher fishing mortality rates than Groups 1 and 2, which is consistent with fishing activities (Ma et al., 2020). F MSY was timevarying and ranged from 0.15 to 1.07 year −1 for Groups 1 and 2, and from 0.15 to 1.32 year −1 for Groups 3, 4, and 5 between 1968 and 2015 (Ma et al., 2020). We took the values in 1968 and 2015 for small yellow croakers with the optimal maximum age and the current maximum age, respectively.
We assumed that small yellow croakers grew up to 1 month older in May for Group 5 and in June for other groups according to the spawning time of each group (Jin, 1996;Cheng et al., 2004;Zhang et al., 2010b). L c was set to 12 cm for all groups, which is the length at 50% retention for the trawl net and the canvas stow net with the minimum allowable mesh size of 54.0 mm (You et al., 2017;Chen et al., 2018;Xu et al., 2019). Heavy fishing pressures have resulted in earlier maturity for small yellow croaker. The proportion of the population that matures at age 1 has been over 90% since 2011 (Yan et al., 2014). The fractions of mature fish at spawning time were assumed to be 1 for all age classes of each group. Currently, the areas where small yellow croakers are distributed are closed to fishing in the summer. The seas to the north of 35.0 • N are closed to fishing from May 1 to August 31 and those to the south of 35.0 • N are closed from May 1 to September 16.

Biological Reference Points
Generally, biological reference points F max , F 0.1 , and F 20% were well defined by credible MCMC samples. F 0.1 existed for all set of credible MCMC samples of each group. However, some sets of MCMC samples could not produce F max and F 20% ( Table 1). Some extremely high values were also derived for F 20% from some sets of parameter estimates (Figure 3). F cur was significantly higher than F MSY for each group because it was greater than the upper limit of 95% credible interval ( Table 2). It was also significantly higher than F 0.1 for Groups 1, 3, 4, and 5 when t max = 3 years. But F cur was not significantly greater than F max and F 20% for each group. It even had a high probability of being lower than F max and F 20% for Group 2 (Figure 4). The probability of F cur being greater than F 0.1 , P (F cur > F 0.1 ), approximated 1 for all groups except Group 2 (92%). In contrast, P (F cur > F max ) and P (F cur > F 20% ) were relatively low and not greater than 74%. However, F cur was greater than F 0.1 , F max , and F 20% with a high probability for all groups (85-97%) when t max = t optmax . Figure 3 shows the distributions of F max , F 0.1 , and F 20% for Group 1 in four scenarios when t max = 3, 5, and 7 years. Reference points of the other four groups had a similar distribution pattern. Uncertainty about each reference point was largest for Scenario 1 at a given maximum age, in which all three types of parameters varied. For Scenarios 2, 3, and 4, the variation in the growth parameters (Scenario 2) led to the largest uncertainty in the estimates of each reference point when t max = 3 years. This indicated that the uncertainty of these FIGURE 3 | Boxplot of biological reference points F max , F 0.1 , and F 20% in four scenarios for Group 1 of small yellow croakers. Scenario 1 denoted the variation in all types of life history parameters. Scenarios 2, 3, and 4 denoted the variation in growth parameters, length-weight relationship parameters, and natural mortality rate, respectively.

Impact of Uncertainty in Parameters on Biological Reference Points
reference points was caused mainly by uncertainty surrounding the growth parameters.

Effect of Maximum Age
The number of the parameter sets that could not define F 20% decreased with an increase in t max ( Table 1). Relative importance of life history parameters also changed as t max increased (Figure 3). The effect of natural mortality rate (Scenario 4) on the uncertainty of reference points became stronger than growth parameters when t max = 5 years. Uncertainty in the natural mortality rate dominated the variation in each reference point when t max = 7 years. The LWR parameters had less influence on the uncertainty of the reference points in all cases.

Yield-Per-Recruit and Spawning Potential Ratio
Both the YPR curve and the SPR curve were variable for each group due to uncertainty in the life history parameters FIGURE 4 | The probability of F cur being greater than the biological reference points F max (green triangle), F 0.1 (red square), and F 20% (blue circle) for each group of small yellow croakers.
( Figures 5, 6). They also varied among groups, depending on parameter estimates of each group and the level of their uncertainty. Maximum YPR ranged from 25.54 g (median value) to 34.20 g when L c = 12.0 cm for small yellow croakers and was significantly different among groups when the exploitation pattern was similar (p < 0.001, Figure 7). The median of maximum YPR was greatest for Group 1, followed by that of Groups 4, 2, 5, and 3. However, the SPR curve of Group 1 was the lowest, whereas the SPR curves of Groups 2, 3, and 5 were relatively high (Figure 6). Median SPR at F cur (SPR cur ) was lower than median SPR MSY for each group. The upper limit of 95% HDI was even lower than 40% for Groups 3, 4, and 5. The distribution of SPR at a fishing mortality rate was different at two maximum ages (Figure 8). SPR cur and SPR MSY had a bimodal distribution or a multimodal distribution when t max = 3 years, whereas they had a unimodal distribution for each group when t max = t optmax . But the pattern of the SSBPR at a fishing mortality was not affected by the maximum age (Figure 9). Median SPR cur was very low (10-15%) for small yellow croakers with the optimal maximum age due to the high fishing mortality rate. A high SPR MSY (55-66%, median value) was needed to attain MSY. Conversely, SPR cur was relatively high when t max = 3 years and a relatively low SPR MSY was needed.
Optimal Length at First Capture SPR MSY was not suitable to be used as the benchmark to derive L copt when t max = 3 years because its distribution had two or more modes (Figure 8). We derived L copt according to the distribution of SSBPR MSY , which ranged from 14.80 to 15.55 cm when t max = 3 years and from 22.19 to 23.30 cm when t max = t optmax (Figure 9). SSBPR at F cur (SSBPR cur ) and SSBPR MSY had a similar distribution when L copt was used for each group. But uncertainty in the SSBPR increased when L c increased to the value of L copt as this increased the distribution of the SSBPR for each group. L copt could also lead to a higher yield when t max = 3 years ( Table 3).

DISCUSSION
This study presented a framework for measuring uncertainty in the outputs of per-recruit analysis for small yellow croaker fisheries. In this framework, the Bayesian approach was used to estimate probability distributions for life history parameters. Credible MCMC samples from Bayesian analysis were used to generate the distributions for a quantity of interest from a YPR model and an SSBPR model. Based on previous work, we checked the status of five groups of small yellow croakers and derived optimal lengths at first capture when SSBPR MSY was used as the benchmark. The results showed that uncertainty in a quantity of interest from per-recruit analysis, such as YPR and SSBPR at a fishing mortality rate, could be well quantified using credible MCMC samples and expressed as a probability distribution (Figures 7, 9). The credible range of a quantity could be further expressed as the 95% HDI (Tables 2, 3). Uncertainty in growth parameters, LWR parameters, and natural mortality rate created variation in the reference points (Figure 3), YPR curve (Figure 5), and SPR curve (Figure 6).
In theory, F 0.1 always exists, but F max is poorly defined in some cases, and even cannot be defined (Rivard and Maguire, 1993). When the natural mortality rate is high relative to growth, the combination of life history paramters will lead to the asymptotic YPR curve, in which F max cannot be found (Jensen, 2000). That F 20% could not be defined was also due to the low growth and a relatively high natural mortality rate. SPR is positively correlated with natural mortality and negatively correlated with various indices of size (Caddy and Mahon, 1995). The SPR curve will be too high to intersect the horizontal line passing through 20% when growth is low relative to the natural mortality rate. Although some estimates of F 20% were extremely high when t max = 3 years (Figure 3), the 95% HDI could exclude them and provide the range of acceptable values ( Table 2).
Overfishing likely occurred for all groups of small yellow croaker because F cur was significantly higher than F MSY for each group ( Table 2). This is consistent with the conclusions made by Ma et al. (2020). A similar conclusion could be drawn if F 0.1 was taken as the benchmark because P(F cur > F 0.1 ) was high for each group when t max = 3 years (Figure 4), indicating that F 0.1 can index the status for small yellow croaker at present. F max is not suited for judging the status of small yellow croaker when t max = 3 years because P(F cur > F max ) was not high, especially for groups 1, 2, and 3. Generally, F max is not a useful conservation standard because it may be extremely high to maximize yield, which will reduce the spawning potential of the stock to near zero (Goodyear, 1993). Recruitment overfishing may not occur for the current stocks (t max = 3 years) because P(F cur > F 20% ) was also not high (Figure 4). But F cur was higher than both F max and F 20% with a high probability for each group when t max = t optmax , indicating that the current fishing mortality rate was too high for small yellow croakers with the optimal maximum age. Recruitment overfishing might have occurred in the past, but the stocks recovered by changing the life history traits. For instance, small yellow croakers now grow faster and mature earlier.
The level of uncertainty in reference points was jointly determined by uncertainty in three types of life history parameters. However, the contribution of each type of parameter was different, depending on its level of uncertainty. Generally, the LWR parameters have a low level of uncertainty for most fished species because length-weight data are easy to collect and can be measured accurately, whereas uncertainty in growth parameters and natural mortality rate is relatively high. For each group of small yellow croaker, the estimated LWR parameter was almost normally distributed and had the lowest uncertainty among the three types of life history parameters (Zhu et al., 2020). Thus, the LWR parameter had limited influence on all reference points in this study (Figure 3). Natural mortality rate is more difficult to estimate accurately in comparison to growth parameters. Estimates of the natural mortality rate have a high degree of uncertainty (Hamel, 2015). Thus, natural mortality rate has been commonly identified as a key factor influencing the stock assessment for a variety of species (e.g., Grabowski and Chen, 2004;Jiao et al., 2005;Chang et al., 2009). However, uncertainty surrounding growth parameters was the main contributor to the uncertainty in all reference points for small yellow croaker stocks when t max = 3 years (Figure 3). But the relative importance of uncertainty in growth parameters and natural mortality rate changed with the assumed maximum age. The impact of uncertainty in natural mortality rate on perrecruit analysis became stronger as the maximum age increased. Natural mortality rate became the key source of uncertainty in stock assessments when t max = 7 years. These results implied that growth parameters may be a key factor influencing the assessment conclusions for a stock with a younger maximum age, unless uncertainty in growth parameters is lower than that of the natural mortality rate.
Maximum age affected the per-recruit analysis. The number of parameter sets that could not generate F 20% reduced with an increase in maximum age ( Table 1). Larger maximum age values resulted in more SSBPR terms according to Eq. 3. Correspondingly, the SPR curve may have a shape that can define F 20% . The estimates of all reference points decreased with an increase in maximum age (Figure 3). That is, fishing mortality rate must be further reduced to let small yellow croakers grow older. A very low fishing mortality rate is needed to maintain the optimal maximum age for small yellow croakers ( Table 2). An appropriate choice of the maximum age is particularly important for SSBPR analysis because the maximum age determines the magnitude of maximum SSBPR and further influences SPR. For instance, SSBPR cur was very close for the two maximum ages (Figure 9), but different SPR cur was produced (Figure 8) due to the effect of maximum SSBPR. Normally, a maximum age exists for an unexploited fish stock. A younger maximum age may be observed from a sample after this stock is exploited. If this younger maximum age is employed in the SSBPR model, maximum SSBPR will be underestimated (Mace and Sissenwine, 1993). As a result, a positive bias will be introduced into the estimate of SPR (Goodyear, 1993). In this situation, the maximum age of unexploited stock must be used in the SSBPR analysis. It is worth discussing the choice of maximum age for an overexploited stock that has a much younger maximum age for many years due to long-term heavy fishing. For example, the reported oldest age for small yellow croaker is 21 years in the Bohai Sea (Guo et al., 2006) and 23 years in the southern Yellow Sea (Mao et al., 1987), but the maximum age of small yellow croakers in each region may have not exceeded 3 years for about 15 years. Therefore, it is not realistic to adopt t max = 21 or 23 years to calculate the YPR and the SSBPR for current stocks. But these optimal maximum ages can be used to calculate management quantities, such as the reference points in Table 2.
There was a significant difference in maximum YPR among groups (Figure 7). YPR at a fishing mortality rate is positively correlated with growth indices and negatively correlated with natural mortality rate according to Eq. 1. Therefore, maximum YPR of Group 1 was greatest mainly because of its highest asymptotic average body size and lowest natural mortality rate among the five groups (Supplementary Tables 2, 4). The lowest SPR curve in Group 1 was also caused by its high growth indices and lowest natural mortality rate, which can be attributed to a negative correlation between SPR and growth, and a positive correlation between SPR and natural mortality. Median SPR MSY was higher than 40% for all groups of small yellow croaker except Group 4 when t max = 3 years (Figure 6). Each group required a higher SPR MSY when t max = t optmax (Figure 8). At present, spawning stock may drop below the level for supporting MSY. Although seas are closed to fishing in summer every year in China, this closed season policy does not prevent overfishing for small yellow croakers; rather, it temporarily enhances the biomass of small yellow croakers Yue et al., 2016). Relative abundance of small yellow croakers increases significantly during the closed months and is the greatest at the end of the closure . However, relative abundance is reduced to a low level again in late December due to high fishing pressure . Additional management measures are needed to ultimately enhance stocks.
Spawning stock biomass-per-recruit is superior to SPR when deriving L copt . It was easy to derive L copt using SSBPR MSY because it had a unimodal distribution (Figure 9). However, a younger maximum age could cause the complex distribution for SPR MSY (Figure 8), which limited the application of SPR. The derived L copt could improve both SPR and YPR for each group when t max = 3 years (Table 3 and Figure 8). But L copt was different among groups. Groups 1 and 2 are commonly managed by one agency, whereas Groups 3-5 are co-managed by another agency. Vessels can operate anywhere within the jurisdiction of an agency. One policy that could be implemented for Groups 1 and 2 is to use a value of 15.05 cm for L copt . This policy could be easily implemented in practice. Similarly, a L copt of 15.55 cm can be used for Groups 3-5. These policies would be more conservative for Groups 2-4. It is not realistic to increase L c to 22.19-23.30 cm (Figure 9). Fishing mortality rate also needs to be reduced to let small yellow croakers reach the optimal maximum age.
The proposed L c from YPR analysis is usually lower than the recommended value in this study. According to the YPR isopleth, Lin et al. (2006) proposed L c = 17.0 cm (L copt range = 15.48-17.43 cm) for Groups 3-5. Zhang et al. (2010a) suggested L c = 14.83 cm (L copt range = 11.89-15.42 cm) for Group 3. Gao et al. (2019) recommended L c = 15.0 cm (L copt range = 15.15-15.91 cm) for Group 5. Although the upper limit of the L copt range estimated by Gao et al. (2019) was higher than 15.55 cm, a relatively small value was chosen. Only the recommended L c (17.0 cm) by Lin et al. (2006) was high enough to conserve Groups 3-5. But the YPR produced by an L c of 17.0 cm at F cur was lower than that produced by an L c of 15.55 cm for these three groups. Zhai and Pauly (2019) derived an L copt of 13.5 cm from YPR analysis and 14.0 cm from utility-per-recruit analysis for Group 1 based on the parameter estimates of Liu et al. (2018). These two values were lower than our result (15.05 cm). L copt might have been underestimated by Zhai and Pauly (2019) due to limited data.
The estimated parameters are often imprecise and different between independent studies because the oberved data are a subsample from the population of interest (Doll et al., 2017). Uncertainty in per-recruit analysis will be dominated by the observed data when parameter estimates and their SDs derived from these data are used to generate random parameters. Uncertainty of the parameters estimated from a subsample may FIGURE 8 | Distribution of spawning potential ratio at F cur (red lines) and at F MSY (blue lines) for each group of small yellow croakers when L c = 12.0 cm (t max = 3 years: solid lines; t max = t optmax : dotted lines).
be narrower than when data from other years are considered (Froese et al., 2014). Thus, the resulting uncertainty in a quantity from per-recruit analysis may be underestimated. When random parameters are obtained by specifying SDs or CVs for a particular set of parameter estimates, the resulting uncertainty in the outcome of a per-recruit model will be influenced by both the specified estimates and SDs or CVs. In contrast, existing information regarding a parameter can be introduced to Bayesian analysis and the resulting uncertainty of this parameter may be more realistic. Therefore, parameter estimates from Bayesian inference may lead to more realistic uncertainty for a quantity of interest from per-recruit analysis. In addition, the correlated parameters from different types of life history can also be easily incorporated into a Bayesian model to estimate the parameter of interest using their correlation (Pulkkinen et al., 2011;Zhu et al., 2020). Moreover, uncertainty in a parameter with less information can be reduced through their correlation. Thus, estimating parameters using Bayesian inference is suited for quantifying uncertainty in per-recruit analysis.
Generating extremely low or high values from Bayesian inference is unavoidable when sampling in the posterior distribution of a parameter using the MCMC technique. The extreme values may have no biological meaning or exceed the reasonable range. We used MCMC samples inside 95% HDI to avoid the impact of extreme values. The total probability of MCMC samples in the 95% HDI is 95% and these samples are the 95% most credible values (Kruschke, 2014). These MCMC samples can generate uncertainty for a quantity for per-recruit analysis. The median of a quantity was similar when all MCMC samples were used. The mean had an increase of 2-5% for F max , 1-4% for F 0.1 , and 1-3% for F 20% when t = 3 years. The 95% HDI became wider due to the increase in the range of MCMC samples; 95% HDI of F max had a largest change among these three reference points. Its lower limit decreased by 3-11% and upper limit increased by 5-10%. The variation in the mean and the median of L copt was no more than 1% and was less than 2% in the limits of 95% HDI.
Markov chain Monte Carlo samples used in per-recruit analysis were derived from the existing estimates of life history parameters. Growth parameters and natural mortality rate were slightly skewed for most groups of small yellow croakers due to variation in parameter estimates and the relatively small size of the dataset (Zhu et al., 2020). As a result, the distribution of a quantity from per-recruit analysis was also not exactly normal (Figures 7, 9). New data are needed to reduce uncertainty in the output of per-recruit analysis by improving parameter estimates. According to the framework shown in Figure 2, when new data are collected for a group of small yellow croakers in the future, the distributions of parameters estimated by Zhu et al. (2020) will serve as the prior values to derive the newest distributions FIGURE 9 | Distribution of spawning stock biomass-per-recruit at F cur and at F MSY for each group of small yellow croakers when L c = 12.0 cm and L c = L cop . L copt was 15. 05, 14.80, 15.15, 15.20, and 15.55 cm for Groups 1-5 when t max = 3, and 23. 30, 22.67, 22.19, 22.76, and 22.35 cm when t max = t optmax .
TABLE 3 | Summary of yield-per-recruit at the current fishing mortality rate for each group of small yellow croakers when t max = 3 years (unit: g).

Statistic
Group 1   for life history parameters. Then the resulting credible MCMC samples are used to conduct per-recruit analysis to update our knowledge on a quantity of interest from per-recruit analysis. This procedure for measuring uncertainty in per-recruit analysis is general and suitable for other species. In this framework, other Bayesian approaches can also be adapted to obtain MCMC samples of life history parameters for quantifying uncertainty in the outcomes of per-recruit analysis. The conclusion regarding the status assessment and the recommended L c in this study were also influenced by the estimates of F cur and F MSY . When new knowledge on F cur and F MSY is available, per-recruit analysis needs to be carried out again to update our understanding of the status of each group. In addition, the derived L copt for small yellow croakers with the optimal maximum age was influenced by parameter values and F MSY . Current parameter values may be different from those when small yellow croakers attained the optimal maximum age. Similarly, F MSY in 1968 may also be biased. The groups of small yellow croakers were divided in terms of conventional investigation region. It is still unclear whether or not Groups 3 and 4 belong to the same subpopulation. Although Group 3 is the main part of Group 4 in this study, their maximum YPRs were significantly different (Figure 7). Group 3 had larger median values of F max , F 0.1 , and F 20% than Group 4 (Table 2), whereas Group 4 had higher YPR and SSBPR (Table 3 and Figure 9). There is a need to discuss the division of groups for small yellow croakers to conduct better resource surveys and stock assessment for this species.

CONCLUSION
Overfishing likely occurred for all groups of small yellow croakers in China, but recruitment overfishing may not occur if the maximum age is maintained at 3 years. High fishing mortality and small length at first capture might cause SPR cur to be lower than SPR MSY . The current closed season policy coupled with the recommended lengths at first capture could recover the biomass to support MSY when t max = 3 years. The required fishing mortality rate was very low for small yellow croakers with the optimal maximum age.
Besides increasing length at first capture, the fishing mortality rate needs to be reduced to attain the optimal maximum age for each group.
The combination of the SSBPR model and F MSY from the surplus production model can lead to a more conservative length at first capture than the YPR model. SSBPR MSY is more suitable to derive the optimal length at first capture than SPR MSY .
The credible parameter values from Bayesian inference could quantify uncertainty well in the output of per-recruit analysis. The presented framework for measuring uncertainty in the outcome of per-recruit models is suited for other fish species, and the Bayesian approach to parameter estimation in this framework is not limited to that utilized in our previous work: other Bayesian methods can also be used to derive input parameters for perrecruit models.

DATA AVAILABILITY STATEMENT
The data analyzed in this study are subject to the following licenses/restrictions: The MCMC samples from Bayesian inference will be shared on reasonable request to the corresponding author. All other datasets generated for this study are included in the article/Supplementary Material. Requests to access these datasets should be directed to ZL, liangzhenlin@sdu.edu.cn.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because our manuscript was based on data from literatures and no live vertebrates or higher invertebrates were involved.

AUTHOR CONTRIBUTIONS
LZ performed the computing work and wrote the manuscript. CG analyzed the data. ZJ and CL drafted the figures. GH collected the data. ZL designed this study and provided financial support. All authors contributed to the article and approved the submitted version.