Evaluation of Biological Reference Points for Conservation and Management of the Bigeye Thresher Shark, Alopias superciliosus, in the Northwest Pacific

Abstract

Full stock assessment of sharks is usually hindered by a lack of long time-series catch and effort data. In these circumstances, demographic and per-recruit analyses may provide alternate approaches to describe population status because these methods can be applied to estimate biological reference points (BRPs) for shark stocks. However, the appropriate level of BRPs for sharks is difficult to determine, given the expected low reproductive rates. To determine which BRPs are most appropriate for the CITES-listed species—bigeye thresher shark, Alopias superciliosus, a stochastic demographic model with Monte Carlo simulations and per-recruit models were used to estimate BRPs in this study. The results indicated that conventional fishing mortality-based BRPs (F_BRPs) derived from per-recruit models may result in a clear population decline. Our analyses also demonstrated that the bigeye thresher population in the Northwest Pacific will stabilize only if demographic-based F_BRP is implemented. The F_BRP estimated based on the stochastic demographic model was 0.079–0.139 y⁻¹, which was equivalent to SPR = 50–70%. The findings strongly suggested that more conservative threshold F_BRPs should be implemented to ensure sustainable utilization of the bigeye thresher stock. The present study provides new and strategically important information on the population dynamics of the bigeye thresher in the Northwest Pacific, which can be used to help fishery managers to adopt more efficient management measures for this stock. It is also suggested that this approach can be applied to other shark species with limited catch and effort data.

1. Introduction

Most pelagic sharks exhibit prolonged life span, late maturity, and low fecundity [1,2,3,4], and are vulnerable to perturbations imposed by anthropogenic factors such as fisheries [5,6]. Sharks are commonly exploited worldwide for their meat, skins, fins, livers, cartilage, jaws, and teeth [7]. Heavy exploitation and largely unregulated trade in shark species, however, are considered to have resulted in the decline of global shark stocks [8]. Accordingly, shark conservation and management have attracted great attention in recent years. Oceanic sharks, although heavily exploited by various fisheries, remain among the least studied and managed fish due to the limited information in their biology and fishery [9]. Although the knowledge about the biology, stock status and population dynamics of some common bycatch shark species have been advanced, evidence of some least productive species is still insufficient and is urgently needed [10].

A common problem for shark stock assessment is that the data required in conventional stock assessment models are rarely available due to the fact of low commercial value, and a lack of regular records in fisheries statistics. In this situation, demographic models can provide valuable insights into the development of management advices for fish species until sufficient fishery statistics data become available to support more complex conventional stock assessments [11,12]. The demographic methods might be relatively simple, and only require some biological information such as survival rate, age at maturity, litter size, longevity, and other reproductive parameters. Therefore, the status of a fish population can be simply described by the primary outputs obtained from demographic analysis (e.g., intrinsic rate of population growth) [11]. Demographic models have several advantages compared to conventional fishery stock assessment models. For example, conventional stock assessment models (such as surplus production models or age structured population models) require large quantities of data (e.g., catch, efforts, abundance indices) to be carried out. Unlike conventional modeling methods, demographic matrix models only require life history information. These models can be applied to estimate biological reference points (BRPs) [13] or used as stock status indicators [14,15,16]. In addition, all life history parameters and characteristics such as age-at-maturity, reproductive cycle, or sexual dimorphism can be taken into account in demographic analyses.

The bigeye thresher shark, Alopias superciliosus, an apex marine predator, is found in temperate and tropical oceans worldwide [17]. It is commonly caught by offshore fisheries and is one of the important by-catch shark species for tuna longline fisheries. This species has been identified as one of the least productive pelagic sharks, and there is increasing concern about its conservation status [9,18]. This species is susceptible to overexploitation due to its life history characteristics of slow growing, late maturity, and few offspring [2,6,9,19,20]. It has been listed as vulnerable (VU) on the International Union for Conservation of Nature (IUCN) Red List [21] and listed on the Appendix A II at the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES) CoP17 meeting due to a decline of abundance in certain waters [22]. Moreover, regional fisheries management organizations (RFMOs), such as the International Commission for the Conservation of Atlantic Tunas (ICCAT) and the Indian Ocean Tuna Commission (IOTC), have prohibited retention of this species on board for commercial use [23,24]. Although it is designated by the Western and Central Pacific Fisheries Commission (WCPFC) as a key shark species, apart from the Pacific-wide sustainability risk assessment [18], no full stock assessment of the bigeye thresher in the North Pacific has been conducted due to lacking of reliable catch, effort and abundance index information. Recently, Tsai et al. (2019) [20] assessed the stock status of bigeye thresher in an area subset of the western North Pacific using a separable virtual population analysis, per-recruit models and age-structured demographic analysis, and concluded that the current stock status is overexploitation. However, the appropriate Biological Reference Points (BRPs) under various harvest strategies were not evaluated in that study.

BRPs are widely used to assess the relative health of fish stocks and the relative intensity of fishing pressure, and they are often further subdivided into the limit (as levels not to be exceeded) and the target reference points (ultimate goal of management measures). BRPs are a key component of how RFMOs formulate fisheries management advice. The most commonly used BRPs are based on many components such as growth, mortality, or the maintenance of appropriate levels of recruitment to the stock. Some of them have been extensively used as target and limit in the management strategy evaluation and harvest control rules [25]. Conventional management reference points are commonly derived from per-recruit, production, or stock-recruitment models. For sharks, biomass-based BRPs are less common compared to fishing mortality-based BRPs because of the lack of long time-series catch and effort data and stock-recruitment relationship for conducting a full stock assessment [26,27]. Conversely, fishing mortality-based BRPs (F_BRPs) based on per-recruit or demographic models provide useful information on data-limited shark populations. Most adopted BRPs are ad hoc and are based on the life history and fishery processes of managed species. However, only a few of them have been rigorously examined [28]. As BRPs can be obtained from various methods, it is very important to evaluate the effectiveness and consistency of BRPs before applying to fishery management.

In the present study, we used the best available life-history parameters from previous studies and estimated fishing mortality while incorporating uncertainty from Bayesian inference to construct demographic age-structured stochastic population matrices. These matrix models were then used to evaluate various biological reference points by using stochastic simulations and to provide useful information regarding fishery management and conservation for the bigeye thresher shark in the Northwest Pacific. This approach can be applied to other shark species that have limited catch and effort data.

2. Materials and Methods

2.1. Source of Data

The pelagic fish, including tuna, billfish, and sharks caught by the Taiwanese small-scale longline vessels (<100 gross tonnage) in the western North Pacific were mainly (>90%) landed at the Nanfangao fish market, eastern Taiwan. These fish were weighed before auction. Thus, the species-specific individual whole weight (W) of these fish can be obtained from sales records. However, the sex of each individual shark was not available in sales records. A sub-sample of 4855 fish (3285 females, 1570 males) collected at the Nanfangao fish market between 2015 and 2019 was used to develop the weight-sex ratio (the proportion of females) relation, which was further used to derive the sex of individual fish. The sex ratios of sharks smaller than 40 kg and greater than 195 kg were set as 0.5 and 1.0, respectively, based on the observation of sub-sample. For fish between 45 and 195 kg, the sex ratio of weight (${\mathsf{\Phi }}_W$) was obtained from the whole weight-sex ratio relation: ${\mathsf{\Phi }}_W=\alpha \times W^{\beta }$, where $\alpha$ and $\beta$ are parameters to be estimated. The sex of each individual in each 5 kg class was randomly assigned sex based on the above equation. Sex-specific weight data were converted into pre-caudal length (PCL) based on the length-weight relationship [19]. Furthermore, catch-at-age was then estimated from the converted PCL by using the sex-specific growth equation [19]. All of the biological parameters used in this study are presented in

Table 1. Life history parameters of the bigeye thresher shark used in this study.

Parameter	Female
Sex ratio (${\mathsf{\Phi }}_W$) 1	${\mathsf{\Phi }}_W$
$weight<40\mathrm{kg}$	0.5
$\alpha$	0.218
$\beta$	0.262
$weight>195\mathrm{kg}$	1
Length–weight relationship 2
$a$	$6.87\times {10}^{-5}$
$b$	2.769
VBGE 3
$L_{\infty }$	224.6
$K$	0.092
$t_0$	−4.21
Maturity fraction 4
$r_m$	−0.747
$a_m$	12

¹ In this study, the sex ratios (the proportion of females) of sharks smaller than 40 kg and greater than 195 kg were set as 0.5 and 1.0, respectively, based on our observations. For fish between 40 and 195 kg, the sex ratio by weight (${\mathsf{\Phi }}_W$) was obtained from the equation: ${\mathsf{\Phi }}_W=\alpha \times W^{\beta }$, where $\alpha$ and $\beta$ are estimated parameters. (R² = 0.979; n = 4855, 5-kg classes, p < 0.0001). ² Liu et al. (1998). where $a$ and $b$ are estimated parameters for length–weight relationship. ³ Liu et al. (1998). where $L_{\infty }$ is the maximum attainable length, k is a Brody growth constant, $t_0$ is a hypothetical age at length of 0. ⁴ Data from Chen et al. (1997). where $r_m$ is the slope and $a_m$ is age at maturity estimated from logistic maturity model.

.

2.2. Mortality Estimation

As the lack of direct natural mortality estimate for bigeye threshers, four empirical formulae (as shown in Cases 1–4 below) were adopted for deriving constant or age-dependent natural mortality ($M_a$) and to account for the possible variations of natural mortality in the model simulations:

Case1 [29]:

$$\ln (M_a)=0.941-0.873\ln (a_{\max })$$

(1)

Case2 [30]:

$$M_a=-\ln 0.01/a_{\max }$$

(2)

Case3 [31]:

$$M_a=1.92y{r}^{-1}\times W_a{}^{-0.250}$$

(3)

Case4 [32]:

$$M_a=3.00y{r}^{-1}\times W_a{}^{-0.288}$$

(4)

where a is age, a_max is longevity, and W_α is the age-specific mean weight. To avoid the possible effect on natural mortality by growth parameters, the above four empirical equations were adopted in this study because the nature mortality was estimated based on body weight or longevity.

The model developed in this study was applied to females only (as no significant difference in growth between sexes was noted in Liu et al., 1998 [19]), and the dynamics of a simulated year class was projected forward using the Ricker’s (1975) [33] exponential survival equation: $N_{a+1}=N_a{e}^{-(M_a+F\times S_a)}$. Here, the gear selectivity ($S_a$) was assumed to exhibit a dome-shaped distribution following Tsai et al. (2011) [27].

The expected catch (${\widehat{C}}_a$) of a fish at age a can be estimated from the catch equation [33]:

$${\widehat{C}}_a=\frac{F\times S_a}{(M_a+F\times S_a)}{N}_a\left(1-{\mathrm{e}}^{-(M_a+F\times S_a)}\right),$$

(5)

where $N_a$ is the initial number of fish of age $a$; $S_a$ is the probability of the bigeye thresher being captured at each age. The estimated value of $F$ was considered to be the current fishing mortality (F_curr). All parameters were estimated by minimizing the sum of squared difference between the observed catch-at-age and model-predicted catch. Following the least-squares optimization approach, the objective function to be minimized is:

$${\displaystyle \sum _a(C_a-{\widehat{C}}_a}{)}^2,$$

(6)

where $C_a$ is the observed catch of fish at age $a$.

2.3. Model Fitting and Convergence

The parameters that minimize the negative log-likelihood function were estimated using the AD Model Builder [34]. In addition to a deterministic estimate of F_curr, the MCMC method based on the Metropolis-Hastings algorithm is used to estimate the Bayesian posterior distributions. The posterior distributions were obtained from samples generated by conducting 12,000,000 cycles of the Markov Chain Monte Carlo (MCMC) algorithm, selecting every 1000th parameter vector thereafter and ignoring the first 2000 cycles as the “burn-in” period. Convergence of the MCMC samples was evaluated by monitoring the density plots, trace plots, and autocorrelation diagnostics of model parameters. All subsequent diagnostic analysis was implemented in the CODA package [35] of the R program [36].

2.4. Biological Reference Points

The yield per recruit (YPR, [37]) and spawning per recruit models (SSB/R, [38]) were adopted in this study to estimate the fishing mortality-based BRPs (F_BRPs) for bigeye thresher sharks.

The Thompson-Bell model was used to calculate yield per recruit curves (Y/R) following the formula:

$$\raisebox{1ex}{$Y$}\!\left/ \!\raisebox{-1ex}{$R$}\right.={\displaystyle \sum _{a=a_c}^{a_{\max }}\left({\overline{W}}_{a,s}\frac{F\times S_a}{F\times S_a+M_a}\left(1-e^{-\left(F\times S_a+M_a\right)}\right)e^{-{\displaystyle \sum _{i=a_c}^{a-1}\left(F\times S_i+M_a\right)}}\right)},$$

(7)

where $a_c$ is the age of a fish at first capture (set as age 1) and $a_{\max }$ is the longevity. The subscript “i” denotes the accumulated survivorship for each age of the cohort.

The spawning potential ratio (SPR) can be calculated as [38]:

$$SPR=\frac{SSB/R}{SSB/R|{}_{F=0}}\times 100\%$$

(8)

where SSB/R is the spawning stock biomass per recruit.

Similarly, assuming a constant year class, SSB/R can be obtained by following equation [38]:

$$\raisebox{1ex}{$SSB$}\!\left/ \!\raisebox{-1ex}{$R$}\right.={\displaystyle \sum _{a=1}^{a_{\max }}\left(m_a·{\overline{W}}_a·e^{-{\displaystyle \sum _{i=a_c}^{a-1}(F\times S_i+M_i})}\right)},$$

(9)

where m_a is the proportion of mature females at age a (further details can be found in [20]).

In this study, a number of biological reference points were estimated including (1) the management targets F_BRP (F_0.1) and threshold F_BRP (F_max) obtained from YPR model; (2) reference points based on SPR model: the threshold (F_SPR30%), and the target (F_SPR35%) reference points that corresponded to SPRs of 30% and 35%, respectively. The above F_BRPs were compared with current fishing mortality rate to evaluate the status of the bigeye thresher population.

2.5. Demographic Model Development

Thresher sharks such as pelagic threshers or bigeye threshers generally exhibit year-round parturition life history characteristics [19,39,40]. Therefore, the birth-flow approximation is likely more appropriate for population analysis of thresher sharks than conventional matrix population model. To account for the continuous reproduction for bigeye threshers, the following age-based matrix model (A) is commonly used for sharks [20,41,42,43]:

$$\mathrm{A}=\left[\begin{array}{ccccc}{f}_0& f_1& f_2& \cdots & f_{a=a_{max}}\\ P_0& 0& 0& 0& 0\\ 0& P_1& 0& 0& 0\\ 0& 0& \cdots & 0& 0\\ 0& 0& 0& P_{a_{max}-1}& 0\end{array}\right]$$

where P_a is the annual natural survivorship for age α, f_α represents the age-specific fecundity. In this case, birth is assumed to have a continuously and uniform distribution throughout the year [39]. More details on parameters estimation for demographic model can be found in Tsai et al. (2019) [20].

Demographic matrix model (A) was then used to estimate finite rate of population increase (λ), intrinsic rate of population growth (r) and the critical fishing mortality (F_crit) at which population is in equilibrium (r = 0 or λ = 1). The following life history parameters were assumed for a deterministic base run:

(1)

Age at maturity = 12 years

(2)

Fecundity = 2 pups

(3)

Sex ratio = 0.5 for embryos

(4)

Selectivity (assumed constant dome-shaped distribution).

(5)

A knife-edge maturity was assumed in this model and age-at-first-reproduction calculated as the mean age at maturity + the gestation period (set as 1 year in this study).

To reflect uncertainties on estimation of BPRs, the sensitivity runs for the three possible longevities (a_max = 35, 30, 25) were examined for the four cases of mortality estimates. In total, 12 runs were conducted:

• Case 1: natural mortality was estimated from Hoeing (1983) [29].
• Case 2: natural mortality was estimated from Campana et al. (2001) [30].
• Case 3: natural mortality was estimated from Peterson and Wroblewski (1984) [31].
• Case 4: natural mortality was estimated from Lorenzen (1996) [32].

The estimations of BRPs described above are deterministic (set as the reference cases).

2.6. Design of the Simulation Study

2.6.1. Biological Reference Points

In addition to deterministic estimates of BPRs, a stochastic method was also applied to include the possible uncertainty in the estimation process. However, for simplicity, the simulations were only conducted for longevity of 35 years, which is the most likely value of female bigeye thresher shark based on von Bertalanffy growth equation of Liu et al. (1998) [19]. For both per-recruit and demographic models, 10,000 replicates of BRPs were estimated by using the posteriors of F_curr and selectivity derived from MCMC. The central tendency and variation for the distributions were quantified by the median and the interquartile range.

2.6.2. Estimates of Population Growth Rates

To deal with the uncertainty regarding life history parameters, we created plausible parameter ranges and propagated these uncertainties through the model to cover the plausible range of the rate of population increase. Three main possible uncertainties in the demographic estimates included the age at maturity, natural mortality and fishing mortality rates. The triangular or lognormal distributions can be applied to represent the uncertainty of life-history parameters that precedes demographic modelling [15,44]. The median age-at-maturity was estimated to be 12 years old for female bigeye threshers, with maturation occurring between 10 and 13 years of age [39]. A triangular distribution was assumed to account for the uncertainty of age at maturity. Age-specific natural mortality ($M_a$) was randomly selected from the estimates derived from the four methods mentioned above. All estimates were given equal weight. A lognormal error structure for $F_a$ can ensure that the generating survival estimates range between 0 and 1. The mean and standard deviation of the age-specific fishing mortality ($F_a$) obtained from the MCMC were used to define a lognormal distribution.

The uncertainty related to age-at-maturity, natural mortality and the fishing mortality rate were then incorporated into the simulations. The BRPs obtained from per-recruit analyses were also set as input values of the demographic model to investigate the possible differences between the per-recruit and demographic models. In total, seven harvest strategies were conducted to assess the population status and to explore the implications of potential management strategies. These harvest strategies were:

• Scenario 1: fishing mortality for all ages set to 0.
• Scenario 2: fishing mortality equal to its current level by age.
• Scenario 3: fishing mortality set to the F_0.1 level.
• Scenario 4: fishing mortality set to the F_max level.
• Scenario 5: fishing mortality set to the F_SPR35% level.
• Scenario 6: fishing mortality set to the F_SPR30% level.
• Scenario 7: fishing mortality set to the F_crit level.

To compute the 95% confidence intervals for both population increase rate (λ) and intrinsic rate of population growth (r): for each scenario, 10,000 replicates of population growth rate were estimated by incorporating parameter uncertainty in Monte Carlo simulation. All demographic and simulation analyses were conducted using CSIRO program-PopTools [45].

3. Results

3.1. Deterministic Estimates

3.1.1. Sex-Specific Catch and Weight Compositions

The relationship between the sex ratio (Φ_W) and weight over the range of 40 to 195 kg was calculated by the following equation: Φ_W= 0.218 × W^0.262 (r² = 0.979; n = 4855, 5-kg classes, p < 0.0001) (Table 1). Based on the weight-specific sex ratio, a total of 20,804 bigeye thresher sharks landed at Nanfangao fish market between January 2015 and December 2019 and were divided into 13,778 females and 7026 males. The major group of the catch fell in the range of 60–80 kg (Figure 1a) corresponding to ages 6–9 years for both sexes (Figure 1b).

3.1.2. Mortality and Selectivity

The range of age-specific natural mortality (M_a), produced by the four indirect methods, was 0.088–0.199 y⁻¹. The lowest estimates of M = 0.107 y⁻¹ by average (calculated as the mean of age-specific M) was obtained using the empirical equation of Peterson and Wroblewski (1984) [31], which relies on the weight at age. The highest estimates of M = 0.184 y⁻¹ were obtained using the method of Campana et al. (2001) [30], based on age at longevity of 25 years (

Table 2. Estimated natural mortality (M) for the bigeye thresher shark using four empirical method based on longevity of 35 years.

Age	Weight	Case 1	Case 2	Case 3	Case 4
1	12.336	0.115	0.132	0.182	0.199
2	18.698	0.115	0.132	0.164	0.177
3	26.049	0.115	0.132	0.151	0.160
4	34.166	0.115	0.132	0.141	0.148
5	42.839	0.115	0.132	0.133	0.139
6	51.872	0.115	0.132	0.127	0.132
7	61.099	0.115	0.132	0.122	0.126
8	70.373	0.115	0.132	0.118	0.121
9	79.575	0.115	0.132	0.114	0.116
10	88.608	0.115	0.132	0.111	0.113
11	97.396	0.115	0.132	0.109	0.110
12	105.880	0.115	0.132	0.106	0.107
13	114.017	0.115	0.132	0.104	0.105
14	121.778	0.115	0.132	0.103	0.103
15	129.144	0.115	0.132	0.101	0.101
16	136.106	0.115	0.132	0.100	0.100
17	142.660	0.115	0.132	0.099	0.098
18	148.812	0.115	0.132	0.098	0.097
19	154.567	0.115	0.132	0.097	0.096
20	159.940	0.115	0.132	0.096	0.095
21	164.943	0.115	0.132	0.095	0.094
22	169.592	0.115	0.132	0.095	0.094
23	173.905	0.115	0.132	0.094	0.093
24	177.899	0.115	0.132	0.093	0.092
25	181.593	0.115	0.132	0.093	0.092
26	185.005	0.115	0.132	0.093	0.091
27	188.152	0.115	0.132	0.092	0.091
28	191.053	0.115	0.132	0.092	0.090
29	193.723	0.115	0.132	0.092	0.090
30	196.179	0.115	0.132	0.091	0.090
31	198.437	0.115	0.132	0.091	0.089
32	200.510	0.115	0.132	0.091	0.089
33	202.413	0.115	0.132	0.091	0.089
34	204.159	0.115	0.132	0.090	0.089
35	205.760	0.115	0.132	0.090	0.088
Mean	132.264	0.115	0.132	0.107	0.109

Values for M2 and M3 were 0.132 y⁻¹ and 0.154 y⁻¹ based on longevity of 30 years and 0.154 y⁻¹ and 0.184 y⁻¹ based on longevity of 25 years.

). Generally, the exponential survival equation [33] fit the observed catch data well for each case (Figure 1c). The estimated mean (μ) and standard deviation (σ) of the dome-shaped component based on different estimates of M were very close (

Table 3. Estimates of fishing mortality (F), selectivity, mean ($\mu$) and standard deviation ($\sigma$) of the dome-shaped component based on four empirical estimators of natural mortality (Cases 1–4) for female bigeye thresher sharks in the Northwest Pacific Ocean.

M	Case 1	Case 2	Case 3	Case 4
F	0.504	0.497	0.507	0.506
$\mu$	8.796	8.795	8.798	8.797
$\sigma$	2.053	2.048	2.049	2.048
Age	Selectivity
1	0.000	0.000	0.000	0.000
2	0.001	0.001	0.001	0.001
3	0.004	0.004	0.004	0.004
4	0.019	0.018	0.018	0.018
5	0.066	0.065	0.065	0.065
6	0.182	0.181	0.180	0.180
7	0.398	0.396	0.396	0.395
8	0.686	0.685	0.684	0.684
9	0.932	0.932	0.932	0.932
10	1.000	1.000	1.000	1.000
11	0.846	0.845	0.846	0.846
12	0.565	0.563	0.564	0.564
13	0.297	0.295	0.297	0.296
14	0.124	0.122	0.123	0.122
15	0.041	0.040	0.040	0.040
16	0.011	0.010	0.010	0.010
17	0.002	0.002	0.002	0.002
18	0.000	0.000	0.000	0.000
19	0.000	0.000	0.000	0.000
20	0.000	0.000	0.000	0.000
21	0.000	0.000	0.000	0.000
22	0.000	0.000	0.000	0.000
23	0.000	0.000	0.000	0.000
24	0.000	0.000	0.000	0.000
25	0.000	0.000	0.000	0.000
26	0.000	0.000	0.000	0.000
27	0.000	0.000	0.000	0.000
28	0.000	0.000	0.000	0.000
29	0.000	0.000	0.000	0.000
30	0.000	0.000	0.000	0.000
31	0.000	0.000	0.000	0.000
32	0.000	0.000	0.000	0.000
33	0.000	0.000	0.000	0.000
34	0.000	0.000	0.000	0.000
35	0.000	0.000	0.000	0.000

). Overall, the dome-shaped selectivity for female bigeye threshers revealed that most of the catch was immature and peaked at ages 8–10 years (Table 3, Figure 1d).

3.1.3. Biological Reference Points

The computation of the BRPs were conducted for different M indirect methods at three values of longevity (35, 30, and 25 years). The results of BRP analysis are summarized in

Table 4. Estimates of current fishing mortality (F_curr) and biological reference points derived from YPR and SPR models based on four empirical estimators of natural mortality for female bigeye thresher sharks in the Northwest Pacific Ocean.

Natural Mortality	Longevity	Reference Points
		Fcurr	YPR	F0.1	Fmax	SPR(%)	FSPR35%	FSPR30%	Fcrit
Case 1	amax = 35	0.504	28.327	0.437	0.975	8.578	0.211	0.243	0.139
Case 2		0.497	25.160	0.455	1.143	9.169	0.213	0.245	0.079
Case 3		0.507	24.506	0.438	1.004	8.355	0.211	0.242	0.116
Case 4		0.506	23.450	0.440	1.031	8.405	0.211	0.242	0.102
Case 1	amax = 30	0.497	25.162	0.455	1.142	9.260	0.214	0.246	0.070
Case 2		0.487	21.514	0.480	1.620	10.078	0.217	0.249	-
Case 3		0.507	24.506	0.438	1.004	8.467	0.211	0.243	0.101
Case 4		0.506	23.450	0.440	1.031	8.520	0.212	0.243	0.087
Case 1	amax = 25	0.487	21.402	0.481	1.648	10.284	0.218	0.251	-
Case 2		0.474	17.293	0.519	4.199	11.493	0.223	0.256	-
Case 3		0.507	24.506	0.438	1.004	8.688	0.213	0.245	0.074
Case 4		0.506	23.450	0.440	1.031	8.744	0.213	0.245	0.060

. All the estimates of YPR, SPR and corresponding BRPs fluctuated largely with M and longevity. The lowest estimates of F_curr was obtained for the M scenario that assumed the lowest value of longevity. This implied that the low longevity contributed to a low fishing mortality. The estimated range of YPR was 17.293–28.327 kg, F_0.1 was 0.437–0.519 y⁻¹ and F_max was 0.975–4.199 y⁻¹. SPR analysis indicated that the current SPR was between 8.405% and 11.493%. The estimated range of F_SPR35%r, F_SPR30%r, and F_crit were 0.211–0.223 y⁻¹, 0.242–0.256 y⁻¹ and 0.060–0.139y⁻¹, respectively. In some cases, however, the F_crit cannot be estimated because of the high value of M, particularly those relying on longevity, particularly in the case of $a_{\max }$ = 25 (Table 4). For the YPR model, F_curr was higher than the corresponding biological reference points F_0.1, but was lower than F_max. However, the results from the SPR analysis showed that the current SPR% was significantly lower than target (SPR35%) and limit (SPR30%) levels (Table 4). To sum up, aside from the case of F_max, current fishing mortality was greater than any level of BRP suggesting that bigeye thresher stock was experiencing overexploitation.

3.1.4. Population Increase Rate

The population increase rate (λ) was estimated ranging from 0.964 to 1.039 y⁻¹. The results based on the longevity of 35 y indicated λs were higher than those of 30 and 25 y (

Table 5. The demographic outputs of each scenario from deterministic models.

Natural Mortality	Longevity	Population Increase Rate
		F = 0	F = Fcurr
Case 1	amax = 35	1.039	0.913
Case 2		1.022	0.900
Case 3		1.031	0.911
Case 4		1.027	0.909
Case 1	amax = 30	1.020	0.892
Case 2		0.998	0.876
Case 3		1.029	0.901
Case 4		1.024	0.898
Case 1	amax = 25	0.993	0.863
Case 2		0.964	0.840
Case 3		1.022	0.885
Case 4		1.018	0.882

). However, even without fishing mortality, some cases still resulted in λ less than 1 (Table 5). In addition, the analyses also indicated that the stock would almost certainly decrease under current fishing conditions (Table 5).

3.2. Estimates with Uncertainty

3.2.1. Model Convergence

The posterior mean and standard deviation of fishing mortality and selectivity obtained from MCMC are listed in Appendix A

Table A1. Estimated fishing mortality (posterior mean and standard deviation, year⁻¹) based on four empirical estimators of natural mortality for the bigeye thresher shark.

	Fishing Mortality
Age/Case	Case 1		Case 2		Case 3		Case 4
	Mean	SD	Mean	SD	Mean	SD	Mean	SD
1	0.0001	0.0000	0.0000	0.0000	0.0001	0.0000	0.0000	0.0000
2	0.0004	0.0000	0.0004	0.0000	0.0004	0.0000	0.0004	0.0000
3	0.0021	0.0000	0.0020	0.0000	0.0021	0.0000	0.0021	0.0000
4	0.0094	0.0000	0.0091	0.0000	0.0093	0.0000	0.0092	0.0000
5	0.0331	0.0000	0.0322	0.0000	0.0329	0.0000	0.0327	0.0000
6	0.0918	0.0000	0.0896	0.0000	0.0915	0.0000	0.0911	0.0000
7	0.2005	0.0001	0.1967	0.0001	0.2006	0.0001	0.2001	0.0001
8	0.3457	0.0002	0.3399	0.0002	0.3468	0.0002	0.3461	0.0002
9	0.4700	0.0004	0.4628	0.0004	0.4723	0.0004	0.4716	0.0004
10	0.5042	0.0007	0.4965	0.0007	0.5071	0.0007	0.5063	0.0007
11	0.4266	0.0009	0.4197	0.0009	0.4291	0.0009	0.4282	0.0009
12	0.2848	0.0009	0.2795	0.0009	0.2861	0.0009	0.2853	0.0009
13	0.1499	0.0007	0.1467	0.0006	0.1504	0.0006	0.1497	0.0006
14	0.0623	0.0004	0.0606	0.0003	0.0623	0.0004	0.0619	0.0004
15	0.0204	0.0002	0.0198	0.0001	0.0203	0.0001	0.0202	0.0001
16	0.0053	0.0000	0.0051	0.0000	0.0052	0.0000	0.0052	0.0000
17	0.0011	0.0000	0.0010	0.0000	0.0011	0.0000	0.0010	0.0000
18	0.0002	0.0000	0.0002	0.0000	0.0002	0.0000	0.0002	0.0000
19	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
20	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
21	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
22	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
23	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
24	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
25	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
26	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
27	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
28	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
29	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
30	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
31	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
32	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
33	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
34	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
35	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000

and

Table A2. Estimated selectivity (posterior mean and standard deviation) based on four empirical estimators of natural mortality for the bigeye thresher shark.

	Selectivity
Age/Case	Case 1		Case 2		Case 3		Case 4
	Mean	SD	Mean	SD	Mean	SD	Mean	SD
1	0.0001	0.0000	0.0001	0.0000	0.0001	0.0000	0.0001	0.0000
2	0.0007	0.0000	0.0007	0.0000	0.0007	0.0000	0.0007	0.0000
3	0.0042	0.0000	0.0041	0.0000	0.0041	0.0000	0.0041	0.0000
4	0.0187	0.0000	0.0184	0.0000	0.0184	0.0000	0.0183	0.0000
5	0.0657	0.0001	0.0649	0.0001	0.0648	0.0001	0.0646	0.0001
6	0.1820	0.0002	0.1806	0.0002	0.1804	0.0002	0.1800	0.0002
7	0.3977	0.0005	0.3961	0.0005	0.3957	0.0005	0.3952	0.0005
8	0.6856	0.0007	0.6845	0.0007	0.6838	0.0007	0.6836	0.0007
9	0.9323	0.0006	0.9321	0.0006	0.9315	0.0005	0.9315	0.0005
10	1.0000	0.0000	1.0000	0.0000	1.0000	0.0000	1.0000	0.0000
11	0.8462	0.0007	0.8453	0.0007	0.8461	0.0007	0.8457	0.0007
12	0.5648	0.0010	0.5630	0.0010	0.5642	0.0010	0.5635	0.0010
13	0.2974	0.0009	0.2954	0.0009	0.2965	0.0009	0.2958	0.0009
14	0.1235	0.0005	0.1221	0.0005	0.1228	0.0005	0.1223	0.0005
15	0.0405	0.0002	0.0398	0.0002	0.0401	0.0002	0.0398	0.0002
16	0.0105	0.0001	0.0102	0.0001	0.0103	0.0001	0.0102	0.0001
17	0.0021	0.0000	0.0021	0.0000	0.0021	0.0000	0.0021	0.0000
18	0.0003	0.0000	0.0003	0.0000	0.0003	0.0000	0.0003	0.0000
19	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
20	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
21	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
22	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
23	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
24	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
25	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
26	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
27	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
28	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
29	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
30	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
31	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
32	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
33	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
34	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
35	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000

. As listing the values for all of the convergence statistics for all of the parameters is not practical, this study presents Appendix A Figure A1, Figure A2, Figure A3, Figure A4 to demonstrate the convergence statistics for major parameters. In Appendix A Figure A1, Figure A2, Figure A3, Figure A4, the convergence statistics suggest the trace of the posterior samples and the posterior density function, which is estimated by using a normal kernel density estimator. The trace and the cumulative patterns do not show any obvious patterns; meanwhile the posterior density functions appear smooth and unimodal. In short, the trace and density plots of the major parameters do not indicate any lack of convergence (Appendix A Figure A1, Figure A2, Figure A3, Figure A4).

3.2.2. Biological Reference Points

The box plots of the F_BRPs derived from YPR and SPR models in 4 cases (4 different M) are shown in Figure 2. Similar to deterministic model, the lowest and highest estimates of BRP obtained from Per-Recruit models were found in Case 1 and Case 2, respectively (Table 4). After considering the uncertainty of fishing mortality in per recruit analysis, the median values of all BRPs, appeared close to those estimated from deterministic models. While the medians of these quantities remained close to those of the deterministic case for most scenarios, there was a large variation in values which highlighted the need for taking uncertainty in estimating BRPs into account. The additional BRPs (F_crit) estimated from demographic model also showed that variation in M will affect the level of fishing mortality that a population can sustain without decline (Figure 3). However, in contrast to per-recruit based BRPs, the lowest and highest estimates of BRPs (F_crit) were produced by Case 2 and Case 1, respectively (Figure 3).

3.2.3. Population Increase Rate

The stochastic estimates of demographic parameters for bigeye threshers obtained from the seven scenarios are shown in

Table 6. The demographic outputs of each scenario from stochastic models.

Scenario	Type of F	$\mathit{\lambda }$	Lower CL	Upper CL	r	Lower CL	Upper CL
1	F = 0	1.023	1.010	1.039	0.023	0.010	0.039
2	Fcur	0.906	0.894	0.915	−0.098	−0.112	−0.089
3	F0.1	0.919	0.905	0.928	−0.085	−0.100	−0.075
4	Fmax	0.815	0.789	0.829	−0.204	−0.237	−0.187
5	FSPR35%	0.969	0.957	0.981	−0.031	−0.044	−0.019
6	FSPR30%	0.962	0.950	0.974	−0.039	−0.052	−0.027
7	Fcrit	0.995	0.974	1.017	−0.006	−0.026	0.017

Values in parentheses are lower and upper 95% confidence intervals calculated as the 2.5th and 97.5th percentiles.

. Figure 4 shows the box plots of λ for Scenarios 1–7. As expected, the lowest λ was derived from the Scenario 4, which suggested that F_max is not a good management reference point for bigeye threshers. The large variation of λ reflects the uncertainty of input parameters. The simulation results clearly indicated that the stock would increase (λ = 1.023, 95% C.I. = 1.010–1.039 y⁻¹) without fishing mortality. Nevertheless, under current fishing conditions (Scenario 2), a λ less than 1 (λ = 0.906, 95% C.I. = 0.894–0.915 y⁻¹) was produced (Table 6). Both of the per-recruit-based BRPs (Scenarios 3–6) management strategies still resulted in clear population declines. Only the demographic-based BRP (Scenario 7) resulted in a relatively stable population (λ = 0.995, 95% C.I. = 0.974–1.017 y⁻¹).

4. Discussion

4.1. Biological Reference Points

Setting fishing mortality at YPR-based BRPs (F_max and F_0.1) indicated that it produced λ less than 1 (Table 6), implying that it may not be a suitable BRP candidate for the management of bigeye thresher sharks. Similarly, SPR-based simulations (F_SPR30% and F_SPR35%) also produced λ less than 1 (Table 6). All of these findings also imply that conservative BRPs are more appropriate for this species. The demographic approach is a preferable alternative stock assessment tools for pelagic shark fishery management because it provides additional information regarding population responses to fishing at levels different from the reference point. Furthermore, when population growth rates are different from zero, equal levels of SSB/R do not result in the same population growth rate for different partial recruitment vectors [13]. Comparing the BRPs derived from demographic models with those from YPR models showed that F_max and F_0.1 were not appropriate BRP candidates for the management of bigeye thresher sharks. The analysis of SPR-based BRPs implied that more conservative BRPs are needed for this species. Clarke and Hoyle (2014) [46] recommended the use of a proxy benchmark (limit F_BRP) of at least F_SPR50% for long-lived and low-productivity shark stocks, and F_SPR60% for shark species having very low compensation after the removals by fishery (e.g., species with a particularly low natural mortality or steepness). Consistently, our analyses also demonstrated that the Northwest Pacific bigeye thresher population will only stabilize if demographic-based BRP is implemented. The most likely estimate of BRP (F_crit) based on the demographic model was 0.079–0.139 y⁻¹, which is equivalent to SPR = 50–70% (

Table 7. Deterministic estimates of critical fishing mortality (based on demographic model) and their corresponding SPR% for the bigeye thresher shark in the Northwest Pacific Ocean.

Natural Mortality	Longevity	Fcrit	Corresponding SPR%
Case 1	amax = 35	0.139	49.907
Case 2		0.079	67.644
Case 3		0.116	55.989
Case 4		0.102	60.044
Case 1	amax = 30	0.070	70.723
Case 2		-	-
Case 3		0.101	60.447
Case 4		0.087	64.896
Case 1	amax = 25	-	-
Case 2		-	-
Case 3		0.074	69.099
Case 4		0.060	74.262

). The findings reported herein strongly suggest that more conservative threshold BRPs should be implemented to ensure sustainable utilization of the bigeye thresher stock.

4.2. Demographic Model

Demographic matrix population models such as age-structured (also known as Leslie Matrix) and stage-structured models are commonly used in the assessment of shark populations. The choice between age- and stage-structured models is basically depending on personal preference. Both approaches will provide similar results if the same life history parameters are used [42]. In some situations, the life history of a shark species can be represented by several discrete stages (e.g., neonate, juvenile, sub-adult, pregnant adults, and resting adults for sandbar sharks) [42,47]. In this case, the stage-based model can be useful if there is only limited age information for a species or complex reproductive physiologies exhibit in the life history (e.g., resting stages and extended gestation periods) [16,48]. In the present study, however, most information is age-based, such as natural mortality, fishing mortality, and age-at-maturity. It would be more consistent and reasonable to use an age-structured matrix model to interpret assessment results. The demographic model adopted in this study was a single-sex model carried out exclusively for females and did not consider the density-dependent compensatory effects for the population. Tsai et al. (2014, 2015) [16,48] demonstrated that the probability of population decline may be underestimated based on single-sex demographic models when life history parameters differ between sexes. As no significant difference in vital parameters between sexes was found for the bigeye thresher (Tsai et al., unpub. data), only the female population dynamic was taken into account in this study. On the other hand, density dependence effect may result in decrease in reproductive output because a consequence of an earlier age at maturity or a decreased asymptotic size because of faster individual growth rate [49]. However, the bigeye thresher is a viviparous species, usually producing two pups at a time and the litter size does not change with maternal size [39]. Therefore, the compensation on reproductive output for this species, if it exists, is likely to be negligible.

4.3. Uncertainty

Natural mortality and longevity may also be factors affecting the results of our analysis. There is currently no direct information to estimate natural mortality for bigeye thresher shark. Unfortunately, estimating natural mortality for shark species is often difficult as they are widely distributed and highly migratory. Many empirical equations have typically been developed to estimate natural mortality. In general, the use of multiple indirect methods has been applied in many demographic studies [15,43,48], which may reduce the bias imposed by any one method. Such methods usually rely on longevity (e.g., [29,30]), age-at-maturity [50] or other growth parameters [51]. However, it is also difficult to accurately estimate longevity for shark species. Although the uncertainties of longevity were not considered in the present study, sensitivity analysis of longevity to examine the possible effects on demographic estimates was adopted (Table 4 and Table 5).

The estimated population increase rate by demographic analysis from the previous studies ranged between 1.008 and 1.046 [20,52,53] with the assumptions of longevity of 35 years [20] and 28 years [52,53] for females in their demographic analysis. The variation of estimated λ may have resulted from different methods used in estimations of M and longevity. The longevity of 35 years was estimated by substituting the maximum observed length of female into the growth equation of Liu et al. (1998) [19]. However, a recent study [54] demonstrated that previous methods used to determine the age of sharks, such as vertebral band counting, have underestimated those ages, particularly in older sharks. Therefore, the longevity of 35 years for female bigeye thresher sharks is believed to be a reasonable estimation.

4.4. Stock Status

Most sharks and their relatives are usually characterized as slow-growing species. Furthermore, Musick (1999) [55] concluded that species with annual intrinsic growth rates less than 10% tend to be particularly vulnerable to increases in fishing mortality. The most likely demographic models for bigeye thresher shark produced a mean λ of 1.023 y⁻¹ (Scenario 1, Table 6) under natural conditions in this study. The low mean rates of λ showed that bigeye thresher sharks have low tolerance of exploitation and will recover slowly from fishery induced mortality and have high risk of extinction [6,56]. The estimated λs for bigeye thresher sharks are extremely low, and any added source of mortality to this population will likely result in a population decline since even under stable conditions the population growth rate was ~2% per year (Table 6). The bigeye thresher shark in the Northwest Pacific was identified as one of the least productive and most vulnerable shark species, with a significantly low population increase rate, low intrinsic rate of population growth of 0.023 y⁻¹, and generation time of 19.63 years. These demographic factors arguably make the bigeye thresher vulnerable to any level of exploitation.

Overestimation of longevity may result in overly optimistic estimates of population growth rate, particularly for long-lived sharks [57]. However, even at the highest longevity (35 years) assumed in this study, the simulations still resulted in clear population declines under current conditions. These finding implies that the Northwest Pacific bigeye thresher stock is declining in population size under current conditions of fisheries, and this conclusion is congruent with the results from per-recruit analyses. Tsai et al. (2019) [20] conducted a risk assessment study of bigeye thresher shark using Bayesian population model in an area subset of the western North Pacific. Their assessment found that the bigeye thresher experienced higher fishing pressure in years 2011–2016 and that current fishing mortality is higher than the target reference point F_0.1 as well as F_SPR35%, suggesting that overfishing is likely occurring for the bigeye thresher shark. This conclusion is consistent with the results obtained from the present study (Scenario 2, Table 6).

5. Conclusions

Our study presents alternative approaches for assessing the population dynamics of pelagic sharks using the bigeye thresher in the Northwest Pacific as an example. The results highlight the high vulnerability of bigeye threshers to fishing pressure and can be used to help fishery managers to adopt more efficient management decisions and conservation measures for this stock. Owing to general lack of catch and effort data, the current fishing pressure of the bigeye thresher shark in the Northwest Pacific Ocean has not yet been tuning with CPUE time series. Better estimates of current fishing level are needed to obtain a more robust estimate of the impact of commercial fishery on the bigeye thresher shark population. Given the increasing trend in global shark catches and landings, the bigeye thresher population should be constantly monitored to ensure their sustainability. It is also suggested that this approach is applicable to other shark species with limited catch and effort data.