Life-history, exploitation and extinction risk of the data-poor Baraka's whipray ( Maculabatis ambigua ) in small-scale tropical fisheries

The Baraka's whipray ( Maculabatis ambigua ) is a major constituent of small-scale fisheries catch in the south-western Indian Ocean. Despite this, little is known of its life-history or exploitation status. We provide the first estimates of crucial life-history parameters and the maximum intrinsic population growth rate r max , using specimens collected from small-scale fisheries landings in Kenya, Zanzibar and Madagascar (with northern Madagascar representing a range extension for this species). We assess the relative risk of overexploitation by combining r max with estimates of total Z , fishing F , and natural M mortality, and an estimate of the exploitation ratio E . The data indicate that Baraka's whipray is a medium-sized, fast-growing, early maturing species, with a relatively long lifespan. This results in a high r max relative to many other elasmobranchs, which when combined with estimates of F suggests that the species is not at imminent risk of extinction. Yet, estimates of exploitation ratio E indicate likely overfishing for the species, with full recruitment to the fishery being post-maturation and exploitation occurring across a broad range of age and size classes. Thus, Bar-aka's whipray is unlikely to be biologically sustainable in the face of current fisheries pressures. This paper makes an important contribution to filling the gap in available data and is a step towards developing evidence-based fisheries management for this species. Further, it demonstrates a simple and widely applicable framework for assessment of data-poor elasmobranch exploitation status and extinction risk.


| INTRODUCTION
Elasmobranchs (sharks and rays) generally display slow growth, late maturity and low fecundity (Compagno, 1990). These life-history traits mean elasmobranchs are intrinsically sensitive to non-natural mortalities and limits their recovery potential (Dulvy et al., 2014a). However, there is considerable variation in the life-history traits both among (Jacobsen & Bennett, 2011;Stevens & McLoughlin, 1991) and within (Jacobsen & Bennett, 2010;Lombardi-Carlson et al., 2003;O'Shea et al., 2013) species. Fisheries are the most prominent source of nonnatural mortalities for elasmobranchs at the global level (Dulvy et al., 2014a;Worm et al., 2013). Understanding species and stock-composition of elasmobranch exploitation in SWIO fisheries is extremely poor. Recent vulnerability assessments based on small-scale fisheries landings suggest that several coastal rays, primarily whiptail stingrays (Family Dasyatidae), may be at risk across the SWIO (Temple et al., 2019). Many of these species have either limited or no regional life-history data available. Rays contribute nearly half of SWIO smallscale fisheries landed elasmobranch catch (by weight and number) and originate from many of the same fisheries as sharks (Temple et al., 2019). Despite this, rays have received little consideration in SWIO elasmobranch management formulation.
Although only recently described (Last et al., 2016b), Baraka's whipray (Maculabatis ambigua) is the most common ray in Kenyan smallscale fisheries catch and a common constituent of small-scale fisheries catch in Zanzibar. Baraka's whipray are primarily caught in bottom-set gillnets (Barrowclift et al., 2017;Temple et al., 2019) and are regularly bycatch in trawl fisheries elsewhere within their range (Last et al., 2016b).
Baraka's whipray was originally thought to be distributed from Zanzibar to the Red Sea and possibly further into the northern Indian Ocean (Last et al., 2016b). Little is known of the life history of this species.
Here, we investigate key life-history parameters of Baraka's whipray through the production of disc width-weight relationships, age-structured growth models, and maturation and longevity estimates. Subsequently, we estimate maximum intrinsic population growth rate, which is combined with estimates of female fisheries and natural mortalities and female exploitation ratio to assess relative extinction and overexploitation risk of this data-poor species.

| Sample collection
Baraka's whipray specimens (n = 48) from small-scale fisheries landings were sampled with the consent of fishers and/or merchants at the village of Mkokotoni and the Darajani Market, Stone Town in Zanzibar between 28 July 2015 and 19 August 2015 (Figure 1). Where possible, we recorded the disc width (cm), weight (kg) and sex of each specimen. Male maturity status (immature or mature) was recorded based on calcification of claspers (Walker, 2005), with only those specimens exhibiting complete calcification considered mature. Maturity status could not be recorded for females because fishers and merchants did not consent to examination of the reproductive tract.
Vertebrae were extracted from the mid-disc of 47 specimens to be used for aging and stored at −20 C. Supplementary data for Baraka's whipray (n = 155) were recorded by trained fisheries observers across Kenya, Zanzibar and northern Madagascar during a 12-month landings monitoring programme between June 2016 and June 2017 (Temple et al., 2019) (Figure 1). Disc width, weight, sex and male maturity status were recorded. Observers also opportunistically recorded female maturity, with only those observed with fertilized eggs considered as mature. Both sets of specimens were caught across a range of gear types (bottom-set gillnets = 78, drift gillnets = 38, handlines = 8, longlines = 6, beach seine = 2, ringnet = 2, and unknown = 69). We assume that all catches originate from the same population or stock. Specimens where disc width was not recorded (n = 29) were disregarded from all analyses. Data analyses and visualizations were carried out and produced using R statistical software, version x64 3.6.0 (R Core Team, 2019).

| Disc width-weight relationship
Disc width and weight data of all specimens that had not been gutted (n = 100) were natural log (ln) transformed and their relationship described using linear models. Cook's distance (4/n) was used to identify data points exerting undue influence on linear models, likely resulting from measurement and/or data entry errors; these points were removed (n = 5) and the models subsequently re-run. The sex of specimens was included as an interaction variable and compared to the null model using ANOVA to determine whether the effect of sex was significant. The 95% confidence interval (95%CI) of the linear models was derived using bootstrapping with replacement for 10,000 iterations.

| Age estimation
Two vertebrae from each sampled specimen were cleaned of excess muscle and connective tissue, and both neural and haemal arches were removed. Subsequently, vertebrae were immersed in a 5% sodium hydrochloride solution for 10-30 min, dependent on the vertebrae size and quantity of remaining tissues. Samples were washed, then towel and air dried. Cleaned vertebrae were embedded in clear epoxy resin (Buehler EpoxiCure; Lake Bluff, Illinois, United States). A single sagittal-plane section was taken from each vertebra using a slow-speed precision saw with a diamond wafering blade (Buehler IsoMet Low Speed Precision Cutter). Several section widths were trialled (600, 450, 300, 200 and 150 μm), with 200 μm producing the highest readability. Sections were mounted permanently onto glass slides (Fisher Chemical DPX Phthalate Free Mounting Media; Hampton, New Hampshire, United States) and photographed using a digital macro-lens camera (Nikon SLR D7200; Tokyo, Japan). Image enhancement for growth band reading was carried out in Adobe Photoshop CS3 (San Jose, California, United States) (Campana, 2014).
The age (to the nearest 0.5 years) for each individual specimen was determined by examination of paired opaque and translucent banding in the vertebral corpus calcareum. The birth band was distinct but did not show a clear change in angle on the corpus calcareum ( Figure 2). Age was estimated from the independent examination of both sagittal sections (images were randomized before reading) by two independent readers. Mean age estimates for each specimen were generated for each reader and compared. Where reader mean estimates differed by <1 year the mean was taken as the best estimate. This allowable difference was more conservative than in other studies in light of the restricted sample size (Jacobsen & Bennett, 2011;Smith et al., 2008). If differences in mean age estimates between readers were >1 year, ages were re-estimated with both readers present. If readers could not agree (n = 0) samples would have been discarded (Goldman, 2005).
Commonly, measures of agreement, precision and bias in age reads are given as percentage agreement (PA), PA ±1 year, the coefficient of variation (CV) and the average percentage error (APE) (Beamish & Fournier, 1981;Chang, 1982). These measures are presented but are commonly recognized as imperfect (Cailliet et al., 2006; This individual was aged 6 years Goldman, 2005). We suggest that the Bland-Altman approach, designed for method comparison, provides improved quantification and visualization of agreement, precision and bias among reads and readers compared to the standard methods used in aging studies (Bland & Altman, 1999. Bias in the relationship between reads (within and among readers) is assessed through linear modelling of the mean age read for each specimen against the difference between reads for each specimen. Precision in age reads (within and among readers) is described by the limits of agreement (LOA) defined by the 95% mean confidence interval of the difference between reads. We display the results of the Bland-Altman method as the primary measure of agreement, precision and bias.
Validation of growth band periodicity was not possible within this study. The short temporal period within which samples were collected and the restricted sample size meant that both marginal incremental analysis and edge analysis, which are the most common validation methods for elasmobranchs (Cailliet et al., 2006), were not possible to conduct. Furthermore, mark-recapture of chemically marked or captive reared individuals was not feasible within the constraints of this study.

| Growth modelling
Age and disc width data were fit to three growth models using nonlinear regression (package nls). Reasonable starting values for growth model parameters were estimated using self-starter functions (packages FSAtools and stats) before fitting. Models were run for males and females both separately and combined. Model selection was made through comparison of Akaike's information criterion corrected for small sample size (AICc). Kimura's likelihood ratio test was used to determine whether sex had a significant effect on the growth model (Haddon, 2010). The 95%CI for growth curves and growth curve coefficients were derived via bootstrapping with replacement for 10,000 iterations. The growth models fit were: the three-parameter von Bertalanffy growth function (Von Bertalanffy 1938), the three-parameter Gompertz growth function (Ricker 1975), and the logistic growth function (Ricker, 1979), where DW t is the disc width at age t, DW ∞ is the asymptotic disc width, DW 0 is disc width at age zero, k is a growth constant and α is the inflection point.

| Estimation of the maximum intrinsic population growth rate
The maximum intrinsic population growth rate r max for female Baraka's whipray, as it is females that constrain population growth, was estimated using the simplified Euler-Loktka equation, accounting for juvenile mortality (Cortés, 2016;Pardo et al., 2016b): where l tmat is survival to maturity, b is the annual reproductive output of female offspring, t mat is age at maturity and Mis the instantaneous rate of natural mortality. We calculate l tmat as: There are no data for the annual reproductive output of Baraka's whipray, thus we provided a plausible estimate of b from related species, as has been done for other data-poor species ( reported as 1-4 and 1-3 pups, respectively (Last et al., 2016a;White et al., 2006). We assume reproduction is annual and that the sex ratio is 1:1. Hence, 0.5-2 female pups per year was considered a plausible range for b for the purpose of estimating r max .
We cannot directly estimate t mat for female Baraka's whipray because the infrequent observations of female maturity (n = 2) were both opportunistic and reliant upon the presence of fertilized eggs.
However, other studies of whiptail stingrays suggest that males and females mature at a similar size (Da Silva et al., 2018;Jacobsen & Bennett, 2010, 2011White, 2007). As such, we assume female size at maturity matches that of males. Disc width at which 50% of males reach maturity (White, 2007) was estimated from logistic regression and bootstrapping with replacement for 10,000 iterations used to define 95%CI. This male size at maturity range was then converted into age using the disc width and age relationship from the relevant selected growth model, providing a range estimate for female t mat .
Lastly, in order to estimate M we first need to estimate longevity t max . Estimates for t max were bound by the theoretical age at which female Baraka's whipray reached 95% DW ∞ and 99% DW ∞ (Fabens, 1965;Ricker, 1979). 95% DW ∞ is calculated as 5ln(2)k −1 and 99% DW ∞ calculated as 7ln(2)k −1 . Subsequently, M was calculated as: The simplified Euler-Loktka equation, solving for r max , was then optimized using the nlminb function (package stats). We made 1,000,000 draws from the estimated ranges of b, t mat and M, assuming a uniform probability distribution within these ranges in order to explore the full range of potential estimates for r max , particularly given the uncertainty in these parameters (Dulvy et al., 2014b).
2.6 | Estimating total mortality, fisheries mortality and the exploitation ratio  (c) Logistic regression with bootstrapped 95% confidence intervals describing the relationship between disc width and maturity status (immature or mature) for male Baraka's whipray, disc width at 50% maturity is displayed beyond age of full recruitment is reasonable given that the samples are drawn from a range of gear types across a large area and are thus unlikely to be strongly size selective. Similarly, the fishery is likely to have reached equilibrium given that the coastal elasmobranch fisheries of this region are well established (Temple et al., 2018). However, predation pressures on this species are unknown. Species which continue to grow substantially after age at full recruitment may show nonconstant M across ages as they outgrow some of their predators (Lorenzen, 2000). Additionally, nothing is known of the population structure of this species. Thus, it is possible that all assumptions of the catch curve may not be fully met with potential resultant biases on Z estimates. The Chapman-Robson method was selected as the regression estimator because traditional catch curves may show negative bias in mortality estimation (Smith et al., 2012).
Age was estimated, using the best fitting growth model, for all recorded catches with disc width data (n = 127) from the 12-month landings monitoring programme. Age of full recruitment to the fishery was assumed as the age class at peak abundance and Z was estimated from 1 year after the age of peak abundance (Smith et al., 2012). We made 1,000,000 draws from the estimated ranges of M and Z, assuming uniform and normal probability distributions, respectively. Ranges for M were subtracted from those of Z to estimate F, which was subsequently compared to r max and the female exploitation rate E. We drew inference on the likely biological sustainability status of Baraka's whipray by comparing F/Z to an optimum value of 0.5 (Gulland, 1971;Pauly, 1983).

| Ethical statement
All samples in this study were derived from specimens originating from fisheries catches. The study did not provide inducements (financial or otherwise) for the provision of specimens. As such, there were no relevant animal welfare implications.

| Age estimation
Age was successfully estimated for all 47 specimens. Bland-Altman analyses showed no evidence of bias within reader 1 (R 2 = 0.038, F = 1.78, P = 0.189; Figure 4a,b) or reader 2 (R 2 = 0.001, F = 0.049, P = 0.825, Figure 4c,d). However, evidence of significant bias between readers was found (R 2 = 0.118, F = 6.04, P = 0.018, Figure 4e,f), with reader 1 producing higher estimates than reader 2 for older specimens. LOAs are presented alongside standard precision metrics (Table 1). Variability in age band counts was not considered to be unusually high when compared to other studies (Baje et al., 2018;Gutteridge et al., 2013;Jacobsen & Bennett, 2010). Consensus was reached between readers for all specimens with initial disparities >1 year, suggesting that the significant bias identified by the Bland-Altman approach was likely to have been overcome. Furthermore, the Bland-Altman approach demonstrated no evidence for an increase or decrease in discrepancies between age-reads with increasing age, indicating consistency in the variability within and among readers, and thus increasing confidence in the validity of band reads across agespectra. Agreement between readers was likely influenced, in part, by the substantial subannuli banding prominent in young specimens. We recommend that future aging studies consider the use of this method when presenting the results of reader agreement, precision and bias in banding counts.

| Growth modelling
Age estimates for specimens ranged from 0 to 12.5 years for males and 0 to 17 years for females. Of the growth models tested, the three-parameter von Bertalanffy yielded the best fit for females, males and combined sex data. Kimura's likelihood ratio test indicated no significant effect of sex on the von Bertalanffy growth model fit

| Estimating maximum intrinsic population growth rate
Disc width at which 50% of males reach maturity was estimated to be 57.9 cm (95%CI 47.5, 62.7) (Figure 3c). These estimates equate to half of males reaching maturity at 2.18 years old (95%CI 1.12, 2.79). The smallest of the female specimens classified as mature had a disc width of 62 cm, equivalent to an age of 2.7 years. For the purpose of estimating r max we assume female t mat to range uniformly from 1.12 to 2.79 years. Estimates for t max were calculated based on age at 95% DW ∞ (Ricker, 1979) and at 99% DW ∞ (Fabens, 1965), yielding 14.4 years and 20.2 years, respectively. However, the maximum age directly observed during this study was 17 years. Thus, for the purpose of estimating r max we assume t max to range uniformly from 17 to 20.2 years. Using the drawn-down estimates of t max and t mat , we calculate median M to be 0.097 (95th percentiles 0.089, 0.107). Using the distributions of b, t max , t mat and M, we calculate the median r max for Baraka's whipray to be 0.446 (95th percentiles 0.220, 0.781; Figure 5a).  Given this and the estimate for M, we estimate median F to be 0.139 (95th percentiles −0.029, 0.307), considerably lower than estimates for r max . Median E is estimated from F and Z at 0.588 (95th percentiles −0.349, 0.764). The proportion of the estimated distribution of E greater than the optimal value of E = 0.5 is 68.7%, hence there is a reasonably high likelihood of overfishing for Baraka's whipray (Figure 5b).

| DISCUSSION
We provide the first estimates for key life-history parameters, extinction risk and exploitation status of the recently described and widely fished Baraka's whipray. Baraka's whipray is a medium-sized, relatively fast-growing and early maturing species of whiptail stingray that exhibits a moderately long lifespan. These life-history traits mean that Baraka's whipray is potentially more resilient to fisheries pressure than many other widely exploited elasmobranchs, reflected in our relatively high estimate for r max . The disparity between r max and female F estimates suggest that Baraka's whipray should not currently be considered at high extinction risk. Yet, the plausible distribution of the  managed effectively there is potential for this species to become one of the biologically sustainable elasmobranch fisheries (Simpfendorfer & Dulvy, 2017). Our study also reveals a larger maximum size (recorded disc width of 112 cm) and wider geographic distribution than previously known (Last et al., 2016b), as evidenced by catches of Baraka's whipray in northern Madagascar. In the following sections we consider (a) the life-history traits of Baraka's whipray relative to other whiprays and (b) the priority data gaps for future research.
Growth models from this study indicate that Baraka's whipray is a fast-growing species with early maturation relative to other similarly sized whiptail stingrays. For example, blackspotted whipray (male: k = 0.03, female: k = 0.03), blue stingray (Dasyatis chrysonota) (male: k = 0.175, female: k = 0.07) and diamond stingray (Hypanus (formerly Dasyatis) dipterura) (male: k = 0.1, female: k = 0.05) (Cowley, 1997;Jacobsen & Bennett, 2011;Smith et al., 2007). However, it must be considered that k may exhibit considerable variability, a product of the constraints of both the study (e.g., sample size and representation of size classes) and the growth-models used (Cailliet & Goldman, 2004;Smith et al., 2007). This may result in wide variability of k even within species, as demonstrated in Kuhl's maskray (Neotrygon kuhlii), where k ranges between 0.08 and 0.38 (Jacobsen & Bennett, 2010;O'Shea et al., 2013), which may reflect locally varying life histories or species complexes. The rapid growth of Baraka's whipray is consistent with the early maturation observed in males and females.
Relatively fast growth rates (k > 0.1) and early maturation are generally associated with higher potential rates of population increase and thus higher rebound potentials (Branstetter, 1990;Frisk et al., 2001;Musick, 1999). However, species with fewer than five female offspring per year tend to have very low estimates of r max (Pardo et al., 2016b). The early maturation of Baraka's whipray, combined with a moderately long lifespan relative to other whiptail stingrays, likely off-sets the apparently low number of female offspring per year and results in a high r max estimate relative to those of many threatened elasmobranchs (Dulvy et al., 2014b). However, uncertainty remains around the fecundity (litter size and inter-birth interval) of this species. The smallest specimens measured in this study and DW 0 values estimated from the joint male and female growth model (23.5 and 33.4 cm disc width, respectively) are large relative to DW ∞ estimates, suggesting high maternal investment in offspring and thus likely few offspring per litter. Verifying our assumptions for sex ratio, given that the degree of overfishing is likely to increase with rising regional fisheries pressure (Temple et al., 2018). In addition to probable overfishing, the age and size class selectivity of the fisheries is of further concern. Restricting elasmobranch fisheries to catches of nonadult age classes can be an effective management strategy for these taxa (Prince, 2002;Simpfendorfer, 1999) but may result in prohibativly low yields. Protection of subadults is likely optimal for maximizing the future reproductive potential of the stock (Kindsvater et al., 2016). Assuming the size-at-maturity estimates presented here are approximately representative, our data suggest that full recruitment to the fishery occurs primarily in young adults. Furthermore, exploitation continues across a broad range of adult age and size classes, indicative of a nonselective fishery. Whilst the nonselective nature of catches is not unexpected given the multigear, multispecies nature of SWIO small-scale fisheries, patterns of pressure across post-maturity life-stages is of great concern for the biological sustainability of both Baraka's whipray and numerous other elasmobranch species caught in these fisheries (Kiszka & van der Elst, 2015;Temple et al., 2019).
Lastly, we must note that our analyses assume that vertebral banding observed in Baraka's whipray are consistent, annual and continual. Such assumptions may lead to misclassification of age in elasmobranchs if violated (Harry, 2018;Kinney et al., 2016;Natanson & Cailliet, 1990), if there is difficulty in counting increasingly small bands in older animals, or if there is variability in band formation among vertebrae within individuals (Natanson et al., 2018). Whilst violation of annual banding assumptions have been seen in elasmobranchs, they have thus far been most common in large and longer-lived species.
Furthermore, annual band deposition has been validated in several other whiptail stingrays (Cowley, 1997;Jacobsen & Bennett, 2010;Pierce & Bennett, 2010;Smith et al., 2007). Thus, we consider our assumption to be reasonable but encourage continued vigilance for and research focus on banding periodicity in this and other elasmobranch species.
The life-history, r max , exploitation and mortality data information presented here for the recently described Baraka's whipray represent an important step towards improving evidence-based fisheries management in the region. However, the limitations of the study also outline the priority next steps required for the sustainable exploitation of this species. Improving our understanding of female maturation and fecundity is crucial in estimation of exploitation resilience. Furthermore, life-history traits may vary significantly between stocks, yet nothing is known of the structure of Baraka's whipray stocks within its range. Stock structure is key for defining management units for fisheries (Pita et al., 2016) and meeting the assumptions of stock assessment approaches, with differing life histories between stocks having implications for their resilience to fisheries exploitation and subsequent management needs. Given the prominence and likely economic importance of elasmobranchs, including Baraka's whipray, in SWIO small-scale fisheries (Temple et al., 2019), understanding stock structure and population dynamics are important in supporting and tailoring accurate and effective fisheries management actions at regional, national and local levels.