Management Implications for Skates and Rays Based on Analysis of Life History Parameters

Introduction

Many batoids, similar to sharks, have the life history characteristics of slow growth, late maturity, and low numbers of offspring (Ebert and Sulikowski, 2009). Excluding manta rays and butterfly rays, most skates and rays inhabit coastal demersal waters and play an important role in the demersal ecosystem (Ebert and Bizzarro, 2007). These skates and rays are vulnerable to anthropogenic pressure and may decline or collapse after experiencing heavy fishing pressure (Hoff and Musick, 1990; Musick, 1999).

The global landings of skates and rays reported to the United Nations Food and Agriculture Organization (FAO) declined almost 20% from 2003 to 2012 (Davidson et al., 2016). Recent assessments of global oceanic rays suggested that several species have been overexploited or have even collapsed (Pacoureau et al., 2021). With the increase in skate and ray catches, it is necessary to develop management plans of these species to ensure sustainable use of these resources (Dulvy and Reynolds, 2002). However, the development of fishery management plans for skates and rays is difficult due to the lack of detailed biological information and species-specific catch data (Stevens et al., 2000) because most of them are of low economic value. The decline of large skates coupling with an increase in the abundance of small skates in the northeast Atlantic resulted in a structural change in the marine ecosystem (Dulvy et al., 2000). These results clearly indicated that species-specific stock status information is urgently needed to ensure the sustainable utilization of skate and ray stocks.

According to the International Union for Conservation of Nature (IUCN) red list criteria, 10% of chondrichthyes fall in the category threatened (including critical endangered, endangered, and vulnerable), of which 40% are sharks and 60% are skates and rays. The Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES) has placed the sawfishes on its Appendix I list (Anon., 2013) and the manta rays and the devil rays Mobula spp. (Anon., 2016) and the guitarfish Rhinobatidae and wedgefish Rhinidae (Anon., 2019) on its Appendix II list. The aforementioned measures highlighted the urgency of conservation and management of skates and rays.

Size and life history parameters of skates and rays vary in wide range. The maximum size ranges from 10-cm total length (TL) for the short-nose electric ray Narcine laticaudus (Froese and Pauly, 2015) to 910 cm disc width (DW) for the manta ray Mobula birostris (White et al., 2006). Growth coefficient (k) ranges from 0.04 year^–1 for the big skate Beringraja binoculata (McFarlane and King, 2006) to 0.454 year^–1 for the starry skate Raja asterias (Serena et al., 2005). The litter size varies remarkably among species even for those falling within the same reproductive mode (oviparity, viviparity, or aplacental viviparity). For example, the fecundity of oviparous stocks ranged from 8 for the roughtail skate Bathyraja trachura (Ebert, 2005) to 360 for the big skate (Ebert and Davis, 2007).

The catch, effort, and bycatch data for most skates and rays are lacking because they are usually discarded or treated as trash fish that have low commercial value. Consequently, conventional models commonly used in stock assessment for teleost fish, such as surplus production and stock-recruitment models, have seldom been used in examining the stock status of skates and rays despite multispecies skate and ray fisheries including three species of Bathyraja spp. and Raja flavirostris in Falkland Islands (Agnew et al., 2000). Instead, yield per recruit and spawning stock biomass per recruit models have been applied to the little skate Leucoraja erinacea (Waring, 1984) and the longnose skate Raja rhina (Gertseva, 2009), respectively. In addition, demographic models, which have been applied to sharks, have been successfully applied to describe the population dynamics of the little skate, the winter skate Leucoraja ocellata, the barndoor skate Dipturus laevis (Frisk et al., 2002), and five deep-water Bering Sea skates Bathyraja spp. (Barnett et al., 2013). However, detailed information on life history parameters including natural mortality, age at maturity, litter size, reproductive cycle, and longevity is needed in this approach. Therefore, these models are difficult to be applied to the species with limited life history information (King and McFarlane, 2003). Hence, it is urgent to develop an alternate way to estimate the finite rates of population increase for skates and rays based on life history parameters to ensure the sustainability of skate and ray stocks.

Multivariate analyses including principal component analysis (PCA), cluster analysis (CA), and regression analysis have been used in marine fish management by various authors (Winemiller and Rose, 1992; Jennings et al., 1999; Cortés, 2000; Frisk et al., 2001; King and McFarlane, 2003). Cortés (2000) identified three life history strategies of sharks based on five life history parameters of 34 species (40 stocks) by using PCA and CA. Liu et al. (2015) also categorized three groups using similar methods and developed empirical equation for each group based on six life history parameters of 38 species (62 stocks) of sharks to estimate their finite rate of population increase. However, none of aforementioned studies has provided an empirical equation to estimate the finite rate of population increase for skates and rays.

Hence, the present study aims to categorize the life history strategies of skates and rays based on their life history parameters using multivariate analysis, to develop empirical equations to estimate the finite rate of population increase, and to propose appropriate management measures for each group. For those skates and rays without detailed life history parameters, they can be classified into one group based on similar species identified in this study, and then provisional management actions can be taken in a precautionary way.

Materials and Methods

Intensive search of the existing literature including scientific journals, reports, and gray literature using keywords of skates, rays, age and growth, and reproduction was used to collect life history parameter data for skates and rays. Only stocks with both age and growth and reproduction information available were included in the analysis.

As conventional demographic analysis assumes that males are not the limiting factors regulating population growth, this study used data only from females. Where sex-specific parameters were not available, sex-combined parameters were used. In total, 12 life history parameters were selected. These included five age and growth parameters [i.e., asymptotic length (L_∞), growth coefficient (k), age at zero length (t₀), maximum age (T_max), and maximum observed length (L_max)] and seven reproduction parameters [i.e., age at maturity (T_m), reproductive strategy (R), size at maturity (L_m), size at birth (L_b), fecundity/litter size (f), gestation period (G_p), and reproductive cycle (R_c)]. Different studies define life history parameters in slightly different ways. To account for this inconsistency, we used the following definitions.

Reproduction Parameters

1. Reproductive strategy (R): oviparity (O), viviparity (V), and aplacental viviparity (OV).

2. Fecundity/litter size (f): the mean litter size of pregnant females, the mean of the maximum and minimum of litter sizes, or the mature eggs assuming all were fertilized.

3. Size at maturity (L_m or DW_m): size (TL or DW) at first maturity, size at 50% maturity, mean size of mature specimens, or the estimated size at maturity by substituting the age at maturity into the growth equation. Several studies used different terminologies; however, based on their description, they were referring to size at 50% maturity.

4. Reproduction cycle (R_c): including gestation and resting periods, if only gestation information was available, R_c was estimated using the value from similar species.

5. Size at birth (L_b or DW_b): size of the smallest free swimmer, the mean of the size range at birth, or the mean of the largest full term embryo and the smallest free swimmer.

Age and Growth Parameters

(1) Maximum observed size (L_max): the maximum size of observed skates and rays.

(2) Maximum age (longevity) (T_max): the maximum age was estimated from Ricker’s (1979) equation, which was 95% of L_∞ as follow: T_max = 2.77/k.

(3) Asymptotic length (L_∞ or DW_∞) and growth coefficient (k) from growth equations.

(4) Age at maturity (T_m): the age at 50% maturity, or the mean age of mature specimens or the mean of the maximum and minimum age at maturity if only the range of age at maturity was available.

Input Parameters

Large variations in L_b, L_m, and L_∞ were found between different species (Supplementary Table 1), and this may affect the results of the analysis. To eliminate the size effect in the analysis, six life history input parameters, namely, the ratio between size at birth and asymptotic length (L_b/L_∞), the ratio between size at maturity and asymptotic length (L_m/L_∞), T_max, T_m, k, and annual fecundity (f/R_c), were used in the analysis.

Habitat Information

To examine the life history strategy associated with the habitat of skates and rays, habitat information including the water depth of habitat (Dep), sea surface temperature (SST), and salinity (Sal) were collected from the literature for the blue stingray Dasyatis chrysonota in the South Africa (Cowley, 1997; Ebert and Cowley, 2008), the diamond stingray Hypanus dipterura in the western Mexico (Mariano-Melendéz, 1997; Smith, 2005; Smith et al., 2007), the western shovelnose stingaree in the western Australia (White et al., 2002), the masked stingaree Trygonoptera personata in the western Australia (White et al., 2002), lobed stingaree Urobatis lobatus in the lower western coast of Australia (White et al., 2001), yellownose skate Dipturus chilensis (Fuentealba and Leible, 1990; Licandeo et al., 2006), and sharpspine skate in the northeastern waters off Taiwan (Chang, 2006; Joung et al., 2011). The habitat information for other species was adopted from the Fish Base (Froese and Pauly, 2015). The habitat information was treated as category variable, as follows: SST, (1) 0–10°C, (2) 10–20°C, and (3) > 20°C; depth, (1) < 50, (2) 50–200, and (3) > 200 m; salinity, (1) 30–33, (2) 33–36, and (3) > 36 psu (Supplementary Table 1).

Demographic Analysis

Conventional demographic analysis requires an input of natural mortality (M) to calculate demographic parameters. Thus, Hoenig (1983) equation was used to estimate the mean M for each stock depending on the longevity, as follows: ln(M) = ln(Z) = 1.46 − 1.01 × ln(t_max) or M = −ln(0.01)/t_max (Cortés, 1998), where Z is total mortality. Natural mortality approaches Z when the fish stock is unfished or at light exploitation levels. Krebs (1985) formula was used to calculate demographic parameters, assuming a sex ratio of 1:1 for the embryos and the population was in equilibrium condition $\left(\sum \frac{1}{2}{m}_x{l}_x{e}^{rx}=1\right)$, and the parameters were calculated as follows: $R_0={\sum }_{x=0}^{t_{max}}\frac{1}{2}{m}_x{l}_x$, $\text{G}=\frac{{\sum }_{x=0}^{t_{max}}\frac{1}{2}x{l}_x{m}_x}{R_0}$, $r=\frac{ln\left(R_0\right)}{G}$, λ=e^r, t_x2 = ln(2)/r.

where R₀ is the net reproductive value per generation, x is the age, t_max is the maximum age, m_x is the fecundity of age x, l_x is the survival rate until age x, G is the generation time, r is the intrinsic rate of population growth, λ is the finite rate of population increase, and t_x2 is the theoretical population doubling or halving time.

Multivariate Analysis

Due to inconsistencies in measurement units, our PCA used correlation matrices, R, rather than variance–covariance matrices. All parameters were log-transformed and normalized, and the eigenvectors and eigenvalues were estimated. A nonparametric multiple dimensional scaling (NMDS) was used to draw the bi-plot. Life history parameters were reduced to several independent principal components, and the scores of principal components were then analyzed using the CA.

The CA with hierarchical Ward’s method was used to estimate the scores of the first to third principal components and to draw the tree plot. Species with similar parameter values were grouped together and named according to their shared life history traits. After grouping, the general linear model (GLM) was used to develop an empirical equation for each group to describe the relationship between λ′ and life history parameters. A variance inflation factor (VIF) (Wooldridge, 2009) was used to examine the multicollinearity of life history parameters in regression models. Multicollinearity exists among life history parameters when VIF ≥ 10, and the parameter is removed from the regression model. The best model was selected using the stepwise Akaike information criterion (AIC; Akaike, 1974) method (Venables and Ripley, 2002).

Correlation of Life History Strategy and Habitat

After life history strategy group was identified by using CA, the scores of the first two principal components were plotted, and the life history strategy associated with the habitat information (Dep, SST, and Sal) for each group was examined.

Results

In total, life history parameters of 35 species of skates and rays (41 stocks), comprising four orders (Myliobatiformes, Rajiformes, Rhinopristiformes, and Torpediniformes) and seven families, were collected (Supplementary Tables 1,2).

Life History Parameters

Age and Growth

The L_max of 35 species (41 stocks) ranged from the minimum of 53.9 and 54.0-cm TL for the sharpspine skate and little skate (Waring, 1984; Joung et al., 2011) to the maximum of 235.0-cm TL for the blue skate Dipturus batis (Du buit, 1977) with a median of 103.4-cm TL. For age and growth parameters, the minimum L_∞ was 52.7-cm TL for little skate (Waring, 1984), and the maximum L_∞ was 293.5-cm TL for the big skate (McFarlane and King, 2006) with a median of 125.8-cm TL. The k value ranged from the minimum of 0.040 year^–1 for the big skate (McFarlane and King, 2006) to the maximum of 0.454 year^–1 for the starry skate (Serena et al., 2005), with a median of 0.112 year^–1 (Supplementary Table 1).

Reproductive Parameters

Among the stocks collected for this study, 23 stocks are oviparity, four stocks are aplacental viviparity, and 14 stocks are viviparity. The fecundity of oviparous stocks ranged from 8 for the roughtail skate and the sandpaper skate Bathyraja interrupta (Ebert, 2005) to 360 for the big skate (Ebert and Davis, 2007) with a median of 40. For viviparous stocks, the fecundity ranged from 1 for the cownose ray Rhinoptera bonasus (Smith and Merriner, 1986), the western shovelnose stingaree, and the masked stingaree (White et al., 2002) to 5 for pelagic stingray (Hemida et al., 2003) with a median of 3. The fecundity of aplacental viviparous stocks ranged from 5 for the southern fiddle ray Trygonorrhina fasciata (Marshall et al., 2007) to 17 for the Pacific electric ray (Neer and Cailliet, 2001). Among those 18 stocks where their gestation periods have been documented, 15 stocks have 1-year gestation period, the Pacific electric ray has 2.5-year gestation period (Neer and Cailliet, 2001), and the diamond stingray and round stingray have 0.5-year gestation period (Babel, 1967; Hemida et al., 2003) (Supplementary Table 1). For those species without gestation period information, gestation period or R_c was estimated using the value from similar species or assumed to be 1 year.

The Ratios of L_b/L_∞ and L_m/L_∞

The L_b/L_∞ ratios of the 35 species (41 stocks) ranged from 0.083 to 0.483, with a median value of 0.165. The maximum value of L_b/L_∞ was 0.483 for the white-spotted stingaree (White and Potter, 2005), and the minimum value was 0.083 for big skate (McFarlane and King, 2006) with a median of 0.165. In total, 26 stocks (63.41%) ranged from 0.102 to 0.200, five stocks (12.20%) ranged from 0.083 to 0.093, and 10 stocks ranged from 0.211 to 0.483 (24.39%) (Supplementary Table 2).

The L_m/L_∞ ratios ranged from 0.307 to 0.855, with a median value of 0.641. The maximum value of L_m/L_∞ was 0.855 for the white-spotted stingaree in southwestern Australia (White and Potter, 2005), and the minimum value was 0.307 for the big skate (McFarlane and King, 2006). In total, 35 stocks (85.37%) ranged from 0.505 to 0.855, and six stocks (14.93%) ranged from 0.307 to 0.50 (Supplementary Table 2).

Maximum Age and Natural Mortality

The maximum ages t_max estimated from Ricker’s (1979) equation ranged from 7 years for the starry skate R. asterias (Serena et al., 2005) to 74 years for the big skate (McFarlane and King, 2006) with a median of 25 years (Table 1).

Finite Rate of Population Increase

The finite rate of population increase estimated from conventional demographic analysis ranged from 0.936 for the white-spotted skate to 1.606 for the thornback ray Raja clavata. Thirty-two of 41 stocks (78.1%) fell in the range 1.002–1.398, six stocks (14.6%) had values in the range 1.417–1.606, and three stocks (7.3%) had values lower than 1.0 (Table 1).

Multivariate Analyses

Correlation Between Life History Parameters

Significant negative relationships were found for the t_max and k (r = −0.973, p < 0.01, n = 41), suggesting that the larger the t_max, the slower the growth. Similar relationships were also found for annual fecundity (litter size) and L_b/L_∞ (r = −0.818, p < 0.01) as well as t_mat and k (r = −0.591, p < 0.01) (Table 2), suggesting that the more the annual fecundity, the less the value of L_b/L_∞ and the earlier the fast growth species mature. Significant positive correlation was also found between t_max and t_mat, suggesting that the species with longer longevity mature late.

Principal Component Analysis

Results of PCA revealed that the top three principal components can explain 96.2% of the variations, of which 52.6, 23.7, and 19.9% of the variations can be explained the first, second, and third components (PC1, PC2, and PC3), respectively. The loadings of PC are the correlations of life history parameters and PC. The life history parameters in PC1 includes L_b/L_∞ (r = 0.732), L_m/L_∞(r = 0.589), k (r = 0.895), T_max (r = −0.878), and T_m (r = −0.635); PC2 includes L_b/L_∞ (r = 0.619), and f/R_c (r = −0.789); and PC3 includes L_m/L_∞ (r = 0.796) and T_m (r = 0.738).

The scatter plot of scores of PC1 and PC2 showed that the positive scores to PC1 represent species with large k and small T_max, such as the starry skate and little skate; the negative scores to PC1 represent species with slow growth and large T_max, such as the big skate and blue skate. For PC2, the positive scores represent species with large L_b, and small litter size (fecundity) such as western shovelnose stingray and the masked stingaree; the negative scores represent species with small L_b, and large litter size (fecundity) such as thornback ray and thorny skate. The scatter plot of scores of PC1 and PC3 showed that the positive scores to PC3 represent late maturing species, such as the barndoor skate and white-spotted skate; the negative scores to PC3 represent early maturing species, such as the common stingray Dasyatis pastinaca and longheaded eagle ray Aetobatus flagellum.

Cluster Analysis and Empirical Equations

Three groups were identified by using the CA with hierarchical Ward’s method (Figure 1), and the box plots of six life history parameters by group are shown in Figure 2. The three groups categorized by the CA were clearly separated by PC1 and PC2 (Figure 3). As there was a high negative correlation between T_max and k, and VIF indicated a multicollinearity for these two parameters, T_max was excluded in empirical equation development.

Group 1

Fourteen stocks with slow growth rates (k < 0.11 year^–1), early maturity (L_m/L_∞ < 0.62), and extended longevity (T_max > 25 years) fell into this group, e.g., the blue stingray D. chrysonota and the big skate. The maximum age ranged from 25 years for the blue stingray to 69 years for the big skate in British Columbia waters, with all stocks being in the range 25–50 years except the big skate. The empirical equation developed for estimating the finite rate of population increase was as follows: λ′ = 1.7127 + 0.3127×ln(L_m/L)−0.272×ln(T_m) + 0.108×ln(f/R_c)(n = 14,p < 0.05) (Table 3A).

Group 2

Fifteen stocks with intermediate growth rates (0.080 year^–1 < k < 0.190 year^–1), intermediate longevity (17 years < T_max < 35 years), and late maturity (L_m/L_∞ > 0.60) fell in this group, e.g., the barndoor skate and the sharpspine skate. The T_m, k, and f/R_c are the significant parameters. The value of T_m ranged from 3.7 years for the common stingray to 14.5 years for the yellownose skate. The value of k for 12 of 15 stocks (80%) fell in the range 0.10–0.19 year^–1, with the remaining three stocks (20%) ranging from 0.08 to 0.09 year^–1. The largest k value was for the blonde skate (k = 0.19 year^–1), while the smallest was for the roughtail skate (k = 0.08 year^–1). The f/R_c ranged from 5.3 for the southern fiddler ray T. fasciata to 140 for the thornback ray R. clavata (Figure 2). The empirical equation for estimating the finite rate of population increase is as follows: λ′ = 1.25230.4038×ln(T_m)0.2302×ln(k) + 0.1166×ln(f/R_c)(n = 15,p < 0.05) (Table 3B).

Group 3

Twelve stocks with fast growth rate (k > 0.160 year^–1), low longevity (T_max < 23 years), and large size at birth (L_b/L_∞ > 0.18), e.g., the pelagic stingray and the little skate, fell in this group. L_m/L_∞ ranged from 0.39 for the pelagic stingray to 0.86 for the white-spotted stingaree Urolophus paucimaculatus, with nine of 12 stocks (75%) in the range 0.63–0.86. A second characteristic of this group was larger values of f/R_c and L_b/L_∞ (Figure 3). The empirical equation for estimating the finite rate of population increase is as follows: λ′ = 0.9205 + 0.0603×ln(L_b/L)0.3744×ln(L_m/L) + 0.1234×ln(f/R_c) (n = 12, p < 0.05) (Table 3C).

The λ′ predicted from the aforementioned empirical equations ranged from 0.9873 to 1.5937 with medians of 1.1084–1.3057. The correlations between λ and λ′ are 0.98, 0.99, and 0.99 for the three groups, and the regression lines are close to the 45° lines (Figures 4A–C). No obvious trends were found for the residual plots of the three regression lines, suggesting that the empirical equations derived from this study can precisely predict the finite rate of population increase for skates and rays.

Habitat Parameters and Grouping

The relationships between habitat parameters and groups related to PC1 and PC2 were as follows. The SST of the habitats of Group 1 and Group 2 species was lower than that of Group 3 species. The species that fall in Groups 1 and 2 are temperate raja in high latitude, but the species that fall in Group 3 are stingaree in tropic waters. The depth of habitat for Group 1 species ranged from 25 to 150 m; that for Group 2 species was deeper than 150 m; and that for Group 3 species was less than 50 m. The grouping had no significant correlation with salinity.

Discussion

The empirical equations based on life history parameters of 41 skate and ray stocks were developed to estimate population increase rates in this study. This approach can be applied to derive useful information in conservation and management for skates and rays particularly for those species with limited data. However, as the life history data used in this study were adopted from the literature, the inconsistence of data quality and design differences among studies may occur. To solve this problem, future work should consider using meta-analyses employing hierarchical models demonstrated by Thorson et al. (2015).

Although the length measurements in the literature for skates and rays were inconsistent, the DWs were not converted to TL because the DW–TL relations were only available for few species. As the values of L_b/L_∞, L_m/L_∞ and DW_b/DW_∞, DW_m/DW_∞ are highly correlated, and only these ratios were used in empirical equations, inconsistencies in length measurement among stocks are not likely to influence the estimate of λ′.

Skates and rays grow slowly compared with teleost fishes; however, the growth coefficients vary largely among species (Musick, 1999; Cailliet and Goldman, 2004). Cailliet and Goldman (2004) mentioned that k ranged at 0.20-0.50 year^–1 for Rhinobatidae and Torpedinidae and Myliobatiformes and that k ranged at 0.05–0.50 year^–1 for Rajiformes. The k values of 41 stocks used in this study ranged from 0.04 to 0.454 year^–1, which fall in the aforementioned ranges. Given that the data used in this study covered a wide range of growth rates, we believe the results derived from this study are robust and can be applied to other skate and ray species.

Winemiller and Rose (1992) and King and McFarlane (2003) identified five groups of life history strategy for fish: opportunistic, periodic, equilibrium, median, and migratory species. Elasmobranchs fall in equilibrium species due to their characteristics of slow growth, late maturity, and prolonged longevity. Liu et al. (2015) identified three life history strategy groups for sharks based on their vital parameters. In this study, three groups were also identified for skates and rays. The species in Group 1 have the characteristics of slow growth and prolonged longevity, which are similar to the characteristics of equilibrium species for sharks, but they mature early. The species in Group 2 have the characteristics between those of opportunistic and equilibrium species for sharks but mature late, which are similar to those of periodical species. The species in Group 3 have the characteristics of fast-growing and short longevity, which are similar to opportunistic species for sharks but with larger L_b/L_∞ values.

The λ values of seven species of skates and rays have been estimated using conventional demographic approach (Simpfendorfer, 2000; Neer and Cailliet, 2001; Frisk et al., 2002; Mollet and Cailliet, 2002; Quiroz and Wiff, 2005; Smith et al., 2008). Slight differences between λ from previous studies and λ′ derived from our empirical equations were found for yellownose skate, diamond stingray, winter skate, little skate, and Pacific electric ray, but the difference was larger for pelagic stingray (Table 4). The differences between λ or λ′ and those values of previous studies (Table 4) were due to the different settings of life history parameters such as t_max, t_mat, M, or f/R_c. For example, the different settings of M, t_max, and f/R_c for pelagic stingray lead the different estimated values (1.174 and 1.433) for λ and λ′, respectively. Nevertheless, the values of λ and λ′ estimated in this study were close even using different approaches, suggesting that the empirical equations derived from this study are robust estimators of λ.

Annual fecundity is important information for demographic analysis. Large variation in fecundity was documented between oviparous and viviparous elasmobranchs. For example, the mean annual fecundity was only 5.5 pups per litter for skates but was 58.9 egg capsules for rays (Musick and Ellis, 2005). Compared with the viviparous species, the oviparous species are more vulnerable to predations or environmental changes; thus, they take a strategy of producing more egg capsules (Cox and Koob, 1993; Lucifora and Garcia, 2004). Among the species analyzed in this study, Myliobatiformes (viviparous species) have the lowest litter size, while Rajiformes (oviparous species) produce the highest number of eggs. Although most skates have high fecundity, those of Bathyraja spp. and Okamejei sp. have relatively low fecundity. The low fecundity is likely due to the relative stable habitats in deeper waters for these species.

Annual fecundity information is very limited particularly for oviparous species. The reproductive cycle was assumed to be 1 year for those species without this information. The big skate can produce 360 egg capsules (three to four embryos per capsule) and 1,440 embryos per year at most, which is the highest productivity species of Elasmobranchs (Ebert and Davis, 2007). The results of CA did not change either embryo per capsule, which was assumed to be one or four for big skate in this analysis. Pardo et al. (2018) examined the effect of uncertainty on the estimating maximum intrinsic rate of population increase and recommended that distributions of litter sizes should be frequently measured and natural mortality estimation should be improved for data–poor elasmobranchs.

The ratio of L_b/L_∞ of most elasmobranchs ranged from 0.150 to 0.350 (Joung, 1993). The ratios of the skates and rays used in this study fall in the range of 0.102–0.200 (n = 26, 63.4%), 0.150–0.350 (n = 20, 48.8%), and <0.150 (n = 17, 41.5%). Only four species (9.76%), i.e., western shovelnose stingaree, masked singaree, lobed stingaree, and white–spotted stingaree, have the L_b/L_∞ > 0.350. The low L_b/L_∞ values for skates and rays may be related to their large fecundity/litter size. Cortés (2000) and Liu et al. (2015) also documented that there is a negative correlation between L_b/L_∞ and annual fecundity. As the annual fecundity of skates and rays was higher than that of sharks, smaller values of L_b/L_∞ were expected.

The L_m/L_max value falls in 0.5–0.9 for most elasmobranchs (Holden, 1974; Compagno, 1984; Pratt and Casey, 1990; Joung, 1993). Similar results were found in this study. Although L_m/L_∞ ranged from 0.307 to 0.855, 35 out of 41 stocks fall in 0.5–0.9. Only six stocks have the L_m/L_∞ smaller than 0.5, and five of them fall in Group 1 (early maturity).

The bias of t_max estimation may influence the estimation of λ. To avoid the underestimates of the maximum age by using the maximum observed age (Skomal and Natanson, 2003) due to lacking of oldest samples, empirical equations such as those from Taylor (1958); Fabens (1965), and Ricker (1979) were commonly applied to estimate the maximum age (Cailliet et al., 2006). In this study, Ricker’s (1979) equation was selected because the values estimated from the other two equations were much higher than the maximum observed age. Consequently, this method was used to estimate t_max for all species. Frisk et al. (2001) developed an empirical equation between t_max and age at maturity for elasmobranchs based on 35 species. However, the t_max values back estimated from this method were smaller than those from other methods for most species, suggesting that the approach may underestimate the maximum age of skates and rays. Thus, this method was not used in this study. Nevertheless, as t_max was not selected in our empirical equations in estimation of λ′, the estimation method of t_max could not affect the results derived from this study.

Hoenig (1983) equations have been commonly used in estimating M of elasmobranchs (Cortés, 1995; Sminkey and Musick, 1995; Au and Smith, 1997) including various skates and rays (Frisk et al., 2002, 2005; Davis et al., 2007). Although additional empirical equations such as Jensen (1996); Mollet and Cailliet (2002), and Hewitt and Hoenig (2005) have been proposed, the estimation of M for elasmobranchs is difficult, as these equations are dependent on life history parameters with uncertainty. The empirical equations derived from this study for estimation of λ are not dependent on M, thus reducing the uncertainty.

The high correlation between predicted λ′ and λ for Groups 1–3 and the randomly distributed residuals suggest that the empirical equations developed in this study can predict λ precisely. To validate the empirical equations derived from this study, one independent data set of Dipturus trachyderma with L_m/L_∞ = 0.8366, f/R_c = 48.7, k = 0.081 year^–1, T_max = 35.6 years, and T_m = 17.4 years (Licandeo et al., 2007) was used. This species falls in Group 2 based on its life history parameters. The predicted values of λ′ (1.1305) from the empirical equation were close to the value of λ derived from conventional demographic method (1.1434), suggesting that the empirical equations derived from this study can predict λ of other skates and rays accurately.

There is closely relation between life history trait and marine habitats. García et al. (2008) proposed three major marine habitats (continental shelves, open sea, and deep sea). The continental shelves include pelagic and benthic waters that have high primary production and large environmental variations, while the open ocean habitat has low primary production and the deep sea habitat is dark and cold with almost no primary production and relies on the organisms falling down for meso-pelagic waters. Compared with pelagic environment, deep sea species grow slowly, mature late, and has prolonged longevity, and their resilience to overexploitation is low.

Rajiformes inhabit high latitude and deep waters, while Myliobatiformes inhabit low latitude and shallow waters (Ebert and Compagno, 2007). King and McFarlane (2003) concluded that opportunistic species stay in shallower waters. while equilibrium species inhabit deeper waters. Similar results were found in this study. The skates and rays in Group 1 or Group 2 identified in this study have habitats in higher latitude with low SST and low primary production, which cause the slow growth of these species. These are the characteristics of equilibrium species for Raja spp. in addition to longheaded eagle ray, blue stingray, diamond stingray, common stingray, bat eagle ray, Myliobatis californica, and cownose ray. The species in the Group 3 inhabit lower latitude with higher SST, shallow waters and have higher primary production. These species have fast-growing characteristics identified in opportunistic species such as starry skate, little skate, common guitarfish, and the species in Urolophidae.

The species in Group 1 have similar life history characteristics as equilibrium species of sharks but with early maturation. These species may be vulnerable or even encounter regional collapsed if catch control or fishing effort monitoring is not implemented such as in diamond stingray (Smith et al., 2008). The species in Group 2 have a general growth rate with late maturity characteristics. It is recommended to apply the management measure of periodical species such as reducing fishing pressure by modifying the fishing gears to these species. In addition, as the exploitation rate of elasmobranchs is related to the size of fish, the large-size species such as yellownose skate (Quiroz and Wiff, 2005) should be monitored with caution. The L_b/L_∞ of the aforementioned two groups was smaller than that of Group 3. Thus, protection of young-of-year is recommended for the species in Groups 1 and 2, as small-size neonates generally mean high mortality. Barnett et al. (2013) applied age-structured demographic models to examine the effect of life history parameters change on population growth rates for five deep sea skate species, which fall in these groups, and they concluded that gear modifications or depth-specific effort controls may be an effective management measure. The species in Group 3 have higher L_b/L_∞ but low fecundity, they are vulnerable to overexploitation, and the recovery time is long. Mollet and Cailliet (2002) proposed to protect adults to ensure the sustainability for stingray in this group. In addition to aforementioned management measures, Enever et al. (2009) mentioned that the survival rate of rays caught by trawl fishery was 55%. Therefore, it is suggested to release the skates and rays to reduce their mortality rate.

It might be too late to take management actions after the species-specific full stock assessment has been made. In this study, an alternate approach was provided to estimate the finite rate of population increase. The empirical equations developed for each group provide accurate predictions of λ by reducing the bias of estimation resulting from parameter uncertainties. The results derived from this study can be used in the implementation of management measures for data-limited skates and rays in a precautionary manner, as the single species stock assessment and management measure are difficult to apply to those skates and rays that were commonly by-caught in trawling fishery. Ecological risk assessment based on the λ′s estimated from our empirical equations can be an alternate approach to evaluate the risk of several species for precautionary purposes. Furthermore, for those skates or rays without detailed life history, parameters can be classified into one of three groups based on similar species identified in this study, and then provisional management actions as aforementioned can be taken until more detailed life history research can be conducted. This study considered only 35 of 574 batoid species (41 stocks) in the world (Ebert and Compagno, 2007). The life history parameters used in this study were mainly from skates (23 stocks), but those for manta rays, devil rays, and guitarfish were lacking. Therefore, the estimates derived from this study may not cover all life history traits of skates and rays. Future studies should focus on collecting updated life history parameters of more species particularly on those not covered by this study to improve the robustness of these empirical equations. These empirical equations should be updated regularly according to the availability of new information of life history parameters.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

K-ML and Y-WH conceived and designed the experiments, performed the experiments, analyzed the data, and wrote the manuscript. H-HH contributed reagents, materials, and analysis tools. All authors contributed to the article and approved the submitted version.

Funding

This study was financially supported by the Ministry of Science and Technology, Taiwan, under Contract No. MOST 105-2313-B-019-005-MY3.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.664611/full#supplementary-material

References