Assessment and Management of Small Yellow Croaker ( Larimichthys polyactis ) Stocks in South Korea

: We aimed to determine the appropriate annual total allowable catch (TAC) levels for the small yellow croaker ( Larimichthys polyactis ). A Bayesian state-space model was used to assess the species stock. This model has been widely used after research conﬁrmed its reliability over other models. However, setting prior distributions for analyzing this model remains controversial. Therefore, a sensitivity analysis was conducted using the model with di ﬀ erent prior distributions and biomass growth functions. Informative and non-informative prior distributions were compared using Schaefer and Fox growth functions. Considering the results of the sensitivity analysis, the assumption of inverse-gamma prior distribution of K , a non-informative distribution, with the Fox function could yield relatively superior estimates than those obtained from other assumptions. Moreover, changing the growth function could have a greater e ﬀ ect on the ﬁtness of the model estimates than changing prior distribution. Therefore, future ﬁshery stock analyses based on this model should consider the e ﬀ ectiveness of various growth functions in addition to the sensitivity analysis for prior distributions. Furthermore, the biomass of small yellow croaker will decrease if the catch increases by 10%. Therefore, the annual TAC levels should be set below the maximum sustainable yield (21,301 tons) for e ﬀ ective small yellow croaker stock management.


Introduction
Fishery stock assessment is an important factor in maintaining the sustainability of fisheries' resources since fish is essential and plays a significant role in the human food supply. Various methods are used for the proper management of fisheries' resources. In addition, diverse analysis models are used to assess the reliable management scheme. However, it is requisite to manage fisheries' resources using an appropriate method with a reliable model. In Korea, because the fish stocks had been estimated to be decreased, a fish stock rebuilding plan (FSRP) has been implemented for many fish species such as sandfish, blue crab, octopus, cod, and filefish since 2005. The results of pilot projects showed that FSRP could contribute to increasing stocks of selected species [1].
The small yellow croaker (Larimichthys polyactis), belonging to the family Sciaenidae and order Perciformes, is mainly distributed in Northeast Asia, including the west coast of South Korea, Bohai Sea, and East China Sea [2]. The small yellow croaker is one of the most commercially important fish in South Korea and China, it has been exploited as a source of food and medicine for a long time [3]. However, despite the importance of the species, it does not have a direct catch regulation such as total allowable catch (TAC) or individual transferable quota (ITQ) and has not been analyzed by its biomass status by a Bayesian state-space model, one of the most powerful stock assessment techniques [4,5].
For analysis using the Bayesian state-space model, prior information on the estimated parameters is required, and care should be exercised when defining the parameters [28][29][30]. However, one of the major challenges in fish stock assessments is inadequate data, in addition to a general lack of prior information to facilitate prior distribution settings [31][32][33][34]. To address this challenge, strategies such as the use of a non-informative prior distribution have been considered [4,29,35]. However, there are debates on the appropriate prior distribution for Bayesian inferences [36]. Furthermore, maximum likelihood estimation is recommended instead of Bayesian inference when appropriate prior information is not available [37]. Therefore, a sensitivity analysis was carried out in the present study to obtain appropriate prior information for small yellow croaker stock analysis using the Bayesian state-space model, as suggested by Punt and Hilborn [36]. A prior distribution that would be effective for assessing small yellow croaker stocks was selected after comparing the results based on uniform distribution and inverse-gamma distribution, which represent informative and non-informative prior distributions, respectively, according to the findings of studies on small yellow croakers [22,35]. We believe this sensitivity analysis would help guide fishery scientists in the lack or absence of reliable prior information for conducting the Bayesian state-space model.
Furthermore, according to Maunder [38], stock assessment results obtained using the Schaefer [39] model could be unreliable. As the Bayesian state-space model suggested by Millar and Meyer [18] only considers the logistic model, we also evaluated the exponential model suggested by Fox [40], which is rarely applied for the Bayesian state-space model. The reference points for the management and recovery of small yellow croaker stocks were set based on the results of small yellow croaker stock assessments, and the appropriate TAC level for use in future stock management policies was explored.

Data
In this study, the available catch and fishing effort data from 1992 to 2018 were used for the small yellow croaker stock assessments. While Sim and Nam [15] used data based on two types of fisheries (stow-net and gill-net) for stock assessments, in this study, we used data from three types of fisheries, including pair-trawl fisheries, which accounted for the largest proportion of small yellow croaker catch in the 1990s. We used the catch and fishing effort (number of vessels) data from stow-net, gill-net, and pair-trawl fisheries, which accounted for approximately 85% of the total small yellow croaker catch from 1998 to 2018. However, in using the number of vessels as fishing effort data based on the three fisheries with varying fishing capabilities, these data had to be standardized into the same unit of fishing capacity [41,42]. Consequently, the fishing effort data of the three fisheries were standardized using the generalized linear model (GLM) [43].
Most of the small yellow croakers were caught by using stow-nets and pair-trawls before 2000; however, after 2000, the small yellow croakers have been mostly caught by using gill-nets ( Figure 1a). This is largely due to changes in fishing grounds and fishing methods across different fisheries, following the enforcement of the South Korea-China Fishery Agreement in 2001. The fishing methods and gears used in pair-trawl fisheries were changed from bottom trawls to middle and surface trawls in the mid-1990s, and the major target species changed from small yellow croakers to Spanish mackerel, chub mackerel, and anchovies [44,45]. In addition, the mesh size used in gill-net fishing has been continuously decreasing since the 1960s, and the catch capacity rapidly increased after the 2000s as the mesh size was reduced to 51 mm [46]. Consequently, the most common type of fishery for small yellow croakers changed from pair-trawl to gill-net fisheries.
The number of vessels associated with the various fisheries has continuously decreased since 1992 due to a vessel buyback program and other initiatives (Figure 1b). In 2018, the number of vessels associated with gill-nets was the highest (395), followed by stow-nets (207), and pair-trawls (70). Furthermore, the engine power (EP) of fishing vessels based on the fishing type has considerably increased from 196 EP in 1992 to 709 EP in 2018 for gill-net fisheries, 355 EP to 820 EP for stow-net fisheries, and 567 EP to 1340 EP for pair-trawl fisheries. The profit per vessel of stow-net fishery recorded at the lowest in 1998 and has increased until 2018 ( Figure 2). Similar to the stow-net fishery, the profit per vessel of gill-net fishery has been increased until 2018 since it recorded the lowest in 1993. However, the profit per vessel of pair-trawl has changed quite extremely since the major change in their target species.

Generalized Linear Model for CPUE Standardization
The GLM method proposed by Gavaris [43] has been extensively used to standardize the fishing capacities of different fisheries for the stock assessments of a single fish species [41,47,48]. The GLM is used to calculate a standardized catch per unit effort (CPUE) using the explanatory variables that influence the CPUE and can account for different fishing capacities across years and fishing method [41]. Previous studies have also used the year and fishing type as the explanatory variables in this model [15]. Lim et al. [49] identified fishing ground as another major variable because the major fishing grounds of small yellow croakers have continuously shifted northeast ( Figure 3).  [49]. Therefore, we selected and analyzed the year, fishing type, and shift in fishing grounds as the explanatory variables in our GLM, which could influence the CPUE of the small yellow croaker, as represented in Equation (1): where I denotes the CPUE and I R denotes the reference CPUE for each explanatory variable, namely the year, fishing type, and shift in fishing grounds. i denotes the explanatory variable, and j denotes the classification criterion in the explanatory variable. P ij denotes the relative fishing capacity of j in the explanatory variable i. Finally, X ij denotes a dummy variable, which is set to 1 if the data fall under the j category, and to 0, if otherwise. denotes a normal distribution where the mean is 0, and the variance is σ 2 . When Equation (1) is rearranged by applying log on both sides, Equation (2) is obtained: When Equation (2) is modified for CPUE standardization based on the fishery type for small yellow croakers, Equation (3) is obtained: where Year i denotes the difference in fishing capacity by year, for which the data from 1992 to 2018 were used. Gear j denotes the difference in fishing capacity by fishing type, including stow-net, gill-net, and pair-trawl. Move k denotes the difference in fishing capacity at the central fishing grounds, for which 1992-1999, 2000-2004, and 2005-2018 were defined as the dummy variables. Finally, the change in fishing capacity based on the shift in the central fishing ground by fishing type was considered as an interaction variable, denoted as Gear × Move jk . Coefficients for standardizing CPUE are inferred by linear regression following Equation (3), and then the standardized CPUE by years, fishing types, and main fishing grounds can be calculated. Standardized efforts can be calculated by catch divided by standardized CPUE of each fishing type. By adding standardized efforts of each fishing type, the total standardized catch efforts for the fisheries' stock assessment can be calculated.

Surplus Production Model
The surplus production model (SPM) can be used to assess fish stocks based on catch and fishing effort time series data. The SPM has been extensively used for the assessment of stocks of various fish species because it can perform stock assessments even with insufficient data [17,19,27,[50][51][52]. The SPM expresses the current biomass based on the relationship between past biomass and the growth and catch of the corresponding stock, as follows: where B y denotes the biomass of a year y, g B y denotes the stock growth, and C y denotes the catch of the year y. Furthermore, based on the g B y assumption, the SPM can be expressed using Equations (5) and (6).
where r denotes the intrinsic growth rate of the fish stocks, and K denotes the carrying capacity. Equation (5) was proposed as a logistic growth function by Schaefer [39], and Equation (6) as an exponential growth function by Fox [40].
In the estimation based on the SPM, the relationship between the catchability coefficient (q) and biomass is generally assumed to be constant [53]. In other words, the CPUE I y can be obtained by dividing the catch C y by the fishing effort E y and the relationship between the catchability coefficient (q) and biomass B y can be expressed using Equation (7).

Bayesian State-Space Model
To estimate the biomass of Korean small yellow croakers using the SPM, we used the Bayesian state-space model suggested by Meyer and Millar [22]. The biomass was reparameterized into a biomass-to-K ratio (P y = B y /K) for speed mixing. Consequently, the surplus production function can be expressed using Equations (8) and (9): Equation (8) was obtained by substituting Equation (5), which is the function of Schaefer (1954), in Equation (4), and Equation (9) was obtained by substituting Equation (6), which is the function of Fox (1970), in Equation (4). Furthermore, Equation (7), which is a function for the observation data, can be reformulated into Equation (10): The Bayesian state-space model can infer the posterior distribution of each parameter by combining prior distribution established from prior information and the likelihood derived from the observation data, according to the Bayes' theorem [54,55]. Furthermore, the parameters of interest required for fish stock assessment can be estimated based on the derived posterior distribution. The joint prior density for the parameters of interest can be expressed as Equation (11): Regarding prior distributions, informative prior distribution was assumed for internal growth (r) and environmental capacity (K), and non-informative prior distribution was assumed for fishing efficiency coefficient (q), according to Millar and Meyer [18]. To set informative prior distributions, we referred to Sim and Nam [15].
In addition, sensitivity analysis was performed for the initial values per the suggestion of Punt and Hilborn [36]. Typically, the r value for fish stocks is predicted to lie within the range of 0.015-1.5 [56]. However, K has a relatively broad range. Therefore, we assumed the r value to be 0.5 for sensitivity analysis and a non-informative condition for K. To assume a non-informative condition for K, the prior distributions of K were assumed to be uniform and to have inverse-gamma distributions [22,35,36]. Table 1 summarizes the prior distribution set for the estimation of Korean small yellow croaker stocks using the Bayesian state-space model. Based on the prior distributions, the likelihood of each parameter of interest can be derived using observation data with Equation (12): p I 1 , · · · , I n K, r, q, σ 2 , τ 2 , P 1 , · · · , P n = n y=1 p I y P y , q, τ 2 (12) According to the Bayes' theorem, the joint posterior distribution of each parameter of interest can be estimated by combining Equations (11) and (12) as follows: p K, r, q, σ 2 , τ 2 , P 1 , · · · , P n , I 1 , · · · , I n = p(K)p(r)p(q)p σ 2 p τ 2 p P 1 σ 2 × n−1 i=1 p P y+1 P y , K, r, σ 2 n y=1 p I y P y , q, τ 2  [53,54]. However, as it is challenging to generate a random sample from the joint probability distribution of multidimensional variables, the variables are sequentially extracted from conditional probability distributions of other variables using Gibbs sampling. Various algorithms are used to estimate the posterior distribution in MCMC simulations by using Gibbs sampling. WinBUGS program derives the posterior distribution by selecting an appropriate algorithm according to the distribution of each parameter of interest [57][58][59][60].
By performing the above process m times, the posterior distribution for each parameter of interest can be estimated using Equation (15):

Model Implementation and Comparison
Six scenarios were set up and analyzed by applying the assumptions of three prior distributions for the Schaefer and Fox functions. A total of 300,000 samples were collected and used to obtain estimates based on the Bayesian state-space model. Among them, the initial 50,000 samples that were not consistent with the posterior distribution were discarded by setting the burn-in period. Furthermore, the twenty fifth sample was extracted to minimize sample correlation that was observed during the extraction process. Finally, the posterior distribution was derived using 10,000 samples.
For comparing model fitness by scenario, we applied the deviance information criterion (DIC), which is useful in the comparison of models using Bayesian reasoning, as shown in Equation (16). The DIC can consider the fitness of the model (model complexity) and that of the analysis results simultaneously [61,62].
where L(data|θ) denotes the likelihood of the estimated parameter, D(θ) denotes the deviation, and p D denotes the posterior mean deviation. In particular, p D and D(θ) are used to evaluate the fitness of the model and that of the analysis results, respectively.

GLM Analysis Results and Model Comparison
The GLM results from 2000 to 2018 are shown in Table A1. The differences in fishing capacity based on the fishing type and shifts in fishing grounds were significant (P < 0.001). When the small yellow croaker stocks were analyzed using the Bayesian state-space model with the standardized CPUE, the Monte Carlo error for the parameters in all scenarios was lower than 5% of the standard error, this indicates convergence of the model [62]. In addition, as illustrated in Figure 4, the shapes of the posterior distributions estimated were not similar to those of the prior distributions in each scenario. Consequently, the data were considered to have sufficient information for model estimation [63]. Most graphs were skewed right, and therefore, the medians were used as the representative values.
The estimates of the model for each scenario have been outlined in Table 2. First, when the Schaefer function was used, the difference between the DICs based on the variation in the prior distributions was lower than two, which indicated no significant difference. In contrast, the Fox function was more suitable when the non-informative prior distribution was assumed than when the informative prior distribution was assumed, as there was a significant difference. The estimates of the model for each scenario have been outlined in Table 2. First, when the Schaefer function was used, the difference between the DICs based on the variation in the prior distributions was lower than two, which indicated no significant difference. In contrast, the Fox function was more suitable when the non-informative prior distribution was assumed than when the informative prior distribution was assumed, as there was a significant difference.  In comparing the two functions, regardless of the prior distribution assumption, the assumption of the Fox function was more suitable than that of the Schaefer function in all scenarios. Therefore, the most suitable analytical model was the Fox model with an inverse-gamma distribution for the K value. In addition, the R 2 value of the estimate was very high, at 0.99. As illustrated in Figure 5, all standardized CPUEs were included in the posterior 95% confidence level of the CPUE estimated using the Bayesian state-space model. Therefore, small yellow croaker stocks were assessed based on the estimates of the Bayesian state-space model, with an inverse-gamma distribution, using the Fox model, which has the lowest DIC. According to the small yellow croaker stock assessments, as illustrated in Figure 6, the fish biomass has decreased since 1992 and reached its lowest level of 16,800 tons in 2001. Subsequently, the fish stocks increased gradually, and the highest value was recorded in 2011; thereafter, it decreased again. The fish biomass was 69,150 tons in 2018, or 121% of the biomass at sustainable yield (B MSY ) (56,985 tons). However, since 1992, its biomass has exhibited considerable variations, and it had dropped below the B MSY between 1993 and 2004. Figure 6. Change in the annual stock biomass trajectories with 95% confidence intervals with biomass at sustainable yield (B MSY ) (a) and catch in comparison with the maximum sustainable yield (b) of small yellow croakers based on the Bayesian state-space model.
Since 2011, the biomass and catch of small yellow croakers have decreased continuously. Small yellow croaker catch has also continuously decreased since 1992, and the lowest value of 6380 tons was recorded in 2001. Since then, the catch has gradually recovered with an increase in biomass, and its highest value of 59,226 tons was recorded in 2011. Subsequently, the catch decreased again, and it was only 17,683 tons in 2018. Comparing it with the profits of gill-net and stow-net fisheries (Figure 2), it seems reasonable that they got relatively low profits in the 1990s when the biomass level of small yellow croaker in most years was lower than B MSY .

Analysis of Appropriate TAC Levels
To analyze the appropriate TAC levels that would facilitate the effective management and recovery of small yellow croaker stocks in the future, scenario analysis was performed based on the stock assessment results. Specifically, future biomass variations were estimated for each scenario under the following annual TAC level settings: a 20% increase (TAC 1.2 ), 10% increase (TAC 1.1 ), 10% decrease (TAC 0.9 ), 20% decrease (TAC 0.8 ), and 20,800 tons (TAC 1.0 ), which was the mean catch of small yellow croakers by gill-net, stow-net, and pair-trawl fisheries in 2014-2018, and the MSY (TAC MSY ) of 21,301 tons.
According to the analysis results, as shown in Figure 7, a 10% increase from the current level of catch could decrease the future biomass of small yellow croakers. However, when the TAC level was set to the current MSY level, their biomass was relatively stable.

Discussion
The results of the small yellow croaker stock assessments using the Bayesian state-space model revealed that it is necessary to manage the biomass of the small yellow croaker to ensure stable catch in the future. In particular, the estimates of future small yellow croaker biomass under different TAC scenarios suggested that their biomass would decrease with even a 10% increase in catch when compared with the current catch levels. Therefore, to prevent the reduction in small yellow croaker catch, as in the past, efforts should be made to maintain the small yellow croaker biomass at the current or MSY levels.
According to the sensitivity analysis results obtained using the Bayesian state-space model, the adoption of a non-informative prior distribution yielded results similar to or better than those observed with the prior information distribution. In other words, even if prior information on the stocks is insufficient, it is not necessary to discard the Bayesian state-space analysis, as demonstrated by Thorson and Cope [37]. Particularly, the estimates based on the Fox function demonstrated that the assumption of non-informative prior distribution could yield superior results when compared with the estimates under the assumption of informative prior distribution. Consequently, even if prior information exists for stock assessment activities based on the Bayesian state-space model in the future, performing an additional sensitivity analysis using non-informative prior information could facilitate more reliable stock assessment results.
Furthermore, in the sensitivity analysis, the DIC difference due to the change in growth function was greater than that due to the change in prior distribution. Specifically, the largest DIC difference based on prior distribution was 3.712 in the Fox model. However, the DIC difference based on the Schaefer and Fox function was 11.727. Therefore, a change in growth function could have a greater influence on the reliability of the estimates of the Bayesian state-space model than a change in prior distribution. Notably, when the Schaefer function was applied, the small yellow croaker biomass of 2018 was estimated at 85% of the B MSY , which was substantially different from the results obtained with the Fox function, where the biomass was greater than the B MSY . The sensitivity analysis for the Bayesian state-space model estimates to date has focused on prior distributions. However, in future Bayesian state-space model estimates, more reliable stock assessments could be obtained if growth function analyses are performed in addition to sensitivity analysis for prior distributions.
In the present study, the Pella-Tomlinson [64] function, which has been extensively used in fish stock assessments in addition to the Schaefer and Fox functions, could not be applied. This is because in the estimation with the Bayesian state-space model, the estimates of the Pella-Tomlinson function under an assumption of a similar prior distribution were not consistent with the parameter of interest. It could be more challenging for the Pella-Tomlinson function to achieve convergence with the model because variables are added based on parameter p in the analysis with the growth function model. However, more reliable results could be obtained in the future if the Pella-Tomlinson function can also be considered. Furthermore, other approaches should be adopted, such as model analyses based on age-structure models. In addition, there has been an increasing impact of climate change on fisheries' resources. However, in this study, only a "shift in fishing grounds" was considered as an environmental factor. In considering other variables, such as sea surface temperature, more reliable results could surface [65,66]. Therefore, there is a need for future studies to explore this variable, among others. Finally, we could not consider maximum economic yield (MEY) and ITQ, which is another catch regulation due to the lack of available economic data. If these are considered with the appropriate data, a more effective management plan for the small yellow croaker could be derived.

Conflicts of Interest:
The authors declare no conflict of interest.