Growth, productivity, and relative extinction risk of a data-sparse devil ray

Devil rays (Mobula spp.) face intensifying fishing pressure to meet the ongoing international demand for gill plates. The paucity of information on growth, mortality, and fishing effort for devil rays make quantifying population growth rates and extinction risk challenging. Furthermore, unlike manta rays (Manta spp.), devil rays have not been listed on CITES. Here, we use a published size-at-age dataset for the Spinetail Devil Ray (Mobula japanica), to estimate somatic growth rates, age at maturity, maximum age, and natural and fishing mortality. We then estimate a plausible distribution of the maximum intrinsic population growth rate (rmax) and compare it to 95 other chondrichthyans. We find evidence that larger devil ray species have low somatic growth rate, low annual reproductive output, and low maximum population growth rates, suggesting they have low productivity. Fishing rates of a small-scale artisanal Mexican fishery were comparable to our estimate of rmax, and therefore probably unsustainable. Devil ray 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 and that a similar degree of protection for both groups is warranted.


Age-structured population modelling
The catch curve analysis makes the assumption that selectivity is constant for rays older than age 3. While it is possible that selectivity is biased by size in this fishery, we wanted to explore how constant selectivity affected the population age structure given the life history of M. japanica and our rage of fishing mortality estimates.
To show the relationship between fishing mortality, fishery selectivity, and population age structure, we simulated the deterministic dynamics of an age-structured population of female M. japanica with overlapping generations as functions of both age a and time t. We adapted the age-structured modelling approach in Kindsvater et al. (2016). Model equations and parameter descriptions are given in Table S1; parameter values are given in Table S2. We assumed that age a and disc width (which we refer to as L) are related according to the von Bertalanffy growth function (see Part 1 results, main text). We can calculate mass-at-age ( ) as a power function of L; this allows us track both age structure and biomass dynamics of our simulated population.
Our model follows the standard conventions of age-structured fisheries models (Mangel 2006). Annual egg production of the population ( ) depends on the spawning stock biomass, which we calculate as number of individuals of age , the probability they are mature at that age ( ), and the annual birth rate of female Spinetail Devil Rays . The age at which 50% of individuals are mature is defined as . We assumed the population dynamics are regulated by density-dependent recruitment, where the number of recruits N0 is determined by the population growth rate from low densities and the asymptotic population size, which depends on the parameters and in the Beverton-Holt stock-recruitment function (see Supplement 2 of Kindsvater et al. 2016). Lacking data on the stock recruitment relationships of mobulid ray populations, we assumed that the slope of the stock-recruitment function near the origin α (i.e., steepness; Mangel et al. 2010) was less than one, but was large enough for our population to be viable, which we determined through sensitivity analyses. We assumed all recruits are the same size, LR (in cm), and live to at most Amax years. Natural mortality is consistent with our assumptions in the main text, though the estimate is applicable to a population growing without density-dependent mortality, and is likely an underestimate here. Table S1. Equations used to model the different biological processes in the age-structured population model.

Biological process Equation Parameter interpretations
Von Bertalanffy growth function ∞ and determine size-at-age Maturation probability ogive is the age at 50% maturity; ogive steepness set by q Recruits[t+1] alpha = 0.9 beta = 4.5e-05 Figure S1. Life history functions that define growth, maturation, fecundity and recruitment for M. japanica (see Table S1 and Table S2). These functions are used in the age-structured population dynamics model used to generate Figs. S2-S3.
We simulated the population for several generations in the absence of fishing until it reached a stable age distribution (equilibrium). The proportion of individuals in each age class is constant at this point. We then added fishing mortality, again assuming constant selectivity for ages 3 and above. We varied fishing mortality according to the range of estimates we found with our catch curve analysis (Part 2, main text).
In Fig. S2 we show the stable age distribution with and without fishing for the lower value of our estimate of fishing mortality . In Fig. S3, we show the same distribution with the upper value of our estimate of fishing mortality . Notice that with constant selectivity for disc widths greater than 800 mm or 2 years of age, we are able to recreate the demographic pattern in our catch curve (Part 2, main text). This is because the population quickly moves into a stable age distribution; fishing the younger age classes thus erodes the number of individuals in larger age classes.  Figure S2. The stable age distribution before and after fishing. Both raw abundance (top row) and the natural log of abundance (bottom row) are plotted. Fishery selectivity was assumed to be constant after age 3, when M. japanica disc widths are greater than 80 cm (upper left panel, Fig. S1). In this case, fishing mortality F was 0.108, which is the lower limit on our estimate of F from the catch curve (main text).  Figure S3. The age distribution before and after fishing. This model run was identical to Figure S2, except fishing mortality was 0.16, which is the upper limit on our estimate of F from the catch curve (main text). This run of the model generates the age distribution most similar to that of the observed catch curve (Part 2, main text).