Natural mortality estimates throughout the life history of the sea cucumber Isostichopus Badionotus (Holothuroidea: Aspidochirotida)

A natural mortality coefficient in sea cucumbers

The natural mortality coefficient (M) is an important input parameter in age-structured models or other stock assessment models that actually uses the M as their primary input parameter; also as sustainability indicators, optimal exploitation rates and other biological reference points (Siegfried & Sansó, 2009; Brodziak et al., 2011). The M coefficient includes the result of all possible causes of natural death other than fishing, and direct estimates can be calculated for exploited and unexploited stocks (McPherson, García-Rubies & Gordoa, 2000; Bacheler et al., 2009). The M coefficient can be calculated using tag-recapture and telemetry methods (Pollock, Jiang & Hightower, 2004; Polacheck et al., 2006; Hewitt et al., 2007; Bacheler et al., 2009), visual population census data, catch-at-age models (McPherson, García-Rubies & Gordoa, 2000; Giménez-Hurtado, Arreguín-Sánchez & Lluch-Cota, 2009), and mortality by predation using Multispecific Virtual Population Analysis (MVPA) (Sparholt, 1990). Other methods such as integral evaluation models or those derived from ecosystem models (Walters, Christensen & Pauly, 1997; Arreguín-Sánchez, 2000) can also be applied when large data sets are available. Given the complexity of direct methods, the use of indirect methods based on life span estimates can provide an adequate solution to infer values of M while validating these values through direct methods (Then et al., 2014).

The natural mortality coefficient M may vary during an organism’s life cycle (Walters, Christensen & Pauly, 1997; Jennings, Kaiser & Reynolds, 2009); especially in species with high fecundity, higher M values are obtained during early life history stages, and these values are gradually reduced to attain an equilibrium value once the organism reaches the adult phase (Caddy, 1996). Natural mortality coefficients can also decrease exponentially until an organism reaches the maximum age (Chen & Watanabe, 1989; Caddy, 1996; Jennings, Kaiser & Reynolds, 2009). The use of a constant M coefficient to evaluate a given population/stock, however, results from the complexity and uncertainty in the estimate of this parameter, as it assumes that all the exploited age groups have the same probability of death (Hewitt & Hoenig, 2005). Therefore, at present, it is necessary to make this assumption when size- or age-specific estimates of M for a particular resource are unavailable (Caddy, 1996; Andrews & Mangel, 2012).

Thus, for some sea cucumber species, there are estimates of M based mainly on the assumption that M remains constant during the life cycle of the species (Ebert, 1978; Herrero-Pérezrul et al., 1999; De la Fuente-Betancourt et al., 2001; Reyes-Bonilla & Herrero-Pérezrul, 2003). Other studies commonly estimate M using empirical equations based on a consideration of growth parameters (Pauly, 1980). These empirical equations use as input information such to the longevity of the organisms (Hoenig, 1983), age-at-maturity (Rikhter & Efanov, 1976), growth (only for the von Bertalanffy model), and mortality as a function of temperature (Beverton & Holt, 1959; Pauly, 1980).

The sea cucumber Isostichopus badionotus

Isostichopus badionotus (Selenka, 1867), also known as the four-sided sea cucumber, can grow to a maximum length of 45 cm (Toral-Granda, Lovatelli & Vasconcellos, 2008). There is at present very limited, information on the status of sea cucumber fisheries in the Caribbean region, only the studies of Buitrago & Boada (1996), De la Fuente-Betancourt et al. (2001), Alfonso et al. (2004) and López-Rocha, Cisneros-Reyes & Arreguín-Sánchez (2012), have addressed in some way the subject of their fisheries assessment. Yet intensive exploitation of Isostichopus species is carried out in several areas of the Caribbean Sea and the Gulf of Mexico. As a result, these populations are declining due to overexploitation and a lack of harvest regulations (Toral-Granda, Lovatelli & Vasconcellos, 2008). Only few countries (e.g., Cuba and Mexico) have regulated fishing activities by implementing catch quotas, a minimum capture size and closed fishing season (Alfonso et al., 2004; DOF, 2015). In some countries (e.g., Colombia), a total ban has been imposed on sea cucumber fishing (Lovatelli & Sarkis, 2011; Rodríguez-Forero & Agudelo-Martínez, 2017).

Different studies on biological and ecological aspects of Isostichopus badionotus populations have been conducted, these studies include aspects dealing with their geographic distribution (Guzmán & Guevara, 2002; Hasbún & Lawrence, 2002; Zetina-Moguel et al., 2003; Rodríguez et al., 2007; Tagliafico, Rangel & Rago, 2011; López-Rocha, Cisneros-Reyes & Arreguín-Sánchez, 2012; Hernández-Flores et al., 2015), reproduction, feeding and growth (Hammond, 1982; Sambrano, Díaz & Conde, 1990; Guzmán, Guevara & Hernández, 2003; Foglietta et al., 2004; Vergara & Rodríguez, 2015; Poot-Salazar, Hernández-Flores & Ardisson, 2015), and more recently immunology and aquaculture (Klanian, 2013; Zacarías-Soto et al., 2013; Zacarias-Soto & Olvera-Novoa, 2015; Martínez-Milián & Olvera-Novoa, 2016). There has been, however, no prior assessment of natural mortality coefficients for this species to date, and no such values have been reported in databases such as http://www.sealifebase.org (Palomares & Pauly, 2010).

Therefore, the aim of the present study was to estimate, for the first time, the natural mortality coefficient, M, using two indirect methods that take the age of the organism into consideration. These are a model based on gnomonic intervals (Caddy, 1996) and a model of the natural mortality coefficient at age developed by Chen & Watanabe (1989).

The gnomonic intervals model

The natural mortality coefficient (M) estimate obtained through application of the gnomonic interval model (Caddy, 1996) assumes that M decreases rapidly with age. Additionally, the life cycle of the species may be divided into intervals such that their duration increases as a function of age (Caddy, 1996). This method that subdivides the life cycle is referred to as “gnomonic”. A constant probability of death is assumed during each gnomonic interval. Estimation of the M parameter requires definition of the number of development stages considered during the life cycle of the species, the duration of the first stage, which corresponds to the first gnomonic time interval, the mean fecundity (annually) and the duration of each subsequent stage of the life cycle.

For implementation of this model, six different stages of development were considered during the life cycle of I. badionotus: four larval stages in addition to one juvenile stage and the adult stage. The first four stages comprised the early auricularia with a duration of two days, a medium auricularia stage lasting three days, and finally a late auricularia and doliolaria-pentactula stage, lasting six and 11 days, respectively (Zacarías-Soto et al., 2013) (Table 1). The last two larval stages were combined to allow convergence in the gnomonic calculation, since the model assumes a series of intervals of increasing duration (Caddy, 1996). The juvenile and adult stages were defined according to the age-at-maturity. It is important to know the effect of the age-at-maturity (t_m) in the estimation of mortality. Two different scenarios varying the age-at-maturity of 1.5 and five years (early and late age-at-maturity respectively) were compared (Table 1). Poot-Salazar (2014) estimated the length-at-maturity based on microscopic observations of the gonads of ∼1,000 individuals caught in the coasts of Yucatan, Mexico in which the highest concentration of mature organisms (25% and 75% percentiles of dorsal length) was found between 20 and 30 cm of dorsal length (equivalent to 18.5 and 32.2 cm Le respectively) corresponding to ∼1.5 and ∼5 years of age (Poot-Salazar, Hernández-Flores & Ardisson, 2015). For this reason, we decided to use both scenarios to avoid discarding the variability in the size and age-at-maturity of this species. Therefore, the age-at-maturity of 1.5 years was defined as early maturity (scenario 1) and age-at-maturity of five years as late maturity (scenario 2) (Table 1).

The age-at-maturity was calculated from the length-at-maturity (23 cm dorsal length equivalent to estimate length (Le) = 21.9 cm) reported for this species by (Poot-Salazar, Hernández-Flores & Ardisson, 2015). The inverse of the von Bertalanffy growth model was used to calculate the age at which the species reaches 23 cm in dorsal length (Le = 21.9 cm). Two sets of parameters of the seasonal von Bertalanffy growth model reported by (Poot-Salazar, Hernández-Flores & Ardisson, 2015) were used. The growth parameters were estimated for 2,347 sea cucumbers collected by previous authors in three zones in Yucatan coast through the analysis of modal progression. Electronic Length Frequency Analysis ELEFAN-1 and Projection Matrix PROJMAT were methods used to estimate growth parameters (Poot-Salazar, Hernández-Flores & Ardisson, 2015).

The duration of the first larval stage (t ₁) corresponded to the first interval of gnomonic time (Δ₁), then Δ₁ = t₁. That of the second gnomonic interval was obtained as Δ₂ = α × t₂₋₁, where α is a proportionality constant (Caddy, 1996). This constant minimizes the residuals until the sum of the duration of all the gnomonic intervals equals the longevity of the species (t_n) using the Newton algorithm (Neter et al., 1996). The successive gnomonic time intervals were calculated as Δ_i = (α × t_i−1) + t_i−1 where i ≥ 3 to the n –th gnomonic interval, and n is the number of intervals. The natural mortality coefficient in each gnomonic interval, M_i, was estimated from the following expression: (1)${\displaystyle M_i=\frac{G}{{\theta }_i-{\theta }_{i-1}}}$ where G is the constant death probability for each gnomonic interval determined by numerical methods using Newton algorithm (Neter et al., 1996). The value of G is the product of M_i by Δ_i and G parameter which allows to calibrate the duration time of the gnomonic intervals. To date, there is no published evidence to determine if G parameter can be considered as an attribute of the species studied. The number of surviving individuals in the last interval is N_i = 2 (Caddy, 1996), θ is the difference of the annual proportional duration of each interval, calculated as θ_i = (Δ_i∕t_n)∕365, and where t_n is the total longevity of the species in days (Martínez-Aguilar et al., 2010). Because the actual longevity of this species has not been proven, we tested simulating M using three longevity (10, 20 and 30 years) in the present study and found that after 10 years the values of M did not change significantly (P > 0.05). Therefore, the results correspond to a longevity of 3,650 days (∼10 years); this longevity is similar to that supported by a maximum predicted age for this species (based on length-age and weight-relationships) from 10 to 12 years (Poot-Salazar, Hernández-Flores & Ardisson, 2015).

The number of survivors (N_i) at the beginning of each interval (i –th) is equal to the number of survivors from the previous interval, except for the first interval for which it is assumed that the total number of individuals at the beginning of the interval equals the individual average fecundity per year of the species ( $\bar{f}$), determined from the following equations: (2)${\displaystyle N_i=\bar{f}{e}^{-M_i{\theta }_i}\mathrm{for}i=1}$ (3)${\displaystyle N_i=N_{i-1}{e}^{-M_i{\theta }_i}\mathrm{for}i>1}$

Variability in M was evaluated considering that the main source of uncertainty in the estimates is contributed by f. A total of 1,000 repeated samples of f were simulated to determine the uncertainty of the M_i estimates; then M_i estimates were obtained from 1,000 simulations per i gnomonic interval. A uniform, non-informative distribution (U) of f was assumed as follows: (4)${\displaystyle f\approx \mathrm{U}\left(f_{min}.f_{max}\right)}$ where f_min and f_max are the minimum and maximum fecundity values reported, respectively (Martínez-Aguilar, Arreguín-Sánchez & Morales-Bojórquez, 2005).

For I. badionotus, the f estimate is uncertain and to date no specific estimate of the annual fecundity is available. Direct counts of the number of eggs released per spawning in captivity were used for this calculation, which yielded an average of 200,000, and a range from 13,500 to 1,066,000 (minimum and maximum) eggs per spawn per female (n = 83 spawnings) (Zacarías-Soto et al., 2013). Additionally, the estimated fecundity using a volumetric method from gonadal microscopic analysis yielded a range between 62,496 and 5,062,490 ova with a mean of 1,417.330 ±151,186 (standard deviation, SD) ova (n = 175 specimens) (Palazon, 2001). This was based on the assumption that the number of hatched larvae is equivalent to the individual fecundity of the species (Caddy, 1996), and includes all estimates until capture, thus representing the total lifespan variability.

Estimates of M were obtained with the Gnomonic-Interval Natural-Mortality Method (GIM) routine for Excel (Morales-Bojórquez et al., 2003), which allows the incorporation of observed data during the specific duration of each developmental stage of the species to calibrate estimates of the M vector (Martínez-Aguilar et al., 2010).

The Chen & Watanabe model

This model is derived from the relationship described by Chen & Watanabe (1989), in which there is a qualitative correspondence between the growth stage and M. Therefore, M can be calculated at any age (t). It is necessary to know the body growth coefficient (k) and the theoretical age when the length is equal to zero (t ₀) of the von Bertalanffy growth model (Von Bertalanffy, 1938). For this model, M is assumed to be inversely proportional to growth (k), using age-at-maturity (t_m) explicitly. The M_t estimate was derived using the following equations: (5)${\displaystyle M_t=\frac{k}{\left(1-e^{-k\left(t-t_0\right)}\right)}\mathrm{for}t\le t_m}$ (6)${\displaystyle M_t=\frac{k}{{\gamma }_0}+{\gamma }_1\left(t-t_m\right)+{\gamma }_2{\left(t-t_m\right)}^2\mathrm{for}t>t_m}$ where M_t is the average natural mortality at age t, γ₀, γ₁ and γ₂ are parameters of the model calculated based on the continuity of growth to the age t_m as follows: (7)${\displaystyle {\gamma }_0=1-e^{-k\left(t_m-t_0\right)}}$ (8)${\displaystyle {\gamma }_1=k{e}^{-k\left(t_m-t_0\right)}}$ (9)${\displaystyle {\gamma }_2=-0.5k^{2e^{\left(-k\left(t_m-t_0\right)\right)}}}$

As in the GIM method above, two different scenarios of age-at-maturity ( t_m) were considered (early and late age-at-maturity).

The number of survivors to the age N_t were estimated from: (10)${\displaystyle N_t=N_0\frac{1-e^{k{t}_0}}{1-e^{-k\left(t-t_0\right)}}{e}^{-kt}\mathrm{for}t\le t_m}$ (11)${\displaystyle N_t=N\left(t_m\right)\left[\frac{2{\gamma }_2\left(t-t_m\right)\left({\gamma }_1-\sqrt{{\gamma }_1^2-4{\gamma }_0{\gamma }_2}\right)}{2{\gamma }_1\left(t-t_m\right)\left({\gamma }_1-\sqrt{{\gamma }_1^2-4{\gamma }_0{\gamma }_1}\right)}\right]\frac{k}{\sqrt{{\gamma }_1^2-4{\gamma }_0{\gamma }_1}}\mathrm{for}t>t_m}$ where: (12)${\displaystyle N\left(t_m\right)=N_0\frac{1-e^{k{t}_0}}{1-e^{-k\left(t_m-t_0\right)}}{e}^{-k{t}_m}}$ and N ₀ is the estimated survival for time t − 1 or the fecundity of the species when t = 1 (Chen & Watanabe, 1989).

It is important to characterize the variability associated with these estimates, because the use of a single, discrete estimate without error measures underestimates the uncertainty associated with the results for future evaluations (Brodziak et al., 2011). In this study, error measures were added to the (Chen & Watanabe, 1989) model, assuming that the main source of uncertainty is generated by the parameters k and t ₀.

The present study was based on the minimum and maximum values of the parameters k and t ₀ of the von Bertalanffy model reported for this species in different locations of the Gulf of Mexico (Poot-Salazar, Hernández-Flores & Ardisson, 2015) (Table 1). The k parameter varied from 0.2 to 0.7 year⁻¹ and the t ₀ parameter ranged between -0.8 and -0.4 years. A total of 1,000 repeated samplings were conducted to generate the k and t ₀ parameters, assuming a uniform, non-informative distribution (U). Finally, M_t values were estimated from the following equations: (13)${\displaystyle k\approx \mathrm{U}\left(k_{min}.k_{max}\right)}$ (14)${\displaystyle t_0\approx \mathrm{U}\left[{\left(t_0\right)}_{min}.{\left(t_0\right)}_{min}\right]}$

Statistical comparison of Natural mortality between scenarios

Generalized Linear Models (GLM) were applied to know the statistical effect of age-at-maturity scenarios on M for GIM and the Chen & Watanabe (1989) methods. M was the response variable while the explanatory variables were the maturity scenarios, the gnomonic stages and the ages. An interaction between gnomonic stages (for GIM), age (for Chen y Watanabe method) and stages of maturity was carried out assuming an identity link function and Gaussian error distribution (η = μ) in order to know if there is a difference between and stages with respect to the maturity scenarios (Nelder & Baker, 1972).

Natural mortality coefficients estimated by the gnomonic intervals model

The estimates of the M coefficient derived using the gnomonic intervals model included the highest values of natural mortality during the larval phase covered by the first four gnomonic intervals (Table 2). A considerable reduction in M values was observed in the juvenile and adult phases (Table 2). The M estimates varied depending on the different scenarios of age-at-maturity considered, such that magnitude of the coefficients in the adult stage was higher in late maturity scenario, whereas in the juvenile stage the magnitude of M decreased (inverse effect) (Table 2). The daily survival rate from one gnomonic interval to another interval is shown in Table 2.

Natural Mortality Coefficients estimated by the Chen & Watanabe (1989) Model

The M estimates obtained using the Chen & Watanabe (1989) model suggested that there was considerable variation among ages. It was generally observed that in the late maturity scenario the (Chen & Watanabe, 1989) model tends to estimate lower natural mortality values (Table 3). The variability of the estimates obtained for the older ages (>5 years) decreased in the late maturity scenario (Table 3).

Statistical comparison of Natural mortality between scenarios

The differences in M among stages were highly significant for both the gnomonic interval method and the Chen & Watanabe method. However, for the GIM method, the maturity scenario was not significant, nor was the interaction between scenarios and stages. This suggests that the age-at-maturity applied by both scenarios (1.5 and five years) does not have a significant effect on natural mortality. For the method of Chen & Watanabe it was found that, both maturity scenarios, the age and the interaction between both variables showed an incidence on the values of natural mortality where the mortality values were lower when the age and the age at maturity are more advanced (Table 4).