Assessing Five Major Exploited Tuna Species in India (Eastern and Western Indian Ocean) Using the Monte Carlo Method (CMSY) and the Bayesian Schaefer Model (BSM)

: The status of data-limited tuna ﬁshery stocks in India has been tested using the latest and most advanced computerized methods, CMSY and BSM. Five tuna ﬁsh stocks from both the Eastern and Western Indian Ocean were assessed using both catch and catch per unit effort (CPUE) details available from 1990 to 2015. Both methods help to calculate the maximum sustainable yield (MSY) and exploitation of MSY relative to biomass (B/B MSY ). The results of maximum intrinsic rate (r) and carrying capacity are also estimated. The results revealed that all tuna stocks in both the regions were overﬁshed, with one, the longtail tuna ( Thunnus tonggol ) in the Western Indian Ocean strongly overﬁshed (B/B MSY = 0.44). Such observations, although still preliminary since the techniques used to produce them are relatively new, often associated with the situation and exploitation of all the stock in question, making the CMSY and BSM methods promising for stock assessment in data-deﬁcit situations. The study concludes that in order to restore the status of these ﬁve tuna stocks in both regions, it would be necessary to reduce the ﬁshing pressure.


Introduction
The data of fisheries obtained from the Food and Agriculture Organization (FAO) revealed that global fish landings have been declining steadily since 1996 [1], with catches declining by 1.2 million tons annually since then due to illegal, unreported fishing and undocumented (IUU) fishing [2,3]. The proportion of fish stocks within biologically sustainable levels, which was 90% in 1990 and 65.8% in 2017, can be used to understand the overfishing status of global fisheries [4]. For Asian fisheries, the situation is even worse [5] due to high fishing pressure on the stock [6]. The highest depletion has been seen for predatory fishes with only 10 percent biomass available of its actual population [7]. Now, in order to reduce the fishing pressure of a particular fish stock, total allowable catch (TAC) should be reduced [8].
Different species of tuna are distributed widely all around the globe but are specifically found from temperate to tropical waters around the equator between 45 degrees North and South. The species are widely divided into neritic, coastal and marine tunas. Taxonomically tuna's are grouped into the Scrombridae family, approximately consisting of fifty species, constitute the third-largest international seafood trade commodity and accounting for almost 10% of total value trade [9]. Tuna catches continued to rise, hitting their highest level at about 7.9 million tons in 2018, which was possibly due to catches from the Central and Western Pacific [9]. Within the species group, the highest catch was dominated by CMSY and BSM computerized methods were used and compared for five major tuna species in India, from the Eastern Indian Ocean and Western Indian Ocean regions. The exploitation status of the stocks from the Eastern Indian Ocean included the regions of Andaman Islands, Lakshadweep Islands, Nicobar Islands, Pamban Islands and St. Mary's Islands. Western Indian Ocean's stock assessment comprised of Vypin Islands, Vallarpadam and Willingdon Islands. The focus of this research is to determine the MSY, CMSY estimated biomass, exploitation rate and other BRPs of five major tuna species in the Indian Ocean using CMSY and BSM models, and provide baseline stock information to the fishery administrator for sustainable tuna management. The significance of this research is that unlike complex stock assessment methods that require fisheries-independent data sets such as research surveys and catch-at-age data, the recently proposed CMSY method quantifies biomass, exploitation rate, Maximum Sustainable Yield (MSY) and related fisheries reference points for a given population using only a time series of catches and ancillary qualitative information.

Materials and Methods
The study is based on assessment of five major tuna fishery stocks from India in the Eastern Indian Ocean and Western Indian Ocean . The catch and effort data for the study were accessed from FAO statistical database and Indian Ocean Tuna Commission (IOTC), respectively. The data source and basic information of the five tuna stocks covered here are given in Table 1. Each the five tuna stocks with catch and catch per unit effort (CPUE) data series of the above-mentioned period was analyzed using the CMSY and BSM methods.
It is possible to see the CMSY approach as an accounting practice, where annual catches are removed from stock biomass by an annual biological product that is included in it and is usually initiated at the expected biomass. This leads to a biomass pathway that is defined by the stock and the growth rate of the stock, as determined by the estimated population growth (r; 'resilience') and the extend of similarity between the biomass and carrying capacity (k). A feature of the Monte Carlo method is then used to produce thousands of such trajectories, each with different values of r and k, and to separate trails with usable solutions, i.e., trajectories can 'crash' the stock and follow different ('priors') methods To determine status of the fishery, the CMSY method was applied to stocks in the Northwest Pacific [25]. Additionally, part of the Bayesian Schaefer State-Space Surplus Production Model (BSM) code for CMSY R was used to respond to a variety of stock assessment (process error) measurements [25]. JAGS program with the Markov chain Monte CMSY and BSM computerized methods were used and compared for five major tuna species in India, from the Eastern Indian Ocean and Western Indian Ocean regions. The exploitation status of the stocks from the Eastern Indian Ocean included the regions of Andaman Islands, Lakshadweep Islands, Nicobar Islands, Pamban Islands and St. Mary's Islands. Western Indian Ocean's stock assessment comprised of Vypin Islands, Vallarpadam and Willingdon Islands. The focus of this research is to determine the MSY, CMSY estimated biomass, exploitation rate and other BRPs of five major tuna species in the Indian Ocean using CMSY and BSM models, and provide baseline stock information to the fishery administrator for sustainable tuna management. The significance of this research is that unlike complex stock assessment methods that require fisheries-independent data sets such as research surveys and catch-at-age data, the recently proposed CMSY method quantifies biomass, exploitation rate, Maximum Sustainable Yield (MSY) and related fisheries reference points for a given population using only a time series of catches and ancillary qualitative information.

Materials and Methods
The study is based on assessment of five major tuna fishery stocks from India in the Eastern Indian Ocean and Western Indian Ocean . The catch and effort data for the study were accessed from FAO statistical database and Indian Ocean Tuna Commission (IOTC), respectively. The data source and basic information of the five tuna stocks covered here are given in Table 1. Each the five tuna stocks with catch and catch per unit effort (CPUE) data series of the above-mentioned period was analyzed using the CMSY and BSM methods.
It is possible to see the CMSY approach as an accounting practice, where annual catches are removed from stock biomass by an annual biological product that is included in it and is usually initiated at the expected biomass. This leads to a biomass pathway that is defined by the stock and the growth rate of the stock, as determined by the estimated population growth (r; 'resilience') and the extend of similarity between the biomass and carrying capacity (k). A feature of the Monte Carlo method is then used to produce thousands of such trajectories, each with different values of r and k, and to separate trails with usable solutions, i.e., trajectories can 'crash' the stock and follow different ('priors') methods To determine status of the fishery, the CMSY method was applied to stocks in the Northwest Pacific [25]. Additionally, part of the Bayesian Schaefer State-Space Surplus Production Model (BSM) code for CMSY R was used to respond to a variety of stock assessment (process error) measurements [25]. JAGS program with the Markov chain Monte Carlo process sampled the probability distributions of the parameters [26,27] [12] CMSY and BSM methods measure parameters that include MSY, biomass that enables a fish stock to deliver the maximum sustainable yield (B MSY ) and the maximum rate of fishing mortality (F MSY ) based on r-k pairs that may be screened for Monte Carlo tests [25].
By using the catch and abundance data, r, k and MSY can be predicted by using the BSM model and are represented by Equation (1).
where B t+1 is the biomass used next year t + 1, B t is the present biomass and C t is the actual catch in year t, r is intrinsic population rate and k is the carrying capacity. In order to understand the reduced recruitment of the highly declining stock, surplus production (linear decline) [28] that is a function of somatic growth, natural mortality and recruitments is applied if biomass declines beyond 0.25 k, as shown in Equation (2).
where the term 4rB t /k assumes that r declines linearly with half of biomass that is actually capable of producing MSY. Table 2 represents the resilience information (FishBase [29] and SeaLifeBase [30] that is translated into prior r ranges. The ratio of higher and lower r and maximum catch (C) were used to enlighten the upper and lower k bounds. Table 2. Suggested resilience categories translated into range of rate of population increase (r) as provided in FishBase [29] and SeaLifeBase [30] of stocks under assessment.  (3) represented stock, which at the end accounted for lower initial biomass, which was later modified for higher biomass.

Resilience Category
Values of k are predicted to have a normal log distribution, a standard deviation that was thought to be 1 4 of the distance between the median value and the lower range of k [31]. Data-scarce stocks might have CPUE (number/hour) as an estimate of stock abundance for some years. The catchability coefficient (q) relates this abundance index to the stock biomass, as shown in Equation (5).
where CPUE t represents the average catch per unit effort, B t is biomass in year t and q represents the catchability coefficient. Catch per unit effort's dynamic abundance is expressed in Equation (6).
Equations (1) and (5) represent the variables and parameters, and the prior catchability coefficient is obtained from Equation (7).
By using the simulations, the multipliers for prior r ranges and high and low biomass are derived. Now, by using Equations (8) and (9), prior catchability coefficient (q) can be derived for fish stocks with high B. q low = 0.25 r pgm CPUE mean /C mean (8) q high = 0.5 r high CPUE mean /C mean (9) where q low is the lower prior catchability coefficient and q high is the upper catchability coefficient, both for the stocks with higher present biomass. The geometric mean of r is represented by r pgm , CPUE mean represents the average CPUE of the last half-decade or full decade and C mean represents the average catch. Pauly [32] suggested changing the multiplier values of q low and q high to 0.5 and 1.0, respectively, for the actual stocks that presently have low biomass. For the species with medium and high resilience, the mean catch and CPUE were applied over last 5 years. For the species with low to very low r values, the same were applied, but for recent decadal years [33].
Froese [25], in his work, depicted that the initial and terminal priors of the relative biomass were represented by B start /k and B end /k of individual time series, and the assumed depletion levels were responsible for estimating their ranges [34]. For every tuna stock under assessment in this study, Table 3 presented the B start /k ranges (suggested). In order to evaluate the actual status of the stock and exploitation level of the species under study, the previous year's biomass (B/B MSY ) was used and further estimated using the CMSY and BSM methods [35]. The results are presented in Table 4. Table 3. Default range of biomass relative to k at the start (B start /k) and the end (B end /k) of the time series of stocks under assessment.

Stock Assessment by IOTC Using CMSY and BSM Models
CMSY is a data-limited Monte Carlo method of assessing fisheries reference points and relative stock size using catch data. The Schaefer model is run multiple times in CMSY to generate annual biomasses for r-K pairs drawn at random from previous distributions. The model decides whether r-K pairs are valid, for example, those that produce a biomass time series that does not (1) result in a stock collapse or (2) allows the stock to exceed carrying capacity. Additionally, r-K pairs that result in a final relative biomass estimate between the values specified in the inputs (the final depletion range) are accepted and used to calculate MSY (rK/4) and biomass over time. The most probable values are the geometric means of the resulting density distributions of r, K and MSY. These values are then used to establish managerial reference points. The CMSY method performs an enhanced Bayesian state-space implementation of the Schaefer surplus production model (BSM) when abundance indices are known. The BSM has the advantage of focusing on informative priors and allowing short (minimum 3 years) and fragmented abundance data [36]. The key difference between our evaluation and the IOTC's stock assessment for neritic tuna is that IOTC assessed tuna stocks across the Indian Ocean, whereas, in our analysis, we distinguished comparative stock assessments in the Eastern (Andaman Islands, Lakshadweep Islands, Nicobar Islands, Pamban Islands and St. Mary's Islands) and Western Indian Oceans (Vypin Islands, Vallarpadam and Willingdon Islands).

Results
The B/B MSY values for all 10 tuna fishery stocks assessed in India from the Eastern Indian Ocean and Western Indian Ocean were less than one, which clearly indicates that all the stocks are overfished (Table 5; Table 6). Active r-K pairs analyzed by CMSY and BSM methods produce triangular-shaped clouds in the plot's space, where the r-K pair most likely (and its 95% confidence intervals) is located at the top of the triangle (Figure 2). The pattern revealed that all the stocks were depleting, although the trends were different for each stock in both the Eastern part of the Indian Ocean (EIO) and the Western part of the Indian Ocean (WIO) (Figures 3 and 4). The estimated biomass of five tuna species for the Eastern Indian Ocean has been depicted in Figure 3. As shown in Figure 3A-E, all the fish stocks were subjected to continuous overfishing. Yellowfin tuna and skipjack tuna showed a similar pattern in their biomass estimation ( Figure 3B,C). It is revealed that the actual relative biomass (B/B MSY ) observed for yellowfin tuna was 0.82 and for skipjack tuna was 0.63. Similarly, longtail tuna and kawakawa ( Figure 3D,E) had related patterns in biomass estimation with better relative biomass reaching up to 0.81-0.83, respectively (Table 5). Of all the stocks, frigate tuna ( Figure 3A) had the least B/B MSY of 0.55 and was overfished, with a substantial decrease in the stock since 1995.
Amongst the studied tuna stocks in the Western Indian Ocean (WIO) (Figure 4), the longtail tuna was highly exploited and overfished with relative biomass (B/B MSY ) of 0.44 ( Table 6). The stocks started depleting from the year 2000, and later, after 2012, the unsustainable exploitation of the species led to a sudden decline. In WIO, the biomass trajectories were also similar for yellowfin tuna and skipjack tuna, with B/B MSY of 0.82 and 0.84, respectively (Table 6).
2). The pattern revealed that all the stocks were depleting, although the trends were different for each stock in both the Eastern part of the Indian Ocean (EIO) and the Western part of the Indian Ocean (WIO) ( Figure 3; Figure 4).        The estimated biomass of five tuna species for the Eastern Indian Ocean has been depicted in Figure 3. As shown in Figure 3A-E, all the fish stocks were subjected to continuous overfishing. Yellowfin tuna and skipjack tuna showed a similar pattern in their biomass estimation ( Figure 3B,C). It is revealed that the actual relative biomass (B/BMSY) observed for yellowfin tuna was 0.82 and for skipjack tuna was 0.63. Similarly, longtail tuna and kawakawa ( Figure 3D,E) had related patterns in biomass estimation with better relative biomass reaching up to 0.81-0.83, respectively (Table 5). Of all the stocks, frigate tuna ( Figure 3A) had the least B/BMSY of 0.55 and was overfished, with a substantial decrease in the stock since 1995.
Amongst the studied tuna stocks in the Western Indian Ocean (WIO) (Figure 4), the longtail tuna was highly exploited and overfished with relative biomass (B/BMSY) of 0.44 The findings of CMSY could be replicated with the BSM method for tuna stock in EIO (Table 7) and WIO (Table 8). For EIO, the fifth-ninety-fifth percentile ranges for MSY, B/B MSY , k and r were thinner for BSM than CMSY (Table 7). Only the estimate r of Auxis thazard (frigate tuna) and k of Katsuwonus pelamis (skipjack tuna), Thunnus tonggol (longtail tuna) and Euthynnus affinis (kawakawa) had broader ranges in Bayesian Model. Similarly, for the tuna stocks in WIO, fifth-ninety-fifth percentile ranges for MSY and B/B MSY were thinner for BSM than CMSY except for r estimate of Auxis thazard (frigate tuna) and k of Thunnus tonggol (longtail tuna) ( Table 8).

Discussion
Sustainable fisheries management is primarily based on stock assessment results, which provide useful information for successful policy formulation of fisheries resources. Many models of stock assessment require a lot of data [37] which ultimately makes their implementation limited to only valuable species and stocks [38,39]; however, other species are less considered [40,41]. Maximum Sustainable Yield is frequently used as a reference point for fish stocks assessment [41], and B/B MSY is used as a context for determining the status of fisheries [19]. Scientists and managers representing fisheries recognize principles of MSY, q, r and k for attaining sustainability in their respective fisheries [25,42,43]. In order to suggest biological reference points (BRPs) estimated by various statistical methods, scientists assess fishery stock status [44,45]. In this study, the biological reference points of five major tuna species have been estimated by two newly developed computerized methods: Monte Carlo Catch-MSY (CMSY) and the Bayesian Schaefer Model (BSM) for estimating the level of exploitation of five major tuna fishery stocks in the Eastern and Western Indian Ocean. An important benefit of using these sophisticated approaches is that they operate in data-scarce circumstances and develop outcomes that can assist national and regional policy and can support key decisions in management. Moreover, overfishing in marine water can disturb the equilibrium of fish species and other marine organisms, and overfishing of specific species may lead to extinction [46,47]. Thus, management strategies for controlling the overfishing status of tuna fisheries should therefore be stringent in order to ensure the sustainability of the species in the Indian Ocean.
The k and r were the necessary prior information for the CMSY and BSM [20] that can be obtained from resilience data obtained from the Fishbase [36]. CMSY shows some inefficiencies when used for very-low-resilience or less-captured fisheries, as per [48] simulation studies. The resilience value for all the five tuna fishery stocks in both the Eastern and Western Indian Ocean was in the range of low to medium. To determine the status of fishery stocks, catch, efforts and CPUE are key indicators. Catchability coefficient (q) is also important to study the fish population dynamics and usually denotes the ratio between fishing mortality and fishing efforts. Hoggarth [49] suggested that q positively affects the spatial distribution of fish stocks and fishing zones.
The ratio of B/B MSY less than 1.0 produced by both CMSY and BSM models determines the over-fished condition and unhealthy biomass of the tuna fishery in both the Eastern and Western Indian oceans. In the EIO, the least relative biomass (B/B MSY ) was observed in frigate tuna (0.55) ( Table 5), and in WIO, the least (B/B MSY ) was for the longtail tuna (0.44) (Table 6), with both having an F/F MSY value > 1, which depicts that both fisheries are overfished but are still sustainable. The r values for the longtail tuna in the entire Indian Ocean assessed by IOTC [50] ranged from 0.96 to 1.48 (mean of 1. 19), whereas in our study, the r values for longtail tuna in the Eastern Indian Ocean ranged from 0.216 to 0.428 (mean of 0.304) and for Western Indian Ocean, ranged from 0.224 to 0.473 (mean of 0.326). Similarly, another study of IOTC for assessment of kawakawa (Euthynnus affinis) in the Indian Ocean using data-limited methods revealed the status of fish as overfished with BMSY of 0.97 but is not subject to overfishing (FMSY of 1.16) [51]. In addition, the catches in the last 5 years have exceeded the MSY estimate.
The strong assumptions imposed on population dynamics and stock productivity, particularly the intrinsic growth rate and carrying capacity parameters, are largely responsible for the consistent estimates among CMSY simulations. The terminal depletion level assumption is subjective, yet it has a significant impact on stock status assessments. Because the model relied on standardized CPUE indices to offer information on abundance trends, it is less reliant on subjective assumptions (particularly those relating to depletion levels, which are critical in stock reduction assessment). The BSM also includes a more robust statistical prediction framework, which allows for an even more precise estimation of key parameters and management quantities. However, the availability of the standardized CPUE as a potential abundance index and its inclusion in the assessment constitutes a significant step forward in the development of more reliable methodologies to assess five major tuna species in both EIO and WIO in the context of data scarcity. Future studies should look towards developing more accurate population models, such as age-structured models that take into account more biological and fisheries data rather than catch series.

Conclusions
This contribution assessed five major exploited tuna stocks of India in the Eastern and Western Indian oceans using CMSY and BSM computer-intensive methods. The results revealed that all the tuna stocks in both the regions were overfished, with specifically the longtail tuna (Thunnus tonggol) fishery strongly overfished in the Western Indian Ocean Region.
The methods of stock assessment are based primarily on catch and catch per unit effort data (such as the CMSY and BSM methods) and are highly precise when stocks are lightly fished, while for over-exploited fisheries, catch variations typically represent differences in the underlying biomass. If the applications of CMSY and BSM methods are to be increased and made more efficient, it is compulsory to have time series catch data of many years, and it should be reliable