Length data analyzed in this study were from bottom-trawl surveys (1960–2015) conducted by the SCS Fisheries Research Institute. The sampling stations (Figure 1) were predetermined before the surveys and consistent from year to year. Each station was investigated once and trawled for 1 h, with an average hauling speed of 3–4 knots in all surveys. The mesh size of the bottom-trawl nets ranged from 120 to 200 mm, with 30–40 mm cod-end mesh size in all the surveys. Surveys were conducted monthly in the 1960s and quarterly in other periods (Table 1). All captured fishery samples were identified to the species level, and biological data including length, weight, sexual maturity and stomach fullness for main commercial species were measured. The individuals were randomly sampled for measurement and laboratory bioassays. For each species, if fewer than 50 individuals were caught in a station, all were cryopreserved for laboratory bioassays; otherwise, 50 individuals were sampled randomly and measured. For each fish, the standard length was measured to the nearest millimeter.

Nine exploited fish stocks, including Japanese scad (Decapterus maruadsi), threadfin porgy (Evynnis cardinalis), yellowbelly threadfin bream (Nemipterus bathybius), golden threadfin bream (Nemipterus virgatus), red bigeye (Priacanthus macracanthus), purple-spotted bigeye (Priacanthus tayenus), brushtooth lizardfish (Saurida undosquamis), Japanese jack mackerel (Trachurus japonicus), and largehead hairtail (Trichiurus japonicus), were selected in this study regarding their high economic values and large catches in the northern SCS (Qiu et al., 2008; Zhang et al., 2017; Su et al., 2021). There were total 53, 652 length records of the nine fish species during different sampling years used in this study (total 30 assessment sequences in Table 1). Anal length was used for T. japonicas, fork length for D. maruadsi and T. japonicas, and body length for the other six fish species.

The growth of the fish stocks was modeled by the von Bertalanffy equation (von Bertalanffy, 1938):

where L_t is length (cm) at age t, L_inf the asymptotic length, K is the von Bertalanffy growth coefficient, and t₀ is the theoretical age at length zero. The ELEFAN I routines incorporated in the FiSAT II (Gayanilo and Pauly, 1997) software were used to fit growth curves to the restructured length-frequency data. Using both the “automatic search routine” and the “response surface analysis” within ELEFAN, it was possible to achieve the best fit for the growth curve (best-fitting combination of L_inf and K) to the length-frequency data.

The parameter t₀ were calculated using the empirical equation (Pauly, 1983):

Total mortality (Z) was estimated by the length-converted catch curve procedure (Pauly, 1983):

where N_i is the number of fish caught in a given length class i, t_i’ is the relative age corresponding to length class i, Δt_i is the time needed for growing through the length class i, and c is the intercept of the linear equation, respectively.

The instantaneous natural mortality (M) was calculated (Pauly, 1983) by:

where T is the mean environmental temperature. Fishing mortality (F) was calculated by subtracting M from Z, and the exploitation ratio (E) was obtained from F/Z.

Growth in body length is also assumed to follow the von Bertalanffy growth function (von Bertalanffy, 1938) in the LBB method (Froese et al., 2018).

Most of commercially exploited fish species grow throughout their lifetime, and their body size would approach the asymptotic length L_inf if mortality were zero, which can be expressed by:

where P_L/Linf is the probability to survive to length L/L_inf, which is solely a function of the M/K ratio.

The LBB method assumes that the selectivity of fishing gear is trawl-like, i.e., small individuals (length < L_x) can not be caught, all individuals will be caught if exceed a certain body size (length > L_start), and part of the individuals are caught when length between L_x and L_start. The gear selectivity can be expressed by the following equation:

where S_L is the fraction of individuals that are retained by the gear at length L, L_c is the length where 50% of the individuals are retained by the gear, and α represents the steepness of the ogive (Quinn and Deriso, 1999).

Combining the equations (1), (5), and (6), and rearranging lead to:

where N_Li and N_Li–1 are the numbers of individuals in length class L_i and the previous length class L_i–1, respectively. To minimize the required parameters, the ratios M/K and F/M are estimated, instead of the absolute values of F, M, and K in the LBB analysis. In other words, the increase in fish body length can be used as a proxy for its life time, and by using ratios instead of absolute values the units of time and biomass cancel out (Froese et al., 2018).

The Bayesian Gibbs sampler JAGS within R statistical language (version 4.0.3) was used to fit the observed proportions at-length to their expected values:

where p_Li is the observed proportions-at-length,p^Li is the mean values for p_Li, N^Li denotes the mean values for N_Li, which has been mentioned in equation (7).

The observed and predicted length distributions were then fitted by assuming Dirichlet-multinomial distribution (Thorson et al., 2017), which was proposed for fitting size and age composition in stock assessment models using a Bayesian framework. Proportions-at-length assume Dirichlet-multinomial distribution with an effective sample size of 1,000, which was chosen based on desirable performance across various simulation-testing trial scenarios (Froese et al., 2018).

The following equations are used to approximate the population status through the estimated quantities L_inf, L_c, M/K, and F/K. First, the length L_opt representing the maximum biomass of unexploited cohort is obtained from:

With a given fishing pressure F/M, the length at first capture L_{c_opt} that maximizes catch and biomass can be obtained from:

An index catch per unit of effort (CPUE’/R) is obtained as dividing relative yield-per-recruit (Y’/R) by F/M, which can be described as:

The relative biomass in the exploited phase of the fish population if no fishing takes place is given by:

where B₀’ > L_c denotes the exploitable fraction (> L_c) of the unfished biomass (B₀).

The ratio of fished to unfished biomass is described as:

A proxy for the relative biomass that can produce B_MSY/B₀ was obtained by re-running Equations (12–14) with F/M = 1 and L_c = L_{c_opt} (Froese et al., 2018).

Hordyk et al. (2019) indicated that the LBB analysis did not correct for the pile-up effect (pile-up of abundance observations in length classes used as bins in length-frequency analyses), and may result in a biased estimate of F and M/K. Therefore, we applied other two modified LBB model (Froese et al., 2019) on the length data of the 9 exploited fish species from Beibu Gulf. The two models, LBB-1 (full correction for the pile-up effect), and LBB-2 (let the Bayesian model determine the degree of correction based on the best fit to the available data) were based on the original LBB equation, and corrected for the pile-up effect.

In this study, we also analyzed the performance of LBB method with and without prior information. The prior information of parameters L_inf and Z/K were from the output of ELEFAN method. All the analysis was implemented using LBB_33a.R, an R-code algorithm presented by Froese et al. (2018, 2019). Fish stocks were classified to three exploitation statuses based on the estimates of B/B_MSY, overexploited status was assigned where B/B_MSY < 0.8, fully exploited status where 0.8 ≤ B/B_MSY ≤ 1.2, and underdeveloped status where B/B_MSY > 1.2 (Amorim et al., 2019). Besides, the stocks are considered as suffering from growth overfishing when the estimated L_c/L_{c_opt} < 1 (Liang et al., 2020; Zhang et al., 2021b).

The estimated asymptotic lengths for all assessment sequences ranged from 22.0 to 70.0 cm using ELEFAN method while ranged from 22.0 to 70.3 cm using LBB method. There were not significant differences in estimated L_inf between ELEFAN method and LBB method, in all assessment sequences (p > 0.05), except D. maruadsi stock of 2006 (t = 5.37, p < 0.05), and T. japonicas stock of 1982 (t = 3.12, p < 0.05). The estimated Z/K for all assessment sequences ranged from 3.36 to 11.19 using ELEFAN method while ranged from 2.3 to 12.0 using LBB method. There were not significant differences in estimated Z/K between ELEFAN method and LBB method in 11 assessment sequences (No. 1, 2, 4, 5, 7, 8, 10, 16, 21, 25, and 27, Table 2).

The estimated exploitation rates for all assessment sequences ranged from 0.34 to 0.87 using ELEFAN method while ranged from 0.26 to 0.86 using LBB method. The lowest value of exploitation rate occurred in N. virgatus stock of 1960 using both the two method. The highest value of exploitation rate occurred in T. japonicus stock of 1999, and D. maruadsi stock of 2010 for ELEFAN method and LBB method, respectively. There were significant differences in estimated E between ELEFAN method and LBB method in 16 assessment sequences (Table 2). Estimated exploitation rates for three assessment sequences, E. cardinalis stock of 1962 and N. virgatus stock of 1960 and 1962 were below 0.5 using LBB method, and only N. virgatus stock of 1960 were below 0.5 using LBB method (Table 2). Therefore, most of the fish stocks faced with overfishing during the assessment years.

There were significant differences in 20 assessment sequences (p < 0.05) in estimated L_inf using LBB method with and without prior information, and 10 assessment sequences (No. 5, 7, 12, 18, 20, 23, 24, 25, 26, and 28) were insensitive to the prior information (Table 3). As for estimated B/B_MSY, there were not significant differences using LBB method with and without prior information in only 2 assessment sequences (No. 6 and 17). In terms of exploitation statuses (B/B_MSY), 9 assessment sequences (No. 3, 4, 8, 13, 14, 16, 18, 23, and 24) showed different exploitation status using LBB method with and without prior information. For example, D. maruadsi stock of 2006 and 2007 were in overexploited status using LBB with prior information, but in fully exploited and underdeveloped status without prior information. Four assessment sequences (No. 12, 13, 15, and 17) showed growth overfishing were not happening using LBB method with prior information, and only 1 assessment sequence (No. 17) showed the same results without prior information (Table 3).

Reference point outputs (Table 4) showed that the three types of LBB methods (original LBB, LBB-1, and LBB-2) produced the same results when detecting the occurrence of growth overfishing, i.e., 86.7% of the 30 target stocks were facing growth overfishing (L_c/L_{c_opt} < 1), except for 4 assessment sequences (No. 12, 13, 15, and 17). The original LBB and LBB-2 models produced the same decisions if the stocks had been overfished (B/B_MSY). LBB-1 model produced similar results with the other two models, except for 5 assessment sequences (No. 3, 8, 18, 23, and 24). The underestimation of estimated B/B_MSY for the 5 assessment sequences made their exploitation status negatively, i.e., the original LBB and LBB-2 models indicated 4 stocks (No. 3, 8, 23, and 24) were fully exploited while LBB-1 showed they were overexploited; the original LBB and LBB-2 models showed N. virgatus stock of 1998 (No. 18) were underdeveloped while LBB-1 indicated it was fully exploited. The estimated B/B_MSY showed 83.4% of the 30 target stocks were in overexploited status or fully exploited status while E. cardinalis stock of 1999, N. bathybius of 1992 and 1997, N. virgatus of 1960, and P. tayenus of 1999 were in underdeveloped status (Table 4). In summary, the original LBB and LBB-2 models produced similar results, indicating Bayesian model can help determine the degree of correction.

There have been burgeoning literatures on the use of relatively simple methods to evaluate data-poor fisheries status. These approaches range from using life history characteristics as a guide to the vulnerability of fishing (Goodwin et al., 2006; Punt et al., 2011; McCully Phillips et al., 2015), to more holistic evaluations for obtaining management reference points and harvest control rules (Cope and Punt, 2009), or extensive assessment-based meta-analysis (Armelloni et al., 2021) using catch data, and auxiliary information (e.g., age-length data, life history parameters). Among the various data, fish growth in body length is discrepant by species and convenient to access.

Fish growth in body length of ELEFAN and LBB method both are assumed to follow von Bertalanffy equation (Pauly and David, 1981; Froese et al., 2018). However, the parameters in the two methods are estimated in different ways: the growth and mortality parameters (i.e., M, F, and K) in ELEFAN method can be directly calculated from empirical formulas or procedures while the ratios M/K and F/M are estimated, instead of the absolute values of F, M, and K in the LBB method. Compared with LBB method, ELEFAN method can provide the estimated growth and mortality parameters which are essential in full stock assessment or ecosystem based assessment, but relatively limited management reference points. Gulland (1983) recommend that 0.5 may be the suitable exploitation rate for fish stocks in temperate water. However, fish stocks in tropical and subtropical areas (e.g., the Beibu Gulf) have short life-cycles and rapid growth, and they can sustain high exploitation rates (Wang et al., 2012). Therefore, the exploitation rate may be insufficient for the fishery management in Beibu Gulf while L_c/L_{c_opt} and B/B_MSY provided by LBB method are commonly used management reference points (Zhang et al., 2021b). Our study showed the prior information from ELEFAN method was effective for LBB method (Supplementary Figure 1) because most of the assessment sequences produced significantly different results with and without prior information.

The parameters of LBB method are calculated by a Bayesian Monte Carlo Markov Chain (MCMC) approach (Cowles and Carlin, 1996), which has been widely used in fishery data analysis (Haddon, 2010). With this Bayesian framework, it is straightforward to calculate credible intervals for multiple parameters. In this study, 95% credible intervals of B/B_MSY and L_c/L_{c_opt} were calculated, and these results can provide alternative information sources which can support decision-making. Besides, we also applied two modified LBB models (LBB-1 and LBB-2) with corrections for the pile-up effect (Froese et al., 2019; Hordyk et al., 2019). The results (Table 4) showed estimations of the original LBB method have been little affected by the pile-up effect for most of the assessment sequences.

In addition, we suggest using both ELEFAN and LBB methods to fit length-frequency data of data-poor fish stocks because they are complementary in estimating management reference points.

The Beibu Gulf has multiple ecosystems with estuaries, mangroves, coral reefs, and shelves, which provide comfortable habitats for spawning, feeding, and nursery areas for abundant fish species (Qiu et al., 2008). As reported to date, 960 fish species inhabit this embayment, which belonging to 162 fish families. The main fishing gears are trawl, purse seine, gill net, hook and set nets, and the trawl fishery accounts for more than 70% of the total catch (Zou et al., 2013). There are not available statistical catch data in the Beibu Gulf because catch data are gathered by administrative districts (e.g., provinces in China), instead of sea areas. Therefore, it was considered to be a data-poor fishing ground, and the length-based methods used in this study may fill a gap in knowledge of biomass levels and exploitation status of main fish stocks in the Beibu Gulf.

Zou et al. (2013) has estimated the fishery catch of the Beibu Gulf to be 85.7 × 10⁴ t in 2012, including 65.7 × 10⁴ t caught by China, and 20.0 × 10⁴ t caught by Vietnam. The numbers of motorized fishing boats and total fishing powers of Guangxi Province were selected to represent the long-term trend of fishing efforts in the Beibu Gulf (Figure 2). Before China’s reform and opening up in 1978, the number of motorized fishing boats in Guangxi Province was less than 1200, and total fishing powers was no more than 8.5 × 10⁴ kW. In 1980s and 1990s, the fishing efforts had been rapidly increasing and some indications of overfishing have appeared, e.g., catch rates trend downward (Qiu et al., 2008), catch composition has changed significantly (Chen et al., 2011), and miniaturization, early sexual maturity and accelerated growth have occurred in main commercial fish species (Zhang et al., 2020a). The indications of overfishing were consistent with our assessment results, which showed the assessed fish stocks had been overexploited in the 1990s.

In the recent two decades, a series of conservative management measures have been implemented by Chinese government in the SCS, including “double control” system, summer fishing moratorium, “zero-growth” and “negative-growth” strategies (Shen and Heino, 2014; Cao et al., 2017). Benefiting from the management measures, the fishing efforts had been stabilized and then decreased (Figure 2). Recent studies have shown that the management measures to reduce fishing pressure had been playing an important role for fishery resources recovery, e.g., the average daily yields of fishing boats and output values have increased after the summer fishing moratorium (Su et al., 2019), and these measures have a positive influence on the biological characteristics of this commercial fish species (Zhang et al., 2021a). However, our results showed the main commercial fish stocks were still in overexploited status in recent years (Table 4). Beibu Gulf is a co-developed area of fishery resources by both China and Vietnam, but until now, only China makes policy efforts to reduce fishing pressure. For example, during the 3.5-months summer fishing moratorium, all fishing gears (except for rod fishing) of China stopped fishing in the Beibu Gulf, but the fishing boats of Vietnam were still active in this sea area. Therefore, we emphasize collaboration mechanism should be established by the two countries for management and sustainability of fishery resources in the Beibu Gulf.