Genetic diversity assessment for the vulnerable migratory cownose ray Rhinoptera bonasus (Myliobatiformes: Rhinopteridae) from the southwestern Atlantic Ocean

Bruno C. Souza¹ , Vanessa P. Cruz¹, Talita R. A. Almeida², João Bráullio L. Sales³, Luís Fernando S. Rodrigues-Filho⁴, Marcelo Vianna⁵, Matheus M. Rotundo⁶, Claudio Oliveira¹ and Fausto Foresti¹

PDF: EN    XML: EN  | Supplementary: S1  S2  S3  S4  S5  S6 | Cite this article

Abstract​

---

Introduction​

---

Sharks and rays are overexploited by overfishing, due to unsustainable fishing activities aimed at obtaining the meat and fins of these animals. Furthermore, due to their K-selected life history strategy, such as late sexual maturity, low fertility, slow growth and prolonged gestation (Dulvy et al., 2014), several species will suffer population declines in the near future. In addition, several species are categorized as under high threat or facing extinction according to the International Union for the Conservation of Nature (IUCN) red list of threatened species (IUCN, 2021).

Myliobatiformes form one of the largest batoid groups, comprising 12 families, six subfamilies and 64 genera, totaling about 242 species (Fricke et al., 2021). This group exhibits a worldwide distribution, inhabiting both marine and freshwater environments (Dulvy et al., 2014). Among this group, cownose rays (Rhinopteridae Family) display benthopelagic habits, inhabiting both continental shelves, the vicinity of offshore islands and bays and estuaries (Fisher et al., 2013). These rays are migratory, presenting significant dispersion and a wide distribution in temperate and tropical regions, usually occupying coastal areas (Collins et al., 2008; Goodman et al., 2011; Ogburn et al., 2018).

The cownose ray Rhinoptera bonasus (Mitchill, 1815) is a highly mobile, batoid ray that occurs in temperate and tropical coastal waters from the New England in the United States to northern Argentina (Last et al., 2016). Females reach sexual maturity between seven and eight years old, and males, between six and seven, with a gestation period lasting from 11 to 12 months, followed by the birth of one embryo a year (Fisher et al., 2013; Fisher et al., 2014). Ogburn et al. (2018) using acoustic telemetry to track tagged in female and male of R. bonasus in US Atlantic coast, detected movements migratory between summer pupping and mating habitats in estuaries, with strong philopatry behavior.

Rhinoptera bonasus is currently categorized as vulnerable (VU) (Carlson et al., 2020a) by the IUCN red list, as it undergoes significant fishing pressures, often captured as bycatch in trawling fisheries (beach seine, simple and double trawling) and gillnets (Fisher et al., 2013), however, due to the high morphological similarity with the vulnerable Brazilian cownose ray Rhinoptera brasiliensis Müller, 1836 (Carlson et al., 2020b), as a result, fishing indices off the coast of Brazil may be flawed.

Notably, despite the wide geographic distribution of R. bonasus, a lack of information concerning the genetic population diversity of this species is noted. Effective conservation of threatened species requires accurate taxonomic classification, information on their life histories, distribution and genetic diversity information (Moyle et al., 2013), which is paramount for the development of regulatory programs assisting in the recovery of endangered species (Sellas et al., 2015).

In this context, considering the significant deficiency of molecular data available for R. bonasus, the present study aimed to describe the genetic diversity of this species, for three Brazilian ecoregions, the Amazon (Pará – PA), Northeastern (Pernambuco – PE) and Southeastern (Rio de Janeiro – RJ; São Paulo – SP and Santa Catarina – SC) ecorregions, through five microsatellite nuclear markers, thereby, testing the (1) panmixia hypothesis for Brazilian ecoregions and (2) test the in order to assess the natal philopatry potential of the species in the southwestern Atlantic Ocean.

Material and methods

---

The specimens analyzed in this study were obtained by professional fishers or donation from other researchers from five states of Brazil, defined for ecoregions based in Spalding et al. (2007) and Bernardino et al. (2015), South Western Atlantic region, from Amazonic ecoregion (Pará – PA), Northeastern ecoregion (Pernambuco – PE), and Southeastern ecoregion (Rio de Janeiro – RJ, São Paulo – SP and Santa Catarina – SC) (Fig. 1). The samples collected by artisanal fishers, a small fragment of muscle tissue was removed from each individual and preserved in of muscle tissue (~1.0 cm³) were extracted and deposited in the collection of the UNESP Laboratory of Fish Biology and Genetics (LBGP) in Botucatu, São Paulo, Brazil and preserved in 96% ethanol at -18°C (Tab. S1).

FIGURE 1 | Map of sampling sites of Rhinoptera bonasus in the Brazilian coast. South Western Atlantic region, from Amazonic ecoregion (Pará – PA), Northeastern ecoregion (Pernambuco – PE), and Southeastern ecoregion (Rio de Janeiro – RJ, São Paulo – SP and Santa Catarina – SC). Numbers in parenthesis indicate number of samples. Red markers illustrate nursery regions described for rays species (Yokota, Lessa, 2006; Feitosa et al., 2017; Rangel et al., 2018; Silva, Vianna, 2018).

Total DNA was extracted from muscle tissues following the protocol proposed by Ivanova et al. (2006). Polymerase Chain Reactions (PCR) were performed in a thermocycler (Applied Biosystems) following the protocol proposed by McDowell, Fisher (2013). Fragment analysis was performed using an ABI 3130X1 Genetic Analyzer sequencer (Applied Biosystems^TM), and genotyping was done with the program Genemapper v. 4.0 (Applied Biosystems, Foster City CA, USA).

The polymorphism information content (PIC) was obtained by the CERVUS v3.0 program (Marshall et al., 1998). Analyzes for estimating the observed heterozygosity (H_O), expected heterozygosity (H_E), number of alleles per locus (A), number of private alleles (Ap), inbreeding coefficient (F_IS), total genetic diversity (F_IT) and the Hardy-Weinberg Equilibrium (HWE) were performed using the GENEALEX v6.5 (Peakall, Smouse, 2012). We calculated the same estimates average (number of alleles and number of private alleles per locus), using rarefaction in order to account for differences in sample size using the HP-Rare program (Kalinowski, 2005).

 Linkage disequilibrium was tested between all pairs of loci with 10,000 permutations, using the GENEPOP 4.2 (Rousset, 2008). Contemporary effective population size (Ne) was estimated with NeESTIMATOR v2 (Do et al., 2014), using the linkage disequilibrium-based and the molecular coancestry methods based on the microsatellite data.

Pairwise F_ST (Weir, Cockerham, 1984) and R_ST (Slatkin, 1995; Gemmell et al., 1997) values were estimated among all populations within species to determine the degree of genetic differentiation using ARLEQUIN 3.5 (Excoffier, Lischer, 2010). The statistical significance index was determined by non-parametric testing with 10,000 permutations (Excoffier, Lischer, 2010). In addition, a hierarchical Analysis of Molecular Variance (AMOVA) (Excoffier et al., 1992) was conducted to examine the genetic diversity within and among populations using Arlequin 3.5 (Excoffier, Lischer, 2010) and the statistical significance was assessed using 10,000 permutations. Eight distinct scenarios of population structure of Rhinoptera bonasus were tested according to different locations analyzed, grouped into scenario one: (SC + SP + RJ + PE) vs. (PA), scenario two: (SC + SP + RJ) vs. (PE + PA), scenario three: (SC + SP + RJ) vs. PE vs. PA, scenario four: (SC + SP) vs. RJ vs. (PE + PA), scenario five: (SC + SP) vs. RJ vs. PE vs. PA, scenario six: SC vs. (SP + RJ) vs. (PE + PA), scenario seven: SC vs. SP vs. RJ vs. (PE + PA) and scenario eight: SC vs. (SP + RJ) vs. PE vs. PA.

To summarize the patterns of variation the principal coordinate analysis (PCA) was also performed using GenAlEx version 6.5 software (Peakall, Smouse, 2012) based on the matrix of pairwise Nei’s genetic distance (Nei, Takezaki, 1983). A multivariate approach, Discriminant Analysis of Principal Components (DAPC) (Jombart et al., 2020) was applied, as implemented in R-package ADEGENET v.2.1 (Jombart et al., 2020). The cross-validation function xval-Dapc was used to set the number of principal components to 50 (Jombart, Ahmed, 2011). Cluster assignments were pre-defined to correspond to a priori defined collect locations. Bayesian clustering for microsatellites loci was also used to assess population relatedness and gene flow, using the program STRUCTURE v 2.3.4 (Pritchard et al., 2000). The Markov chain Monte Carlo (MCMC) was run for 1 million generations, with an initial burn-in of 10% steps discarded, and 20 interactions of each K and the admixture model. The K values were selected using the delta K (Evanno et al., 2005) method described by (Earl, vonHoldt, 2012) in Structure Harvester (https://taylor0.biolo gy.ucla.edu/structureHarvester/).

Gene flow was estimated by calculating the direction relative migration rate, assuming asymmetric bidirectional gene flow, using the divMigrate function (Sundqvist et al., 2016) of the package diveRsity (Keenan et al., 2013), based on effective number of migrants (Nm) (Alcala et al., 2014). The statistical significance of directional migration was calculated on the basis of 1000 bootstraps.

Contemporary effective population size (Ne) was estimated with NeESTIMATOR v2 (Do et al., 2014), using the linkage disequilibrium-based and the molecular coancestry methods based on the microsatellite data. In order to detect recent bottlenecks or expansions in each population we used the program BOTTLENECK 1.2.02 (Piry et al., 1999), using two tests: “Sign test” (Cornuet, Luikart, 1996) and the Wilcoxon sign-rank test” (Luikart et al., 1998), both based on the Infinite Alleles Model (IAM), Stepwise Mutation Model (SMM) and Two-Phase Model (TPM), with a P-value < 0.05. Isolation by distance (IBD) by a Mantel test, was tested by computing the regression of F_ST/1-F_ST on geographic distances and the level of significance determined by performing a test with ISOLDE in GENEPOP 4.2 (Rousset, 2008) based on 1000 randomization.

Results​

---

A total of 69 Rhinoptera bonasus individuals were genotyped at five microsatellite loci. No evidence of null alleles was observed. All comparisons indicate that the investigated loci were not in linkage disequilibrium after global and population-specific analyses (p> 0.05). The five microsatellites were proven suitable for the identification of R. bonasus polymorphisms, comprising four to 12 alleles, with polymorphic information content (PIC) ranging from 0.108 to 0.711 and the highest value associated with the Rbon_56 locus (Tab. S2). Additionally, observed (H_O) and expected (H_E) heterozygosities presented average values of 0.336 and 0.382, respectively (Tab. S2), exhibiting the lowest value in Rbon_41 (0.110, 0.137) and the highest in Rbon_56 (0.680, 0.733). Furthermore, the inbreeding coefficients (F_IS) for four of the analyzed loci exhibited a heterozygosity deficit, with the exception of Rbon_52 (-0.048). The mean value for all loci was of 0.144, p<0.008. The F_IT ranged from 0.129 (Rbon_56) to 0.507 (Rbon_30), and the F_ST for all loci was significant, with an average F_ST = 0.200 (Tab. S2). In addition, no loci deviated from Hardy–Weinberg equilibrium expectations.

The diversity estimators indicate that the number of alleles per locus (A) are similar among localities, except for São Paulo, which presented the highest value (A = 4.2) (Tab. 1). In addition, 17 private alleles (Ap) were identified for all sampling areas, most for the state of São Paulo (Ap = 8), though the estimates average based on rarefaction were extremely low, the values for a estimate ranged from 1.85 in Pará to 2.22 in Pernambuco, for the estimative Ap, ranging from 0.08 in Santa Catarina to 0.43 in Pernambuco.

H_E was slightly higher than H_O in most localities, indicating a heterozygosity deficit among the investigated samples. A positive inbreeding coefficient (F_IS) was observed, except for Pará, with H_O = 0.3916, H_E = 0.3830 and F_IS -0.0465 (Tab. 1). The overall effective numbers (Ne) estimated for R. bonasus for all sampling sites combined was of 302 (35–∞), with populations ranging from 1.8 (Santa Catarina) to ∞ (São Paulo) (Tab. 1).

TABLE 1 | Genetic diversity indices of the five microsatellite loci for Rhinoptera bonasus. n: number of individuals analyzed, Na: number of alleles. A: number of alleles per locus and Ap: number of private alleles, in parentheses is the same estimates average based on rarefaction; H_O: observed heterozygosity, H_E: expected heterozygosity, F_IS: inbreeding coeficiente, Ne: effective population size based on the molecular coancestry method; 95% CI: Single standard deviation with 95% confidence interval, % PL: percentage of polymorphic loci.

The pairwise F_ST indicated significant values for both São Paulo and Pará (0.1536) and Rio de Janeiro and Pará (0.2788), whereas the pairwise R_ST did not result in significant values (Tab. S3). The pairwise F_ST and R_ST indices demonstrated no consistent genetic structuring concordance among sampling sites. Interestingly, when analyzed F_ST and R_ST pairwise by ecoregions, the results of F_ST were significant between Amazonic ecoregion (Pará – PA) and Southeastern ecoregion (Rio de Janeiro – RJ, São Paulo – SP and Santa Catarina – SC) (Tab. S3).

A global AMOVA analysis of the 69 R. bonasus individuals revealed significant genetic heterogeneity (Φ_ST = 0.12; p < 0.05, Tab. S4). A hierarchical AMOVA was performed considering a hypothetical population grouping, with the variance components associated with each scenario as the result of the variations observed within populations (86.64–93.76%), and simulations were conducted with populations grouped in eight scenarios, revealing a significant genetic variance Φ_ST for six scenarios (Tab. S4).

The clustering pattern of the principal component analysis (PCA) based on Nei’s genetic distance (DA) indicated that the first axis was responsible for 74% of the molecular data variance, followed by 13% explained by the second axis, where the Pará distinctly separated from the other geographic populations (Fig. 2). At least two main clusters were identified by the DAPC analysis (Fig. 3) and a subdivision was observed. Individuals from Santa Catarina, São Paulo and Pernambuco were most closely grouped, while the Rio de Janeiro group was established a little more to one side and the Pará group was more isolated from the other clusters. The structure analysis conducted under the admixture model and K = 1–6 populations indicated a highest likelihood (ln(P)D) in a population structure of K = 2 (-2113.46±0.079) (Fig. 4).

FIGURE 2 | Principal component analysis (PCA) analysis based on Nei distances using five Rhinoptera bonasus microsatellite markers. Collection areas: Santa Catarina (SC, red), São Paulo (SP, green), Rio de Janeiro (RJ, blue), Pernambuco (PE, orange) and Pará (PA, black).

FIGURE 3 | Discriminant Analysis of Principal Components (DAPC) scatterplots for the five Rhinoptera bonasus microsatellite genotypes. Dots represent individuals, whereas coloured ellipses correspond to geographical populations. Collection areas: Santa Catarina (1, red), São Paulo (2, green), Rio de Janeiro (3, blue), Pernambuco (4, orange) and Pará (5, black).

FIGURE 4 | Bayesian clustering evidenced by STRUCTURE with K = 2. A vertical bar represents each individual, and the length of each bar indicates the probability of membership in each cluster. Collection areas: Santa Catarina (SC), São Paulo (SP), Rio de Janeiro (RJ), Pernambuco (PE) and Pará (PA).

Although no significant asymmetries were detected, relative pairwise migration and gene flow demonstrated a directional connectivity pattern among Santa Catarina, São Paulo, Rio de Janeiro and Pernambuco, except Pará which showed the lowest values of gene flow, with rates of immigration ranging from 0.026 (Rio de Janeiro) to 0.179 (São Paulo), and rates of emigration ranging from 0.061 (Rio de Janeiro) to 0.185 (Pernambuco) (Fig. 5; Tab. S5).

The recent population size reduction (bottleneck effect) test demonstrated a significant excess heterozygosity applying the infinite alleles model (IAM), two-phase model (TPM) and the stepwise mutation model (SMM), in the Sign test for Pará, indicating the occurrence of past bottlenecks. São Paulo and Pernambuco also presented significant excess of heterozygosity based on the SMM (Sign test) (Tab. S6). The Mantel test using the microsatellite data demonstrated a correlation between geographic distance and fixation index values (R² = 0.1472) with P = 0.001, respectively.

FIGURE 5 | Contemporary relative migration rates of Rhinoptera bonasus based on number effective of migrants (Nm) estimated from 5 microsatellite loci, estimated by divMigrate package. The arrows were weighted according to number of migrant values presented in S5, and arrowheads show the estimated direction of gene flow. Collection areas: Santa Catarina (1), São Paulo (2), Rio de Janeiro (3), Pernambuco (4) and Pará (5).

Discussion​

---

The present study assessed the genetic diversity and population genetic structure of cownose ray Rhinoptera bonasus populations sampled off the Brazilian coast using five microsatellite markers. The results show a moderate genetic diversity and a significant genetic structure was detected, demonstrating the presence of at least two regional management units, the first comprising Santa Catarina, São Paulo, Rio de Janeiro and Pernambuco and the other, Pará, thus rejecting the panmixia hypothesis.

Genetic diversity. The allelic variability level of the evaluated microsatellites in R. bonasus was of four, at 12 allelles per locus, considered moderate when compared to another marine ray, Dipturus trachyderma (Krefft & Stehmann, 1975), which presented from two to six (Vargas-Caro et al., 2017), and similar to Bathytoshia brevicaudata (Hutton, 1875), with two, at 15 alleles per locus (Le Port et al., 2016) and Dipturus chilensis (Guichenot, 1848), with four, at 10 allelles per locus (Vargas-Caro et al., 2017). Observed and expected heterozygosity estimates were similar across the four sampling localities, and all presented lower than expected values. However, the rarefaction is an approach used to produce comparable estimates of the number of alleles and number of private alleles in populations with varying sample sizes, if we compare the same rarefied estimates, the data found in this study are the lowest ever reported for rays. A proportion of alleles were identified as private alleles (17) across all localities, with the highest number of private alleles were observed in São Paulo (8), thereby, after rarefaction, the highest value was in Pernambuco.

Compared to other endangered rays species, the overall effective population size (Ne) for R. bonasus was modest (302), similar to that reported for the Sawfish smalltooth sawfish Pristis pectinata Latham, 1794 from the Gulf of Mexico, classified as Critically Endangered (Carlson et al., 2013), ranging from 250 to 350 (Chapman et al., 2011). In contrast, Aetobatus narinari (Euphrasen, 1790) from the Gulf of Mexico and coastal Atlantic waters off the United States, classified as near threatened (Kyne et al., 2006), presented a high value of 1.893 (Newby et al., 2014). Notwithstanding, the results of this study should be carefully evaluated, as the sample size at some locations was low and could affect these estimates. Even so, estimating the Ne of endangered populations and species is paramount, as loss of genetic diversity is unavoidable in small populations (Newby et al., 2014).

In this study, although we did not observe a heterozygosity deficit for the Pará locality (H_O = 0.3916, H_E = 0.3830), the results on population reduction indicate genetic bottleneck effects for the Pará locality in both models used (IAM, TPM and SMM) (Tab. 1; Tabs. S2 and S6). However, we carefully emphasize these identified results, given the low sample size analyzed for this region. Although the genetic variability of a population is usually measured in terms of average heterozygosity, it is important to know the number of alleles per locus, as if this number is drastically reduced after a population goes through a bottleneck, the adaptability of this population may be limited, even if heterozygosity remains high (Nei, 2005; Poirier et al., 2019).

Information that corroborates the small number of alleles per locus and the reduction of heterozygosity by loci observed in this study for the Pará locality (Tab. 1; Tabs. S2 and S6), are also found by Schultz et al. (2008), who identified recent bottleneck effects for a population of the species Negaprion acutidens (Rüppell, 1837) in French Polynesia, even though there was no deficit in heterozygosity (H_O = 0.45, H_E = 0.44), but low allelic diversity. The locations of São Paulo and Pernambuco showed significance for a possible bottleneck effect according to the SMM model, where we observed a deficit of heterozygosity, number of alleles and polymorphism per loci (Tab. 1; Tabs. S2 and S6). The excess heterozygosity observed for the Pará locality can be explained by the possible retention of genetic variability mediated by evolutionary forces, such as balancing selection, strong enough to contain the effects of drift (Aguilar et al., 2004; Fisher et al., 2021).

Population structure. Different studies have been carried out in the last decade involving biological cownose ray aspects, such as age, growth, reproductive and migration patterns, along the US Atlantic coast (Fisher et al., 2013; Ogburn et al., 2018). Several studies report migrations between summer habitats in estuaries south of Long Island for breeding and winter habitats along the coast of Florida for feeding, indicating concentrated ecological interactions in the spring and fall migrations (Ogburn et al., 2018).

Our results indicate a population structure with significant F_ST in two pairwise analyses (p <0.001), between São Paulo and Pará and between Rio de Janeiro and Pará, with a non-significant pairwise fixation index R_ST between all sampling areas.

Interestingly, the F-statistics (Φ_ST estimates) generated by the AMOVA revealed significant genetic heterogeneity (p < 0.05), with a significant genetic variation partition into six group within the eight tested scenarios. The microsatellite population-level based differentiation and PCA, DAPC, STRUCTURE analyses indicate that a moderate restricted genetic connectivity is present in the Amazon ecoregion (Pará) when compared to the genetic connectivity observed between the Northeastern (Pernambuco) and Southeastern (Rio de Janeiro, São Paulo and Santa Catarina) ecoregions, indeed, these data were corroborated by the analysis of gene flow, we observed high levels of gene flow between Santa Catarina, São Paulo, Rio de Janeiro and Pernambuco, but with low gene flow between Pará and others areas. Furthermore, reinforcing all the analyzes performed, the F_ST pairwise analyzed by ecoregion was significant between Amazonic ecoregion (Pará – PA) and Southeastern ecoregion (Rio de Janeiro – RJ, São Paulo – SP and Santa Catarina – SC.

In general, the combined population-level analyses outcomes suggest the possibility of ecological, oceanographic, or behavior barrier dispersion strong enough to restrict the connectivities between the sampling areas. The Brazilian coastline comprises approximately 8,000 km, featuring different marine ecoregions and boasting of highly productive and dynamic areas (e.g., coral reefs, estuaries, mangrove forests) (Moura et al., 2016), which vary from place to place, and may include isolation, upwelling events conditions, nutrient inputs, freshwater influxes, topography, temperature and sediment conditions, all of which may exert direct and indirect pressure on marine species populations.

The values detected by the fixation indices between the analyzed areas can reflect the different environments that these populations occupy. For example, in the Amazon ecoregion, where the Amazon River joins the saline waters of the continental shelf, a plume comprising the low salinity and turbid water masses that persists both temporally and spatially and extends for many kilometers is noted (Grodsky et al., 2014). The Northeastern ecoregion is characterized by a semi-arid climate undergping a discreet freshwater input from local rivers, and exhibiting a topographic barrier close to the Archipelago of Fernando de Noronha ecoregion and the São Pedro and São Paulo islands ecoregion (Castro Filho, Miranda, 1998). The Southeastern ecoregion presents a heterogeneous environment, with lagoons and a minor archipelago off the coast of Rio de Janeiro state (Gonzalez-Silveira et al., 2004), presenting less saline waters due to the mixing of the Brazil Current with coastal waters, as well as the Cabo Frio Magmatic Lineament, located in the Cabo Frio area (Zalán, Oliveira, 2005).

Another hypothesis that could explain the findings reported herein is the fact the cownose rays exhibit strong behavior philopatry with onsite reproductive fidelity (Ogburn et al., 2018), as reported in studies concerning for different taxonomic ray groups (Flowers et al., 2016), Pristis pristis (Linnaeus, 1758) (Feutry et al., 2015), A. narinari (Sellas et al., 2015), B. brevicaudata (Le Port et al., 2016) and Neotrygon kuhlii (Müller & Henle, 1841) (Borsa et al., 2012).

Concerning the Brazilian coast, four nursery regions are reported for rays, one for P. pectinata in the northern region of the Maranhão Amazon coast (Feitosa et al., 2017), another for R. bonasus and 11 other elasmobranch species in the state of Rio Grande do Norte, in the Brazilian northeast (Yokota, Lessa, 2006), and two nursery regions in southeastern Brazil, one in Rio de Janeiro, for Gymnura altavela (Linnaeus, 1758) (Silva, Vianna, 2018) and one in Bertioga, São Paulo, for R. bonasus (Rangel et al., 2018). Therefore, we believe that individuals from Pará may use the Maranhão nursery, while individuals from Pernambuco may use the Rio Grande do Norte nursery, and individuals from Santa Catarina, São Paulo and Rio de Janeiro, tend to return to the São Paulo nursery or other still unidentified nursery regions for reproduction.

In addition to a restricted genetic connectivity among R. bonasus populations from the Brazilian coast, the IBD test results suggest that individuals do not disperse over long distances, similar as observed to some shark species from the Atlantic Ocean, such as the lemon shark Negaprion brevirostris (Poey, 1868) (Ashe et al., 2015) and the scalloped hammerhead Sphyrna lewini (Griffith & Smith, 1834) (Pinhal et al., 2020), in which reproductive philopatric behavior is used to explain IBD results. Therefore, reproductive philopatric behavior is probably one of the main factors driving the observed R. bonasus population subdivision.

The present study indicates that the coastal ray R. bonasus exhibits a significant genetic differentiation between the sampled locations throughout the Brazilian coast through microsatellite marker analyses. Knowledge in this regard is important not only to understand the population structure of the species, but also to determine management unit conservation activities and the priorities of actions that must be taken for all R. bonasus geographic distribution regions. Therefore, future studies should include a broad sampling effort along the northern coast of Brazil, mainly in the states of Amapá and Maranhão, in order to shed new light on the philopatry of the Pará R. bonasus population. Furthermore, the evidence of new genetic stocks is critical to reinforce conservation polices for this vulnerable species.

Acknowledgments​

---

This study received support from Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), grants #2017/01288-9 (BCS), #2016/16483-9 (VPC). The authors thank all fishers from the “Pró-Pesca Project: fishing the knowledge” for the samples.

References​

---

Ashe JL, Feldheim KA, Fields AT, Reyier EA, Brooks EJ, O’Connell MT et al. Local population structure and context-dependent isolation by distance in a large coastal shark. Mar Ecol Prog Ser. 2015; 520:203–16. https://doi.org/10.3354/meps11069

Aguilar A, Roemer G, Debenham S, Binns M, Garcelon D, Wayne RK. High MHC diversity maintained by balancing selection in an otherwise genetically monomorphic mammal. Proc Natl Acad Sci USA. 2004; 101(10):3490–94. https://doi.org/10.1073/pnas.0306582101

Alcala N, Goudet J, Vuilleumier S. On the transition of genetic differentiation from isolation to panmixia: What we can learn from G_ST and D. Theor Popul Biol. 2014; 93:75–84. https://doi.org/10.1016/j.tpb.2014.02.003

Bernardino AF, Netto SA, Pagliosa PR, Barros F, Christofoletti RA, Rosa Filho JS et al. Predicting ecological changes on benthic estuarine assemblages through decadal climate trends along Brazilian marine ecoregions. Estuar Coast Shelf Sci. 2015; 166(Part A):74–82. https://doi.org/10.1016/j.ecss.2015.05.021

Borsa P, Arlyza IS, Laporte M, Berrebi P. Population genetic structure of blue-spotted maskray Neotrygon kuhlii and two other Indo-West Pacific stingray species (Myliobatiformes: Dasyatidae), inferred from size-polymorphic intron markers. J Exp Mar Biol Ecol. 2012; 438:32–40. https://doi.org/10.1016/j.jembe.2012.09.010

Carlson J, Charvet P, Avalos C, Blanco-Parra MP, Briones Bell-Iloch A, Cardenosa D et al. Rhinoptera bonasus. The IUCN Red List of Threatened Species Version 2021-2 [Internet]. IUCN; 2020a. https://dx.doi.org/10.2305/IUCN.UK.2020-3.RLTS.T60128A3088381.en

Carlson J, Charvet P, Avalos C, Blanco-Parra MP, Briones Bell-lloch A, Cardenosa D et al. Rhinoptera brasiliensis. The IUCN Red List of Threatened Species Version 2021-2 [Internet]. IUCN; 2020b. https://dx.doi.org/10.2305/IUCN.UK.2020-3.RLTS.T44595A2997621.en

Carlson J, Wiley T, Smith K. Pristis pectinata (errata version published in 2019). The IUCN Red List of Threatened Species Version 2021-2 [Internet]. IUCN; 2013. https://dx.doi.org/10.2305/IUCN.UK.2013-1.RLTS.T18175A141791261.en

Castro Filho BM, Miranda LB. Physical oceanography of the western Atlantic continental shelf located between 4 N and 34 S. The sea. 1998; 11(1):209–51.

Chapman DD, Simpfendorfer CA, Wiley TR, Poulakis GR, Curtis C, Tringali M et al. Genetic diversity despite population collapse in a critically endangered marine fish: The smalltooth sawfish (Pristis pectinata). J Hered. 2011; 102(6):643–52. https://doi.org/10.1093/jhered/esr098

Collins AB, Heupel MR, Simpfendorfer CA. Spatial distribution and long-term movement patterns of cownose rays Rhinoptera bonasus within an estuarine river. Estuaries Coast. 2008; 31(6):1174–83. https://doi.org/10.1007/s12237-008-9100-5

Cornuet JM, Luikart G. Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics. 1996; 144(4):2001–14. https://doi.org/10.1093/genetics/144.4.2001

Do C, Waples RS, Peel D, Macbeth GM, Tillett BJ, Ovenden JR. NeEstimator v2: Re-implementation of software for the estimation of contemporary effective population size (N_e) from genetic data. Mol Ecol Resour. 2014; 14(1):209–14. https://doi.org/10.1111/1755-0998.12157

Dulvy NK, Fowler SL, Musick JA, Cavanagh RD, Kyne PM, Harrison LR et al. Extinction risk and conservation of the world’s sharks and rays. eLife. 2014; 3:e00590. https://doi.org/10.7554/eLife.00590

Earl DA, vonHoldt BM. STRUCTURE HARVESTER: A website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv Genet Resour. 2012; 4(2):359–61. https://doi.org/10.1007/s12686–011-9548-7

Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study. Mol Ecol. 2005; 14(8):2611–20. https://doi.org/10.1111/j.1365–294X.2005.02553.x

Excoffier L, Lischer HEL. Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010; 10(3):564–67. https://doi.org/10.1111/j.1755-0998.2010.02847.x

Excoffier L, Smouse PE, Quattro JM. Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data. Genetics. 1992; 131(2):479–91. https://doi.org/10.1093/genetics/131.2.479

Feitosa LM, Martins APB, Nunes JLS. Sawfish (Pristidae) records along the Eastern Amazon coast. Endanger Species Res. 2017; 34:229–34. https://doi.org/10.3354/esr00852

Feutry P, Kyne PM, Pillans RD, Chen X, Marthick JR, Morgan DL et al. Whole mitogenome sequencing refines population structure of the Critically Endangered sawfish Pristis pristis. Mar Ecol Prog Ser. 2015; 533:237–44. https://doi.org/10.3354/meps11354

Fisher RA, Call GC, Grubbs RD. Age, growth, and reproductive biology of cownose rays in Chesapeake Bay. Mar Coast Fish. 2013; 5(1):224–35. https://doi.org/10.1080/19425120.2013.812587

Fisher RA, Call GC, McDowell JR. Reproductive variations in cownose rays (Rhinoptera bonasus) from Chesapeake Bay. Environ Biol Fishes. 2014; 97(9):1031–38. https://doi.org/10.1007/s10641-014-0297-9

Fisher KJ, Vignogna RC, Lang GI. Overdominant mutations restrict adaptive loss of heterozygosity at linked loci. bioRxiv. 2021. https://doi.org/10.1101/2021.05.13.444075

Flowers KI, Ajemian MJ, Bassos-Hull K, Feldheim KA, Hueter RE, Papastamatiou YP et al. A review of batoid philopatry, with implications for future research and population management. Mar Ecol Prog Ser. 2016; 562:251–61. https://doi.org/10.3354/meps11963

Fricke R, Eschmeyer WN, Fong JD. Eschmeyer’s catalog of fishes: Genera/Species by Family/Subfamily. [Internet]. San Francisco: California Academy of Sciences; 2021. Available from: http://researcharchive.calacademy.org/research/ichthyology/catalog/fishcatmain.asp

Gemmell NJ, Allen PJ, Goodman SJ, Reed JZ. Interspecific microsatellite markers for the study of pinniped populations. Mol Ecol. 1997; 6(7):661–66. https://doi.org/10.1046/j.1365-294X.1997.00235.x

Goodman MA, Conn PB, Fitzpatrick E. Seasonal occurrence of cownose rays (Rhinoptera bonasus) in North Carolina’s estuarine and coastal waters. Estuaries Coast. 2011; 34(3):640–51. https://doi.org/10.1007/s12237-010-9355-5

Gonzalez-Silvera A, Santamaria-del-Angel E, Garcia VMT, Garcia CAE, Millán-Nuñez R, Muller-Karger F. Biogeographical regions of the tropical and subtropical Atlantic Ocean off South America: Classification based on pigment (CZCS) and chlorophyll-a (SeaWiFS) variability. Cont Shelf Res. 2004; 24(9):983–1000. https://doi.org/10.1016/j.csr.2004.03.002

Grodsky SA, Reverdin G, Carton JA, Coles VJ. Year-to-year salinity changes in the Amazon plume: Contrasting 2011 and 2012 Aquarius/SACD and SMOS satellite data. Remote Sens Environ. 2014; 140:14–22. https://doi.org/10.1016/j.rse.2013.08.033

International Union for Conservation of Nature (IUCN). IUCN Red List of Threatened Species. The IUCN Red List of Threatened Species Version 2021-2 [Internet]. IUCN; 2021. Available from: https://www.iucnredlist.org

Ivanova NV, Dewaard JR, Hebert PDN. An inexpensive, automation-friendly protocol for recovering high-quality DNA: Technical Note. Mol Ecol Notes. 2006; 6(4):998–1002. https://doi.org/10.1111/j.1471-8286.2006.01428.x

Jombart T, Ahmed I. Adegenet 1.3-1: New tools for the analysis of genome-wide SNP data. Bioinformatics. 2011; 27(21):3070–71. https://doi.org/10.1093/bioinformatics/btr521

Jombart T, Kamvar ZN, Collins C, Luštrik R, Beugin M-P, Knaus BJ et al. Adegenet: Exploratory analysis of genetic and genomic data. Version 2.1 [Internet]. CRAN; 2020. Available from: https://cran.r-project.org/web/packages/adegenet/index.html

Kalinowski ST. HP-RARE 1.0: A computer program for performing rarefaction on measures of allelic richness. Mol Ecol Notes. 2005, 5(1):187–89. https://doi.org/10.1111/j.1471-8286.2004.00845.x

Keenan K, McGinnity P, Cross TF, Crozier WW, Prodöhl PA. DiveRsity: An R package for the estimation and exploration of population genetics parameters and their associated errors. Methods Ecol Evol. 2013; 4(8):782–88. https://doi.org/10.1111/2041-210X.12067

Kyne PM, Ishihara H, Dudley SFJ, White WT. Aetobatus narinari. The IUCN Red List of Threatened Species Version 2012.2 [Internet]. IUCN; 2006. Available from: https://www.iucnredlist.org/

Last PR, White WT, Carvalho MR, Séret B, Stehmann MFW, Naylor GJP, editors. Rays of the world. Australia: CSIRO Publishing; 2016.

Le Port A, Roycroft EJ, Thakur V, Lavery SD. Characterisation of eleven new polymorphic microsatellite markers for the coastal stingray Dasyatis brevicaudata (Dasyatidae Hutton, 1875), and cross-amplification in seven dasyatid species. Biochem Syst Ecol. 2016; 65:234–37. https://doi.org/10.1016/j.bse.2016.02.009

Luikart G, Allendorf FW, Cornuet J-M, Sherwin WB. Distortion of allele frequency distributions provides a test for recent population bottlenecks. J Hered. 1998; 89(3):238–47. https://doi.org/10.1093/jhered/89.3.238

Marshall TC, Slate J, Kruuk LEB, Pemberton JM. Statistical confidence for likelihood-based paternity inference in natural populations. Mol Ecol. 1998; 7(5):639–55. https://doi.org/10.1046/j.1365–294x.1998.00374.x

McDowell JR, Fisher RA. Final report discrimination of cownose ray, Rhinoptera bonasus, stocks based on microsatellite DNA markers. Marine Resource Report No. 2013-8 [Internet]. Virginia: Virginia Institute of Marine Science; 2013. http://dx.doi.org/doi:10.21220/m2-7jky-2s80

Moyle PB, Kiernan JD, Crain PK, Quiñones RM. Climate change vulnerability of native and alien freshwater fishes of California: A systematic assessment approach. PLoS ONE. 2013; 8(5):e63883. https://doi.org/10.1371/journal.pone.0063883

Moura RL, Amado-Filho GM, Moraes FC, Brasileiro PS, Salomon PS, Mahiques MM et al. An extensive reef system at the Amazon River mouth. Sci Adv. 2016; 2(4):e1501252. https://doi.org/10.1126/sciadv.1501252

Nei M. Bottlenecks, genetic polymorphism and speciation. Genetics. 2005; 170(1):01–04. https://doi.org/10.1093/genetics/170.1.1

Nei M, Takezaki N. Estimation of genetic distances and phylogenetic trees from DNA analysis. Proceeding 5th World Congress on Genetics Applied to Livstock Production. 1983; 21(21):405–12.

Newby J, Darden T, Shedlock AM. Population genetic structure of spotted eagle rays, Aetobatus narinari, off Sarasota, Florida and the Southeastern United States. Copeia. 2014; 2014(3):503–12. https://doi.org/10.1643/CG-13-122

Ogburn MB, Bangley CW, Aguilar R, Fisher RA, Curran MC, Webb SF et al. Migratory connectivity and philopatry of cownose rays Rhinoptera bonasus along the Atlantic coast, USA. Mar Ecol Prog Ser. 2018; 602:197–211. https://doi.org/10.3354/meps12686

Peakall R, Smouse PE. GenAlEx 6.5: Genetic analysis in Excel. Population genetic software for teaching and research—an update. Bioinformatics. 2012; 28(19):2537–39. https://doi.org/10.1093/bioinformatics/bts460

Pinhal D, Domingues RR, Bruels CC, Ferrette BLS, Gadig OBF, Shivji MS et al. Restricted connectivity and population genetic fragility in a globally endangered Hammerhead Shark. Rev Fish Biol Fish. 2020; 30(3):501–17. https://doi.org/10.1007/s11160-020-09607-x

Piry S, Luikart G, Cornuet J-M. Computer note. BOTTLENECK: A computer program for detecting recent reductions in the effective population size using allele frequency data. J Hered. 1999; 90(4):502–03. https://doi.org/10.1093/jhered/90.4.502

Poirier M-A, Coltman DW, Pelletier F, Jorgenson J, Festa–Bianchet M. Genetic decline, restoration and rescue of an isolated ungulate population. Evol Appl. 2019; 12(7):1318–28. https://doi.org/10.1111/eva.12706

Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000; 155(2):945–59. https://doi.org/10.1093/genetics/155.2.945

Rangel BS, Rodrigues A, Moreira RG. Use of a nursery area by cownose rays (Rhinopteridae) in southeastern Brazil. Neotrop Ichthyol. 2018; 16(1):e170089. https://doi.org/10.1590/1982-0224-20170089

Rousset F. GENEPOP’007: A complete re-implementation of the GENEPOP software for Windows and Linux. Mol Ecol Resour. 2008; 8(1):103–06. https://doi.org/10.1111/j.1471-8286.2007.01931.x

Schultz JK, Feldheim KA, Gruber SH, Ashley MV, McGovern TM, Bowen BW. Global phylogeography and seascape genetics of the lemon sharks (genus Negaprion). Mol Ecol. 2008; 17(24):5336–48. https://doi.org/10.1111/j.1365-294X.2008.04000.x

Sellas AB, Bassos-Hull K, Pérez-Jiménez JC, Angulo-Valdés JA, Bernal MA, Hueter RE. Population structure and seasonal migration of the spotted eagle ray, Aetobatus narinari. J Hered. 2015; 106(3):266–75. https://doi.org/10.1093/jhered/esv011

Silva FG, Vianna M. Diet and reproductive aspects of the endangered butterfly ray Gymnura altavela raising the discussion of a possible nursery area in a highly impacted environment. Braz J Oceanogr. 2018; 66(3):315–24. https://doi.org/10.1590/S1679-8759201801906603

Slatkin M. A measure of population subdivision based on microsatellite allele frequencies. Genetics. 1995; 139(1):457–62. https://doi.org/10.1093/genetics/139.1.457

Spalding MD, Fox HE, Allen GR, Davidson N, Ferdaña ZA, Finlayson M et al. Marine ecoregions of the world: A bioregionalization of coastal and shelf areas. BioScience. 2007; 57(7):573–83. https://doi.org/10.1641/B570707

Sundqvist L, Keenan K, Zackrisson M, Prodöhl P, Kleinhans D. Directional genetic differentiation and relative migration. Ecol Evol. 2016; 6(11):3461–75.

Vargas-Caro C, Bustamante C, Bennett MB, Ovenden JR. Towards sustainable fishery management for skates in South America: The genetic population structure of Zearaja chilensis and Dipturus trachyderma (Chondrichthyes, Rajiformes) in the south-east Pacific Ocean. PLoS ONE. 2017; 12(2):e0172255. https://doi.org/10.1371/journal.pone.0172255

Weir BS, Cockerham CC. Estimating f-statistics for the analysis of population structure. Evolution. 1984; 38(6):1358–70. https://doi.org/10.2307/2408641

Yokota L, Lessa RP. A nursery area for sharks and rays in northeastern Brazil. Environ Biol Fishes. 2006; 75(3):349–60. https://doi.org/10.1007/s10641-006-0038-9

Zalán PV, Oliveira JAB. Origem e evolução estrutural do Sistema de Riftes Cenozóicos do Sudeste do Brasil. B Geoci Petrobras. 2005; 13(2):269–300.

Authors

---

Bruno C. Souza¹ , Vanessa P. Cruz¹, Talita R. A. Almeida², João Bráullio L. Sales³, Luís Fernando S. Rodrigues-Filho⁴, Marcelo Vianna⁵, Matheus M. Rotundo⁶, Claudio Oliveira¹ and Fausto Foresti¹

^[1]    Departamento de Biologia Estrutural e Funcional, Instituto de Biociências, Universidade Estadual Paulista Júlio de Mesquita Filho – UNESP, 18618-689 Botucatu, SP, Brazil. (BCS) brunocampos.ibb@gmail.com (corresponding author); (VPC)
cruzvp@outlook.com; (CO) claudio.oliveira@unesp.br; (FF) f.foresti@unesp.br.

^[2]    Departamento de Ciências Químicas e Biológicas, Instituto de Biociências, Universidade Estadual Paulista Júlio de Mesquita Filho – UNESP, 18618-689 Botucatu, SP, Brazil. (TRAA) aleixo.talita@yahoo.com.br.

^[3]    Universidade Federal do Pará, Grupo de Investigação Biológica Integrada, Centro de Estudos Avançados da Biodiversidade, 66075-750 Belém, PA, Brazil. (JBLS) braullio@ufpa.br.

^[4]    Universidade Federal Rural da Amazônia, Campus de Capanema, 68600-030 Capanema, PA, Brazil. (LFSRF) luis.filho@ufra.edu.br.

^[5]    Laboratório de Biologia e Tecnologia Pesqueira, Universidade Federal do Rio de Janeiro – UFRJ, 21941-901 Rio de Janeiro, RJ, Brazil. (MV) mvianna@biologia.ufrj.br.

^[6]    Universidade Santa Cecília – UNISANTA, Acervo Zoológico – AZUSC, 11045907 Santos, SP, Brazil. (MMR) mmrotundo@unisanta.br.

Authors Contribution

---

Bruno C. Souza: Conceptualization, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing-original draft, Writing-review and editing.
Vanessa P. Cruz: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing-original draft, Writing-review and editing.
Talita R. A. Almeida: Conceptualization, Methodology, Software, Visualization, Writing-original draft, Writing-review and editing.
João Bráullio L. Sales: Data curation, Methodology, Visualization, Writing-original draft, Writing-review and editing.
Luís Fernando S. Rodrigues-Filho: Data curation, Methodology, Visualization, Writing-original draft, Writing-review and editing.
Marcelo Vianna: Data curation, Methodology, Visualization, Writing-original draft, Writing-review and editing.
Matheus M. Rotundo: Data curation, Visualization, Writing-original draft, Writing-review and editing.
Claudio Oliveira: Data curation, Funding acquisition, Project administration, Resources, Supervision, Writing-original draft, Writing-review and editing.
Fausto Foresti: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Writing-original draft, Writing-review and editing.

Ethical Statement​

---

All samples and specimens of R. bonasus were collected in strict accordance with the regulations of the Brazilian Federal Animal Ethics Committee (SISBIO 13843–1) and the analyzes followed the International Guidelines for Animal Experiments, as authorized by CEEAA-IBB/UNESP, protocol number 556.

Competing Interests

---

The authors declare no competing interests.

How to cite this article

---

Souza BC, Cruz VP, Almeida TRA, Sales JBL, Rodrigues-Filho LFS, Vianna M, Rotundo MM, Oliveira C, Foresti F. Genetic diversity assessment for the vulnerable migratory cownose ray Rhinoptera bonasus (Myliobatiformes: Rhinopteridae) from the southwestern Atlantic Ocean. Neotrop Ichthyol. 2021; 19(4):e210077. https://doi.org/10.1590/1982-0224-2021-0077

Copyright​

---

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

Distributed under

Creative Commons CC-BY 4.0

© 2021 The Authors.

Diversity and Distributions Published by SBI

Accepted XXXX

Submitted April 6, 2021

Epub December 13, 2021