Evaluation of traditional and bootstrapped methods for assessing data-poor fisheries: a case study on tropical seabob shrimp (Xiphopenaeus kroyeri) with an improved length-based mortality estimation method

The improved length-converted catch curve (iLCCC)

The standard method to estimate total mortality (Z) from size-frequency data is the length-converted catch curve (LCCC, Baranov, 1918; Pauly, 1983; Pauly, 1984a; Pauly, 1984b). This method transforms length classes into relative ages according to input growth parameters. A considerable limitation of the LCCC is its inapplicability when the asymptotic length (L_∞) is lower than the maximum size in the sample. This limitation essentially excludes populations with lower mortality/growth (Z/K) ratios, as these typically exhibit higher frequencies of large-sized individuals due to lower relative harvest pressure (Schwamborn, 2018). To address this issue, we developed a new algorithm, the iLCCC (“improved length-converted catch-curve”, de Barros, 2022), which automatically excludes size classes that are larger than the estimated L_∞ and performs bootstrapping to propagate uncertainty from growth parameters into the resulting mortality estimate. The iLCCC is based on the following framework:

First, we assume numbers-at-age (N_t) decrease over time according to an instantaneous rate of total mortality (Z): (1)${\displaystyle N_t=N_{t-1}{e}^{-Z_t}}$ (2)${\displaystyle ln\left[N_t\right]=ln\left[N_{t-1}\right]-Z_t}$

We then transform the length composition to relative ages (a_L) using growth parameter estimates from the von Bertalanffy growth function (Essington, Kitchell & Walters , 2001; Von Bertalanffy, 1957) as follows: (3)${\displaystyle a_L=t_0-\left(\frac{1}{K}\right)ln\left[1-\frac{L}{L_{\infty }}\right]}$ where a_L is the relative age (the estimated age for the size class L), t ₀ (year⁻¹) is the theoretical age when length equals zero, which determines the horizontal position of the growth curve (Kirkwood, 1983; Von Bertalanffy, 1957), K (year⁻¹) is the instantaneous growth rate, and L_∞(mm) is the theoretical asymptotic length. The instantaneous rate of mortality (Z, year⁻¹) is then estimated as the slope of a simple linear regression fitted to the decline in numbers-at-age through time, using relative age (a_L) as a proxy (Eqs. (1) and (2)).

The major advantage of the iLCCC lies in the use of bootstrapping to resample relative age estimates based on the uncertainty in growth parameters. Specifically, we employ a “two-step” bootstrap (Schwamborn & Schwamborn, 2021), where Z is calculated multiple times by iterating through randomly generated samples of K, L_∞, and t₀, as well as by resampling the input regression data. Despite its relative simplicity and similarity to the traditional length-converted catch-curve, we argue this method represents a substantial improvement because it essentially propagates the uncertainty in life-history traits to the resulting mortality estimate. This is often not incorporated in length-based assessments, possibly resulting in underestimates in the underlying uncertainty in mortality rates.

Simulation study

We first conducted a simulation study to (1) determine the reliability of the iLCCC method in estimating total mortality rates from length composition data, and (2) compare growth and mortality estimates between the fishboot package and the traditional PW approach, using known “true” parameter values as baseline. The synthetic population was simulated using the “fishdynr” package (Taylor & Mildenberger, 2017) with pre-determined growth and mortality rates (Table 1). A year of length frequency data was extracted from the model. This data was then used to estimate growth and mortality rates using bootstrapped ELEFAN and the PW method, as detailed in the following sections.

Real world scenario—study area and sampling

The seabob shrimp Xiphopenaeus kroyeri is a widely distributed penaeid inhabiting shallow coastal waters in the tropical and subtropical Western Atlantic. It is a highly popular target species for artisanal small-scale fisheries (Musiello-Fernandes, Zappes & Hostim-Silva, 2018; de Barros et al., 2022). Shrimp were sampled monthly from October 2018 to September 2019 at coastal sampling stations of four important shrimp fishing grounds in the Eastern Brazilian Marine Ecoregion (sensu Spalding et al., 2007) (Fig. 1). These areas included: Caravelas (CV), São Mateus (SM), Ipiranga (IP), and Rio Doce (RD). The Caravelas estuary encompasses the delta of the Peruípe River, which is 58 km long with a watershed of 4.780 km², and includes the Cassurubá Extractive Reserve, a Marine Protected Area of extractive use. The São Mateus River, approximately 76 km long, has a basin that drains about 8.000 km², formed by the confluence of the Cotaxé and Cricaré rivers (Condini et al., 2022). The Ipiranga River extends about 138 km long, as the southern extension of the Barra Seca River (Silva et al., 2018). The Rio Doce River runs for 850 km from the Mantiqueira and Espinhaço mountains to the Atlantic Ocean, with a watershed basin covering about 86.715 km² (Schettini & Hatje, 2020).

Artisanal fisheries are an important source of income for local communities in the study areas. The coastal zones near the mouths of the Doce and Caravelas rivers are major fishing grounds, where seabob shrimp are exploited (Oliveira et al., 2020). Seabob shrimp samples were collected using trawl nets (13 mm mesh size in the net body and 15 mm in the cod end) for five minutes. Three replicate tows were conducted at four stations in each area at the beginning of each month, resulting in 48 (3 × 4 × 4) samples per campaign (Condini et al., 2022; Oliveira-Filho et al., 2023). After 12 monthly campaigns, this resulted in a total of 576 samples. Seabob shrimp were immediately euthanized in ice slurry and kept frozen laboratory processing, where a subsample of 30 individuals was randomly selected from each haul for measuring total length (TL), carapace length (CL) and wet weight (WW) (Fig. 2).

Assessing growth and mortality with traditional length-based methods

For both synthetic (simulated) and real-world scenarios, we initially estimated L_∞ and Z/K for each sample using a “traditional” length-based approach, specifically the modified version of the PW-plot method (Pauly, 1986). This method relies on the partition of the size frequency distribution into size classes of equal intervals using a series of arbitrary cutoff lengths (L_c). The mean lengths of all fish within each i-th size class (L_mean[i]) are then calculated and used to fit a simple linear regression of the difference between the midpoints and cutoff lengths (L_mean - L_c) against L_c. The intercept (α) and slope (β) estimates from this regression are then used to calculate L_∞ as L_∞ = α/- β and Z/K as Z/K = β/(1- β).

In this framework, data points were automatically selected using standardized “gamma” selection (initial point: “mode + 10%”; last point: 80% of the maximum length, R code available at http://rpubs.com/rschwamborn/267465). Although this selection was arbitrary, we previously verified that other values for the last point yield similar results. The estimated L_∞ values from the PW method were then used as fixed inputs for subsequent estimation of the instantaneous growth rate K through the “ELEFAN_GA_boot” function (Schwamborn, Mildenberger & Taylor, 2019) Then, total instantaneous mortality (Z) was estimated by the iLCCC method using these K and L_∞ values as inputs.

Length-frequency analysis with fishboot

VBGF parameters (L_∞, K) were estimated simultaneously using modern bootstrapped algorithms within the fishboot toolbox (Schwamborn, Mildenberger & Taylor, 2019). First, LFDs for each region were organized by sampling months with a bin size of four mm following the rule of thumb by Wang et al. (2020) where optimal bin size = 0.23 × $L_{max}^{0.6}$. These LFDs were then used as input data for the bootstrapped ELEFAN with a genetic search algorithm (ELEFAN_GA_boot) implemented in the “fishboot” R package (Schwamborn, Mildenberger & Taylor, 2018; Schwamborn, Mildenberger & Taylor, 2019), which builds on functions from the “TropFishR” R package (Mildenberger, Taylor & Wolff, 2017). This algorithm employs a complex genetic framework based on natural selection (Xiang et al., 2013) to improve the goodness-of-fit in the ELEFAN curve fitting process, initially applied to LFA within the “TropFishR” package (Mildenberger, Taylor & Wolff, 2017). It searches for optimal von Bertalanffy growth equation parameters by testing several combinations through multiple iterations until finding a local maximum (i.e., the curve that best performs with the given LFDs). For each study area, LFD data was resampled 500 times in the bootstrap routine using the ELEFAN_GA_boot algorithm in fishboot (Schwamborn, Mildenberger & Taylor, 2018). A very wide and unconstrained search space was used (K: 0.05 to 3 yr⁻¹, L_∞: 50 to 200 mm, t_anchor: 0–1), following standard procedures for unbiased parameter search (Schwamborn, Mildenberger & Taylor, 2019). This approach was considered more accurate as it is intrinsically unrestricted, and therefore assumed to represent the most likely population parameter distribution without potential artifacts from a priori defined narrow search space limits (Zhou et al., 2022). The moving average (MA) was set to a value of 11 since higher values improve modal detection when using intermediate bin sizes (Taylor & Mildenberger, 2017).

Natural mortality

Natural mortality rates (M) for seabob shrimp in the real-world scenario were estimated by the growth-based empirical method of Then et al. (2015) as follows: ${\displaystyle \mathit{M}=4.118{\mathit{K}}^{0.73}{\mathit{L}}_{\infty }^{-0.33}}$ where K is the instantaneous growth rate (yr⁻¹), and L_∞ is the theoretical asymptotic length in cm.

Comparing total mortality (Z) estimates fromtraditional and bootstrapped methods

Relative bias of the PW method (RB, in %) was calculated relative to the fishboot posterior distributions as following: ${\displaystyle \text{RB}\left(\text{\%}\right)=\left[\left(\text{Observation}-\text{Estimate}\right)/\text{Observation}\right]\times 100.}$ For the simulation study, relative bias was calculated based on the true input parameters from the synthetic population (Table 1). For the seabob shrimp real-world scenario, fishboot outputs were considered to be “observations” and outputs from the PW model were considered as “estimates”. Subsequent mentions of “bias” in the results refer to relative bias outlined in this section.

For each parameter (L_∞, K, Z, etc.), 95% confidence intervals (CI) were calculated from bootstrapped posterior distributions as the 0.25 and 97.5% quantiles.

Simulation study

After 500 bootstrap iterations, the genetic optimization algorithm consistently converged to the “true” parameter values (K, L∞) used to simulate the population, with less than 5% of bias in the mean posterior distributions relative to the known parameters (Table 1), and relatively narrow uncertainties. Subsequent mortality estimates from the iLCCC using the estimated growth parameters to calculate relative age were also accurate, with less than 5% of bias (Fig. 3). Contrastingly, growth parameter and subsequent mortality estimates from the PW method showed significant bias, exceeding 10% for asympptotic length L∞, and almost 30% and 50% of bias for the instantaneous growth rate K and total mortality Z, respectively (Fig. 3).

Real-world scenario

The median best-fit parameter estimates obtained from the fishboot analysis of seabob shrimp length compositions resulted in L_∞ values considerably below the overall L_max for all regions, with reasonably small confidence intervals (CI = 8%). Growth coefficients (K) showed a rather distinct behavior between samples, with high intra-cohort variability and wide distributions (ranges between 0.44 and 1.85, CI = 25%) for all regions. Most importantly, the best-fit L_∞ was much smaller than L_max for all samples, with only values larger than the third quantile of the posterior distributions of L_∞ being larger than L_max. This discrepancy precluded the use of the traditional LCCC method and led to the development of the iLCCC method (see methods section). Similarly to the growth parameter estimates, Z estimates displayed relatively large values, with wide 95% confidence intervals (Table 2).

Growth/mortality (Z/K) ratios obtained using fishboot were very close to 1 and showed persistent multimodal distributions (Fig. 4 and Table 2). Further comparisons between the outputs of both methods suggested that the PW method is inappropriate to estimate L_∞ and subsequent analyses, at least in populations with a low Z/K ratio. Mortality estimates were remarkably different (by more than four-fold) between the fishboot (ELEFAN_GA_boot and subsequent iLCCC) and traditional (PW-plot, ELEFAN_GA_boot and subsequent iLCCC) methods, with considerable between-region variability (Fig. 5). This resulted in an astonishing amount of bias in the mortality estimates from the traditional PW approach, likely due to a multiplicative effect from the bias in L_∞. In fact, a bias of only 20% in L_∞, which might be considered as “moderate”, can lead to more than 200% of bias on total mortality Z. The PW method, with such a huge amount of bias, may be considered an extremely biased or “erroneous” approach. (Fig. 6). Additional analysis revealed that the total mortality estimates resulting from the PW method are highly correlated with the bias in L_∞. Consequently, bias in both total mortality and Z/K values appear to be highly dependent on the bias in L_∞ (Fig. 7). Interestingly, natural mortality (M) for the four different areas were remarkably similar to or even greater than total mortality values (Table 2). While this would imply minimal fishing at the four fishing grounds, it is more likely that such estimates are invalid for the studied populations.

Pitfalls in parameter estimation with traditional methods

According to the simulations conducted by Schwamborn (2018), the PW-plot routine considerably overestimates L_∞ for type “A” populations (low mortality/growth relationships and maximum length larger than L_∞). This was also observed in our study, where the simulated population showed very low Z/K values using fishboot, which appear to be strongly influenced by L_∞ (Fig. 3 and Table 1). For our seabob shrimp data, the PW method consistently produced L_∞ values that were much larger than the maximum length of the sample. Contrastingly, the best-fit L_∞ values obtained using the fishboot method were much smaller than the maximum length for all samples (Table 1). This culminated in deviations larger than 50% between the results of both methods in most runs (see Fig. 5). Ideally, mathematical and statistical models attempting to describe the dynamics of natural processes should strive to predict system behaviors with a minimal number of parameters and reduced complexity. However, unrealistic and oversimplifying assumptions can lead to misleading overconfidence on the model’s ability to describe a process or estimate a crucial parameter (Collie et al., 2016; Haddon, 2011). This is the case with the PW routine: as stated by Schwamborn (2018), its assumptions—such as steady-state equilibrium conditions, no intra-cohort/seasonal variability in growth parameters, and no gear avoidance (perfectly-sampled populations)—are major sources of errors in estimating asymptotic size. Another oversimplifying assumption, present in widely referenced fisheries textbooks (e.g., Sparre & Venema, 1992), is that L_∞ should always be above the maximum observed size. Traditional methods to estimate growth parameters through LFA such as the von Bertalanffy plot do not accept lengths that are larger than L_∞ (e.g., “The argument of the logarithm…must be positive as the logarithm would otherwise not be defined. Thus, the von Bertalanffy plot cannot accept a length greater than L_∞”, Sparre & Venema, 1992). This assumption (L_t >L_∞) holds true only when there is little and homogenous variability in growth (not realistic), or when there is a high mortality/growth ratio (as in “dwarfed” populations, see Schwamborn & Moraes-Costa, 2019).

Overestimation of L_∞ from the PW method can lead to considerable errors in estimating total mortality with the widely used length-converted catch-curve (LCCC) method. This occurs because relative age estimates, which are used as a solution to convert the observed sample length composition into ages (the X-axis in the catch-curve plot), will be negatively biased with a larger asymptotic length. Consequently, a much steeper slope (higher mortality estimate) will be produced in the linear regression model (Pauly, 1983). This is simply because it produces an effect where the largest size/age classes of the population, known to always exert a high influence in parameter estimates (Wang et al., 2020), wrongly appear to be missing from the sample and are then assumed by the model to be depleted, while in fact those would not exist. According to our study, errors when converting the length composition of the sample into relative ages derived from erroneous L_∞ estimates have a multiplicative effect on the mortality estimation. As such, this can create an undesirable and misleading situation where even an “acceptable” 10% bias in the asymptotic length from the PW method can lead to extreme (>100%) deviations from the “true” total mortality and Z/K values (Fig. 5). Bias in L∞ estimation can also result from sampling-related issues such as selection bias, where, for example, when larger individuals are not well sampled. This could occur in cases where distinct size classes or life stages use different habitats (e.g., larger individuals occupying different depths, Castrejón, Pérez-Castañeda & Defeo, 2005), which would also impact subsequent mortality estimates.

Errors in estimating L_∞ can also lead to erroneous estimates of natural mortality (M). This is one of the most important parameters in stock evaluations, but extremely difficult to directly estimate in the absence of pre-exploitation data (Punt et al., 2021). Therefore, stock assessments of data-limited stocks often rely on empirical estimates that either directly incorporate asymptotic length in its equation or use parameters such as the theoretical maximum attainable age, often calculated using the inverse von Bertalanffy equation which uses L_∞ (Dureuil et al., 2021; Hewitt & Hoenig, 2005; Jørgensen & Holt, 2013; Kenchington, 2014; Pauly, 1980; Then et al., 2015). It is likely that the natural mortality of the seabob shrimp population in this study, if estimated using outputs from the PW method, would likely be significantly overestimated. According to Punt et al. (2021), positive bias in M can lead to considerable overestimation of stock biomass and a subsequent stumble in its ability to withstand fishing pressure. Ultimately, errors in natural mortality will lead to unreal estimates of fishing mortality if the catch curve method is used, which are crucial input values for models that provide biological reference points for management such as virtual population analysis (VPA) and yield per recruit (YPR) models (Pauly & Soriano, 1986).

Given the absence of information about the “real” parameter distribution for each sample, we acknowledge and emphasize that a central assumption in our real-world study was that outputs from fishboot and length-converted catch-curve methods can be considered reliable reference data for measuring relative bias from the PW-plot outputs. Despite using an unconstrained, robust approach for parameter search in the fishboot, it should be noted that this indirect method (since it does not rely directly rely on length-at-age data to estimate growth) has potential shortcomings that could affect accuracy in parameter estimation such as its sensitivity to the reduced abundance of large individuals in the sample and overall sample size (Schwamborn, Mildenberger & Taylor, 2019). Unconstrained search ranges are recommended for the optimization routines used in TropFishR and fishboot, but current research has demonstrated that the choice of search ranges from wide to narrow impacts parameter estimates (Zhou et al., 2022), possibly because restricting the search space induces the optimization algorithm to get trapped in local maxima (Mildenberger, Taylor & Wolff, 2017). The length-converted catch-curve, used in this study to derive total mortality estimates, has a set of simplifying assumptions that are unlikely to hold on dynamic fished populations such as constant recruitment (steady-state) and mortality rates during the analyzed period, and no variation in selectivity-at-age (Chapman & Robson, 1960). Currently, it is known that mortality greatly varies within a population, usually decreasing with body size due to reduced predation risks (Ramírez-Rodríguez & Arreguín-Sánchez, 2003).

Estimation of natural mortality for invertebrate stocks

In a priori trials and data exploration aimed at conducting a data-limited style stock assessment with the seabob shrimp data, estimation of natural mortality (M) using empirical equations (Then et al., 2015) yielded estimates that were very close (or even above) our estimates for total mortality (Z), which was unexpected. This could indicate, at first sight, a negligible fishing mortality F (assuming Z = M + F). However, considering the intensive trawl fisheries in the study area, a more likely explanation might be a serious inadequacy in the methods used to estimate M, or Z, or both. Since the overall robustness and reliability of the LCCC was verified to be acceptable under most conditions to approximate total mortality rates (Huynh et al., 2018), it is more likely that the estimate of M obtained from empirical equations inferred by relationships between mortality and life-history parameters in fish stocks may not be suitable for penaeid shrimps. Generally, benthic shrimp species spend most of their time buried in sediment, have cryptic, nocturnal habits, and possess long, hard spines (especially X. kroyeri, which has an unusually long rostral spine). Therefore, they may suffer lower predation mortalities than small pelagic fish of similar size (“small pelagics don’t die of old age”). This relationship between morphological features and M has already been demonstrated for fish, where reliable estimates of M are available (Griffiths & Harrod, 2007). In a reanalysis of an updated data set based on the one used for Pauly’s (1980) empirical equation, Griffiths & Harrod (2007) showed that 15 of the 18 species belonging to the 11 families with the lowest mortality rates had morphological defense structures (e.g., venomous spines or defensive armor) or burrowing habits (e.g., in flatfish).

Thus, it is likely that popular empirical equations for estimating M (e.g., Pauly, 1980; Then et al., 2015), which are constructed from fish mortality data and designed for (mostly pelagic and unarmored) fish, are not suitable for well-armored shrimp with distinct behaviors. This is not a widely accepted rationale, except within the study of Mamie et al. (2008). Therefore, we conclude that F cannot be reliably calculated from the difference between Z and empirical M for our data set. This is mainly because our results show that M estimated from empirical equations is virtually equal to Z, which would absurdly imply that there are no fisheries in any of the four study areas, where in fact there are important fisheries targeting seabob shrimp. However, numerous authors have used such empirical equations for many penaeid shrimp stocks (Mehanna, 2000; Leite Jr & Petrere Jr, 2006; Amin et al., 2008; Fernandes, Keunecke & Di Beneditto, 2014; de Barros et al., 2021; Razek et al., 2022), and this approach remains standard for estimating M in data-poor shrimp fisheries. Failures by previous authors to detect issues in M estimation for shrimp stocks might be related to biased Z estimates, which are strongly influenced by L_∞ and K during LCCC procedures. In such cases, errors in total mortality estimation could be due to the assumption that L_∞ >L_max, which also leads to biased K estimates if the search space for L_∞ is constrained. If we assume that K estimates have been severely biased in previous studies with fixed or constrained L_∞ (Schwamborn, 2018, this study), and Z estimates (e.g., from LCCC) are strongly affected by the input K value (Schwamborn, 2018), it is obvious that the Z estimates in these previous studies could have been biased, affecting all subsequent analyses.

Reliable in situ estimates of M for penaeids are scarce (Pérez-Castañeda & Defeo, 2005; Aranceta-Garza et al., 2016). For the brown shrimp Crangon crangon (Caridea: Crangonidae), a small-sized species without a long rostral spine, detailed studies on its mortality were conducted by Oh, Hartnoll & Nash (1999) and Hufnagl et al. (2010), both of which showed a substantially high M/K ratio. Oh, Hartnoll & Nash (1999) found very close estimates for Z (3.96 yr⁻¹) and M (3.60 yr⁻¹), with a Z/K ratio of approximately 4 for the same species in an area with negligible fisheries. Conversely, estimates of M for other penaeids in the adult phase seem to be lower, with a M of about 0.96 year⁻¹ for banana prawn Penaeus merguiensis in tropical Australia (Wang & Haywood, 1999), and around 1.63, 2.52, and 2.89 yr⁻¹ using gnomonic intervals for Farfantepenaeus californiensis, Penaeus stylirostris, and Penaeus vannamei, respectively (Aranceta-Garza et al., 2016). While our data suggests that “empirical” equations used to estimate M may not be appropriate for some crustacean stocks, we recognize our discussion on this issue remains largely speculative, and further studies should examine the relationships between field-based M estimates and life-history/morphological traits.

Estimating Z for data-poor stocks is far from trivial. Severe bias in mortality estimation using common length-based methods has been demonstrated by Hufnagl, Huebert & Temming (2013) and Schwamborn (2018). Obtaining reliable, robust estimates for the mortality components M and F is even more challenging. Ideally, F should be derived from extensive catch and biomass data (Hart, 2012) or tagging (mark and recapture studies) for reliable population modeling and stock assessment. However, tagging shrimp has been proved to be costly, time-consuming, and often problematic, especially due to tagging-related mortality in these more fragile, small organisms (Jewett, 1986). Similarly, catch and biomass time series are rarely available for small-scale artisanal shrimp fisheries. In published shrimp stock assessments, M has generally been generally estimated from empirical equations and used to estimate F from Z. In a few exceptions, M was fixed at an arbitrary value (e.g., 0.3 y⁻¹ or 0.5 y⁻¹) and F was estimated in situ from shrimp catch and biomass data (Mamie et al., 2008; Hart, 2012) using statistical catch-at-age models (Methot Jr & Wetzel, 2013).

In this study, the improved length-converted catch-curve (iLCCC) is of particular novelty. In the old, commonly used LCCC method (Pauly, 1984a; Pauly, 1984b), L_∞ necessarily must be above the maximum length. This is not always a realistic assumption. In many cases (e.g., when growth is highly variable and the Z/K ratio is low, as in simulated populations from Schwamborn, 2018), L_∞ can be below L_max, as shown in the present study (also see Schwamborn, 2018; de Barros et al., 2021; Hossain Yeamin & Ohtomi, 2010 for examples of populations with asymptotic length below maximum length). The implicit, un-circumventable rule of the common LCCC that L_∞ must be above L_max can lead to limitations and bias in the choice of L_∞ search spaces, thus impeding an unconstrained search. Surprisingly, the serious drawback of the traditional LCCC method has not been adequately highlighted and discussed so far in the published literature. In the present study, we highlight this drawback and present a new, improved iLCCC method (de Barros, 2022), that has no intrinsic limitations regarding the amplitude of a priori defined search spaces for L_∞. While our data suggested the iLCCC is adequate to estimate mortality rates, we also highlight the need for more studies to further explore the accuracy of this method against known values for other populations, particularly those with different mortality/growth relationships.