Evidence of overfishing of geoduck clam Panopea globosa from a length-based stock assessment approach

Shell length structure and catch data

Biological data were obtained from Puerto Peñasco, Sonora, in the upper Gulf of California, Mexico (Fig. 1), during two time periods: 2010 to 2012 and 2014 to 2016. A total of 19,445 individuals of P. globosa were collected from two sources of information. The first source was a sampling design based on the estimation of the fishing ground. To this goal, a survey was conducted to identify the patches and beds with high abundance; the fishing ground was estimated based on the boundaries of sampling stations, and the area was expressed as km².

The second source of biological data was obtained from commercial landings, which were endorsed by fishing licenses that establish basic conditions for the extraction activities. During the predevelopment phase of geoduck fishery, a legal requirement for the fishers was annual data collection and analysis of the harvest population; however, when the analysis of the fishing ground indicates the availability of sufficient biomass for exploitation, then the fishing ground is classified as being in growth phase. Consequently, the data collection is not mandatory for the fishers, and the biological information may be limited due to the cost of the geoduck monitoring program, which was not implemented during 2013. According to the bathymetry around Puerto Peñasco (Ramírez-Mendoza & Álvarez, 2009), the commercial harvest was performed between 10 and 30 m deep; greater depths are restricted by the Mexican government to avoid the decompression sickness of divers. Geoduck harvest was performed by hookah diving during low tides, with small boats equipped with an outboard motor and a low-pressure compressor connected to a hose. Geoducks were located through their siphon holes on the substrate, and they were removed using a water jet to loosen the sediment. Shell length (mm) and total weight (g) data were recorded monthly; the specimens were grouped by shell length frequency distributions by year, and the size classes varied from 78 to 200 mm shell length.

Data origin

The annual catch data were obtained from the Subdelegación Federal de Pesca at Puerto Peñasco, Sonora, Mexico (Comisión Nacional de Acuacultura y Pesca, CONAPESCA) (https://www.conapesca.gob.mx/wb/cona/estadisticas_de_produccion_pesquera).

Scientific collecting and the use of biological and fishery data was endorsed by the Mexican federal government (CONAPESCA) through the commercial fishery permit numbers DGOPA.05338.050710.3326 (2010), DGOPA.09450.221111.3442 (2011), 126054025018 (2012), PPF/DGOPA-186/14 (2014), 126070025019 (2015), and 126047025017 (2016) with fishermen of the SCPP Buzos de Puerto Peñasco, SCPAP Islas de Sonora, SCPP Jaiberos y Escameros, and SCPP Mar y Tierra del Golfo de Cortez.

Population dynamics

The population dynamics of P. globosa was analyzed using an integrated catch-at-size assessment (ICSA) model (Morales-Bojórquez, López-Martínez & Beléndez-Moreno, 2013). To provide a better description of the model, a summary of the symbols and descriptions of the parameters and variables are defined in Table 1. The ICSA model describes the population changes over time for multiple shell length classes according to the following: (a) the total mortality of geoducks in shell length class l at time t, representing the sum of fishing and natural mortality (Z_l,t = F_l,t + M_l,t), and (b) the relationship between N_l,t and N_l,t′, which denotes the number of geoducks at shell length surviving and growing during the next time period according to the equation N_l,t′ = N_l,te^−Z_l,t. A natural mortality value of 0.046 was used for modeling the population dynamics of geoduck clam (González-Peláez et al., 2015a); this fixed parameter was based on the number of gnomonic intervals determined by dividing the life history of the geoduck into a predetermined number of biological units (life stages). The advantage of this method is that uses the biological information of the lifespan for knowledge of the ontogenic changes in natural mortality; this procedure is better than estimations based on meta-analysis because these approaches were developed for specific taxa, such as fishes (Pauly, 1980), empirical approaches based on biological parameters (Hoenig, 1983) or based on life strategy (r-K selection) (Gunderson, 1980). A gnomonic interval is a systematic strategy for the unequal subdivision of the lifespan of an individual into time intervals, which increase in duration and proportion for each time interval. Therefore, in terms of elapsed time, two gnomonic intervals in a subdivided life history can be considered equivalent if they each form the same constant proportion of the time elapsed (Caddy, 1996; Ramírez-Rodríguez & Arreguín-Sanchez, 2003; Martínez-Aguilar, Arreguín-Sánchez & Morales-Bojórquez, 2005; Aranceta-Garza et al., 2016; Romero-Gallardo et al., 2018). Thus, the natural mortality of P. globosa was modeled as follows: (1) egg to trochophore larvae (24 h); (2) early larvae (6.5 days); (3) late larvae (11 days); (4) early juveniles (35 days); (5) juveniles (3–9 months); (6) late juveniles (1–2 years); and (7) preadult to adults (3–47 years). Basically, the fishery maintains high fishing pressure on biological unit 7; such that M = 0.046 is a plausible value; additionally, the comparison between the gnomonic method and alternative procedures yielded similar values for geoducks larger than 130 mm; therefore, the natural mortality was relatively stable for individuals from preadult to adult (Calderon-Aguilera et al., 2010; González-Peláez et al., 2015a). The catch ${\bar{C}}_{l,t}$ was calculated using the Baranov catch equation, ${\bar{C}}_{l,t}=N_{l,t}{\mu }_{l,t}$, where μ_l,t is the harvest rate defined as the proportion of individuals that die from fishing mortality compared to total mortality ${\mu }_{l,t}=F_{l,t}∕Z_{l,t}\left(1-e^{-Z_{l,t}}\right)$.

Initial conditions

The number of individuals of the geoduck clam population at the beginning of the time period denotes the initial stock abundance, which was estimated as N₀ = ∑_lN_l,0, where N_l,0 was estimated as $N_{l,0}=\frac{C_{l,0}\left(Z_{l,0}\right)}{F_{l,0}\left(1-e^{-Z_{l,0}}\right)}$. In this way, the proportion of geoducks in shell length class l was calculated as $P_l=\frac{C_{l,0}{Z}_{l,0}∕F_{l,0}\left(1-e^{-Z_{l,0}}\right)}{\sum C_{l,0}{Z}_{l,0}∕F_{l,0}\left(1-e^{-Z_{l,0}}\right)}$. Thus, the shell length distribution of abundance was computed as ${\hat{N}}_{l,0}={\hat{P}}_l{N}_0$. This procedure simplifies the parameter space for estimating only one N₀ instead of estimating a vector of N_l,0, increasing the performance of the objective function.

Recruitment

Usually, estimates of recruitment are based on a stock-recruits relationship, as in the studies by Beverton & Holt (1957), Ricker (1954) and Ricker (1975). However, establishing a functional stock-recruits relationship for the geoduck clam is very difficult; therefore, recruitment in the Panopea genus is only measured in terms of recruitment rates (Zhang & Hand, 2006). Sullivan, Lai & Gallucci (1990), Fisch et al. (2019) and Amezcua-Castro et al. (2019) used length-structured models and proposed that recruitment can be estimated from a gamma probabilistic density function, assuming that it can represent the variability in recruitment. Biologically, recruitment represents the number of individuals at specific age or length classes that are added to the exploitable stock in the fishery (Myers, 2002). Recruitment is modeled as the product of a time-dependent variable (R_t) corresponding to the annual recruitment and a length-dependent variable (φ_l), representing the proportion of the annual recruitment into each length class l. This proportion φ_l of recruits in each length class was estimated following a gamma distribution with recruitment parameters α_r and β_r. Thus, the equation for estimating recruitment-at-shell length was R_l, = R_tφ_l (Fisch et al., 2019); this procedure of separating variables allows R_t to be compared with recruitment estimates from standard procedures. According to Sullivan, Lai & Gallucci (1990), recruitment specified in this way represents the type of recruitment observed in nature, where variation in growth, behavior, or food supply can result in individuals entering the population at several length classes. Given that the recruitment to the fishery occurs over a range of shell length classes, we grouped the number of individuals with shell length between 78 and 130 mm (which were caught from the fishery) to represent the recruits added to fishable stock.

Fishing mortality and selectivity

The assumption in the ICSA model is that the fishing mortality vector can be partitioned into two components: (i) a shell length-specific component that does not vary over time (e.g., a constant exploitation pattern) and (ii) an annual multiplier commonly specified as catchability, vulnerability or selectivity (Megrey, 1989). Thus, in this study, the fishing mortality was estimated as a separable product of the shell length selectivity and the full-recruitment fishing mortality rate: F_l = s_lf_t. The separability assumption allows the estimation of only a fishing mortality value, which is distributed proportionally in the shell length classes, improving the performance of the model parameterization (Fisch et al., 2019). To analyze s_l, a gamma probabilistic function was applied; the advantage of this approach is that it provides great flexibility to the model, allowing different selectivity patterns (Carlson & Cortés, 2003). Thus, s_l was estimated as $s_l={\left(l∕{\alpha }_s{\beta }_s\right)}^{{\alpha }_s}exp\left({\alpha }_s-l∕{\beta }_s\right)$. The parameterization was based on the residual sum of squares (RSS_s) function: (1)${\displaystyle RS{S}_s=\sum _{i=1}^n{\left(s_l-s_{o,l}\right)}^2}$

where s_o,l is defined in Table 1.

Individual growth and survival

Individual growth and survival were modeled by using a transition matrix, which represents the combination of a growth matrix and a survival matrix (Table 1). According to Cao, Chen & Richards (2016), the growth increments assume shell length variability, which can be described as Δ_l = l_t+1 − l_t. Therefore, the mean growth increments were expressed by the Fabens model (1965): ${\bar{\Delta }}_l=\left(L_{\infty }-l_{\ast }\right)\left(1-e^{-k}\right)$. The growth matrix was based on a gamma distribution expressed as $g\left({\Delta }_l|{\alpha }_l{\beta }_g\right)=\frac{1}{{\beta }_g^{{\alpha }_l}\Gamma \left({\alpha }_l\right)}{\Delta }_l^{{\alpha }_l-1}{e}^{-{\Delta }_l∕{\beta }_g}$. The expected proportion of individuals growing from shell length class l to shell length class l + 1 was found by integration of the shell length ranges, l + 1₁, l + 1₂, representing the lower and upper ends of shell length classes, respectively: (2)${\displaystyle G_{l,l+1}={\int }_{l+1_1}^{l+1_2}\mathrm{g}\left(\mathrm{x}|{\alpha }_l,{\beta }_{\mathrm{g}}\right)\mathrm{dx}.}$

The survival matrix was estimated as e^−Z_l,t, and the transition matrix was calculated as follows: (3)${\displaystyle T_{l,t}=\sum _l{G}_{l,l+1}{S}_{l,t}.}$

Parameter estimation

The ICSA model was fitted through the residual sum of squares (RSS) criterion: (4)${\displaystyle RSS=\sum _{l,t}{\left({\bar{C}}_{l,t}-C_{l,t}\right)}^2.}$

The parameters were estimated by minimizing the objective function (RSS) with a nonlinear fit using the Generalized Reduced Gradient (GRG) algorithm (Neter et al., 1996). The optimization was performed in phases; this procedure consists of estimating only a subset of parameters through the optimization of the objective function and adding more parameters in each sequential phase until all parameters are estimated. An advantage is that to the extent that new parameters are included into a new phase, the previously estimated parameters are still estimated, and they are not fixed as prespecified values, allowing the progressive improvement of the goodness of fit of the RSS (Legault & Restrepo, 1998; Luquin-Covarrubias et al., 2016a; Luquin-Covarrubias et al., 2016b). A non-linear fit in phases is an statistical procedure accepted when there is high uncertainty in seed values and where the number of parameters to be estimated is greater than 20 (Fournier et al., 2012; Punt, Huang & Maunder, 2013; Canales, Company & Arana, 2016; Cao, Chen & Richards, 2016; Fisch et al., 2019). According to Fournier et al. (2012), this procedure is also called “model sculpting”. Thus, the F_l,t vector was known when f_t was included in the numerical optimization process in the first group of parameters; the second group of parameters included N_l,0, the third group integrated R_t, α_r and β_r, and the fourth group incorporated β_g. The α_s and β_s parameters were independently estimated and later integrated into the RSS as a fifth group. For ICSA model if any parameters are known, then they can be prespecified instead of estimated (Sullivan, Lai & Gallucci, 1990); thus, the parameter L_∞ = 200 mm was assigned according to the largest observed shell length, k = 0.25 was assigned from the von Bertalanffy growth model, reported by Aragón-Noriega, Calderon-Aguilera & Pérez-Valencia (2015), and M = 0.046 was assigned based on the report by González-Peláez et al. (2015a).

Management quantities

In this study, the total abundance-at-shell length of the geoduck clam population was estimated using T_l,t and ${\hat{N}}_{l,0}$. The recruitment-at-shell length was added to represent the population dynamics of the geoduck clam as follows: (5)${\displaystyle N_{l^{{}^{\prime }},t^{{}^{\prime }}}=\sum _l{T}_{l,t}{\hat{N}}_{l,0}+R_{l,t}}$

The total and vulnerable biomass-at-shell length were estimated as follows:

(6)${\displaystyle T{B}_l=\sum _{i=1}^n{N}_{l^{{}^{\prime }},t^{{}^{\prime }}}{\omega }_l}$ (7)${\displaystyle V{B}_l=\sum _{i=1}^n{N}_{l^{{}^{\prime }},t^{{}^{\prime }}}{\omega }_l{s}_l}$

where ω_l was estimated from the power equation ω_l = αl^β (González-Peláez et al., 2015b).

Shell length structure

The time series showed that from 2010 to 2012, the shell length frequency distributions were similar; the most harvested shell length classes were between 126 and 174 mm, while the least harvested shell length classes were between 178 and 198 mm. From 2014 to 2016, a marked change was observed in the total shell length frequency distribution. We noted an increase in harvest of individuals with a shell length less than 110 mm. Although the shell length classes mostly harvested (118 and 178 mm) and less predominant (182 and 198 mm) were similar to those caught during the first three years, we observed that shell length frequency of geoduck clam was fragmented from 2014–2016, and the bell shape of the shell length frequency observed during 2010–2012 changed to an irregular shape during 2014–2016. During the first three years, the total number of individuals caught per shell length class decreased. In contrast, during the second period, the catches increased along the shell length frequency distribution; thus, two periods were clearly distinguishable (Fig. 2).

Recruitment

The recruitment estimated for P. globosa showed a negative trend from 2010 to 2012, with values of 2.78 × 10⁶, 1.93 × 10⁶ and 3.75 × 10⁵, respectively. This sudden decrease to the lowest level represented a change by one order of magnitude in 2012 and a loss of 86% of the recruitment abundance, indicating a decrease over a very short period of time. Subsequently, recruitment increased from 8.85 × 10⁵(2014) to 1.76 × 10⁶ (2015) and 3.02  × 10⁶ (2016), reaching recruitment levels slightly greater than those estimated in 2010. Figure 3 shows the changes in number of recruits to the fishery (from 78 to 130 mm) estimated in our study.

Fishing mortality and selectivity

The fishing mortality estimates calculated by the ICSA model varied over time. The geoducks of 166 mm shell length experienced the main fishing pressure along the time series, exhibiting the highest fishing mortality during 2016 (0.79), followed by 2011 (0.63), 2010 (0.62) and 2012 (0.42). In 2014, fishing mortality showed two peaks at 0.36 in the 158 and 166 mm shell length classes. Finally, during 2015, a peak of 0.22 was observed for the 158 mm shell length class. Comparatively, the lowest fishing mortality (<0.05) occurred in individuals in the shell length classes of 78–130 mm (2010–2016) and those larger than 190 mm (2011–2015). Fishing mortality increased in the 194 mm shell length class only during 2010 and 2016 (Figs. 4A, 4B).

The selectivity showed a decrease in shell length from 150 (2012) to 142 mm (2016). This change was observed in the shell length selectivity of geoducks that had a 50% probability of being caught by a diver; this value was defined as S_l50% (Fig. 5). The selectivity estimates for the analyzed years showed the lowest values (<0.07) for the shell length classes between 78 and 118 mm. Changes in selectivity from 2010 to 2015 showed an accelerated increase from 0.2 to 0.99 for individuals with shell lengths between 122 and 174 mm. Comparatively, during 2016, the changes in selectivity included smaller individuals and a wider shell length range (110–174 mm), with values between 0.03 and 0.99. Finally, the geoducks in the shell length classes larger than 178 mm attained the maximum selectivity.

Management quantities

Harvest rate

The tendency of the harvest rates of geoducks between 78 and 115 mm varied slightly throughout the time series, exhibiting values less than 0.05; subsequently, the harvest rate estimates gradually increased for geoducks from 118 to 130 mm, reaching a maximum value of 0.16 (2014) for geoducks of 130 mm (Figs. 6A, 6B). For individuals larger than the minimum legal size, an increasing trend in the harvest rate-at-shell length was also observed. Thus, the individuals between 154 and 178 mm reached maximal values of 0.48 (2010), 0.49 (2011), 0.36 (2012), 0.37 (2014), 0.26 (2015) and 0.56 (2016). A decreasing harvest rate-at-shell length occurred from 2011 to 2015 for the shell length classes between 182 and 198 mm; however, for 2010 and 2016, a patch of large geoducks was located and harvested; consequently, the harvest rate of the 194 mm shell length class suddenly increased. Although the fishing pressure on geoducks, expressed as the harvest rate-at-shell length, was highly variable among the years analyzed, we observed a general pattern over time: the shell length classes between 142 and 174 mm accounted for more than 50% of the harvest (Figs. 7A, 7B).

Total and vulnerable biomass

The estimates of total biomass-at-shell length showed changes along the time series; from 2010 to 2012, the total biomass-at-shell length decreased; subsequently, an increase was observed during 2014 and 2015, with a slight decrease in 2016 (Table 2). In the first period, the maximum estimates of total biomass were 204.54 t (2010) and 161.10 t (2011) for the shell length class of 122 mm; then, the total biomass-at-shell length rapidly decreased to 65.41 t (2012) for geoducks in the shell length class of 146 mm. During the second period (2014–2016), the highest values of total biomass-at-shell length were 266.66 t (2014), 265.75 t (2015) and 224.36 t (2016) for the shell length classes of 158 mm, 146 and 126 mm, respectively (Figs. 8A, 8B). The shell length ranges most representative of geoduck clams associated with the total biomass were different for each time period analyzed. For the first period, the range of shell length changed from 110-138 mm (2010–2011) to 130–170 mm (2012). This change indicated that the total biomass of younger individuals decreased, and during the third year, although the individuals were larger, their total biomass was diminished, representing the year with the lowest biomass. During the second period, the range of shell length progressively diminished, and the ranges among years were 138-170 mm (2014), 130–138 mm (2015), and 110–142 mm (2016) (Figs. 8A, 8B).

The time series showed that the vulnerable biomass-at-shell length decreased from 2010 to 2012; subsequently, during 2014 and 2015, the vulnerable biomass exhibited an increase, diminishing in 2016 (Table 2). During the first period, the peaks of vulnerable biomass were 72.93 t (2010) and 52.77 t (2011) for the 154 mm shell length class and 43.62 t during 2012 for individuals with a shell length of 162 mm. Comparatively, during the second period, the highest values of vulnerable biomass were 237.90 t (2014), 170.66 t (2015), and 99.72 t (2016) for the 178, 162 and 150 mm shell length classes, respectively. In the time series, vulnerable biomasses less than 25 t were observed for individuals between 78 and 110 mm and larger than 186 mm, excluding 2014, when the vulnerable biomass for individuals between 190 and 198 mm was not estimated because these shell length classes were not harvested. Conversely, during 2015, the vulnerable biomass for the individuals in the 198 mm shell length class showed a slight increase (Figs. 9A, 9B).

Parameter estimation

The parameters associated with the population dynamics of the geoduck clam for each year of fishing analyzed are shown in Table 3. The interaction between the different modules of the ICSA model exhibited convergence in the objective function. Thus, the ICSA model showed high performance in fitting the observed catch-at-shell length data, exhibiting low residual values (RSS) (Table 3).