That’s All I Know: Inferring the Status of Extremely Data-Limited Stocks

There is a growing number of methods to assess data-limited stocks. However, most of these methods require at least some basic data, such as commercial catches and life history information. Meanwhile, there are many commercial stocks with an even higher level of data limitation, for which the inference of stock status and the formulation of advice remain challenging. Here, we present a stepwise approach to achieve the best possible understanding of extremely data-limited stocks and facilitate their management. As a case study we use a stock of the shrimp Plesionika edwardsii (Decapoda, Pandalidae) from the eastern Mediterranean Sea, where the only available data was a sub-optimal sample of length frequencies coming from a small-scale trap fishery. We use a suite of different methods to explore and process the data, estimate the growth parameters, estimate the natural and fishing mortalities, and approximate the reference points, in order to provide a preliminary evaluation of stock status. We implement multiple methods for each step of this process, highlighting the strong and weak points of each one of them. Our approach illustrates the better insights that can be gained by applying ensembles of models, rather than a single ‘best’ model when working with limited data of poor quality. The stepwise approach we propose here is transferable to other extremely data-limited stocks to elucidate their status and inform their management.


INTRODUCTION
Depending on the amount of available information, fish stocks can be characterized as datarich or data-limited. Data-rich stocks contain enough information to carry out analytical stock assessments, while data-limited ones do not. However, there are several levels of data-limitation. The International Council for the Exploration of the Sea (ICES) identifies six different stock categories with regards to data availability (ICES, 2012). Categories 1 and 2 include data-rich stocks with age-structured catch and survey data allowing quantitative assessments. These assessments are considered to describe adequately the true trends of stock size and exploitation levels; as such, trends of category 1 and 2 stocks are used to monitor the effectiveness of fisheries regulations (STECF, 2020). Categories 3-6 include stocks with progressively increasing data limitations. In category 3, survey data are available which can indicate trends of mortality rates, recruitment, and biomass. In category 4, a time-series of catch data is available which allows an approximation of maximum sustainable yield (MSY). In category 5, only landings data are available. Finally, category 6 includes negligible landings stocks and stocks caught as bycatch. This latter category of extremely data-limited stocks is the focus of the current study.
There is an ever-growing number of data-limited assessment methods focusing on stocks falling primarily within data categories 3 to 5. For example, for category 3 stocks, the surveybased assessment method (SURBA) (Beare et al., 2005) or the time-series analysis assessment (TSA) (Fryer et al., 1998;ICES, 2008) can be used to estimate population numbers and fishing mortality rates based on survey data. For category 3/4 stocks, surplus production models can be used to estimate biomass and exploitation level for commercial stocks when their age and size data is absent (Punt, 2003). These models are suitable for stocks with data from commercial catches along with indices of exploitable biomass (from catch-per-unit-effort, or survey data) (Polacheck et al., 1993). For example, Pedersen and Berg (2017) presented a stochastic surplus production model in continuous time (SPiCT) which combines dynamics of biomass and fisheries with remarked error of catches and biomass indices. For category 4/5 stocks, where only time-series of catches or landings are available, methods such as the CMSY have been proposed (Martell and Froese, 2013) to estimate extracted yields in relation to MSY. Froese et al. (2017) updated the CMSY method by using catch and productivity to assess biomass. In addition, this method can approximate exploitation rate, MSY and fishing reference points. Froese et al. (2017) also used a Bayesian state-space estimation model (BSM) (Meyer and Millar, 1999) to verify and evaluate the CMSY model.
The examples of data-limited methods mentioned earlier are not exhaustive, but they are indicative of the fact that most data-limited methods require at least some information from surveys, commercial catches and productivity in order to estimate stock status. However, in the case of extremely data-poor stocks (category 6), it is not possible to apply such methods. In such cases, the starting point is the estimation of life history traits, such as growth and maturity. These life history traits can be used order to infer sustainable harvesting strategies, even if the exact stock status is unknown (Froese, 2004;Froese et al., 2008;Prince and Hordyk, 2019). Growth parameter estimates can also be used to estimate mortality rates, approximate reference points, and infer the stock status; a suite of different methods exists for every step of this way (Gayanilo and Pauly, 1997). Recently, a new method for the analysis of extremely data-limited stocks of fish and invertebrate species was proposed: the Length-based Bayesian Biomass (LBB) method (Froese et al., 2018). LBB requires only length frequency distributions (LFDs) as an input, as it makes a series of assumptions for the estimation of the missing life history information. LBB's key outputs are the current exploited biomass relative to unexploited biomass (B/B0) and the fishing mortality relative to natural mortality (F/M) (Froese et al., 2018).
Typically, in extremely data-poor situations it is difficult to identify a single method that produces the 'best' estimate of a given variable. In that case, an ensemble modeling approach combining the outputs from multiple methods can help produce more robust estimates (Dormann et al., 2018). This process can in turn inform fisheries management more effectively and facilitate measures to promote fisheries sustainability.
In this study, we present a stepwise methodological framework to estimate the stock status of extremely data-limited stocks. We illustrate the use of this framework by estimating the stock status of an extremely data-limited shrimp stock (Plesionika edwardsii, Decapoda, Pandalidae), which is a bycatch of a small-scale trap fishery in the eastern Mediterranean Sea. Our proposed framework includes various methods for estimating growth parameters, calculating mortality rates, approximating reference points and eventually characterizing stock status. We synthesize the outputs from different model combinations, illustrating the advantages from using an ensemble modeling approach, and compare our findings with the outputs from a relevant LBB. This way, we elucidate the stock status of the studied shrimp stock and deduce implications for its management, in a way that is reproducible to other extremely data-limited stocks.

Sampling
Samples were collected in the Dodecanese archipelago (southeastern Aegean Sea) (Figure 1) from October 2014 to October 2015 on a monthly basis (except of February 2015 due to adverse weather conditions), under the framework of PLESIONIKA MANAGE project. PLESIONIKA MANAGE studied the biology and exploitation of Plesionika narval, a valuable fishery resource in the area (Kalogirou et al., 2017;Maravelias et al., 2018;Vasilakopoulos et al., 2019). P. edwardsii was a commercial bycatch of the fishery for P. narval, far less valuable than the targeted congeneric species. Sampling depth ranged from 0 to 280 m and was divided in strata A (0-45 m; depth of thermocline in the summer), B (46-100 m; to the end of the continental shelf), C (>100 m). Circular traps covered with nylon-based net of a mesh size of 12 mm (knot to knot) and with an upper trap opening of 13 cm were used . Fishing took place during night hours (20:30-06:00); after 9.5 h all traps were hauled, and the catch was separated in two categories: target (P. narval) and by-catch species, the latter including P. edwardsii. A random sample of approximately 100 shrimps, mainly P. narval, was collected from each stratum. The number, size (carapace length and body weight), sex and maturity stage of P. edwardsii shrimps within these samples was also documented (Figure 2).

Exploring and Processing Data
Carapace length (CL) was rounded at the nearest millimeter and LFDs were estimated by sex, maturation stage, month and depth strata, using length classes of 1 mm Due to the great variability of the sample size across months, standardized LFDs by month were estimated, by dividing the numbers of individuals within each length class by the total number of monthly samples (Figure 2). These standardized LFDs were used to: (1) identify the main spawning periods; (2) identify the depth distribution of males and females; (3) estimate the sex ratio within each length class; (4) create a length frequency object (LFQ) to be used for the estimation of the growth parameters. For this, the sexcombined LFD by month was restructured as a list object containing a catch frequency matrix, a vector of midlengths corresponding to rows of the catch matrix and a vector of sample dates corresponding to the columns of the catch matrix.

Estimation of Growth Parameters
This analysis was conducted using R programming language (R Core Team, 2020) and the TropFish R package Taylor and Mildenberger, 2017). To estimate growth parameters (Figure 2), the ELEFAN (Electronic Length Frequencies Analysis) program (available as tool in the TropFishR package) was used. ELEFAN calculates a moving average (MA) over the LFDs bin, and then compares the observed frequency with this average; values much above the average indicate a "true" mode (Pauly and David, 1981;Pauly, 1985). The LFQ file was prepared for running the ELEFAN by posing a MA to generate the Von Bertalanffy Growth Function (hereafter VBGF) estimations. For this, we set ranges for the infinite length (L inf ) and growth coefficient (K) values, and a theoretical time zero (t 0 ) at which individuals of this species hatch. Preliminary estimations of L inf were done based on three different approaches: • the Powell and Wetherall method (Powell, 1979;Wetherall et al., 1987): a linearizing transformation of length classes to estimate L inf by plotting L mean -L' and L'. L mean is the mean length of all individuals greater than L' and L' is the smallest length of fully represented individuals in catches. • the empirical formula from Froese and Binohlan (2000):L inf =e 0.44+0.984log(Lmax) ? with L inf being infinite length and L max being maximum observed length. • the maximum carapace length observed in the samples.
The Powell and Wetherall method is very sensitive to intra-cohort variability in growth and to changes in the occurrence of large individuals in the sample, resulting often in underestimation of the L inf value (Schwamborn, 2018). Exclusion of the largest size classes during the regression procedure or weighing by abundance does not resolve these issues (Schwamborn, 2018). By contrast, the Froese and Binohlan formula tends to overestimate L inf values. Accordingly, two different length ranges for L inf were chosen to run ELEFAN. The first length range had a lower limit of L inf derived from the Powell and Wetherall method and the second one had a lower limit of L inf equal to the maximum length observed in the samples. The L inf value derived by the Froese and Binohlan formula was chosen as the upper limit for the length range of L inf in both cases.
The range of K values was set between 0.4 to 0.9 y −1 , based on previous publications on the growth of this species (Santana et al., 1997;Company and Sardà, 2000;García-Rodriguez et al., 2000;Colloca, 2002).
Four different scenarios of the month when length is equal to zero (t anchor in ELEFAN, conceptually similar to t 0 of the VBGF) were explored: (i) February, (ii) May, (iii) August, and (iv) November. The one resulting in the best fit and agreeing with the spawning information was chosen as the optimal t 0 (Supplementary Figures 6, 7).
The MA value used in the ELEFAN analysis was set based on two different scenarios. The first scenario used the default setting in FISAT II (Gayanilo and Pauly, 1997), i.e., a width of 5 bins (1 bin = 1 mm) for each cohort. Indeed, the smallest cohort width was up to 5 bins so it was plausible to compare each bin with the average across 5 consecutive bins (i.e., ±2 bins to either side). However, taking into consideration that the MA settings can significantly affect the scoring of the growth curve  a second scenario was also tested, with a 7-bin cohort width.
The estimation of the best fit was based on searching for the VBGF parameters with the maximum score value (Rn) as a measure of relative fit: where the Estimated Sum of Peaks (ESP) is the sum of peak values crossed by the growth curves, with the caveat that positively crossed bins are only counted once, while negatively crossed bins are counted every time they are encountered (Pauly, 1985). The Available Sum of Peaks (ASP) is the sum of all positive peaks, which represents a maximum possible score (if negative bins are crossed). Rn can attain a maximum value of 1. Fitting scores across the whole range of L inf and K combinations was visualized by a Response Surface Analysis (RSA).

Estimating Total, Natural and Fishing Mortality
The VBGF parameters were used both to compute the length at which 50% of the cumulative catch is captured (L 50 ) and to estimate the total mortality (Z) (Figure 2) Z was computed using two methods: • the Length Converted Catch Curve (LCCC) (Pauly, 1990) • the relevant B&H formula (Beverton and Holt, 1956) The Length Converted Catch Curve (Pauly, 1990) is a way to estimate Z plotting the natural logarithm (loge) of the number of fishes in the sample (N) against the relative age corresponding to the midrange of the length class in question [ t is the time needed to grow from the lower (t 1 ) to the upper (t 2 ) limit of a given length class]: where a and Z are the regression parameters.
In the second case (Beverton and Holt, 1956), Z value is calculated as: where L inf and K are parameters from VBGF, L mean is the mean length in the catches and L' is the smallest length of animals that are fully represented in the catch samples. A suite of different methods and formulas were used to estimate natural mortality (M) values (Figure 2 and Supplementary Table 1). This great number of methods reflects the fact that M is notoriously difficult to estimate (Kenchington, 2014), and it has a big effect on our perception of stock status (Mannini et al., in review). Four methods (and their variants) provided empirical scalar values (Alverson and Carney, 1975;Pauly, 1980;Hoenig, 1983;Hewitt and Hoenig, 2005;Then et al., 2014). Pauly's (1980) equation was computed using three different bottom temperature values (14, 16, 18 • C), based on the seasonal difference in bottom temperature. Seven more methods (and their variants) produced natural mortality vectors by age (Gulland, 1965;Chen and Watanabe, 1989;Caddy, 1991;Abella et al., 1997;Lorenzen, 2000;Gislason et al., 2010;Brodziak et al., 2011;Martiradonna, 2012) (Supplementary Figure 10). The estimated VBGF parameters were used to convert the maximum CL observed in the catch into age. For each of the M vectors by age a mean over the range between age 0 and maximum observed age was computed to obtain the corresponding scalar value. This process resulted in sixteen different M scalar values. These values were subtracted from the two Z values calculated earlier, to provide a set of 32 different values of fishing mortality (F) (Figure 2).

Reference Points and Advice
In the final stage of the analysis, reference points for management (Caddy and Mahon, 1995) were computed according to a Yield per Recruit (YpR) model (Beverton and Holt, 1957) (Figure 2). In the Mediterranean, two reference points are being used for exploitation levels: F 0.1 and E 0.4 (STECF (17-15), 2017; STECF (19-16), 2019). The fishing mortality level F 0.1 is the F rate at which the slope of the yield per recruit curve as a function of F is 10% of its value at the origin (Gulland and Boerema, 1973). E 0.4 comes from Patterson (1992) who suggested as reference point in terms of exploitation rate (E), in particular for pelagic stocks, a value of } E = F/Z = 0.4.
Thirty two different YpR analyses were run, corresponding to the 32 estimates of F, to extract the respective F 0.1 values. For all scenarios, L 50 of the selectivity was set as the length of 50% of the cumulative catch. The status of exploitation was estimated according to two ratios: F F0.1 and E E0.4 , for which values over 1 indicate overfishing.
A kobeplot was used to visualize the 32 scenarios results (statuses of exploitation), as well as the mean and the median of them (Figure 2). Four areas with different colors based of E status were plotted in a Cartesian system in which x-y intersection is set at point equal to 1, 1. In the green area fell points having ratios below or equal to 1 for both E E0.4 and F F0.1 . In the red area fell points having both ratios over 1 and, in the yellow areas points for which one of the two ratios was over 1.

Applying the LBB Method
LBB (Froese et al., 2018) requires only LFDs and is able to estimate L inf , length at first capture, relative natural mortality (M/K) and relative fishing mortality (F/K). M/K and F/K can be combined to give fishing mortality relative to natural mortality (F/M). LBB also estimates an approximation of current exploited biomass relative to unexploited biomass (B/B 0 ).
To apply this method, the monthly LFDs we had available were aggregated to a yearly sample. The L inf prior was set equal to our best estimation from ELEFAN. In assigning M/K priors, K was set equal to our best estimation from ELEFAN while M values changed according to the 16 methods presented in chapter "Estimating total, natural and fishing mortality." All other inputs were set to the default values suggested by Froese et al. (2018).

Length Frequency Distributions
In total, 1993 P. edwardsii individuals were sampled. Plotting the LFD using the raw data showed strong monthly variability in the number of individuals (unbalanced data), ranging from 4 in October 2014 to 693 in November 2014, and irregular population distributions (Tables 1, 2). The length at which 50% of the cumulative catch is captured (L 50 ) was estimated as 21.70 mm. (Supplementary Figure 8).  After standardizing the LFDs and examining separately mature and immature individuals, a dominance of mature females emerged between spring and mid-summer indicating a spawning period (Tables 1, 2).
Depth was found to be influential, with more individuals found in the deeper strata (Supplementary Figure 1). A dominance of females was evident at all depths with a sex ratio of 0.82 (Supplementary Figure 2), complemented by larger females (up to 30 mm) at increased depths (Supplementary  Figures 3, 4). Mean size of males (19.83 mm CL) was lower than that of females (21.40 mm CL), with a maximum size of 26 mm.

Growth Parameters
Visualizing the raw and standardized LFQ data (Tables 1, 2 and Figure 3) showed that recruits to the fishery appear in March at a CL of ∼12 mm. After exploring possible months for t anchor (Supplementary Figure 4), May was selected (t anchor = 5/12) (Supplementary Figure 5). This month provided the best fit, and it was also when the higher fraction of mature females was observed.
L inf estimations varied between the different methods used. L inf was estimated at 26.64 mm by the Powell and Wetherall method (Supplementary Figure 5), at 44.13 mm by the Froese and Binohlan formula, while the maximum CL observed in the samples was 30.00 mm ( Table 3). Although the Powell and Whetherall method estimation with MA = 7 got the highest Rn value (Table 3), its resulting VBGF parameter estimates were not retained. That was because the relevant L inf estimation (27.24 mm) was underestimating the population's true L inf , being lower than the larger individuals sampled (30.00 mm) ( Table 3). By contrast, the L inf estimation by the Froese and Binohlan formula was found to be greater than the larger individuals sampled. Therefore, to select the optimal VBGF parameter estimates we focused on the two runs using the maximum length as a lower limit of the L inf range ( Table 3). Among these two runs, the one with a MA width of 5 bins had the highest Rn value and a t anchor value closer to the assumed spawning period ( Table 3). The final VBGF parameter estimates were: L inf = 31.94 mm, K = 0.56 y −1 , t anchor = 0.31 y (Figure 4). Five age cohorts were estimated by the VBGF model ( Figure 5).

Total, Natural and Fishing Mortality
Z was estimated as 2.25 y −1 by the LCCC method (Supplementary Figure 9) and 2.36 y −1 by the Beverton and Holt formula. The sixteen different methods used for M, produced values ranging between 0.44 y −1 (Alverson_Carney) and 2.01 y −1 (Lorenzen) ( Table 4). The combination of the two values of Z with 16 values of M produced 32 different FIGURE 4 | Response surface analysis. The best fit for the combination between Linf (in mm) and K-values is indicated with a red asterisk on the plot. The color scale corresponds to the Rn values from ELEFAN analysis. K ranges between 0.4 and 0.9 and L inf ranges between 30 mm (maximum observed length) and 44.13 mm (estimated by Froese and Binohlan). This analysis corresponds to the optimal combination from Table 1.  values for F ranging between 0.24 y −1 (Lorenzen) and 1.92 y −1 (Alverson_Carney) ( Table 4).

Stock Status
In total, only seven out of the 32 estimates of the exploitation state indicated sustainable levels of fishing both in terms of F and in terms of E (Figure 6). These included all six cases where the Lorenzen, ProdBiom and Gulland methods were used and one case where the Gislason method was used for the estimation of M. By contrast, 22 estimates of the exploitation state indicated overfishing, both in terms of F and in terms of E (Figure 5). These included all cases where the Pauly_T1, Pauly_T2, Pauly_T3, Alverson_Carney, Then_1, Then_2, Hoenig_1, Hoenig_2, Then_scaled, Brodziak_Tmax, and Brodziak_K methods were used for the estimation of M. In three cases, two using ChenWatanabe and one using Gislason for the estimation of M, the stock was found to be overfished in terms of E but non-overfished in terms of F. Both the mean and median stock status was estimated as overfished (Figure 6).  Table 5 summarizes the main results from the LBB. L inf ranged between 31.92 mm (Lorenzen) and 32.56 mm (Hoenig); M/K ranged from 0.96 (Alverson_Carney) to 3.60 (Lorenzen) and F/M ranged between 0.41 (Prodbiom) to 5.08 (Alverson_Carney).

Results From LBB
All these values were similar to ones obtained from our original analysis. The exploitation status of P. edwarsii was estimated as non-depleted in three M scenarios (Gulland, Lorenzen and Prodbiom), close to equilibrium in one M scenario (Gislason) and as in moderate or severe depletion in all other scenarios ( Table 5).

DISCUSSION
This study proposes a stepwise methodological framework to assess stock status of an extremely data-limited exploited stock. Lack of data constitutes a common restrictive factor for fisheries management and the applied methods can be decisive for the outcomes. Our proposed framework is a sequence entailing estimations of basic life traits, mortality rates and biological reference points to infer stock status. The combination of high data uncertainty and multiple method availability (Gayanilo and Pauly, 1997) often lead to tradeoffs in the method(s) to be used. By contrast, the numerous estimates for the same variable offer the opportunity to calculate the average (weighted or not) of values instead of presenting a single-method prediction (Dormann et al., 2018). This is especially relevant when the analyzed data come from small-scale fisheries with low data availability, which are prevalent in the Mediterranean Sea.
The range of outputs estimated by this study underlines the necessity of a multi-method approach. For each of the steps presented in this study, a single method selection would be misleading and would constitute a sub-optimal methodological path selection. Previous studies on the biology of P. edwardsii stock (e.g., Colloca, 2002;González et al., 2016) used a single method approach but with much larger and balanced datasets and thus the relevant results were more reliable; however, these studies did not touch upon stock status. In our study, we use a multi-method approach through a stepwise process; the most pronounced example being the estimation of natural mortality rate (M) for which we used 16 different methods. Such a parallel implementation of several methods for M is not usually applied (Kenchington, 2014), and the M value is typically calculated using only one of a handful of M equations (e.g., Pauly, 1980;Hoenig, 1983). The wide range of our results demonstrates the great variability that exists across these 16 approaches used and the complexity of pinpointing the most suitable one. Using various M estimators is often recommended for reducing bias, errors, underestimations and uncertainties of the methods applied (Gunderson et al., 2003;Simpfendorfer et al., 2005). Notably, the range of outputs produced when using our proposed framework was similar to that observed when implementing the novel LBB method (Froese et al., 2018). This highlights how in extremely data-limited situations the choices made with regards to key input parameters (such as M) have a greater impact on the outputs than the analytical method used.
Ensemble modeling use several options and generates more robust outputs. Selection of the most appropriate method may often prove to be more difficult than initially expected because the theoretical and/or empirical background is not precise enough. Applying ensemble modellin is recommended in extremely datalimited situations, because various models can be used for estimating a value and all results can be statistically tested (Kuhn and Johnson, 2013). Using model averaging to infer stock status could be a sufficient solution for decreasing the predicted error. When models are unbiased and with high variance, using an increasing number of different models could minimize the error (Dormann et al., 2018). An alternative to using a simple model mean or median, as used in this study, would be to use a weighted average. Not all estimators may be equally reliable; e.g., some of the M estimators used here may be more suitable for fish than crustaceans. In such cases, a weighted estimation could provide more accurate results (Kenchington, 2014;Dormann et al., 2018).
Extremely data-limited stocks often have a high importance for both fisheries and the marine ecosystem; hence one should strive for analytical examinations to infer stock status. However, the use of a multi-method approach alone does not provide certainty about the real situation in extremely data-limited situations. The analytical framework proposed here, when applied to extremely data-limited stocks, it should be viewed as the first stage of exploration; a starting point to reveal stock status. Next steps should involve a more extensive data collection so that the initial insights could be corroborated or rejected. In the meanwhile, controlling size selectivity and safeguarding stock productivity constitute a sound strategy for the management of extremely data-limited stocks (Prince and Hordyk, 2019), especially in areas such as the Mediterranean Sea which are known to be severely overfished (Vasilakopoulos et al., 2014).
The same principles of better monitoring and precautionary measures also apply to the Dodecanese P. edwardsii stock used as a case study here. Our study suggests that the P. edwardsii stock is likely overfished, in line with its sympatric P. narval stock (Maravelias et al., 2018) and most other assessed crustacean stocks in the Mediterranean Sea  (Vasilakopoulos and Maravelias, 2016). In this context, the managers of the P. narval fishery need to strive for better monitoring of the bycatch P. edwardsii to obtain a better understanding of its status; e.g., by carrying out targeted sampling in deeper strata where P. edwardsii is known to be more abundant than P. narval. If the overexploited state of P. edwardsii is confirmed, management plans for P. narval should take into consideration the state of P. edwardsii as well and adjust the fishing activities accordingly. This study provides a stepwise analytical methodology to be applied to any extremely data-limited stock. The range of methods that we have used in each step is not exhaustive and fisheries scientists are encouraged to use more and/or different methods according to the specific characteristics of the stock at hand. Inferring the stock status of extremely data-limited stocks usually involves unique challenges in every individual case; nevertheless, we are confident that following part or all of the methodological steps proposed here can prove extremely helpful.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
SK, VP, PV, and PM conceived the study. SK and KK contributed to the data collection. AM and VP did the analysis. VP, AM, PV, and SK interpreted the results. VP led the writing of the manuscript with input from all authors. All the authors contributed to the article and approved the submitted version.