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

Abstract

The anadromous tropical Hilsa shad formed the largest single-species fishery in Bangladesh, making the highest contribution to the country’s total fish production (14%) and nearly 83% of the global Hilsa catch in 2018. However, increased fishing pressure made the fishery 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 fishing mortality and exploitation, Froese’s length-based indicators for fishing sustainability, and a surplus production-based Monte Carlo method-CMSY, for fisheries reference points estimation) for the best possible estimates of the Hilsa stock status in the water of Bangladesh. The present findings revealed that the stock is likely to be overfished due to over-exploitation. 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 fishery in Bangladesh.

1. 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 river and tributary along the coast during the spawning seasons [12,13,14,15]. This species has significant market demand, especially in South and Southeast Asia, and the global annual production is more than 0.60 million tons [12,16]. Though Hilsa shows a greater distribution, the northern Bay of Bengal regions (India, Myanmar, and Bangladesh) mainly support this fishery contributing approximately 99 percent of the world production, and Bangladesh contributes nearly 83 percent of reported global Hilsa catches in 2018. This fishery is recognized as the largest single-species fishery in Bangladesh, making the highest contribution, about 14 percent, of the country’s total fish production and 47 percent of the total marine catch [17]. Thus, this fishery became the “social and economic engine”, providing food and livelihoods of about three million people, directly and indirectly [18,19]. The average annual Hilsa production in Bangladesh is about 0.45 million tons, representing more than one percent of its total Gross Domestic Product (GDP) [11,17]. The growing market demand of Hilsa shad in South and Southeast Asia, the Middle East, Japan, North America, Australia, and Europe has stimulated the increased fishing effort and makes the fishery vulnerable. Considering this, knowledge about the stock condition and its response to the current degree of removal is essential for both users and managers to explore the future potential for sustainable exploitation of the fishery [2,20,21,22,23,24,25].

Overexploitation and dominance of first-spawners in the catch have been reported by several studies that carried out the stock assessment of Hilsa shad from the waters of Bangladesh [25,26,27,28,29,30,31,32,33,34,35], which could be the warning sign of growth and recruitment overfishing of the species [2,36,37,38]. These studies were carried out using only length-frequency data in FAO-ICLARM stock assessment tools (FiSAT), aiming to explore the stock exploitation status, the spawning stock biomass was neglected. Furthermore, this fishery is regulated mainly by mesh size regulation of nets, closure of spawning ground, and banning fishing during spawning seasons [2,39,40,41,42]. However, banning fishing only in the spawning seasons, not before or after the spawning seasons, did not prove to be effective [43]. Therefore, there is no annual fish landing limit, an essential tool for sustainable fisheries management.

Considering the limitations of previous studies and the fishery’s data-limited condition, this research carried out an assessment based on three sets of indicators to provide the best possible estimates of Hilsa stock status in Bangladesh’s water, with fisheries representative length-frequency and catch and resilience data. For the first set of indicators (i.e., growth, mortality, and fisheries reference points), FAO’s traditional stock assessment methods [44] were used. However, to enhance the robustness of the growth parameters estimation, new procedures have been incorporated using a recently developed R package, “TropFishR” [38,45,46]. The second set of indicators based on catch proportions related to length-reference points [2,21] was used to monitor the species’ population status related to exploitation. These indicators provide the basis for avoiding growth and recruitment overfishing and distinguishing selectivity patterns to develop harvest control rules [2]. Finally, a surplus production-based Monte Carlo method (CMSY) has been used to evaluate the stock biomass with the exploitation and estimate fisheries reference points based on catch series and resilience data [47]. This method provides the most reliable and robust estimates for a fishery with increasing catch and catch per unit effort (CPUE) [48]. This study concludes by highlighting the knowledge gaps and presenting viable management recommendations based on the results and lessons learned from evaluating existing management plans.

From the genetic and otolith microchemistry analysis, previous studies [49,50] showed that Hilsa in the northern Bay of Bengal region belongs in a single population and should be managed as a single stock [51]. Therefore, this study’s recommendations based on the length-frequency analysis (first two methods) could also be applicable for managing the fishery in other countries in this region (South-East India and Myanmar).

2. Materials and Methods

2.1. 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., 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°59′59″ N, 90°39′52″ E]; Chandpur [23°13′48″ N, 90°38′32″ E]; Ramgati [22°34′55″ N, 90°59′59″ E] and Charfassion, Bhola [22°11′09″ N, 90°45′48″ E]), estuarine (Jahajmara, Hatiya [22°09′14″ N, 91°02′57″ E]; Ramgati and Charfassion, Bhola), and marine habitats (BFDC Fisheries Ghat, Cox’s Bazar [21°27′06″ N, 91°58′05″ E]; Badarkhali, Chakaria [21°43′09″ N, 91°57′00″ E]; Baro Ghop, Kutubdia [21°49′00″ N, 91°50′58″ E]; Fisheries Ghat, Chittagong [22°19′42″ N, 91°50′48″ E]; Bangla Bazar, Sandwip [22°24′38″ N, 91°28′48″ E]; Kuakata [21°48′53″ N, 90°07′23″ E] and BFDC Fisheries Ghat, Barguna [22°03′07″ N, 89°58′18″ E]) 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 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.

2.2. 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).

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 average catch per boat per annum in marine water. Data interpolation technique was employed to avoid the anomalies in effort data.

2.3. 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].

2.3.1. Estimation of Growth Parameters

The von Bertalanffy (1938) growth parameters (K, L_∞, and t₀) 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 = L_∞(1 − e − (K(t − t₀) + S(t) − S(t₀)))

(1)

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⁻¹, t₀ 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₀)” 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]:

Φ′ = log(K) + 2log(L_∞)

(2)

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]:

L_∞ = L_max/0.95

(3)

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 Taylor & Mildenberger (2017) [59].

2.3.2. 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:

M = 4.118K^0.73L_∞^−0.33

(4)

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]:

W = aL^b

(5)

2.3.3. 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:

$$P_{mat}={\displaystyle \sum _{L_{mat}}^{L_{\max }}{P}_L}$$

(6)

$$P_{opt}={\displaystyle \sum _{0.9L_{opt}}^{1.1L_{opt}}{P}_L}$$

(7)

$$P_{mega}={\displaystyle \sum _{1.1L_{opt}}^{L_{\max }}{P}_L}$$

(8)

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].

2.3.4. 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.

CMSY Analysis

(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]:

k_low = max(C)/r_high,  k_high = 4max(C)/r_low

(9)

k_low = 2max(C)/r_high,  k_high = 12max(C)/r_low

(10)

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} = B_t + rB_t(1 − B_t/k) − C_t

(11)

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 ¼ k, a linear decline of surplus production also incorporated in the prediction formula [47]:

B_{t+ 1} = B_t + 4B_t/kr(1 − B_t/k)B_t − C_t, when B_t/k < 0.25

(12)

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].

Bayesian Schaefer Model (BSM)

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:

CPUE_t = qB_t

(13)

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]:

CPUE_t+1 = CPUE_t + r(1 − CPUE_t/qk)CPUE_t − qC_t

(14)

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).

3. 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² = 0.94).

3.1. Growth Parameters

von Bertalanffy’s growth parameters (VBGP) were obtained through the Electronic Length Frequency Analysis (ELEFAN) within the R package TropFishR (

Table 1. Growth parameters, mortality, and exploitation rates of Hilsa shad estimated by TropFishR.

Parameters	Value	Comments
Growth parameters
Asymptotic Length (L∞)	57.84 cm	If there was no fishing pressure/predation/starvation/disease, Hilsa would reach this length.
Growth Coefficient (K)	0.94 year−1	Moderately higher growth rate.
tanchor	0.40	Yearly repeating growth curves cross length equal to zero on 07 June.
Growth Performance Index (Φ)	3.50	High
Mortality and Exploitation
Natural Mortality (M)	1.03 year−1	Natural mortality is moderately high.
Total Mortality (Z)	3.05 year−1	Comparatively high.
Fishing Mortality (Fcurr)	2.02 year−1	Significantly higher than natural mortality.
Fmax	1.70 year−1	Fishing mortality at maximum yield per recruit.
F0.1	1.00 year−1	Fishing mortality at which the marginal gain in yield per recruit decreases to an arbitrary 10% from that at F = 0.
F0.5	0.70 year−1	Fishing mortality reduces the population to 50% of unfished spawning biomass.
Current Exploitation (Ecurr)	0.66	32% higher than the threshold level, indicating overexploitation.
Emax	0.62	Maximum exploitation level to obtain optimum yield.
E0.1	0.49	Maximum exploitation level to obtain biologically optimum yield.
E0.5	0.40	Maximum exploitation level to obtain 50% of the biomass as annual yield.
Length at First Capture (Lc)	20.48 cm	The length class in the population has a 50% probability of being captured.

). Estimated growth parameters were 57.84 cm for L_∞ and 0.94 year⁻¹ 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.

3.2. 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 Figure 5, Figure 6 and Figure 7. The estimated fishing mortality (F = 2.02 year⁻¹) 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.

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 F_max and E_max values are at 1.70 year⁻¹ and 0.62, respectively. Similarly, 50% of the biomass could be obtained as the annual yield when F_0.5 and E_0.5 are at 0.70 year⁻¹ and 0.40. The biologically optimum yield could be obtained when F_0.1 = 1.00 year⁻¹ and E_0.1 = 0.49, respectively.

3.3. 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. 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) respectively (See Methods)].

Lm	Lopt	Pmat	Popt	Pmega	Pobj	Stock Condition Interpretation	Probability
26	27.13	0.56	0.41	0.33	1.30	SB < RP	100%
27	28.23	0.50	0.39	0.29	1.18	SB < RP	100%
28	29.33	0.45	0.37	0.27	1.09	SB < RP	100%
29	30.44	0.40	0.32	0.25	1.00	SB < RP	100%
30	31.54	0.32	0.31	0.18	0.83	SB < RP	7%
31	32.65	0.29	0.29	0.15	0.73	SB < RP	7%
32	33.76	0.26	0.26	0.13	0.65	SB < RP	19%
33	34.87	0.25	0.23	0.09	0.57	SB < RP	19%

and Figure 8).

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 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 (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).

3.4. 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. Estimates and the 95% confident intervals of r, k, and MSY of Tenualosa ilisha in CMSY.

Models Output in CMSY	r (1/Year)	k/103 MT	MSY/103 MT
CMSY	0.566 0.407–0.785	2230 1357–3665	315 226–439
BSM	0.464 0.326–0.662	2264 1650–3107	263 227–305

.

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³t in CMSY and 2264 × 10³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³t (227 × 10³t – 305 × 10³t at 95% CI) from BSM and 315 × 10³t (226 × 10³t – 439 × 10³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. Result for management based on BSM model analysis in the CMSY method (Units for yield and biomass are ×10³ MT).

Parameters	Value	95% CI
FMSY	0.232	0.163–0.331
MSY	263	227–305
BMSY	1132	825–1554
		2.5th–97.5th percentile
Biomass in the last year	1088	834–1354
B/BMSY in the last year	0.961	0.736–1.20
Fishing mortality in the last year	0.476	0.382–0.620
Exploitation F/FMSY	2.05	1.95–2.67

. 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.

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 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.

4. 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. 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].

4.1. 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].

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].

Table 6. Estimated population parameters of Hilsa in different water bodies of the world.

Regions	L∞ (cm)	K	Φ	M	Z	F	E	Emax	Lc	Year	Reference
Bangladesh waters	58.70	0.90	3.50	1.36	4.19	2.83	0.68	0.66	27.00	2015	[11]
	53.70	0.86	3.40	1.36	3.51	2.16	0.61	0.58	19.87	2002	[20]
	61.10	0.74	-	1.16	2.41	1.25	0.52	-	35.00	1992	[25]
	58.30	0.74	3.40	1.18	2.61	1.43	0.55	-	30.00	1995
	61.50	0.83	3.50	1.28	3.29	2.01	0.61	0.69	30.30	1997
	66.00	0.67	3.46	1.25	3.43	2.18	0.63	0.60	27.06	1998
	60.00	0.82	3.47	1.28	3.77	2.49	0.66	0.59	22.80	1999
	62.50	0.72	3.45	1.17	2.79	1.62	0.58	0.46	13.12	2000
	52.00	0.71	3.28	1.22	2.61	1.39	0.53	-	-	2006	[31]
	53.55	0.61	3.24	1.10	2.83	1.73	0.61	-	-	2003	[78]
	57.84	0.94	3.50	1.03	3.05	2.02	0.66	0.62	20.48	2018	Present study
North Arabian Gulf, Kuwait	52.50	0.36	3.00	0.40	1.20	0.80	0.67	-	-	1992	[8]
Bay of Bengal, India	47.80	1.90	3.64	1.25	1.98	0.73	0.37	0.56	26.56	2010	[32]
Indus River, Pakistan	31.50	1.50	2.13	2.21	2.89	0.67	0.23	1.00	-	2012	[79]
North-west Arabian Gulf, Iraq	61.5	0.28	3.02	0.55	1.66	1.11	0.67	0.72	27.80	2013	[80]
North-west Arabian Gulf, Iran	42.20	0.78	3.16	1.29	4.53	3.24	0.72	-	22.30	2008	[81]

Table 6. Estimated population parameters of Hilsa in different water bodies of the world.

Regions	L∞ (cm)	K	Φ	M	Z	F	E	Emax	Lc	Year	Reference
Bangladesh waters	58.70	0.90	3.50	1.36	4.19	2.83	0.68	0.66	27.00	2015	[11]
	53.70	0.86	3.40	1.36	3.51	2.16	0.61	0.58	19.87	2002	[20]
	61.10	0.74	-	1.16	2.41	1.25	0.52	-	35.00	1992	[25]
	58.30	0.74	3.40	1.18	2.61	1.43	0.55	-	30.00	1995
	61.50	0.83	3.50	1.28	3.29	2.01	0.61	0.69	30.30	1997
	66.00	0.67	3.46	1.25	3.43	2.18	0.63	0.60	27.06	1998
	60.00	0.82	3.47	1.28	3.77	2.49	0.66	0.59	22.80	1999
	62.50	0.72	3.45	1.17	2.79	1.62	0.58	0.46	13.12	2000
	52.00	0.71	3.28	1.22	2.61	1.39	0.53	-	-	2006	[31]
	53.55	0.61	3.24	1.10	2.83	1.73	0.61	-	-	2003	[78]
	57.84	0.94	3.50	1.03	3.05	2.02	0.66	0.62	20.48	2018	Present study
North Arabian Gulf, Kuwait	52.50	0.36	3.00	0.40	1.20	0.80	0.67	-	-	1992	[8]
Bay of Bengal, India	47.80	1.90	3.64	1.25	1.98	0.73	0.37	0.56	26.56	2010	[32]
Indus River, Pakistan	31.50	1.50	2.13	2.21	2.89	0.67	0.23	1.00	-	2012	[79]
North-west Arabian Gulf, Iraq	61.5	0.28	3.02	0.55	1.66	1.11	0.67	0.72	27.80	2013	[80]
North-west Arabian Gulf, Iran	42.20	0.78	3.16	1.29	4.53	3.24	0.72	-	22.30	2008	[81]

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⁻¹) and a high exploitation rate (E = 0.66 year⁻¹) 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.

4.2. 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].

4.3. 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. Summary of the MSY of the Hilsa shad in the waters Bangladesh estimated by different authors.

Studies	Methods	Year	Data Type	MSY (Tons)
Rahman et al. (2018) [11]	Gulland method	2015–2016	Length-frequency	526,000
Amin et al. (2002) [29]	Gulland method †	1999	Length-frequency	162,396
Halder & Amin (2005) [84]	Gulland method	2002	Length-frequency	235,130
Amin et al. (2008) [85]	Gulland method	1980–2000	Length-frequency	210,125
Present Study	CMSY	2017–2018	Catch and resilience	315,000
	BSM		Catch and effort	263,000

† Gulland method is MSY = 0.5 × Z_t × B_t, where Z_t and B_t are the current total mortality rate and biomass, respectively [63].

), 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.

4.4. 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 mega-spawners. 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.

4.5. 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.

5. 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.