Genetic variability and population structure of two sympatric cownose rays Rhinoptera (Myliobatiformes, Rhinopteridae) in the Western Atlantic Ocean

Sample collection and DNA extraction

We acquired tissue samples from specimens collected (as bycatch by artisanal fisheries and research sampling efforts) at different localities from the western Atlantic Ocean coast to southeastern Brazil. The specimens analyzed in this study were collected by personnel from several research institutions: samples from waters in the USA were obtained through collaborations with the Florida Fish and Wildlife Conservation Commission through the Fisheries-Independent Monitoring program. Samples from Mexico and Colombia were obtained and identified by author Palacios-Barreto from artisanal fisheries, and those from Brazil in Bertioga were gathered by researchers collaborators of this article, under authorization from IBB-UNESP CEUA No. 1662170524. Morphological dental features such as the number of rows of tooth plates, were used to distinguish between the two species: R. bonasus typically exhibits seven rows in each jaw, whereas R. brasiliensis typically has nine rows in each jaw (Last et al., 2016). Additionally, DNA-based species identification methods were employed. A similarity test was conducted between our sequences and reference libraries in GenBank^® (Sayers et al., 2018) using similarity scores (e.g., BLAST) (Altschul et al., 1990), with the species name assigned to the sequence showing the highest similarity.

We analyzed samples from a total of 208 cownose ray individuals, of which 103 were identified as R. bonasus (Fig. 1), collected from the Tropical Northwest Atlantic (TNA): IR, Indian River Lagoon (n = 13); CH, Charlotte Harbor (n = 6); TB, Tampa Bay (n = 21); Colombia: Co, Manaure (n = 13); warm temperate Northwest Atlantic (WTNA): CK, Cedar Key (n = 14); AB, Apalachicola Bay (n = 8), and warm temperate Southwestern Atlantic (WTSA): SP, Bertioga (n = 28). The remaining 105 specimens were representative of R. brasiliensis (Fig. 1), which were obtained from the warm temperate Northwest Atlantic (WTNA): CK, Cedar Key (n = 1); tropical Northwest Atlantic (TNA): CC, Cape Canaveral (n = 10); IR, Indian River Lagoon (n = 2); Tampa Bay (n = 1); TT, Tampico (n = 2); TV, Tamiahua (n = 8); CV, Chachalacas (n = 29); CS, Seybaplaya (n = 32); Q, Chiquilá (n = 9); Colombia: Co, Ciénaga (n = 2), and warm temperate Southwestern Atlantic (WTSA): SP, Bertioga (n = 10). All tissue samples were preserved in 95% ethanol until extraction. Additionally, the mtDNA COI dataset was complemented with 100 sequences of R. bonasus from five locations (IR, CH, TB, CK, and AB) and 12 sequences of R. brasiliensis from two locations (IR and CK), obtained in previous studies (Florida Fish and Wildlife Conservation Commission, Fish and Wildlife Research Institute). In addition, we included mitochondrial COI, and Cytb genes available in GenBank (Table S1). Whole-genome DNA (gDNA) was isolated using the traditional phenol:chloroform: isoamyl alcohol (25:24:1 v/v) protocol (Sambrook, Fritsch & Maniatis, 1989). The DNA was resuspended in 50 μL TE buffer 1X (pH 8) and stored at 4 °C for posterior procedures. Field species identification was confirmed using mitochondrial sequence data of mitochondrial cytochrome c oxidase subunit I (COI) subunit I (Ward et al., 2005).

Mitochondrial DNA sequencing

We amplified partial sequences of two mitochondrial (mtDNA) genes: cytochrome oxidase subunit I (COI) and cytochrome b (Cytb) were amplified by Polymerase Chain Reactions (PCR) using the primers: FishF1 (5′-TCAACCAACCACAAAGACATTGGCAC-3′), FishR1 (5′-TAGACTTCTGGGTGGCCAAAGAATCA-3′) (Ward et al., 2005) for COI whereas the primers described by McDowell & Fisher (2013); the RbCytbF (5′-GGCCTHTTYCTRGCTATACACTACACAC-3′), RbCytbR (5′-AGGGRTGGAATGGRATTTT-3′), were for Cytb. The PCR occurred in 15 μL of the final volume containing 1.5 mM of MgCl2, 0.2 mM of each dNTP, 10 pm/μL of each primer, and 1 U Taq DNA polymerase and 1 μL of DNA template (20 ng/μL). For the mitochondrial COI gene, we set thermal profiles as follows: initial denaturation of 95 °C for 2 min, 35 cycles at 94 °C for 30 s, 54 °C for 1 min, 72 °C for 1 min, and a final extension at 72 °C for 10 min. For the Cytb gene, the conditions were: initial denaturation of 94 °C for 3 min, 35 cycles at 94 °C for 1 min, 56 °C for 1 min, 72 °C for 1 min, and a final extension at 72 °C for 7 min. We visualized all PCR products by 1.5% agarose gel electrophoresis, then were sequenced (including amplicon purification) by MacroGen (MacroGen Inc., Seoul, South Korea). GenBank accession numbers for each COI and Cytb haplotypes are listed in supporting information (Table S2; accession numbers OR710910–OR710922). All sequences were analyzed and edited in BioEdit v7.0.5.3 (Hall, 1999) and then aligned using the MUSCLE algorithm (Edgar, 2004) in Mega 6.0.

Genetic clustering

We use a Bayesian Analysis of Population Structure approach implemented in BAPS v. 5.4 (Corander et al., 2008), to identify distinct mitochondrial genetic clusters by assigning individuals to populations. For that, we used a clustering mixture analysis with linked loci, K = 1–5 with 20 replicates and 100 iterations for each run (Corander & Marttinen, 2006; Corander et al., 2008).

Genetic diversity and spatial distribution

The indices of genetic diversity, such as the number of polymorphic sites (S), number of haplotypes (H), haplotype diversity (hd), and nucleotide diversity (π), were estimated for each locality separately, then all were estimated considering the grouping determined by BAPS (Table S3), and subsequently for the entire set of individuals of each species, using DnaSP v.5.1 (Librado & Rozas, 2009). For each molecular marker, we constructed the haplotype networks through the median-joining method (Bandelt, Forster & Röhl, 1999) included in PopART (Leigh & Bryant, 2015). To better predict the pattern of genetic diversity distribution of each genetic group identified with BAPS, we applied the hypothesis of continuous population spreading by interpolating haplotype diversity (hd) data from each locality to adjacent areas along their known geographic distribution. The hd values of each genetic group were interpolated by the Kriging method (geostatistical interpolation method that creates a smooth surface even when sampling is spatially unequally, Oliver & Webster, 1990) in QGIS^® v.3.22 software; then, surface maps of the interpolated hd for each group were made, representing the spatial patterns of the genetic groups inferred by BAPS.

Population structure analyses

The population differentiation through an analog of Wright’s pairwise F_ST fixation index (Φ_ST), and the Analysis of Molecular Variance (AMOVA) (Excoffier, Smouse & Quattro, 1992) were conducted to examine the null hypothesis of panmixia, these analyses were evaluated among all sampling sites and BAPS clustering results using Arlequin 3.5 (Excoffier & Lischer, 2010) with 10,000 permutations and α = 0.05. For R. bonasus three grouping scenarios were considered: (1) by biogeographic bioregionalization of Coastal and Shelf Areas proposed by Spalding et al. (2007) (Fig. 1) which, five of six provinces were represented (cold temperate Northwest Atlantic (CTNA), warm temperate Northwest Atlantic (WTNA), tropical Northwestern Atlantic (TNA), North Brazil Shelf (NBS) and warm temperate Southwestern Atlantic (WTSA)). (2) Previous genetic groups identified in the Northwest Atlantic × Northern Gulf of Mexico by Carney et al. (2017) and (3) considering the clustering determined by BAPS. For R. brasiliensis only scenarios corresponding to Spalding bioregionalization and results from BAPS were considered (scenarios 1 and 3, respectively). All P-values were corrected by the FDR method (Benjamini & Hochberg, 1995).

Environmental association analyses

To assess the correlation between oceanic landscape features with the genetic diversity of cownose rays, the environmental variables, maximum and minimum benthic (B) and surface (S) values of temperature (BTmax, BTmin, STmax, STmin), salinity (BSmax, BSmin, SSmax, SSmin), dissolved oxygen (BOmax, BOmin, SOmax, SOmin), current velocity (BCVmax, BCVmin, SCVmax, SCVmin) and bathymetry (BAT) were downloaded from the Bio-ORACLE database (https://www.bio-oracle.org). Spearman correlations were run to avoid autocorrelation among them (Table S4). We used ecological modeling based on canonical correspondence analysis (CCA) to combine genetic and oceanographic data (Jombart, Pontier & Dufour, 2009) with ‘cca’ function in a vegan package in R (Oksanen et al., 2022). For that, we use scaled haplotypic (hd) and nucleotide diversity (π) considering the genetic groups detected by BAPS for R. bonasus and R. brasiliensis, and non-correlated scaled oceanic variables.

Historical demography analyses

The demographic history of cownose rays was inferred from several approaches: mismatch distribution analysis, neutrality tests Fu’s (1996) Fs and Ramos-Onsins & Rozas’s (2002) R2 implemented in DNAsp v.5.1 (Librado & Rozas, 2009). The goodness of fit between the observed data and the expected mismatch distribution under each demographic scenario was evaluated with Harpending’s (1994) raggedness (r) index in ARLEQUIN (Rogers & Harpending, 1992). Each demographic scenario was run with 10,000 bootstrap replicates.

Mitochondrial DNA sequencing

We analyzed COI sequences of 551 base pairs (bp) from a total of 230 sequences of Rhinoptera bonasus and 181 sequences of R. brasiliensis (including GenBank accession numbers OR710910–OR710922). Additionally, we obtained a fragment of 423 bp from the Cytb locus, including 108 sequences from R. bonasus and 105 sequences from R. brasiliensis. The concatenated databases (COI+Cytb) reduced the number of individuals and the number of locations of R. bonasus and R. brasiliensis, rendering the information inconclusive. Consequently, the results of subsequent analyses are not presented.

Genetic clustering

The clustering analysis conducted by BAPS revealed several coexisting genetic groups, indicative of well-differentiated groups within the same locality, for both genes and species (Figs. 2A, 2B, 2E, 2F). For R. bonasus, five groups were identified with COI: three overlapping groups among the cold, warm, and tropical biogeographic provinces of the Northwestern Atlantic, including Chesapeake Bay (BCh), Indian River Lagoon (IR), Charlotte Harbor (CH), Tampa Bay (TB), Cedar Key (CK), and Apalachicola Bay (AB); the fourth group encompassed localities in Colombia and northern Brazil (Tropical Northwestern Atlantic and North Brazil Shelf), while the fifth group corresponded to the São Paulo locality in the warm temperate Southwestern Atlantic biogeographic province (Figs. 2A, 3, Table S3). Cytb results indicated a similar pattern to COI: two geographically overlapping genetic groups in the Northwest Atlantic, from Chesapeake Bay to Apalachicola Bay, and one group in the tropical Southwest Atlantic encompassing localities in Colombia and São Paulo (Figs. 2B, 4, S1).

For R. brasiliensis, coexisting genetic group patterns were observed at broader scales than for R. bonasus. With COI, a group was identified spanning from Cape Canaveral (northernmost locality) to São Paulo (southernmost locality), the second extending from Indian River Lagoon (Florida) to Pará (Brazil), the third comprising only individuals from Pará (North Brazil Shelf biogeographic province), and two groups from São Paulo (Figs. 2E, 4, Table S3). Clustering obtained with Cytb did not exhibit a clear geographic pattern, as several coexisting groups were found, primarily in Gulf of Mexico localities; however, one group from the São Paulo locality was well-differentiated, although some individuals exhibited 100% similarity with individuals (n = 7) from the north (from Cape Canaveral and the Gulf of Mexico) (Fig. 2F, Table S3).

Genetic diversity

Haplotype diversity (hd) and nucleotide diversity (π) values were higher in R. bonasus as compared to R. brasiliensis for both genes.

Rhinoptera bonasus. The COI dataset exhibited 19 polymorphic sites, resulting in 11 haplotypes whereas the Cytb gene exhibited 11 polymorphic sites and 12 haplotypes (Table 1). Haplotype diversity (hd) was highest when assessed among localities, whereas, within the genetic groups generated by BAPS (Table S3), lower variation was observed among overlapping localities. The overall haplotype diversity (hd) among localities was 0.780 ± 0.016 for the COI gene and 0.735 ± 0.026 for the Cytb gene. Regarding nucleotide diversity (π), the overall value was 0.0140 ± 0.016 for the COI gene and 0.00848 ± 0.00030 for the Cytb gene (Table 1).

The spatial relation of mitochondrial haplotypes from both loci revealed two primary phylogroups. Notably, there were exclusive haplotypes for Colombia and Brazil, which were separated from the CTNA, WTNA, and TNA biogeographic provinces by eight mutational steps in the COI gene and interestingly, there were shared haplotypes between Colombia and Brazil (Fig. 2A). In contrast, the separation occurred with only three mutational steps in the Cytb gene (Fig. 2B).

The interpolation of haplotype diversity (ihd) allowed us to visualize the geographical patterns of the genetic groups detected by BAPS. For COI, the genetic group with the highest ihd was group 4 (Colombia and northern Brazil, ihd = 0.28, Fig. 3B), for the overlapping groups in the northwest Atlantic (from Chesapeake to Apalachicola Bays) it was group 3 (ihd = 0.22, Fig. 3), while group 5 (Brazil) obtained low ihd values (ihd = 0.0005, Fig. 3D). For Cytb, the two overlapping genetic groups in the Northwest Atlantic showed similar values of ihd (ihd = 0,22 and ihd = 0.23, respectively, Fig. S1) and the group that includes Colombia and São Paulo had the highest ihd (ihd = 0.687, Fig. S1).

Rhinoptera brasiliensis. The COI sequences from 181 individuals revealed the presence of eight polymorphic sites, resulting in 10 haplotypes. In contrast, the fragment from the Cytb locus displayed six polymorphic sites, corresponding to eight haplotypes (Table 1).

Like R. bonasus, in R. brasiliensis, a higher haplotype diversity (hd) was observed across localities (Table S3), compared to the genetic groups generated in BAPS when clustering with overlapping localities. Considering the overall dataset, both markers exhibited low haplotype (hd COI = 0.518 ± 0.031; hd Cytb = 0.404 ± 0.057) and nucleotide diversities (π COI = 0.0011 ± 0.0000; π Cytb = 0.0010 ± 0.0001).

In general, the mtDNA haplotype network showed limited haplotype exchange between distant locations and a close relationship among haplotypes within the same localities or areas. Analysis of haplotype networks based on the COI gene (Fig. 2G) revealed the presence of two major haplotypes, encompassing all sampled localities. Conversely, the haplotype network constructed using the Cytb gene exhibited a dominant haplotype encompassing most localities, except for the Tampa Bay and Tampico localities (Fig. 2H).

The ihd of COI, for overlapped groups on a broad scale (from the US to Brazil), ranged from ihd = 0.0001 (group 1, Fig. 4B), to ihd = 0.247 (group 2, Fig. 4C). In São Paulo were identified two overlapped groups whose ihd was contrasting (group 4: ihd = 0.33 and group 5: ihd= 0.00003, Figs. 4C, 4D). Finally, the third group located in Pará, (Brazil) had the lowest idh (ihd= 0.00003, Fig. 4D). For Cytb, the broadly overlapped genetic groups range from ihd = 0.0001 (groups 1 and 4, Fig. S2) to ihd= 0.22 (group 2, Fig. S2). However, the presence of one group mainly inside the Gulf of Mexico (group 3) had the highest ihd (ihd = 0.98, Fig. S1) while the most southern group (group 5 in São Paulo) had ihd = 0.22 (Fig. S2C).

Population structure

The global population structure reveals the presence of significant population genetic structure in both species. Bayesian clusters in BAPS were consistent with the segregation pattern observed for haplotype networks based on COI and Cytb markers for both species.

The optimal scenario from this Bayesian clustering analysis defined K = 5 genetic groups for the COI gene, while for the Cytb gene, the best scenario was K = 3 for R. bonasus (Figs. 2A, 2B). In contrast, in R. brasiliensis, it defined a K = 5 for both genes (Figs. 2E, 2F).

For R. bonasus, the pairwise Φ_ST analysis revealed significant genetic differentiation primarily between the southeastern ecoregion, with the various groups from the Northern Gulf of Mexico and the Western North Atlantic region. However, a high genetic differentiation can be observed even between closely sampled regions (Table 2). When assessing divergence between regions, the pairwise Φ_ST estimates resulted in highly significant differences among these five biogeographic provinces (Table S4).