Fish Stock Assessment for Data-Poor Fisheries, with a Case Study of Tropical Hilsa Shad ( Tenualosa ilisha ) in the Water of Bangladesh

: The anadromous tropical Hilsa shad formed the largest single-species ﬁshery in Bangladesh, making the highest contribution to the country’s total ﬁsh production (14%) and nearly 83% of the global Hilsa catch in 2018. However, increased ﬁshing pressure made the ﬁshery vulnerable, and hence, information on the stock condition and its response to the current degree of removal is essential to explore the future potential for sustainable exploitation. This study carried out a rigorous assessment based on three different methodological approaches (traditional length-frequency based stock assessment method for ﬁshing mortality and exploitation, Froese’s length-based indicators for ﬁshing sustainability, and a surplus production-based Monte Carlo method-CMSY, for ﬁsheries reference points estimation) for the best possible estimates of the Hilsa stock status in the water of Bangladesh. The present ﬁndings revealed that the stock is likely to be overﬁshed due to overexploitation. Depending on the outputs, this study recommended a lower length limit for the catch (>33 cm), distinguished a selectivity pattern (mesh size limit ≥ 8 cm), and proposed a yearly landing limit (within the range of 263,000–315,000 tons) for the sustainable management of the Hilsa ﬁshery in Bangladesh.


Introduction
The sustainability of marine fisheries became a global concern because of the rapid increase in fishing pressure and arbitrarily exploitation, which may cause negative consequences for the ecosystems and societies [1]. This situation is worse, particularly in developing countries where proper management tools and political will lack and pervasive illegal fishing of juvenile and brood fishes is the fisheries' governing features [2][3][4]. This worst-case scenario of the stocks provides the impetus for articulating effective management tools, focusing on all concerned stakeholders' participation to promote conservation of the ecosystems and achieve the optimum yields to sustain the livelihoods and sustainable food supply [4,5]. Again, an effective management plan encompasses various management tools, including effort control, minimum mesh size regulations for nets, and maximum allowable catch limits [6,7]. Before the setup of those management tools, stock status, especially the current level of exploitation and stock spawning biomass status, needs to address first.
The Hilsa shad, Tenualosa ilisha (Hamilton, 1822), is a member of the Clupeidae family shows a wide range of distribution from the northern part of the Bay of Bengal to the Arabian Sea regions and some parts of coastal Southeast Asia [8][9][10][11][12]. Hilsa spends most of their adult lives at sea and migrates towards the freshwaters, virtually in every accessible Bangladesh is located in the northeast corner of the Bay of Bengal with a diverse and dynamic 710 km long coastline. This coastline is supported by several estuaries, i.e., Meghna river estuary, Karnaphully river estuary, and vast open water bodies over 118,813 square kilometers [52], which is characterized by various fisheries resources, including Hilsa (Tenualosa ilisha) fishery. 12 representative sampling stations (Figure 1) along the coast, including riverine (Lakshmipur [22 • [22 • 03 07 N, 89 • 58 18 E]) were selected for lengthfrequency data collection. Different industrial trawlers and artisanal boats (mechanized and non-mechanized) with varying capacities of storage were found to be operated in this area for Hilsa catch. Trawl nets with cod-end mesh sizes ranging from 2 to 6 cm and large mesh drift and set gill nets with different mesh sizes (2-10 cm) were the predominant Hilsa fishing gears. Different Seine nets (mesh size: 0.5 to 6 cm) and fixed purse nets (mesh size: 1 to 5 cm) were also commonly used by the fishermen.

Study Area
Bangladesh is located in the northeast corner of the Bay of Bengal with a diverse and dynamic 710 km long coastline. This coastline is supported by several estuaries, i.e., Me ghna river estuary, Karnaphully river estuary, and vast open water bodies over 118,813 square kilometers [52], which is characterized by various fisheries resources, including Hilsa (Tenualosa ilisha) fishery. 12 representative sampling stations (Figure 1)  were selected for length-frequency data collection. Different industrial trawlers and artisanal boats (mechanized and non mechanized) with varying capacities of storage were found to be operated in this area fo Hilsa catch. Trawl nets with cod-end mesh sizes ranging from 2 to 6 cm and large mesh drift and set gill nets with different mesh sizes (2-10 cm) were the predominant Hilsa fishing gears. Different Seine nets (mesh size: 0.5 to 6 cm) and fixed purse nets (mesh size 1 to 5 cm) were also commonly used by the fishermen.

Data Collection
Based on sampling carried out at 12 fish landing sites from January 2017 to Decembe 2018, length-frequency data of 6132 individuals of Hilsa (both sex) were collected monthly from the fishermen who used drift and set gill nets for their catches. Lengths were meas ured using digital slide calipers. In some instances, data collection was prevented due to

Data Collection
Based on sampling carried out at 12 fish landing sites from January 2017 to December 2018, length-frequency data of 6132 individuals of Hilsa (both sex) were collected monthly from the fishermen who used drift and set gill nets for their catches. Lengths were measured using digital slide calipers. In some instances, data collection was prevented due to the adverse weather condition. In those cases, data was collected immediately when the bad weather was over. In the landing centers, fishers categorized their catches mostly in three categories according to the fishes' individual weight (fish ≥ 1 Kg, fish ≥ 800 gm < 1 Kg, fish < 800 gm). We collected 10% of well-mixed fish from each category. Landings were very poor in January, June, July, August, and December in Badarkhali, Chakaria; Baro Ghop, Kutubdia; Lakshmipur, Chandpur, Ramgati, and Jahajmara, Hatiya. We collected 30-100% of the landed fish from the respective landing center during the low landing. The total length of all collected samples was measured to the nearest 0.1 cm. A total of 427 individuals was randomly selected from the three categories (77 individuals from fish ≥ 1 Kg, 132 individuals from fish ≥ 800 gm < 1 Kg and 218 individuals from <800 gm.) and weight to the nearest 0.1 gm using a digital balance to estimate the length-weight relationship of the Hilsa shad.
To analysis the stock status from catch and resilience, 31 years of catch and effort data of Hilsa were extracted from the Yearbook of Fisheries Statistics of Bangladesh (1988-2018) published by Fisheries Resources Survey System, Department of Fisheries (DoF), Bangladesh, and from the Food and Agriculture Organization of United Nations (FAO) fisheries and aquaculture-online query panel [52,53]. With few exceptions, the annual landing of Hilsa shad shows an increasing trend, and in 2018, the total landing was 0.52 million tons, which is 52 percent higher than the landing of 2014 and more than 78 percent higher than 2009 ( Figure 2).  [52,53]. With few exceptions, the an landing of Hilsa shad shows an increasing trend, and in 2018, the total landing was million tons, which is 52 percent higher than the landing of 2014 and more than 78 per higher than 2009 ( Figure 2). Due to the unavailability of effort data in inland water bodies, the number of b was calculated from the total Hilsa production in the inland water bodies per annum average catch per boat per annum in marine water. Data interpolation technique was ployed to avoid the anomalies in effort data. Due to the unavailability of effort data in inland water bodies, the number of boats was calculated from the total Hilsa production in the inland water bodies per annum and Sustainability 2021, 13, 3604 5 of 23 average catch per boat per annum in marine water. Data interpolation technique was employed to avoid the anomalies in effort data.

Stock Assessment Indicators
For the evaluation of the stock status of Hilsa, the following three sets of indicators were used in this study: (a) Two biological reference points, fishing mortality (F) and exploitation rate (E), estimated from the linearized catch curves analysis and the yield per recruit (YPR) model, as described by Sparre and Venema (1998) [44]. (b) The length-based indicators (LBIs) were proposed by Froese (2004) for sustainable fishing [2]. (c) Fisheries reference points from catch and resilience, as describe by Froese et al. (2017), using an R-package CMSY [47].

Estimation of Growth Parameters
The von Bertalanffy (1938) growth parameters (K, L ∞ , and t 0 ) were estimated applying the following seasonalized von Bertalanffy's growth function to the length-frequency data using the TropFishR package [38,45,54]: L t is the length of a fish at age t, L ∞ is the asymptotic length in cm, K is the growth coefficient in year −1 , t 0 is the theoretical age of a fish at which its length is zero.
For simplicity, the seasonal growth was not studied here; thus, "+ S(t) − S(t 0 )" in the above equation was set to zero. In TropFishR, Electronic Length Frequency Analysis (ELEFAN) is used for the fitting process [55][56][57][58]. The parameter t anchor in ELEFAN indicates the fraction of a year where the von Bertalanffy growth curve crosses length = 0 for a given cohort [59]. The growth performance index (Φ ) was calculated from the formula given by Pauly & Munro (1984) [60]: L max was estimated from the mean of the 1% largest fish present in the sample, and then an initial value for L ∞ was calculated using the following formula proposed by Pauly (1984) [61]: Then the initial seed value for L ∞ was determined as L ∞ ± 20% and fixed the K range between 0.01 and 2. The range of initial seed values for t anchor and C was 0 to 1 [38]. The most suitable moving average (MA) was determined by trial and error basis with different MA values and rule of thumb described in   [59].

Fishing Mortality and Exploitation
For the estimation of morality, a linearized length-converted catch curve was produced using length-frequency data and previously estimated von Bertalanffy's growth parameters [38]. The instantaneous mortality rate (Z) was estimated from the regression line's slope of the catch curve's descending part [38]. For the estimation of the natural mortality (M), the following empirical formula proposed by Then et al. (2015) [62] was used: This formula was selected because of its better prediction power when precise estimations of maximum age are unavailable [38]. Fishing mortality (F) and exploitation rates were then estimated from F = Z − M and E = F/Z, respectively. The estimated value of exploitation rate (E) was then compared to the Gulland (1971) [63] proposed a reference value of 0.5, as the upper level of sustainable exploitation. The estimated fishing mortality rate (F) was also compared to the following reference points obtained from the Beverton & Holt's (1956) [64] Yield per Recruit (YPR) model: (a) fishing mortality at maximum yield per recruit (F max ), (b) fishing mortality reduces the population to 50% of unfished spawning biomass (F 0.5 ), and (c) fishing mortality at which the marginal gain in yield per recruit decrease to an arbitrary 10% (F 0.1 ) from that at F = 0.
The parameters of the length-weight relationship ("a" and "b") are the input parameters for the YPR model analysis and estimated from the formula proposed by Le Cren (1951) [65]:

Length-Based Indicators (LBIs)
For the assessment of the stock status relative to exploitation, and to avoid growth and recruitment overfishing, Froese (2004) proposed three simple length-based indicators (LBIs) based on the principle "let them spawn, let them grow and let the mega-spawners live" [2]. These indicators are denoted as P mat , P opt , and P mega . P mat refers to the proportion of mature fish present in the catch with 100% as a target and can be calculated as: Percentage of fish in the catch having a length greater than the length at sexual maturity (L m ) [38]. Froese (2004) set the target reference point for this indicator as; let all (100%) the fishes in the stock spawn at least once in their life before they are caught, which will help rebuild and keep the spawning biomass of the stock healthy [2].
The 2nd indicator (P opt ) is the percentage of fish in the catch within a 10% range around the optimum length (L opt ) (between 0.9L opt and 1.1L opt ) with 100% as the target reference point. Froese (2004) recommended that the size of all (100%) fishes should be within this range in the catch. P opt can be expressed as: P opt = % of fish in the catch between 0.9L opt and 1.1L opt [2]. where, log(L opt ) = 1.053 × log(L m ) − 0.0565 [66]. P mega refers to the percentage of fish in the catch having a length greater than optimum length plus 10% (≥1.1L opt ) [2]. The target of the management strategy is not to exclude the mega-spawners from the catch. If there is no such strategy, about 30-40% of mega-spawner in the catch should be allowed [2]. The three LBIs are further summarized below: where PL is the percentage of fish in the catch in the length interval L. Though the Froese recommended indicators are best fitted for sustainable fishing, knowledge of selectivity patterns, especially in multi-gear fisheries, is essential for interpreting the values. Considering this, Cope and Punt (2009) introduced a new indicator P obj by summing up Froese's indicators, to follow a decision tree established based on a deterministic age-structured population dynamics model to explore the impact of different fishery selectivity patterns on the spawning biomass (SB) of the stock and to examine whether the current SB is at or above the target reference point [21].
The decision tree determines the fishery's selectivity patterns by evaluating Froese indicators, P obj , and the ratio of L m /L opt . Once the fishery's selectivity pattern is confirmed, P mat or P opt is then compared with an empirically established reference point to reveal whether the stock is above or below the SB reference point corresponds to the overfished status of the fishery [21,36].

Fisheries Reference Points from Catch and Resilience
This study also employed a surplus production-based Monte Carlo Method-CMSY to estimate fisheries reference points of the Tenualosa ilisha from catch and prior range of resilience and compared to the above methods [23,47]. Typical surplus production model (e.g., Schaefer model) is required to catch and abundance data to estimate productivity, while CMSY estimates biomass using catch data and productivity [47]. Besides, a Bayesian state-space implementation of the Schaefer production model (BSM) that uses catch and biomass or CPUE data is also included in this method.

(i) Prior r-k Range Determination
The proxy for the resilience of Hilsa in the FishBase is "medium" and translated into prior r-ranges 0.2 to 0.8 [67].
Next, the prior range for 'k' has been estimated using the following equations [49]: where k low and k high are the lower and upper bounds of the prior range of k, max(C) is the highest catch in the time series, r low and r high are the lower and upper bounds of r-values that the CMSY will find out. Equation (9) is applied if the stock's prior biomass is low at the end of the time series, while Equation (10) is applied if the biomass is high [47].

(ii) Estimation of Prior Biomass Ranges
The prior range of relative biomass at the beginning and end of the time series along with an intermediate year estimated automatically by the default rules of this method from the three biomass ranges; low (0.01-0.4), medium (0.2-0.6), and high (0.5-0.9) [47].

(iii) Estimation of Probable Reference Points from Viable r-k Pairs
CMSY selects a random r-k pair from within the prior ranges of r and k for the viable r-k pairs detection. After that, initial biomass is selected from the prior biomass range for the first year and then calculates the predicted biomass for subsequent years from Equation (11) or (12).
B t + 1 is the exploited biomass in the year t + 1, k is the carrying capacity, r is the resilience, and B t and C t are the biomass and catches in the year t. For a severely depleted stock where the biomass falls below 1 4 k, a linear decline of surplus production also incorporated in the prediction formula [47]: where 4B t /k simulates a linear decline of recruitment below half of the biomass capable of yielding MSY [47]. The r-k and predicted biomass trajectory is considered viable if the predicted biomass is not smaller than 0.01 k and falls within the range of prior biomass of the intermediate and final year [47]. From the ranges of "viable" r-k pair, CMSY estimated the most probable values of r and k by the method's default rules. Then the MSY and related fisheries reference points are calculated as MSY = rk/4, fishing mortality corresponding to MSY (F MSY ) = 0.5r, biomass corresponding to MSY (B MSY ) = 0.5k [68,69]. For Bayesian Schaefer Model (BSM) analysis, the uniform r and k ranges used by CMSY are converted into prior densities with reasonable central values [47]. Then, the catchability coefficient q is estimated from the available CPUE data by the following formula: where CPUE t and B t are the mean catch per unit effort and available biomass in year t. For CPUE, the basic dynamic of the corresponding Schaefer model can be shown as [47]: The R code for the CMSY and implementation of the Bayesian Schaefer model was downloaded from http://oceanrep.geomar.de/33076/ (accessed on 3 February 2019).

Results
The total lengths (TL) of 6127 Hilsa specimens, ranging from 12 to 54 cm, were measured. A unimodal length-frequency distribution has been obtained, which is asymmetrical and highly skewed towards the left, in the region of the smaller lengths ( Figure 3). Though capturing Hilsa smaller than 25 cm is illegal (November to June), about 25% of our samples were of length lower than 25 cm. The dominant size range was 21 to 31 cm (Figure 3), and the estimated mean length of the collected specimens was 31.5 cm. There were no significant monthly size differences observed during the sampling period. The parameters of the length-weight relationship, a and b, were obtained from total length (TL) and weight (W) data using the Equation (5), and the relationship was W = 0.023L 2.81 (R 2 = 0.94). catchability coefficient q is estimated from the available CPUE data by the following formula: where CPUEt and Bt are the mean catch per unit effort and available biomass in year t. For CPUE, the basic dynamic of the corresponding Schaefer model can be shown as [47]: The R code for the CMSY and implementation of the Bayesian Schaefer model was downloaded from http://oceanrep.geomar.de/33076/ (accessed on 3 February 2019).

Results
The total lengths (TL) of 6127 Hilsa specimens, ranging from 12 to 54 cm, were measured. A unimodal length-frequency distribution has been obtained, which is asymmetrical and highly skewed towards the left, in the region of the smaller lengths ( Figure 3). Though capturing Hilsa smaller than 25 cm is illegal (November to June), about 25% of our samples were of length lower than 25 cm. The dominant size range was 21 to 31 cm (Figure 3), and the estimated mean length of the collected specimens was 31.5 cm. There were no significant monthly size differences observed during the sampling period. The parameters of the length-weight relationship, a and b, were obtained from total length (TL) and weight (W) data using the Equation (5), and the relationship was W = 0.023L 2.81 (R 2 = 0.94).

Growth Parameters
von Bertalanffy's growth parameters (VBGP) were obtained through the Electronic Length Frequency Analysis (ELEFAN) within the R package TropFishR (Table 1). Estimated growth parameters were 57.84 cm for L∞ and 0.94 year −1 for the growth coefficient K, indicate that Hilsa would reach this length with a moderate growth rate if they were to grow indefinitely. The growth performance index (Φ) was 3.50, with an Rn score value of 0.270.
The automated procedure of TropFishR returned the value for tanchor was 0.40 that refers to 26 May when yearly curves cross length equal to zero. The super-imposed growth curves and the length-frequency data of the Hilsa population of the sampling areas are shown in Figure 4.

Growth Parameters
von Bertalanffy's growth parameters (VBGP) were obtained through the Electronic Length Frequency Analysis (ELEFAN) within the R package TropFishR (Table 1). Estimated growth parameters were 57.84 cm for L ∞ and 0.94 year −1 for the growth coefficient K, indicate that Hilsa would reach this length with a moderate growth rate if they were to grow indefinitely. The growth performance index (Φ) was 3.50, with an Rn score value of 0.270. The automated procedure of TropFishR returned the value for t anchor was 0.40 that refers to 26 May when yearly curves cross length equal to zero. The super-imposed growth curves and the length-frequency data of the Hilsa population of the sampling areas are shown in Figure 4.

Fishing Mortality and Exploitation Rate
For the linearized length-converted catch curve and yield per recruit (YPR) analysis, estimated growth parameters of Hilsa shad were used as input parameters. The estimated values of M, Z, F, E, L c , and biological reference levels of fishing mortality and exploitation rate (F max , F 0.1 , F 0.5 , E max , E 0.1 , and E 0.5 ) are summarized in Table 1. Graphical outputs of the catch curve and YPR model are presented in Figures 5-7. The estimated fishing mortality (F = 2.02 year −1 ) was significantly higher than the natural mortality (M = 1.03) and the current exploitation rate (E) obtained from the F and Z values was at 0.66, which is 32% above the threshold of E opt = 0.5 [65]. More importantly, current exploitation is 52% higher than the required level of exploitation for removing 50% of stock biomass as annual yield (E 0.5 ) ( Table 1), suggesting an over-exploitation that destroys the stock biomass of the fishery rapidly.
Sustainability 2021, 13, x FOR PEER REVIEW For the linearized length-converted catch curve and yield per recruit (YPR) a estimated growth parameters of Hilsa shad were used as input parameters. The es values of M, Z, F, E, Lc, and biological reference levels of fishing mortality and expl rate (Fmax, F0.1, F0.5, Emax, E0.1, and E0.5) are summarized in Table 1. Graphical output catch curve and YPR model are presented in Figures 5-7. The estimated fishing m (F = 2.02 year −1 ) was significantly higher than the natural mortality (M = 1.03) a current exploitation rate (E) obtained from the F and Z values was at 0.66, which above the threshold of Eopt = 0.5 [65]. More importantly, current exploitation is 52% than the required level of exploitation for removing 50% of stock biomass as annu (E0.5) ( Table 1), suggesting an over-exploitation that destroys the stock biomass of t ery rapidly. The probability of capture for different length classes was analyzed, and th FishR package returned an Lc value = 20.48 cm, which indicates that this length the Hilsa population has a 50% probability of being captured in the fishery when t mortality (Z) is 3.05 per year.

Length-Based Indicators (LBIs)
A simple second set of indicators, LBIs, was also used to assess the current stock status of the Hilsa fishery in the water of Bangladesh. Estimation of LBIs required the value of length at first sexual maturity (Lm). Through the gonado-somatic index (GSI) analysis, Islam et al. (1987) reported that the males attain their maturity at the length range of 26-29 cm and the females at 31-33 cm [70]. Because the length-frequency data used in this

Length-Based Indicators (LBIs)
A simple second set of indicators, LBIs, was also used to assess the current s tus of the Hilsa fishery in the water of Bangladesh. Estimation of LBIs required t of length at first sexual maturity (Lm). Through the gonado-somatic index (GSI) a Islam et al. (1987) reported that the males attain their maturity at the length rang 29 cm and the females at 31-33 cm [70]. Because the length-frequency data use The probability of capture for different length classes was analyzed, and the TropFishR package returned an L c value = 20.48 cm, which indicates that this length class in the Hilsa population has a 50% probability of being captured in the fishery when the total mortality (Z) is 3.05 per year.
The YPR model outputs show that maximum yield could be obtained when the

Length-Based Indicators (LBIs)
A simple second set of indicators, LBIs, was also used to assess the current stock status of the Hilsa fishery in the water of Bangladesh. Estimation of LBIs required the value of length at first sexual maturity (L m ). Through the gonado-somatic index (GSI) analysis, Islam et al. (1987) reported that the males attain their maturity at the length range of 26-29 cm and the females at 31-33 cm [70]. Because the length-frequency data used in this study represent both sexes, the length range of 26-33 cm was used as L m for the LBIs analysis to get a comprehensive idea about the species' stock status (Table 2 and Figure 8). Table 2. The proportion of mature fish (P mat ), optimum-sized fish (P opt ), larger than optimum size fish (P mega ), and P obj for collected samples of Tenualosa ilisha and current stock condition interpretation [based on the indicators, and a decision tree proposed by Froese (2004), and Cope and Punt (2009) Figure 8).    [2]. The percentages of mature fish in the catch were range from 25% to 56% (Mean 38%), which means that about two-thirds of the Hilsa were excluded from reproduction. Subsequently, the percentage of optimum-sized fish at which the total biomass of a year class reaches the maximum value ranges from 23% to 41% (mean 33%), demonstrating that target length-classes for catch contribute only one-third of the total catch [66]. In the absence of any size limit for the catch, relatively low percentages (mean 21%) of larger fish (Pmega < 30-40%) could be interpreted as the significant portion of the older or mega spawners has already been excluded from the stock [1,71].
In the decision tree, the sum of all the three proportions (Pmat, Popt, and Pmega), termed as Pobj, provides comprehensive information to distinguish the selectivity patterns. However, after the application of the decision tree for the estimated length-based indicators, the result shows that the stock of the Tenualosa ilisha is likely to be overfished (probabilities  Table 2 illustrates the estimated values of P mat , P opt , and P mega for the catch-based composition of Tenualosa ilisha using different L m values proposed by Islam et al. (1987) [70].
The estimated values of catch-based length proportions (P mat , P opt , and P mega ) are insufficient to support the individual reference target values recommended by Froese (2004) [2]. The percentages of mature fish in the catch were range from 25% to 56% (Mean 38%), which means that about two-thirds of the Hilsa were excluded from reproduction. Subsequently, the percentage of optimum-sized fish at which the total biomass of a year class reaches the maximum value ranges from 23% to 41% (mean 33%), demonstrating that target lengthclasses for catch contribute only one-third of the total catch [66]. In the absence of any size limit for the catch, relatively low percentages (mean 21%) of larger fish (P mega < 30-40%) could be interpreted as the significant portion of the older or mega spawners has already been excluded from the stock [1,71].
In the decision tree, the sum of all the three proportions (P mat , P opt , and P mega ), termed as P obj , provides comprehensive information to distinguish the selectivity patterns. However, after the application of the decision tree for the estimated length-based indicators, the result shows that the stock of the Tenualosa ilisha is likely to be overfished (probabilities of being overfished are ranging from 7-100% for different L m values) with stocking biomass below the target reference point of RP = 0.4 SB for all L m values ( Table 2).

Stock Status and Exploitation from Catch and Resilience
The outputs for the r, k, and biological reference point "MSY" from the CMSY applied to the Hilsa fishery are provided in Table 3. The difference between the values of r produced by CMSY and BSM (0.566 and 0.464) was low (0.1), suggesting that this method can provide a good fit. Moreover, the estimated values of k (2230 × 10 3 t in CMSY and 2264 × 10 3 t in BSM) were also much closer, given that r and k are strongly interrelated. The MSY estimate for the Hilsa shad was estimated at 263 × 10 3 t (227 × 10 3 t -305 × 10 3 t at 95% CI) from BSM and 315 × 10 3 t (226 × 10 3 t -439 × 10 3 t at 95% CI) from CMSY (Table 3).
Based on the above-estimates, BSM also provided the management information for the Hilsa fishery, summarized in Table 4. The predicted trajectories of MSY for Hilsa indicated that sustainable harvest rates were already exceeded by 2004 and peaked in 2018 at around 98% higher than the expected level (Table 4 and Figure 9). The result shows that fishing mortality in 2018 (0.476) was more than two times higher than the mortality required to produce maximum sustainable yield (F MSY = 0.232) ( Table 4). In F/F MSY , a rapid upturn trend was observed, which crossed the sustainable margin in 2004 and peaked in 2018. The exploitation (F/F MSY ) was 2.05 in the last year of the data series, which is more than 100% higher than the maximum allowable limit, indicates the severe overfishing status of the stock. Here, the F value is different from that in Table 1, which may be due to the catchability coefficient "q" (F = qE). However, the ratios of F/F max from length-frequency distribution data and F/F MSY from catch-resilience/catch-effort data both are larger than 1, showing similar results for the overexploitation status of the fishery. set and moved below the threshold level of BMSY in 2018. This indicates the stock's transition phase that will end up with stock's depletion in the future years if current fishing mortality does not go back to FMSY. The B/BMSY trajectory in Figure 9 shows a sharp decline since 1988 and reached the minimum in 2004, followed by a substantial recovery till 2006. After that, the trend is reversed again, and in 2018 it was below the target reference point (B/BMSY = 0.961). The Kobe plot illustrates the simultaneous development of the B/BMSY and F/FMSY (Figure 10). The plot shows the gradual stock depletion from 1988 and moved from a healthy stock with sustainable fishing pressure to a stock that is already depleted by over-fishing. In 2018, the stock moved in the danger red zone where continuous over-fishing will further reduce the stock size that will not produce maximum sustainable yield (MSY). From the Kobe plot, it can be concluded that the stock will collapse shortly if the current overfishing continues. The biomass trend influenced by both fishing mortality and recruitment shows the highest biomass estimate in 1988, followed by a steady decrease until the end of the data set and moved below the threshold level of B MSY in 2018. This indicates the stock's transition phase that will end up with stock's depletion in the future years if current fishing mortality does not go back to F MSY . The B/B MSY trajectory in Figure 9 shows a sharp decline since 1988 and reached the minimum in 2004, followed by a substantial recovery till 2006. After that, the trend is reversed again, and in 2018 it was below the target reference point (B/B MSY = 0.961).
The Kobe plot illustrates the simultaneous development of the B/B MSY and F/F MSY (Figure 10). The plot shows the gradual stock depletion from 1988 and moved from a healthy stock with sustainable fishing pressure to a stock that is already depleted by overfishing. In 2018, the stock moved in the danger red zone where continuous over-fishing will further reduce the stock size that will not produce maximum sustainable yield (MSY). From the Kobe plot, it can be concluded that the stock will collapse shortly if the current overfishing continues.

Discussion
An explicit diagnosis of a fish stock needs information on growth and mortali isting and allowable fishing pressure and exploitation, biological reference point stock biomass and maximum sustainable yield (MSY). This study combined three ent methodological approaches to provide the best estimates of those reference poin the Hilsa stock in the water of Bangladesh (Table 5).

Methods Application
TropFishR Estimation of growth and mortality from length-frequency distrib data.

LBIs
Examine the current status of stock biomass (SB) in relation to target and limit reference points (TRPs and LRPs).

CMSY
Estimation of exploitation and biological reference points from c and resilience.
Before running the TropFishR, the length-weight parameters "a" and "b" wer mated, as these are the input parameters for YPR analysis within the package. The nent "b" value confirms the species' positive isometric growth (b = 2.81). The max . The green area indicates the fishery's safe zone, producing MSY with sustainable fishing pressure and healthy biomass. In the orange zone, the stock biomass is about to be overfished due to high fishing pressure. The red area is the depleted stock biomass that cannot produce MSY due to continuous over-exploitation. And the yellow zone indicates the recovery phase of the stock with reduced fishing pressure.

Discussion
An explicit diagnosis of a fish stock needs information on growth and mortality, existing and allowable fishing pressure and exploitation, biological reference points, i.e., stock biomass and maximum sustainable yield (MSY). This study combined three different methodological approaches to provide the best estimates of those reference points for the Hilsa stock in the water of Bangladesh (Table 5). Table 5. Methodological approaches used in this study.

Methods Application
TropFishR Estimation of growth and mortality from length-frequency distribution data.

LBIs
Examine the current status of stock biomass (SB) in relation to the target and limit reference points (TRPs and LRPs).

CMSY
Estimation of exploitation and biological reference points from catch and resilience.
Before running the TropFishR, the length-weight parameters "a" and "b" were estimated, as these are the input parameters for YPR analysis within the package. The exponent "b" value confirms the species' positive isometric growth (b = 2.81). The maximum and minimum lengths recorded in this study were 54 and 12 cm. Though existing law strictly prohibited the catch of juvenile Hilsa with a length <25 cm [12], the size distribution shows that 25% of the collected specimens are below this length, recognizing the gravity of growth overfishing. Furthermore, the estimated mean length found to be declined from 35 cm in the 1990s [72] to 31.5 cm (in the present study), a manifestation of increased efforts, indicating overfishing. On the other hand, the reported larger lengths (44-54 cm) in the length-frequency data represent only 5% of the total catch, suggesting that this fishery's current selectivity pattern has already excluded the stock's largest size classes. Since the old and larger fishes can produce more offspring and buffer the environment's unfavorable condition [2], their presence in a low quantity in the size structure can intensify fluctuations of the abundance [73].

Stock Condition Analysis Based on TropFishR
Effective management of exploited fish stock requires a clear understanding of the species' life-history parameters [74,75]. This study provides the first assessment of growth parameters of Hilsa shad using TropFishR in the water of Bangladesh. Many stock assessments of Hilsa have been conducted and showed a remarkable variation among the growth parameters (Table 6). Except for the growth coefficient (K), which is much higher than the previous estimates, this study's findings are well within the range of earlier literature estimates (Table 6). Based on the results of different investigations done in the past, Figure 11 shows a consistent increasing trend in growth coefficient (K) and a decreasing trend in asymptotic length (L ∞ ) of Hilsa shad in Bangladesh waters over the last three decades (though they are from different authors but used the similar types of data and similar methods). Although it is not possible to assert a causal link based on those observations, it does appear likely that an increasing trend in K, followed by a decreasing trend in L ∞ , is at least partially responsible for the species' early maturation and thereby reduces the life span [76]. that this fishery's current selectivity pattern has already excluded the stock's largest s classes. Since the old and larger fishes can produce more offspring and buffer the en ronment's unfavorable condition [2], their presence in a low quantity in the size struct can intensify fluctuations of the abundance [73].

Stock Condition Analysis Based on TropFishR
Effective management of exploited fish stock requires a clear understanding of species' life-history parameters [74,75]. This study provides the first assessment of grow parameters of Hilsa shad using TropFishR in the water of Bangladesh. Many stock asse ments of Hilsa have been conducted and showed a remarkable variation among growth parameters (Table 6). Except for the growth coefficient (K), which is much hig than the previous estimates, this study's findings are well within the range of earlier li ature estimates (Table 6). Based on the results of different investigations done in the p Figure 11 shows a consistent increasing trend in growth coefficient (K) and a decreas trend in asymptotic length (L∞) of Hilsa shad in Bangladesh waters over the last th decades (though they are from different authors but used the similar types of data a similar methods). Although it is not possible to assert a causal link based on those ob vations, it does appear likely that an increasing trend in K, followed by a decreasing tre in L∞, is at least partially responsible for the species' early maturation and thereby redu the life span [76].
The estimated growth performance index (Φ) also complies with the previous e mates (Table 6), established the reliability of the growth parameters estimation for Te alosa ilisha from the length-frequency data using TropFishR. The estimates of growth rameters for the Hilsa population from other studies carried out in other countries var widely (Table 6). This is because length-derived growth parameters can be highly in enced by the selectivity of the fisheries [38] and the productivity of the ecosystem [63, The genetic structure difference may also cause variation in the estimation of growth rameters in different regions [63].  The estimated growth performance index (Φ) also complies with the previous estimates (Table 6), established the reliability of the growth parameters estimation for Tenualosa ilisha from the length-frequency data using TropFishR. The estimates of growth parameters for the Hilsa population from other studies carried out in other countries varied widely (Table 6). This is because length-derived growth parameters can be highly influenced by the selectivity of the fisheries [38] and the productivity of the ecosystem [63,77]. The genetic structure difference may also cause variation in the estimation of growth parameters in different regions [63]. Comparative assessment of the current and desired level of fishing mortality and exploitation is necessary to develop a sustainable management strategy for an exploited fish stock. Though many matrices regulate the fishery's status for the operational purpose, these two indices have been widely used to formulate management action for obtaining the maximum sustainable yield (MSY) and keep the fishing pressure at the sustainable level to maintain the productivity of the stock [82]. It is also necessary to evaluate the gears' performance and associate selectivity and how the fishery responds to the size limitation regulations. The prediction model (YPR) used in this analysis allows us to evaluate the status of the stock in relation to different reference levels [45] and to infer control measures such as controlling fishing effort and regulating gear selectivity (Figure 7). The outputs of the model suggest that the Hilsa stock in the water of Bangladesh has been disrupted with high fishing pressure (F = 2.02 year −1 ) and a high exploitation rate (E = 0.66 year −1 ) that removes the stock biomass more than one and a half time higher than the sustainable level (E 0.5 ) ( Table 1). The estimated size at first capture (L c = 20.48 cm) is much smaller than the minimum allowable length for the catch (25 cm) set by the authority. Though some management measures have been taken [12], previous studies show that this fishery has been facing continuous overfishing and overexploitation since 1992 ( Table 6). The first capture (L c ) length was also decreasing since 1997 and crossed the allowable length in 1999 (Table 6).
Though the impact of high fishing pressure could be reduced by letting the new spawners produce at least one replacement spawners [2], an exploitation level higher than the theoretical optimum level (E = 0.5) with a reduced size in the first capture indicates exclusion of juvenile fish from the stock. Islam et al. (1987) [70] reported the length at sexual maturity of Hilsa shad is between 26 cm to 29 cm for males and 31 cm to 33 cm for females. However, unfortunately, present and previous studies show that the first spawners of the stock have been targeted by the fishery since 1999 (Table 6). In return, this will deplete the stock spawning biomass and leads the fishery towards recruitment failure and thereby collapse.

Stock Condition Analysis Based on Length Based Indicators (LBIs)
Although theoretically informative, the previous method based on fitting population dynamics models to data is less capable of producing the stock's real scenario for a data-limited fishery like Hilsa [21]. Froese (2004) proposed an alternative assessment method based on length-based indicators to determine whether the stock is being harvested sustainably [2]. Estimating these indicators is straightforward using length-frequency data and easily understandable to all concerned stakeholders, including the general public [2]; therefore, their participation can be ensured for effective fisheries management.
According to the indicators, the proportion of mature fish (P mat ) in the catch should be as high as possible to avoid recruitment overfishing, and all catches should be maintained within the 10% range of optimum length (L opt ) to prevent growth overfishing. However, the findings of this study show a different picture of the stock. The mean value of P mat (38%) and P opt (32%) suggesting that the Hilsa stock is experiencing both recruitment and growth overfishing. On the other hand, older and larger fishes (mega-spawners) are essential for the fishery because of their high fecund rate and ability to extend longevity and prolong the reproductive cycle. Those, the mega-spawners, reduce the probability of recruitment failure [83,84] and play a significant role in the stock's long-term survivability [85]. Assuming that there is no upper limit to the size of fish caught in the stock, Froese recommended not to capture more than 30-40% of the mega-spawner to support his manifesto "Let the mega-spawner live" [2]. The estimated mean value of P mega in this study is 21 percent, and in the absence of any size regulations, this low percentage of mega-spawners indicates that this fishery is already excluded a significant portion of mega-spawners from the stock. Furthermore, personal communication with the fishermen confirms that they never released the larger fishes they caught. From this point of view, if the removal of mega-spawners from the stock will continue with the current rate, the stock's resilience will be reduced, which poses a threat to the stock to be depleted.
The application of the Cope and Punt decision tree based on Froese indicators shows that the spawning stock biomass of the Hilsa fishery is below the target reference points. This confirms the fishery's overfished status under the current pattern of selectivity [2,21].

Stock Status Analysis Based on CMSY
Though TropFishR and LBIs have given a comprehensive idea on the current exploitation and stock biomass status of the fishery, a Monte Carlo method-CMSY has also been employed in this study to evaluate the stock status using catch-resilience and catch-effort data and to confirm the robustness of the estimations of previous methods [23,47]. Giving the limited requirement of input data, the estimates of CMSY surprisingly comply with the estimates estimated by the previous two different methods. From the result produced in the CMSY and BSM analysis, the intrinsic growth rate of the Hilsa population in Bangladesh water shows the stock's ability to replenish its 47 to 57 percent every year. Hence, the total mortality of the stock should not exceed this limit. However, the present high fishing mortality and exploitation rate estimated by this method show the alarming condition of the stock.
As a consequence of overfishing and overexploitation rate, the biomass of Hilsa is also found to go below the reference point (B/B MSY = 0.961) in recent years and shows a continuous downward trend. The Kobe plot ( Figure 10) presents the stock biomass movement from the safe zone (green) to the highly exploited zone (red). Because of the moderately high resilience, the stock biomass remains close to the reference point, but if the overfishing and overexploitation are let to continue, this fishery will not produce MSY in the future years.
Although several studies have proposed the MSY for Hilsa stock in Bangladesh's water using Gulland methods (Table 7), a discrepancy exists between those and the present study. The method used in the previous studies using length-frequency data may not produce robust estimates of MSY when the fishery is characterized by a simultaneous increasing trend in the catch and CPUE indices [48]. On the other hand, the Monte Carlo catch-only method (CMSY) has shown its robustness in previous studies [47,48] and reviews [86,87]. The estimated MSY in CMSY is about 64% lower than the fish landed in 2018, which shows how extremely the stock is being disrupted. Therefore, to save the stock, immediate protective management measures are required.

Management Recommendations
From the analysis of TropFishR and CMSY, the present fishing mortality rate of Hilsa appears to be high, which is almost double for CMSY than the expected level. Again, LBIs and CMSY analysis show that the stock biomass of the Hilsa stock in the water of Bangladesh is already overfished. Therefore, reducing fishing mortality is advisable, but save the stock spawning biomass from depletion should be the highest management priority. Imposing size limits for the catch would help maintain the stock biomass and reduce fishing mortality. Though the catches of juveniles smaller than 25 cm in length are prohibited [12], this length is still smaller than the length of first sexual maturity (L m ). Furthermore, the estimated length at first capture in the present and previous studies ( Table 6) are smaller than 25 cm, which indicates the absence of proper management.
To maintain the stock biomass and the fishing mortality at a sustainable margin, a smaller length limit of 27.50(±1.3) cm for males and 32(±1) cm for the female Hilsa species was calculated in this study. These are the mean of the lengths at first sexual maturity (L m ) for both sexes [70]. However, considering the difficulties maintaining different selectivity of the same gear for a fishery, particularly where resources are limited for proper monitoring, this study further recommends a standard smaller length limit for the catch, which is 33 cm. This length is the highest reported length at first maturity for females, and the ban of catching fish smaller or equal to this length will allow all immature fish, both male, and female, to spawn and grow to their optimum size.
Again, to recommend the upper length limit for the catch, this study calculated the length, 39 cm ("L opt + 10% of L opt " when L m is 33 cm). However, the proportion of fish greater than this length in the catch is only 9 percent, presenting a frightful stock scenario. Considering this, the mega-spawner catch should ban immediately. However, if the first recommendation can be executed, then instead of a ban, a motivational campaign can be launched to motivate the fishermen not to catch the mega-spawners. Allowing the immature fish to reach their optimum length will compensate for the catches of megaspawners. After implementing this strategy, the annual yield would have been lower in the first few years but would reach its maximum after those periods [2]. The existing policy of incentivizing the fishing communities during the fishing ban seasons [12] can be extended for those few years to compensate for the loss of fishers' earnings.
To reduce the fishing pressure and bring the fishing mortality to the F MSY level, regulation of effort controls with the gear's nature specification is essential. In Bangladesh, fishing in rivers and marine water is regulated by the "The Protection and Conservation of Fish Act of 1950" and "Marine Fisheries Act, 2020 (previously "Marine fisheries ordinance, 1983"). These rules do not allow a gill net with a mesh size less than 4.5 cm. However, gill net with a mesh size even larger than 4.5 cm would catch fish smaller than 33 cm [28]. According to Rahman and Cowx (2008), the optimum mesh size to catch fishes with a length of 33 cm or bigger is 8.0 cm [28]. So, gill nets with a mesh size smaller than 8.0 cm should be strictly prohibited for Hilsa catch, and implementation of this rule immediately is urgent to save the stock from further depletion. Another issue is licensing. Anyone fishing for Hilsa with nets and boats must have a license, and the number of licenses should be limited. Though licensing of boats is already executing for the marine fisheries but absent in the inland water bodies. From the analysis of this study, it is recommended that the existing number of boats should be reduced to its half to get the fishery free from high fishing pressure and maintain the yearly landing limit within the range of 263,000 to 315,000 tons.
Due to the weak monitoring and management, illegal use of nets (mesh size < recommended size) and unlawful fishing during the ban periods cause the depletion of this fishery rapidly. A formal monitoring and evaluation mechanism should be established for regular monitoring and evaluation of the management process.

Recommendations for Future Research
The length at first sexual maturity (L m ) used in this study was estimated long before (1987), so further research to estimate L m would significantly improve the stock status estimation.
Continuous collection of length-frequency and marine and inland catch-effort data is necessary. This will help to understand the changes in the stock using the methods used in this study.
Finally, an investigation of the impacts of climate change on the biology and ecology of the Hilsa shad in Bangladesh's water could be added to a new dimension in the process of stock status evaluation.

Conclusions
In this study, a complete and rigorous stock assessment of Hilsa shad, the most important fishery of Bangladesh, has been done using three different methods. Along with a maximum sustainable yield (MSY), the outputs indicated that the fishery is in an overfished state due to overfishing. After analyzing the results, this study concluded with the following recommendations: a. Maintain the length at first capture beyond the first maturity level and close to or within the optimum length where maximum potential yield could be obtained.
b. The recommended minimum length for the catch is >33 cm. For the mega-spawners' protection, fish with a size >39 cm should not include in the catch. And c. Immediate control measures should be taken to bring the fishing mortality back to the F MSY and establish the yearly catch landing limit of Hilsa within the range of 263,000 to 315,000 tons.

Informed Consent Statement: Not applicable.
Data Availability Statement: The Datasets generated during this study are available from the corresponding author on reasonable request.