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

: Full stock assessment of sharks is usually hindered by a lack of long time-series catch and e ﬀ ort 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 di ﬃ cult 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 ﬁshing 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 Paciﬁc 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 − 1 , which was equivalent to SPR = 50–70%. The ﬁndings 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 Paciﬁc, which can be used to help ﬁshery managers to adopt more e ﬃ cient management measures for this stock. It is also suggested that this approach can be applied to other shark species with limited catch and e ﬀ ort data.


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

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 (Φ W ) was obtained from the whole weight-sex ratio relation: Φ W = α × W β , where α and β 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.  4 r m −0.747 a m 12 1 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 (Φ W ) was obtained from the equation: Φ W = α × W β , where α and β are estimated parameters. (R 2 = 0.979; n = 4855, 5-kg classes, p < 0.0001). 2 Liu et al. (1998). where a and b are estimated parameters for length-weight relationship. 3 Liu et al. (1998). where L ∞ is the maximum attainable length, k is a Brody growth constant, t 0 is a hypothetical age at length of 0. 4 Data from Chen et al. (1997). where r m is the slope and a m is age at maturity estimated from logistic maturity model.

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]: Case2 [30]: Case3 [31]: Case4 [32]: 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×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 (Ĉ a ) of a fish at age a can be estimated from the catch equation [33]: 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: where C a is the observed catch of fish at age a.

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

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

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]: 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).
Sustainability 2020, 12, 8646 6 of 29 (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).
The estimations of BRPs described above are deterministic (set as the reference cases).

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.

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: 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 Sustainability 2020, 12, 8646 7 of 29 by incorporating parameter uncertainty in Monte Carlo simulation. All demographic and simulation analyses were conducted using CSIRO program-PopTools [45].

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

Mortality and Selectivity
The range of age-specific natural mortality (M a ), produced by the four indirect methods, was 0.088-0.199 y −1 . The lowest estimates of M = 0.107 y −1 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 −1 were obtained using the method of Campana et al. (2001) [30], based on age at longevity of 25 years (Table 2). 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). 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).

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. 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 −1 and F max was 0.975-4.199 y −1 . 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 −1 , 0.242-0.256 y −1 and 0.060-0.139y −1 , 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.

Population Increase Rate
The population increase rate (λ) was estimated ranging from 0.964 to 1.039 y −1 . The results based on the longevity of 35 y indicated λs were higher than those of 30 and 25 y (Table 5). 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).

Population Increase Rate
The population increase rate (λ) was estimated ranging from 0.964 to 1.039 y −1 . The results based on the longevity of 35 y indicated λs were higher than those of 30 and 25 y (Table 5). 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).

Model Convergence
The posterior mean and standard deviation of fishing mortality and selectivity obtained from MCMC are listed in Appendix A Tables A1 and A2. As listing the values for all of the convergence statistics for all of the parameters is not practical, this study presents Figures A1-A4 to demonstrate the convergence statistics for major parameters. In Appendix A Figures A1-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 Figures A1-A4).

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). parameters. The simulation results clearly indicated that the stock would increase (λ = 1.023, 95% C.I. = 1.010-1.039 y −1 ) without fishing mortality. Nevertheless, under current fishing conditions (Scenario 2), a λ less than 1 (λ = 0.906, 95% C.I. = 0.894-0.915 y −1 ) was produced ( Table 6). Both of the per-recruitbased 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 −1 ).

Population Increase Rate
The stochastic estimates of demographic parameters for bigeye threshers obtained from the seven scenarios are shown in Table 6. 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 −1 ) without fishing mortality. Nevertheless, under current fishing conditions (Scenario 2), a λ less than 1 (λ = 0.906, 95% C.I. = 0.894-0.915 y −1 ) was produced ( Table 6)

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 −1 , which is equivalent to SPR = 50-70% (Table 7). The findings reported herein strongly suggest that more conservative threshold BRPs should be implemented to ensure sustainable utilization of the bigeye thresher stock.

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. (2014Tsai et al. ( , 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.

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 (Tables 4 and 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.

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 −1 (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 −1 , 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).

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.