Int. J. Aquat. Biol. (2021) 9(2): 97-104 ISSN: 2322-5270; P-ISSN: 2383-0956 Journal homepage: www.ij-aquaticbiology.com © 2021 Iranian Society of Ichthyology Original Article Optimized exploitation of Pharaoh Cuttlefish (Sepia pharaonis Ehrenberg, 1831) stocks in the Iranian part of Persian Gulf and Oman Sea Seyed Ahmadreza Hashemi* 1, Mastooreh Doustdar2 1Offshore Fisheries Research Center, Iranian Fisheries Science Research Institute, Agricultural Research, Education and Extension Organization, Chabahar, Iran. 2Iranian Fisheries Science Research Institute, Agricultural Research, Education and Extension Organization, Tehran, Iran. s Article history: Received 29 February 2020 Accepted 22 March 2021 Available online 2 5 April 2021 Keywords: Maximum sustainable yield Monte Carlo simulation method ARIMA model Abstract: The purpose of the present study was to investigate the trends in Pharaoh Cuttlefish (Sepia pharaonis) capture fisheries and determine the suitable range for optimized exploitation of S. pharaonis resources in the Iranian part of Persian Gulf and Oman Sea using catch data. The data on Pharaoh Cuttlefish capture fisheries in Iranian southern waters for the twenty-three years was collected and the suitable range for optimized exploitation of S. pharaonis was estimated using the R Software. The average values (95% confidence interval) using the Monte Carlo simulation method for intrinsic population growth rate (r), maximum sustainable yield (MSY), the biomass of maximum sustainable yield (Bmsy) and maximum fishing mortality rate of maximum sustainable yield (Fmsy) were 0.92 (0.73-1.17) per year, 5100 (4200-6200) tons, 1100 (8670-13900) tons, 0.46 (0.36-0.58) per year, respectively. The results showed that the annual catch of S. pharaonis exceeded the maximum sustainable yields and measures should be taken to reduce the number of capture fisheries and fishing effort. With results of the prediction model was observed moving average analysis (MAPE=2.85, MAD=0.10, MSD=0.02) and ARIMA (0, 0, 1) (AIC=9.79, BIC=6.38), are better than other models for a period of five years for modeling annual this species landing. It seems that reducing fishing permits and fishing effort will put the S. pharaonis stock situation in a more favorable condition in the long term and will further benefit the exploiters and the fishing community. Introduction Global total capture fisheries production was about 91 million tons in 2016, of which 87% were in seawater (79.3 million tons) and 13% from inland waters (11.6 million tons) (FAO, 2018). In recent years, there are remarkable signs of over exploitation of major fish stocks and other aquatic resources. The proportion of stocks fished at Biologically Sustainable Levels (BSLs) to Biologically Unsustainable Levels (BULs) in 1974 was about 90%, whereas this proportion in 2016 was about 67%. Thus, BULs has increased and requires immediate management measures (FAO, 2018). The entire amount of fisheries production in Iran was near 800,000 tons, of which 720,000 tons (more than 90 percent) were fisheries production within the seawaters of southern Iran (Fisheries Statistical Yearbook, Iran area, 2020). *Correspondence: Seyed Ahmadreza Hashemi DOI: https://doi.org/10.22034/ijab.v9i2.736 E-mail: seyedahmad91@gmail.com The quantity and quality of knowledge from many of the world’s fisheries are inadequate to enable traditional methods of assessment to be applied. The management of "data-rich" stock approaches is usually focused on complex models of stock evaluation and involves variety of knowledge sources. Stock management is currently faced with numerous fish stocks that have little knowledge that do not help these with data-rich approaches (Dick and MacCall, 2011; Ghaitaranpour et al., 2019). Today, Length-based models, like Length Based Spawning Potential Ratio (LBSPR), Length-Based Integrated Mixed Effects (LIME), and Length-Based Bayesian (LBB), and also as catch-based methods, like Catch-Maximum Sustainable Yield (Catch-MSY), Depletion Based Stock Reduction Analysis (DBSRA), Simple Stock Synthesis (SSS), and Catch-MSY (CMSY) in many fishery scenarios and 97 Int. J. Aquat. Biol. (2021) 9(1): 97-104 in several countries are developed (known as “data- poor” or “data-limited” fisheries) (Wetzel and Punt, 2015). With their unique ecological conditions, the Persian Gulf and Oman Sea host a good range of aquatic species that provide settlers with livelihoods, jobs and various economic activities (Rajaei et al., 2014; Taghavimotlagh and Shojaei, 2017; Eagderi et al., 2019). Iran has quite 120,000 fishermen, mainly engaged in fishing, and fishing has played a serious role in generating jobs in coastal areas and in post-harvest economic activities (Taghavimotlagh and Shojaei, 2017). Mollusks are the largest and most diverse group of invertebrates after arthropods with more than 80,000 extant species (Bruska et al., 2016). Among the economically important mollusks, Pharaoh Cuttlefish (Sepia pharaonis) is a neritic demersal species inhabiting coastal waters to the depth of 130 meters. It is more likely to be found at the depth of 10 to 40 m (Jereb and Roper, 2005). This species is a predominant species in the Persian Gulf and Oman Sea (Valinasab, 1993, 1997; Roper et al., 1984). In addition, S. pharaonis is found in the Indo-Pacific region, the Red Sea, the Arabian Sea to the South China Sea, the East China Sea and the northern Australian waters (Jereb and Roper, 2005). Estimation of optimal fishing based on the time series of various aquatic species fishery has been carried out in some fish species (Martell and Froese, 2013; Froese et al., 2016; Zhou et al., 2017). Hence, the present study aimed to investigated the trends in Pharaoh Cuttlefish capture fisheries in the Iranian part of Persian Gulf and Oman Sea with the aim of determining the suitable range for its optimized exploitation and proper and principled exploitation management. Materials and Methods The data on Pharaoh Cuttlefish capture fisheries landings (based on metric ton) in the Iranian part of Persian Gulf and Oman Sea (Fig. 1) for the twenty- three years was collected from the Iran Fisheries Organization (from 1997 to 2019). Catch-MSY (CMSY): The Catch-MSY model has the same characteristics as the Graham-Shaefer surplus production model. These models rely on only a catch time series dataset and prior ranges of r and k and possible ranges of stock sizes in the first and final years of the time series. The CMSY is a method for estimating maximum sustainable yield (MSY) and related fisheries reference points (Bmsy, Fmsy) from catch data and information on resilience (Froese et al., 2017). This model requires a prior distribution on r and K as well as priors on the relative proportion of biomass at the beginning (Martell and Froese, 2013). The biomass in subsequent years was then generated Figure 1. Map of the study area and main Iranian fishing ports (Green color denoted the fishing ports). 98 Hashemi and Doustdar / Optimized exploitation of Pharaoh Cuttlefish stocks in the Persian Gulf and Oman Sea from a Schaefer model according to equation of By+1=By+rBy (1-By/k) es1-Ct es2, where By is Biomass in the year y+1, r = population instantaneous growth rate, K = carrying capacity, Cy = catch in the time series. In this method, the values of population instantaneous growth rate and carrying capacity are calculated with depletion formula (d) and storage saturation (S): d=1-S=1-By/Ky. The maximum steady-state mortality rate with the formula of Fmsy=r/2 and the maximum sustainable yield is calculated from MSY=rk/4 and Bmsy=K/2 (Zhou et al., 2017). A prior range was set for r based on the resilience of the stock as proposed by Martell and Froese (2013), where stocks with a high resiliency were allocated a r value from 0.6–1.5 (Zhou et al., 2017). The intrinsic population growth rate was calculated based on the inverse range factor and the formula of irf=3/r high–r low (Froese et al., 2016). The es1 and es2 remove the bias in the equation and simulation. The es1 is related to the process error, whereas the es2 is related to the observation error. The desired pattern of the model is achieved when the standard deviation of the process error is set to 0.2 and the observation error is set to 0.1 (Froese et al., 2016). Forecasting Methods: Univariate time series models: Moving average (MA) method is the arithmetic mean of observations of the full data set and uses the arithmetic mean as the predictor of the future period (Karmaker et al., 2017). Exponential Smoothing (ES) method is a kind of weighted averaging method which estimates the future value based on previous forecast plus a percentage of the forecasted error. Trend Analysis (TD) may be a general model to multiple statistic data having trend pattern and provides idea to traders about what is going to happen within the future supported historical data (Karmaker et al., 2017). Winters Method (WM) is employed to smooth data employing A level component, a trend component, and a seasonal component at each period and provides short to medium range forecasting. Decomposition method (DM) is employed to separate the statistic into linear trend and seasonal components and error (Karmaker et al., 2017). The autoregressive integrated moving average (ARIMA) models demonstrated a good performance in terms of explained variability and predicting power. The autocorrelation (ACF) and partial autocorrelation functions (PACF) were estimated, which led to the identification and construction of ARIMA models (Tsitsika et al., 2007). The ARIMA were implemented on the data and the best model was chosen using the test of Akaike coefficient data autocorrelation functions (Lawer, 2016). The general form of the ARIMA models is referred to as ARIMA (p, d, q), where p is the order of the autoregressive term (AR term), d is the degree of differencing involved to achieve stationarity, and q is the order of the moving average term (MA term). The best forecasting model: The values of Mean Absolute Percentage Error (MAPE), Mean absolute deviation (MAD) and mean square deviation (MSD) are used to identify the best model (Karmaker et al., 2017). 𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 = � �𝑒𝑒𝑒𝑒𝐷𝐷𝑒𝑒� 𝑛𝑛 ∗ 100 MAD=∑ |𝐷𝐷𝐷𝐷−𝐹𝐹𝐷𝐷| 𝑛𝑛 MSE=∑(𝐷𝐷𝑒𝑒 − 𝐹𝐹𝑒𝑒)2/ 𝑛𝑛 − 1 Where Dt is actual demand for time period t, Ft is forecast demand for time period t, n is specified number of time periods, and et is forecast error = (Dt − Ft). Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC)) were estimated as follows: AIC = -2 ln (maximum likelihood) + 2m BIC = -2 ln (maximum likelihood) + m ln(n) Where m is that the number of the estimated parameters and n is that the number of the observations. Statistical analyses were performed with R studio (1.1.45), SPSS (22) and Minitab (16) software package and a significance level of 0.05 was adopted. Results The average catch (Yi) of the studied period was 3536 tones with 95% confidence intervals of 4284-2919 tones and the average catch was significantly 99 Int. J. Aquat. Biol. (2021) 9(1): 97-104 decreased for the twenty-three years (R=-0.14, P<0.05) (Fig. 2). The amounts of trend catch of this species showed a decreasing trend (Fig. 2). The CMSY model: Based on the annual catch catches and initial Pharaoh Cuttlefish information (initial growth 0.6-1.5 per year), the software used the initial value to start of modeling based on the CMSY and the Monte Carlo simulation method. The initial relative biomass as 0.5-0.9 and the final relative biomass as 0.2-0.6 were considered. The output values Table 1. Average values (95% confidence interval) of CMSY model parameters for Pharaoh Cuttlefish in the Persian Gulf and Oman Sea (Iranian coastal waters). Indices / models CMSY Average (Maximum-minimum) Biomass (1000 tonnes) 11.9 ( 7.1-13.2) MSY (1000 tonnes) 5.1 (4.2-6.2) Bmsy (1000 tonnes) 11 (8.67-13.9) Fmsy 0.46 (0.36-0.58) F 0.46 (0.41-0.77) B/Bmsy 1.08 (0.64-1.2) F/Fmsy 1 (0.93-1.32) K (1000 tonnes) 22 (17.3-27.9) r 0.92 (0.73-1.17) Bt/K 0.54 Bmsy/K 0.50 Figure 2. The changes trend in the annual (yr) and catch (total) (A) and trend catch of this period with 95% confidence intervals (B) for Sepia pharaonis in the Persian Gulf and Oman Sea (Iranian coastal waters). 100 Hashemi and Doustdar / Optimized exploitation of Pharaoh Cuttlefish stocks in the Persian Gulf and Oman Sea of the Monte Carlo simulation method with 30,000 times the modeling repeat result the average values (95% confidence interval) of state-space surplus production model for intrinsic population growth rate (r), carrying capacity (k), maximum sustainable yield (MSY), biomass of maximum sustainable yield (Bmsy), current or existing biomass (B), maximum fishing mortality rate of maximum sustainable yield (Fmsy) and current or existing fishing mortality (F) that are shown in Table 1. The ratios of the current biomass to the biomass of maximum sustainable yield (B/Bmsy) and the ratio of the current fishing mortality to the maximum fishing mortality rate of maximum sustainable yield (F/Fmsy) were 1.08 (0.64-1.2) and 1 (0.91-1.2), respectively (Fig. 3). According to the kobe plot, which represents the graphic of the values of B/Bmsy and F/Fmsy, it can be concluded that the increase in the fishing Table 2. Forecasting and error calculations of different methods (Mean Absolute Percentage Error (MAPE), Mean absolute deviation (MAD), mean square deviation (MSD) for Pharaoh Cuttlefish. Method MAPE MAD MSD Multiplicative decomposition (Yt=3.431+0.0375*t) 4.43 0.15 0.035 Additive decomposition (Yt=3.429+0.0377*t) 4.43 0.15 0.035 Moving average 2.85 0.1 0.02 Single exponential smoothing 3.88 0.13 0.025 Double exponential smoothing 3.87 0.13 0.027 Trend analysis (Linear) (Yt=3.56+0.0017*t) 4.43 0.15 0.035 Trend analysis (Exponential) (Yt=3.56*(0.99*t) 4.43 0.15 0.035 Trend analysis (Quadratic) (Yt=3.90-0.083*t+0.0033*t2) 3.03 0.10 0.016 Winters multiplicative 4.13 0.14 0.036 Winters additive 4.13 0.14 0.036 Figure 3. Changes in biomass of maximum sustainable yield, intrinsic population growth rate and carrying capacity of Pharaoh Cuttlefish in the Persian Gulf and Oman Sea (Iran) based on 30,000 times repeat of modeling using the Monte Carlo simulation method for CMSY Model. 101 Int. J. Aquat. Biol. (2021) 9(1): 97-104 mortality and the reduction of the current biomass has begun and is still ongoing. The mortality fishing (F) to mortality fishing of the maximum sustainable yield (Fmsy) ratio (F/Fmsy), the biomass (B) to biomass of the maximum sustainable yield (Bmsy) ratio (B/Bmsy) in the CMSY model were 1 (0.93-1.32) and 1.08 (0.64-1.20) (in 2019 year), respectively. Forecasting methods: Performance of sixteen methods was evaluated based on forecasting accuracy. The ten different forecasting techniques, including moving average and using of six methods identified orders of ARIMA (p, d, q) based on the AIC and BIC. With results of the prediction model was observed moving average analysis (MAPE=2.85, MAD=0.10, MSD=0.02) and ARIMA (0, 0, 1) (AIC=9.79, BIC=6.38), are better than other models for a period of five years for modeling annual this species landing (Fig. 4). Various models were then fitted and compared as presented in Tables 2 and 3, using identified orders of ARIMA (p, d, q) based on the AIC and BIC. However, ARIMA (0, 0, 1) with drift was suitable for modeling annual S. pharaonis landings Table 3. Coefficients and summary statistics of multivariate ARIMA modeling for Pharaoh Cuttlefish in the Persian Gulf and Oman Sea (Iran) in the Persian Gulf and Oman Sea (Iran). Method B (coefficient) Constant Standard error of B Log likelihood AIC BIC AR1 MA1 AR1 MA1 ARIMA(1,0,0) 0.58 - 3.58 0.18 - 9.52 13.03 9.62 ARIMA(0,1,0) - - - - - 7.07 12.14 11.05 ARIMA(0,0,1) - 0.37 3.55 - 0.16 7.78 9.79 6.38 ARIMA(1,0,1) 0.77 -0.29 3.59 0.19 0.27 9.95 11.90 7.35 ARIMA(0,1,1) -0.38 - - 0.19 - 8.59 13.05 10.86 ARIMA(1,1,1) -0.16 -0.24 - 0.60 0.59 8.56 11.13 7.85 Figure 4. The moving average trend fitted models (A), the ARIMA (0,0,1) trend fitted (B) models (forecast plot is blue line) for of Sepia pharaonis in the Persian Gulf and Oman Sea (Iranian coastal waters) 102 Hashemi and Doustdar / Optimized exploitation of Pharaoh Cuttlefish stocks in the Persian Gulf and Oman Sea (Fig. 4) based on the selection criteria (AIC=9.79, BIC=6.38). Discussions In recent years, the catching of Pharaoh Cuttlefish in the waters of southern Iran has been decreased. Over the past two decades, its catch rate in the waters of southern Iran has dropped from about 9,000 tons in 1997 to about 5,000 tons in 2019, indicating a sharp decline (Fisheries Statistical Yearbook, Iran, 2020). In this regard, the Sistan and Baluchestan Province had the highest percentage of Pharaoh Cuttlefish capture fisheries (63%), followed by Bushehr (19%) and Hormozgan (10%) provinces. Khuzestan Province had the smallest percentage of Pharaoh Cuttlefish capture fisheries (8%) and the Hormozgan Province showed the highest reduction. The intrinsic population growth rate (r) value of this species is high flexibility (0.6-1.5) indicating its high potential to stand up to fishing pressure and recovery i.e. it demonstrates the ability to withstand fishing pressure and the recovery of declining fish stocks (Martell and Froese, 2013; Froese et al., 2016; Zhou et al., 2017). The earlier work showed S. pharaonis as the most flexible species (Froese and Pauly, 2015). A strong correlation exists between the intrinsic population growth rate (r) and other parameters of life history, especially natural mortality (M). The best model is r=1.73 M for teleosts and r=0.76 M for elasmobranches (Zhou et al., 2017). Froese and Pauly (2015) confirmed that the population intrinsic growth rate (r) is approximate twice the maximum fishing mortality rate of maximum sustainable yield (Fmsy), 2 times the natural mortality (M), 3 times the growth rate coefficient of the von Bertalanffy curve (K), 3 divided by the generation time (tgen) and 9 divided through the most age (tmax) (r ≈ 2FMSY ≈ 2 M ≈ 3 K ≈ 3 / tgen ≈ 9 / tmax). In the present study, the decreasing trend in the amount of current biomass to biomass of maximum sustainable yield (B/Bmsy) showed that Pharaoh Cuttlefish had been in full exploitation condition for catching fisheries in the Iranian part of Persian Gulf and Oman Sea. In addition, the level of current fishing mortality at the fishing mortality rate to be the maximum sustainable yield (F/Fmsy), is shown with an increasing trend and toward overfishing pattern (Arrizabalaga et al., 2012). The fishery status is usually evaluated based on B/Bmsy and is divided into three general categories: B/Bmsy≥1/5 means stocks under exploited, 0.5