Devil in the details: growth, productivity, and extinction risk of a data-sparse devil ray

Devil rays (Mobula spp.) face rapidly intensifying fishing pressure to meet the ongoing international trade and demand for their gill plates. This has been exacerbated by trade regulation of manta ray gill plates following their 2014 CITES listing. Furthermore, the paucity of information on growth, mortality, and fishing effort for devil rays make quantifying population growth rates and extinction risk challenging. Here, we use a published size-at-age dataset for a large-bodied devil ray species, the Spinetail Devil Ray (Mobula japanica), to estimate somatic growth rates, age at maturity, maximum age and natural and fishing mortality. From these estimates, we go on to calculate a plausible distribution of the maximum intrinsic population growth rate (rmax) and place the productivity of this large devil ray in context by comparing it to 95 other chondrichthyan species. We find evidence that larger devil rays have low somatic growth rate, low annual reproductive output, and low maximum population growth rates, suggesting they have low productivity. Devil ray maximum intrinsic population growth rate (rmax) is very similar to that of manta rays, indicating devil rays can potentially be driven to local extinction at low levels of fishing mortality. We show that fishing rates of a small-scale artisanal Mexican fishery were up to three times greater than the natural mortality rate, and twice as high as our estimate of rmax, and therefore unsustainable. Our approach can be applied to assess the limits of fishing and extinction risk of any species with indeterminate growth, even with sparse size-at-age data.


20
Understanding the sustainability and extinction risk of data-sparse species is a pressing problem 21 for policy-makers and managers.
is challenge can be compounded by economic, social and . 66 First, we estimate growth parameters using a Bayesian approach that incorporates prior knowl-67 edge of maximum size and size at birth of this species, using the length-at-age data presented in tail Devil Ray. 78 We analyse a unique set of length-at-age data for a single population of M. japanica caught in a 79 Mexican artisanal shery. Individuals in this sample were limited to 110 and 240 cm disc width 80 (DW), which falls short (77%) of the maximum disc width reported elsewhere [14]. erefore, 81 we use a Bayesian approach to re t growth curves to this length-at-age dataset [15]. We use 82 published estimates of maximum size and size at birth to set informative priors. 83 We t the three-parameter von Bertalan y equation to the length-at-age data, combining 84 sexes: where L t is length at age t, and the growth coe cient k, size-at-age zero L 0 , and asymptotic 86 size L ∞ are the von Bertalan y growth parameters. ese parameters are conventionally pre-87 sented in terms of length; here we use length synonymously with disc width, such that L t , L 0 , 88 L ∞ , and L max are synonymous with DW t , DW 0 , DW ∞ , and DW max , respectively.

89
In order to account for multiplicative error, we log-transformed the von Bertalan y growth 90 equation and added an error term: is can be wri en as: We also constrain our prior for L 0 around size at birth, and use a beta distribution to constrain 106 our prior for growth coe cient k between zero and one, with a probability distribution that is 107 slightly higher closer to a value of 0.1: We compared the e ect of our informative priors on our posteriors with parameter estimates 109 with weaker priors, in which we maintained the mean of the distributions but increased their We also considered a scenario with uninformative priors, where all prior distributions are 112 uniform: In all models we set an weakly informative prior for the variance σ 2 , such that: A summary of the priors used can be seen in selective across each age-or size-class. We also assume that there is limited migration in and out 120 of this population. With these assumptions, counting the number of individuals captured in each 121 age-class represents the population age structure, which can be used to construct a catch curve.  is information is very valuable when inferring whether shing mortality F is unsustainable. 128 We calculate Z as the slope of the regression of the catch curve, including only those ages or sizes 129 that are vulnerable to the shery. 130 We removed age-classes that had zero individuals in our sample to be able to take the natural 131 logarithm of the count. Because there is uncertainty associated with the dataset (due to its rel-132 atively small size), we resampled a subset of the dataset 20,000 times, a er randomly removing 133 20% of the points. is allowed us to quantify uncertainty in our estimate of Z. For each subset, 134 we computed the age-class with the maximum number of samples, and removed all age-classes 135 younger than this peak. With the remaining age-classes, we t a linear regression to estimate the 136 slope which is equivalent to −Z. is method for estimating mortality relies on two assumptions 137 of the selectivity of the shery. First, catch is not size-selective once individuals are vulnerable 138 to the shery. Second, if young age-classes are less abundant than older age-classes, they are 139 assumed to have lower catchability. is is why we removed the younger age-classes before the 140 "peak" abundance of each sample, as this will a ect the steepness of the slope. year. us we draw b from a uniform distribution bound between 0.25 and 0.5.

158
Age at maturity (α mat ) ere are no direct estimates of age at maturity for any mobulid ray,   Fig. 1). e asymptotic size in the model with strong priors was closest to the maximum 180 observed size for this species (Fig. 2). Estimates of k were lowest in the model with strong priors 181 and highest in the model with uninformative priors (  (Fig. 3). ere were no devil rays aged 12 or 13, 185 and therefore these age-classes were removed from the catch curve analysis before ing each 186 regression. Because these age classes are some of the oldest, removing these points is likely to 187 provide more conservative estimates of total mortality. As a reference, we also ran the models Devil and manta rays have low intrinsic rate of population increase relative to other chondrichthyan 203 species (Fig. 4). Among species with similar somatic growth rates, the Spinetail Devil Ray has 204 the lowest r max value (black diamond in Fig. 4a). is contrast is strongest when excluding deep-205 water chondrichthyans (white circles in Fig. 4), which tend to have much lower rates of population  In this study, we examined multiple lines of evidence that suggest devil rays have relatively low 212 productivity, and hence high risk of extinction compared to other chondrichthyans. e r max of 213 the Spinetail Devil Ray is comparable to that of manta rays, and much lower than that of other 214 large planktivorous shallow-water chondrichthyans such as the Whale Shark and the Basking 215 Shark (Fig. 4). We conclude the comparable extinction risks of devil and manta rays, coupled with 216 the ongoing demand for their gill plates, suggest that conferring a similar degree of protection to 217 all mobulids is warranted.

218
Next we consider three key questions that arise from our analyses: (1) Why do mobulid rays Can small-scale sheries cause local declines in devil ray populations? 221 We found that the Spinetail Devil Ray has similar productivity to manta rays, despite this species 223 being half the size of the manta rays. is result suggests that smaller devil ray species are also 224 likely to have very low productivity, probably due to their very low reproductive rates. Mobulid year -1 ) is twice as high than our estimate of r max , which also represents the shing mortality F 259 expected to drive this species to extinction (F ext = 0.077 year -1 ) [28]. Even though our estimate of 260 shing mortality is highly uncertain (Fig. 3b)       (a) Catch curve of natural log abundance at age. e regression lines represent the estimated slopes when omi ing different age-class subsets (as shown by the horizontal extent of each line), and resampling 80% of the data. Note that individual estimates of Z di er in the number of age-classes included for its computation, resulting in regression lines of di erent lengths. (b) Violin plot of estimated total mortality (Z) values calculated using the bootstrap resampling method. Estimates from di erent age-classes suggest an estimate of Z ≈ 0.25 year -1 . e median is shown by the dark grey line.