Recruitment-driven fish production in two regions where fish biomass has drastically declined

Catches have remained relatively high in the Gulf of Thailand and the Bohai Sea, China, despite severe biomass declines (around 95%) evidenced by fishery-independent surveys. Such high production at very low stock sizes is not predicted by simple-surplus production theory, but can be explained by age-structured models that predict high recruitment rates even when biomass per recruit (BPR) has been drastically reduced. Re-cruitment rates can be reconstructed by estimating changes over time in biomass and BPR, for alternative assumptions about survey catchability, growth, and mortality rates. These reconstructions indicate that likely severe decreases in BPR, due to high fishing mortality rates, imply that total recruitment rates have likely been relatively stable over time, with catch consisting largely of new recruits making up most of the low biomass. These considerations may apply to numerous areas in east and southeast Asia where most of the catch of coastal trawlers is used to produce animal feeds, notably for aquaculture.


Introduction
Globally, wild-caught seafood provides high quality and quantity of animal protein to people's diets and is more available and affordable than other sources of animal protein (FAO, 2016;Mohanty et al., 2019).Hence, in the past decades, indiscriminate and intense fishing has occurred in many marine ecosystems around the world (Anticamara et al., 2011).After a long history of overexploitation, major fish stocks, especially large predatory ones, have witnessed widespread decline and collapse (Myers and Worm, 2003;Christensen et al., 2014;Meissa and Gascuel, 2015;Palomares et al., 2020).This overfishing has, in many east and southeast Asian seas, led to a "miniaturization" of the fish caught (Yang et al., 2022), which are considered too small to be marketed for human consumption, but sought after as animal feed, notably for the aquaculture sector (Ahmed et al., 2007;Cao et al., 2015).
Classical surplus production models have typically assumed maximum sustainable yields to occur at biomasses between 30% and 50% of natural unfished biomass (Schaefer, 1954(Schaefer, , 1957)).But there are areas in east and southeast Asia where survey data indicate biomass declines of 90% or more, while surplus production as measured by catch has remained near historical peaks.Although the accuracy of catch reporting in these areas has been questioned, a noticeable change in catch composition occurred.In contrary to the past, when large, long-lived species were abundant, catches in such areas are now dominated by short-lived species with very high growth and natural mortality rates, and very small, newly-recruited individuals make up the majority of the catch (Chotiyaputta et al., 2002;Qiu et al., 2010;Cao et al., 2015;Szuwalski et al., 2017).This hints at the possibility that there was a compensatory increase in recruitment rate per spawner at the community level.
Fisheries scientists often assume that recruitment dynamics remain stable, and recruitment is positively related to the spawning biomass.However, the mismatch between spawning biomass and recruitment has been reported by several studies (Myers, 2001).By analysing 224 stocks in the RAM Legacy Stock Assessment Database, Szuwalski et al. (2015) revealed that more than half of these stocks do not show a positive relationship between recruitment and the biomass of spawning adults over the observed range of stock sizes, among which 85% exhibit shifts in average recruitment.Multiple factors, including density-dependent changes in survival and environmental conditions, can lead to changes in pre-recruit mortality rate (Szuwalski et al., 2015;Hilborn and Walters, 2021).Pre-recruit survival rates at low stock sizes are often seen to increase by factors of 20 or more when compared with that at unfished stock sizes, which leads to improvement in the recruits per spawner as the abundance of spawning stock declines, and sets the common ecological basis for high fishery production despite the decline of community biomass (Barrowman and Myers, 1996;Hilborn and Walters, 2021).
The pattern of very high sustained catch at very low community biomass levels has been particularly striking in two areas (Gulf of Thailand, Bohai Sea, China), where long-term fisheries-independent survey data can be used to document the severity of the biomass decline.Due to the absence of data on maturity, these analyses work with total biomass instead of spawning stock biomass.Here, we use simple population dynamics theory to estimate changes in biomass from the survey data and changes in biomass per recruit (BPR) from estimates of growth, natural mortality, and fishing mortality rates.Dividing biomass by BPR provides a direct estimate of recruitment rate (Pauly, 1982), and we show for the two regions that likely severe decreases in BPR, due to high fishing mortality rates, imply that total (over species contributing to the catch) recruitment rates have likely been relatively stable over time, i.e. the ratio of biomass to BPR has remained high despite very low biomasses.We discuss possible mechanisms for this surprising stability, ranging from very strong recruitment compensation to size-area refuges for egg production to decreases in pre-recruit mortality rates as predator abundances have declined.

Methods
We apply calculations usually done for single species assessments to total community catches and biomasses.The total community biomass consists of the sum of species specific recruitments times species specific BPRs.This sum can be expressed as the total community recruitment (which we estimate) times the sum of products of proportions of total recruitment by species times species BPRs.Thus, this sum of products is a recruitment-weighted community average BPR, and we assume that it can be estimated from community average growth, mortality, and fishing estimated from proportional contributions of species (or functional groups with similar parameters) to the total catch.We recognize that there can be error in this average calculation, but not enough information is available to make the calculations on a species specific basis, then add up the resulting totals.

Data source
Catch  and survey catch-per-unit-effort (CPUE;1959, 1982, 1992, 1998, 2010) time series data for the Bohai Sea came from China Fishery Statistical Yearbook (Ministry of Agriculture, 1979Agriculture, -2020) ) and Jin (2004Jin ( , 2014)), respectively.The fisheries catch  and survey CPUE (1961CPUE ( -2014) ) data for the Gulf of Thailand were originated from Department of Fisheries (2015).Given that these catches came from multiple species with widely varying natural mortality (M) and growth parameters, we used an average M, and corresponding von Bertalanffy growth parameter K value in the calculations, obtained as a catch-weighted average over catches by functional grouping (e.g.small pelagic fish; large demersal fish) for which M estimates were assembled by Christensen et al. (2009Christensen et al. ( , 2015) ) from FishBase (www.fishbase.org),and catch composition data from the Sea Around Us (Pauly and Zeller, 2016;www.seaaroundus.org)catch reconstructions.

Estimation of biomass from survey data given assumed fishing mortality rate
Suppose that for some recent reference period when near-peak catches were being taken, we simply assume a reference fishing mortality rate F ref .Then we can estimate average biomass for that period from the fact that catch is biomass times fishing rate, i.e.
If we have a mean survey catch rate I ref for that period, we can calculate a survey catchability coefficient q for the model I = q•B (index I proportional to biomass B) as simply q = I ref /B ref .Using this q, we can then "invert" the index time series I t to obtain biomass estimates B t for all years when surveys were conducted, as B t = I t /q.We can then obtain fishing mortality rates for all such years, as Estimation of BPR and recruitment from surveys, catches, and fishing mortality rate One of the most basic calculations or predictions in fish population dynamics is biomass-per-recruit, or BPR.At equilibrium, we obtain this prediction by summing or integrating the product of survivorship to age (probability of surviving from recruitment to the age) times predicted weight at age, over ages from the time of recruitment until all fish have died.To do this calculation, we require an assumed fishing mortality pattern over age, an assumed natural mortality rate M, and a body growth curve (Beverton and Holt, 1957).Under the simplest or "knife-edge" assumption that fish become vulnerable at some minimum age or size than are equally vulnerable for all older ages, and the von Bertalanffy growth curve, one can derive expressions for BPR as a function of F, M, and the growth model parameters.For example, as noted by Pauly (1980) the classical Beverton and Holt (1957) yield per recruit model predicts BPR from the equation: where Z is the total mortality rate (Z = F + M), with F the fishing and M the natural mortality, W ∞ is the asymptotic body weight (assumed proportional to the cube of asymptotic length), K is the parameter expressing how fast the asymptotic body weight approaches, and τ is the age at recruitment, and all rates are per year.A simpler prediction for the continuous time delay-differential model (Walters, 2020) is where W r is the body weight at the age of fishery recruitment and K is a weight-growth parameter typically near 0.5 K.Note that BPR is likely to vary over time not only with changes in fishing rate F, but also with recruitment variation, changes in M and growth rate, and departures from equilibrium during periods of rapid fishery development.In particular, during periods of increasing fishing mortality rate, older animals are typically more abundant (and BPR is hence larger) than predicted by equilibrium equations such as (1) and ( 2).
BPR predictions from models like Equations ( 1) and ( 2) are typically very sensitive to changes in fishing rate F, with large declines in BPR with even modest increases in fishing rate.This sensitivity is present even when the basic biological turnover rate of the stock (as expressed by M and K or K ) is very high [e.g.M and K often exceed 2 year −1 for small creatures like shrimp (Pauly, 1982)].As a rough rule, if L r /L ∞ (the ratio between body length at the age of recruitment and asymptotic length) is in the range 0.3-0.5 (implying probable growth overfishing) and F/M is near 2, BPR is predicted to drop by 90%.
Dynamic models can be used to reconstruct changes in BPR with estimated past changes in recruitment and growthmortality parameters, but a simpler approach is just to compare predicted BPR for two or more reference periods during which fishing rates and BPR have likely been stable.For example, we can compare average BPR predicted for a nearunfished situation to BPR for the most recent decade over 1045 which surveys indicate low but stable total biomass.For any such period, average estimated biomass B t (based on assumed F for the period and survey index I t ) has to have been B t = R t •BPR t , i.e. recruitment rate R t times BPR, and from this we can "back-calculate" what recruitment rate R t must have been (see, e.g.Pauly, 1982), i.e. just R t = B t /BPR t . (3) Note that this estimate depends on the assumed peakperiod fishing mortality rate F t used to estimate B t , as well as the period and catch-dependent For cases where annual B t estimates (conditional on assumed F ref ) are available for successive years, annual recruitments can be approximated using the delay-differential population dynamics model: with Z t = M+(C t /B t ).This approximation approaches the simpler prediction from Equation (3) when B t+1 is close to B t and depends on vulnerable numbers of fish remaining near equilibrium with respect to the recruitment rate (which typically occurs with time simulations of vulnerable numbers for developing fisheries using the delay-differential model).
Change in recruitment between the unfished and reference periods can be expressed (conditional on F ref and the growthnatural mortality parameters) in terms of the recruitment ratio R ref /R 0 , where R ref is mean estimated recruitment over the reference period and R 0 is estimated unfished recruitment.This ratio can be predicted just from the estimated stock depletion (5) That is, estimated relative recruitment is negatively related to depletion D t but positively related to the relative change in BPR as measured by the ratio BPR 0 /BPR ref .

Estimation of recruitment compensation
Recruitment compensation refers to the increase in recruitment per unit spawning biomass as spawning biomass decreases.A simple measure of compensation is the Goodyear compensation ratio (CR) (Goodyear, 1977), which is the ratio of the maximum recruitment per biomass at low stock size to the recruitment per biomass for the unfished population.Comparative studies of this ratio using stock-recruitment data for temperate fish stocks have found widely varying values, from as low as 3-5 for species like Pacific salmon that have relatively low fecundity to 50 or more for highly fecund demersal species (Goodwin et al., 2006;Dorner et al., 2008).The average CR across species with relatively high fecundity has been found to be around 20-25.
For the two-period comparison described above, a minimum estimate of the recruitment CR would be This minimum value is based on assuming that recruitment would decline linearly towards zero if stock size were reduced to below B ref , i.e. all of the "compensatory reserve" has been used.A higher value of CR would result from assuming that recruitment has followed a Beverton-Holt relationship of the form R = aB/(1 + bB).For this form, the CR implied by any Such relatively optimistic estimates of CR can also be obtained for other stock-recruitment equations, such as the Ricker model.Note that CRBH from Equation ( 7) can even take negative values, for parameter scenarios that would require R ref > R 0 , i.e. increasing recruitment over time.In case there is interest in predicting future recruitment using the Beverton-Holt recruitment model, the model's parameters are given by a = CRBH/BPR 0 and b = (CRBH-1)/B 0 .

Results
Figure 1 presents a time series of the catch and survey biomass from the Bohai Sea and the Gulf of Thailand.Catches in 1950s were relatively low in both areas, and gradually increased since.Catch in the Bohai Sea experienced a rapid increase from 1980s to 1990s, and reached a peak by the end of 1990s, from which it declined in the 2000s.Catch in the Gulf of Thailand strongly increased starting in 1960s and reached a peak around 1990, then declined somewhat and remained stable thereafter.The survey index data (Figure 1) indicate that biomass declined dramatically with the increase of catch, and current biomass depletion (D ref for years 1998 and later) is very similar for the Gulf of Thailand and the Bohai Sea, i.e. 0.076 for the Gulf of Thailand and 0.074 for the Bohai Sea.That is, current total harvestable biomass is only ∼7.5% of what it was before rapid fishery development occurred.
The interpretation of how much BPR and recruitment has changed in both cases is strongly dependent on what assumptions we make about recent fishing mortality rate F ref and average natural mortality rate M (and corresponding growth parameters K or K and W ∞ ).Low assumed F ref leads to much lower estimates of R ref /R 0 from Equation (5) and minimum CR from Equation ( 6), but the basic shape of the predicted pattern of change with increasing F ref is virtually independent of assumed M (Figure 2).For the M = 0.5 year −1 scenario in Figure 2, a current fishing mortality rate (F ref ) of around 1.75 year −1 would have to be assumed in order to conclude that recruitment has not declined over time (R ref /R 0 > 1.0), and F ref would need to be around 3.5 year −1 in order for estimated recruitment not to have declined if M = 1 year −1 (Figure 2).Put simply, in order to explain the severe observed depletions purely as an effect of F on BPR rather than recruitment, we would need to assume a recent F near 3.5 year −1 if M has averaged 1.0 year −1 , and around 1.75 year −1 even if M has averaged only 0.5 year −1 .
Catch-weighted average M values for the main functional fish groups caught in both areas suggest that an average M near 1.0 year −1 is more realistic (Figure 3).In fact, average M has very likely increased over time in both regions due to shifts in catch composition from larger, longer-lived species towards much smaller, short-lived species that can withstand higher fishing mortality rates (Figure 3).Decreases in mean trophic level of the Bohai Sea catches (Liang and Pauly, 2020) suggest a similar trend for M. Note that M was assumed to be constant for each functional group, but the smaller individuals from more long-lived functional groups now being caught may well have higher M values (Pauly et al., 2001).In any case, it does not seem reasonable to assume average M for  recent years <1.0 year −1 for either of the regions.Note that using M = 1 year −1 to calculate unfished BPR, when mean M was likely closer to 0.8 year −1 before fishery development, has likely caused an underestimate of unfished biomass per recruit (BPR 0 ) and hence an overestimate of early recruitment, of ∼25%.
Assuming that BPR remained near equilibrium over time for the two regions, changes over time in relative recruitment rates are shown in Figure 4 for two extreme F ref assumptions (low, 1.5 and high, 3.5) and with time-varying M esti-mates from Figure 3.For the Beverton-Holt BPR calculation (1) each year, the M for that year was multiplied by 0.8 to give an estimate of K for the year, and W ∞ for the year was assumed to be proportional to K −2/3 based on the auximetric relationship of Pauly (1980).To obtain annual estimates for every year, biomasses for non-surveyed years were interpolated linearly between estimates for surveyed years and were assumed stable at the first surveyed biomass for years before the first survey.For the Bohai Sea, calculated recruitments decrease as the fishery developed (and biomass declined) under the low F ref assumption, but are stable or even increasing under the high F ref assumption (Figure 4a and b).For the Gulf of Thailand, increases in mean M lead to increases in predicted recruitment during the fishery development, though less pronounced for the low F ref assumption (Figure 4c and d).For both cases, recruitment is estimated to have remained relatively high over the whole fishery development under the high F ref assumption, i.e. there is no clear evidence of severe recruitment overfishing for the overall biomass production system.
The trends of number of fishing vessels and fishing power/total standardized fishing effort in the Bohai Sea and Gulf of Thailand implied that fishing mortality increased in both areas (Figure 5).Estimates of trawl swept area do suggest very high potential reference F for both areas (Table 1).The trawl characteristics (speed, net size) reported by Kongprom et al. (2003) imply that the swept area ratio (the ratio of area swept by gear to area of the fishing ground) could easily exceed 10 year −1 in the Gulf of Thailand, suggesting F near 3 even if the trawls have a lower capture efficiency than assumed by Pauly (1980) and Kongprom et al. (2003).This is also the case in the Bohai Sea.Based on the trawl characteristics reported by Jin (2004), the swept area ratio in the Bohai Sea could be around 35 year −1 , which indicates that assuming high F ref for recent years is credible.

Discussion
Is it plausible to assume very high recent fishing mortality rates for the two regions, e.g.rates exceeding 3.0 per year that would imply stable or increasing recruitments, but very low current BPR?In both areas, the size distribution of fish has shifted dramatically towards smaller individuals (Figure 3b; Supongpan, 2001;Jin, 2004), which does indicate a severe drop in BPR, and imply that much of the observed biomass depletion has been due to decreases in BPR rather than decreases in recruitment alone, no matter what we assume about changes in mean M for the fished community.Recruitment    4) are not that sensitive to alternative assumptions about the M/K ratio, but are quite sensitive to predicted changes with K and the asymptotic weight W ∞ .The inverse relationship between K and W ∞ implies a substantial decrease in BPR as mean M and K increased over time (and mean W ∞ decreased), considerably more severe than would be expected from the predicted changes in BPR with F shown in Figure 2 for constant K.
An F of 3.0 year −1 or higher may seem very extreme to scientists working on temperate fish stocks with relatively low Ms, but in fact such high Fs just mean that the biomass is being turned over (by fishing) about three times per year or once per quarter.If fishing mortality rates are indeed as high as we suspect, then since average catch C = F•R•BPR, very low BPR implies that catch C is driven mainly by the recruitment rate R, or in biomass terms, by the input rate R•W r to the biomass pools.Importantly, the minimum recruitment CR estimates needed to explain current catches without any decline in recruitment (but at very high current F ref ) are quite reasonable (Figure 2), on the order of 10-20 and hence below the average estimated for temperate species using stockrecruitment data.Only if we assume particular recruitment equations like the Beverton-Holt would we estimate unreasonably high CRs for reference period recruitments R ref approaching the unfished R 0 estimates.Similar patterns of severe decline in survey biomass and apparently very high ratios of fishing to the natural mortality rate have been estimated for the Gulf of Thailand and two other areas in southeast Asia, the Philippines, and Malaysia (Stobutzki et al., 2006).Such observations suggest that recruitment-driven fisheries based on very small fish may be the general rule in east and southeast Asia as a whole.
At least some of the larger species still being caught in both regions are very likely recruiting to at least trawl gear at sizes well below their size at maturity.This means that their spawning biomass per recruit (SBR) has likely declined even more severely than indicated by BPR changes, and their apparent recruitment CRs must be even higher than indicated by R/B ratios if these species have indeed been subject to knife-edge vulnerability to the fishing gear.
There are at least four possible factors that could explain why high recruitments have been maintained despite very severe biomass declines: (1) Larger species may have dome-shaped vulnerability, with a size "refuge" for larger individuals due to ability to escape at least trawl gear and/or ontogenetic movement into areas with relatively low fishing effort (e.g. to offshore, deeper water).(2) There may be spatial refugia where biomass density remains high despite widespread depletion that act as recruitment sources (e.g.water too shallow for trawlers, areas surrounding dangerous obstructions to trawling like sunken vessels and reefs, far offshore areas, etc.); besides, unaccounted biomass (e.g.seasonally migrating stocks with some of their life history located in other areas) might also contribute to recruitment.(3) Recruitment CRs are likely to be higher in small, lower trophic level species (often r-strategists), and thus the community-averaged CR are expected to increase under depletion of larger species.(4) Natural mortality on recruits, particularly for smaller species, may have decreased considerably due to reduced predation mortality as larger fish have become less abundant, as suggested by Pauly (1982, Figure 4), demonstrated by Walters et al. (2008, Figure 3) and predicted by the size-spectrum model used by Szuwalski et al. (2017).
The factors involving size and spatial refugia cannot be assumed to be stable or persistent over time, since they involve technological and economic factors that may be subject to rapid change.Thus, it should not be assumed for management planning that recruitment rates will "take care of themselves", and research efforts should be initiated to determine which of the factors is/are most important.
Maintenance of high recruitment rates for the fish communities as a whole does not imply that there has not been recruitment overfishing, especially of larger species, but also smaller ones.For example, Pauly (1980) presented evidence of recruitment overfishing for Lactarius lactarius (a species with a relatively high K value and a small maximum body size) in the Gulf of Thailand, with its stock recruitment reconstruction indicating a recruitment CR of only ∼6.0 for that species (below the CR needed for persistence given F ref near 3.5 in Figure 2, M = 1 case), and the stock did collapse.Other species, such as rays and sawfish, have disappeared entirely from the catch data for the two regions, providing clear evidence of recruitment overfishing.
This contribution presents a novel BPR analysis applied at the community level, rather than the typical stock level, to analyse changes in community recruitment over time.This offers a new perspective on the system when detailed species specific information is missing.However, this approach pays less attention to species shifts and biodiversity loss.In fact, large species composition changes of fishery resources occurred in both areas.These changes consisted of the virtual disappearance of very large, long-lived fish, a massive decrease of some of the previously abundant groups, such as croakers (Sciaenidae) and slipmouths (Leiognathidae), and a marked increase of snappers (Lutjanidae) and squid (Pauly, 1988).These two regions also show different developments in catches over time.Compared with the Gulf of Thailand, the total catch of large species in the Bohai Sea exhibited increasing trend, and the proportion of small fish in the Gulf of Thailand has increased significantly.The difference was likely caused by differences in temporal development of F and hence on impact on the whole suite of species.
Species disappearances and other biodiversity issues are usually of little concern to companies deploying fleets of trawlers whose catch is to be mainly converted to animal feed.However, this is different for governments, which could also consider the loss of proteins and micronutrients caused by used "miniaturized" fish as animal feed (Hick et al., 2019) and which could instead be allowed to grow and then be caught, allowing them to be directly consumed by people.

R
ref , B ref , and relative depletion D ref estimate is given by CRBH

Figure 1 .
Figure 1.Trends in total catch and survey CPUE in two heavily exploited water bodies in northeast Asia (the Bohai Sea) and southeast Asia (the Gulf of Thailand).

Figure 2 .
Figure 2. Predicted effect of alternative assumptions about current fishing rate F ref and natural mortality rate M on predicted relative BPR, recruitment, and estimated CRs.Note difference in X-axis scale, with higher assumed M resulting in higher F ref needed to cause severe changes in BPR and unrealistically high estimates of the CR-CRBH [extreme CRBH values (higher than 50 or lower than 0) are not shown in Figure 2].

Figure 3 .
Figure 3. Trends in species (a) and size (b) composition of the catch and catch-weighted mean natural mortality rate M over major functional groups (c) in the two regions.Note that the catches used here for the Bohai Sea include catches from both the Bohai Sea and Yellow Sea, since the two seas belong to the same large marine ecosystem and report catch composition together.

Figure 4 .
Figure 4. Predicted trends in relative recruitment and its relationship with stock biomass for the two regions, under two extreme assumptions about average fishing mortality rates after the year 2000: F ref = 1.5 year −1 is likely below the fishing rate actually experienced, while F ref = 3.5 year −1 could be an overestimate.The top row displays regimes in recruitment and stock biomass (solid line: recruitment; dotted line: biomass; coloured boxes: recruitment changes in different biomass regimes; green box: less exploited period; pink box: transition period; purple box: reference period), where box height represents the within regime, among-year variance, and length represents regime length.The bottom row demonstrates the relationship between recruitment and stock biomass, and the dots are coloured to match recruitment regime.

Figure 5 .
Figure 5. Trends of number of fishing vessels and fishing power/total standardized fishing effort in the two regions.

Table 1 .
Swept area ratio calculation in the Bohai Sea and Gulf of Thailand.