Introduction
Butterfly rays, Gymnura van Hasselt, 1823, are marine and demersal batoids displaying a worldwide geographic distribution, occurring in shallow waters with sandy and muddy bottoms in tropical and temperate climates (Muktha et al., Reference Muktha, Akhilesh, Sandhya, Jasmin, Jishnudev and Kizhakudan2016). Although Gymnura genus rays are easily distinguished from other ray species, significant morphological similarities are noted among congeneric species, making their taxonomy complex. Three nominal species have been described in the Atlantic, by 2017, namely G. altavela (Linnaeus Reference Linnaeus1758), G. natalensis (Gilchrist & Thompson 1911) and G. micrura (Bloch & Schneider 1801). Yokota & de Carvalho (Reference Yokota and de Carvalho2017) presented an extensive taxonomic revision of Gymnura micrura based on external and internal morphological features and considering specimens from the entire Atlantic Ocean geographic distribution. These authors redescribed G. micrura and two newly described species (G. lessae and G. sereti) previously included in G. micrura, highlighting the possibility of the occurrence of cryptic species within the genus. The possible occurrence of cryptic diversity is especially worrisome in species exploited by fisheries with a marked reduction in population size, as in the case of various Gymnura species.
The spiny butterfly ray, Gymnura altavela (Linnaeus Reference Linnaeus1758), is considered as presenting an amphiatlantic distribution, at depths ranging from 5–100 metres, targeted by intense commercial fisheries along its distribution (Yokota et al., Reference Yokota, White, de Carvalho, Last, White, de Carvalho, Séret, Stehmann and Naylor2016; ICMBio, 2018). In the western Atlantic Ocean, the species is distributed from New England (USA) to Argentina (Alkusairy et al., Reference Alkusairy, Ali, Saad, Reynaud and Capapé2014) (Figure 1). In Brazil, it is classified as a resident breeder and its presence is confirmed only for the south–south-east region. Furthermore, the only currently known nursery area for the species on the Brazilian coast is Guanabara Bay (Rio de Janeiro) (Gonçalves-Silva & Vianna, Reference Gonçalves-Silva and Vianna2018a, Reference Gonçalves-Silva and Vianna2018b). In the eastern Atlantic, this species has been recorded from Portugal to Angola (ICMBio, 2018) (Figure 1), and in the Mediterranean Sea, spread from Gibraltar to Lebanon (McEachran & Capapé, Reference McEachran, Capapé, Whitehead, Bauchot, Hureau, Nielsen and Tortonese1984). Spiny butterfly rays are economically important where they occur (Alkusairy et al., Reference Alkusairy, Ali, Saad, Reynaud and Capapé2014), and landings of the species in the Mediterranean have fluctuated spatially and in time (Capapé et al., Reference Capapé, Zaouali, Tomasini and Bouchereau1992).
Distribution for Gymnura altavela (black) and G. natalensis (grey) and G. altavela (black circle) and G. natalensis (grey circle) sampling sites.
Gymnura altavela is currently classified as globally Endangered (criteria A2d) by the International Union for Nature Conservation (IUCN) (Dulvy et al., Reference Dulvy, Charvet, Carlson, Badji, Blanco-Parra, Chartrain, De Bruyne, Derrick, Dia, Doherty, Dossa, Ducrocq, Leurs, Notarbartolo di Sciara, Pérez Jiménez, Pires, Seidu, Serena, Soares, Tamo, Vacchi, Walls and Williams2021) and locally considered Critically Endangered in both the Mediterranean Sea (Walls et al., Reference Walls, Vacchi, Notarbartolo di Sciara, Serena and Dulvy2016) and in Brazil (ICMBio, 2018).
On a time scale of only three generations (33 years, from 1982–2015), estimates based on landing data indicate that G. altavela populations have declined 42.5% in Morocco, 54% in Senegal and 98.8% in southern Brazil, with an overall population decline estimated as 50–79% during this time frame (Dulvy et al., Reference Dulvy, Charvet, Carlson, Badji, Blanco-Parra, Chartrain, De Bruyne, Derrick, Dia, Doherty, Dossa, Ducrocq, Leurs, Notarbartolo di Sciara, Pérez Jiménez, Pires, Seidu, Serena, Soares, Tamo, Vacchi, Walls and Williams2021). This dramatic reduction in population size is due to strong fisheries pressure, both as targeted and as bycatch due to the use of multiple types of fishing gear.
Despite its status, few studies on the species population biology and genetics are available, undermining conservation efforts, although it has increasingly attracted attention in the last decade. For example, several studies on G. altavela have been recently published in Brazil, on the following subjects: length-weight relationships (Silva-Junior et al., Reference Silva-Junior, Andrade and Vianna2011), embryo descriptions (Paiva et al., Reference Paiva, Julio, Marques and Vianna2018), distribution and density in an estuarine zone (Gonçalves-Silva et al., Reference Gonçalves-Silva and Vianna2018a), reinforcing its diet, reproductive aspects and a probable nursery area (Gonçalves-Silva et al., Reference Gonçalves-Silva and Vianna2018b), highlighting its fishing importance and molecular and morphometric relationships (Marques et al., Reference Marques, Guimarães, Sole-Cava and Vianna2019), revealing differences in in situ and ex situ bacterial communities associated with the skin and the stinger areas (Gonçalves-Silva et al., Reference Gonçalves-Silva, Dos Santos, de Assis Leite, Lutfi, Vianna and Rosado2020), and quantifying high Persistent Organic Pollutant contamination levels (Rosenfelder et al., Reference Rosenfelder, Lehnert, Kaffarnik, Torres, Vianna and Vetter2012; Paiva et al., Reference Paiva, Vannuci-Silva, Correa, Santos-Neto, Vianna and Lailson-Brito2021).
Given the taxonomic uncertainties associated with the species (Yokota & de Carvalho, Reference Yokota and de Carvalho2017), the present study aims to assess genetic relationships among Gymnura altavela specimens sampled along both sides of the Atlantic and Mediterranean. Comparisons between cytochrome oxidase subunit 1(COI) gene DNA sequences were performed in a phylogeographic framework. DNA sequence-based approaches are widely employed in studies that aim to inventory taxa diversity in groups affected by extraordinary morphological stasis and ecological traits (Cariani et al., Reference Cariani, Messinetti, Ferrari, Arculeo, Bonello, Bonnici, Cannas, Carbonara, Cau, Charilaou, El Ouamari, Fiorentino, Follesa, Garofalo, Golani, Guarniero, Hanner, Hemida, Kada, Lo Brutto, Mancusi, Morey, Schembri, Serena, Sion, Stagioni, Tursi, Vrgoc, Steinke and Tinti2017; Crobe et al., Reference Crobe, Ferrari, Hanner, Leslie, Steinke, Tinti and Cariani2021), with important contributions to the taxonomy of Gymnura species (Smith et al., Reference Smith, Bizzarro, Richards, Nielsen, Márquez-Flarías and Shivji2009; Shen et al., Reference Shen, Ma, Ni, Xu and Ma2012).
Materials and methods
A total of 95 Cytochrome oxidase I (COI) sequences from nine nominal Gymnura genus species were obtained from GenBank (Table 1; accession numbers in Table S2). Furthermore, we generated 37 sequences (GenBank numbers MW321984–MW322020) from G. altavela samples sampled off Rio de Janeiro, Brazil (Table S2) and identified according to McEachran & Carvalho (Reference McEachran, Carvalho and Carpenter2003). All sequences and information about specimens are recorded in the public BOLD Project ‘Gymnura_altavela_OTUs’ (project code: GYMNU). DNA extraction, PCR and sequencing procedures followed the methodology described in Marques et al. (Reference Marques, Guimarães, Sole-Cava and Vianna2019). Sequences for Himantura uarnacoides (N = 1), Mobula mobular (N = 2), Myliobatis chilensis (N = 1) and Myliobatis longirostris (N = 1) were used as outgroups (accession numbers available in Table S2).
Number of individual specimens per species
The sequences were aligned using the MEGA 7.0 software (Kumar et al., Reference Kumar, Stecher and Tamura2016) through the Clustal W algorithm (Thompson et al., Reference Thompson, Higgins and Gibson1994) and carefully checked visually. The HKY + G + I (Hasegawa et al., Reference Hasegawa, Kishino and Yano1995) nucleotide substitution model was selected using the Bayesian Information criterion as implemented in the aforementioned software (Kumar et al., Reference Kumar, Stecher and Tamura2016). Tree topologies were generated using Maximum likelihood (ML) and Neighbour joining (NJ) in the same software with 1000 bootstrap replicates. The MrBayes 3.2.0 (Ronquist et al., Reference Ronquist, Teslenko, van der Mark, Ayres, Darling, Höhna, Larget, Liu, Suchard and Huelsenbeck2012) was applied for the Bayesian Inference (BI), as implemented in NGPhylogeny.fr (Lemoine et al., Reference Lemoine, Correia, Lefort, Doppelt-Azeroual, Mareuil, Cohen-Boulakia and Gascuel2019), with two independent analyses of four concomitant Markov Chain Monte Carlo (MCMC) runs for 15 million generations and sampling parameters every 1000 generations. The first 25% of the trees were discarded as burn-in and a 50% majority-rule consensus tree was estimated. Intra- and interspecific pairwise genetic distances employing the Kimura 2-parameter distance model (K2P) (Kimura, Reference Kimura1980) were estimated using the MEGA 7.0 software. We used the Automatic Barcoding Gap Discovery (ABGD) species delimitation method (Puillandre et al., Reference Puillandre, Lambert, Brouillet and Achaz2012) to support the genetic validity of the species within our Gymnura dataset. This analysis was performed at https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.html with the P parameter ranging from 0.001–0.1 and a value of 2.0 for relative gap width.
A sub-dataset composed of 70 G. altavela and G. natalensis sequences (Gilchrist & Thompson 1911) (G. altavela Mediterranean N = 4; G. altavela North-west Atlantic USA N = 4; G. altavela South-west Atlantic Brazil N = 57; G. natalensis N = 5) was used to construct a median joining haplotype network (Bandelt et al., Reference Bandelt, Forster and Röhl1999) in PopArt (Leigh & Bryant, Reference Leigh and Bryant2015).
Results
A fragment of 559 COI gene base pairs was used for the phylogenetic analyses. All three methods (NJ, ML and BI) recovered similar topologies (Figures 1, S1, S2) and clearly indicated a subdivision between G. altavela from the West Atlantic (WA; Brazil + USA) and those from the East Atlantic (EA; Mediterranean). Both groups displayed a strict monophyletism, with 2.6% divergence. Furthermore, no haplotype is shared between both groups (Figure 3), suggesting two unique and divergent evolutionary mitochondrial lineages. Surprisingly, G. altavela from the Mediterranean and G. natalensis from South Africa were clustered together, exhibiting low divergence, indicating that these two nominal species belong to the same evolutionary lineage.
The analysis of intra- and interspecific K2P distance ranges also supports the divergence between G. altavela lineages from both sides of the Atlantic (Figure S3). Intraspecific genetic divergence levels for all Gymnura species were extremely low, ranging from 0.000–0.007, whereas values ranged from 0.019–0.265 between species. Two exceptions were noted, between G. altavela from the Mediterranean and G. natalensis (K2P distance = 0.003), and between G. crebipunctata and G. marmorata (K2P distance = 0.001) (Figure S3). The ABGD method also supported the splitting of G. altavela into two distinct lineages, ‘G. altavela’ from the West Atlantic (Brazil + USA) and ‘G. altavela’ from the East Atlantic (Mediterranean + G. natalensis in South Africa) (Figure 2).
Maximum likelihood tree for the Gymnura taxa, comprising two groups (WA lineage in grey and EA lineage in black) as supported by ABGD method.
The phylogeographic analysis limited to G. altavela and G. natalensis sequences indicated a low level of haplotype variation (Table S1) and low divergences between North-west Atlantic (i.e. EUA) and South-west Atlantic (i.e. Brazil) G. altavela (K2P distance = 0.001) specimens and between East Atlantic specimens (Mediterranean) and G. natalensis in South Africa (K2p distance = 0.001). On the contrary, the median joining haplotype network revealed a deep divergence between the two groups of haplotypes from both sides of the Atlantic (Figure 3).
Median joining haplotype network from the COI gene, indicating two distinct groups.
Discussion
The present study indicates the existence of two genetically differentiated Gymnura altavela lineages on both sides of the Atlantic Ocean. Both lineages exhibited strict monophyletism, with a genetic distance of 2.6%, a higher divergence than the intrapopulational divergence of any species belonging to the Gymnura genus, and within the range found for pairwise comparisons between other Gymnura species (0.019–0.265; Figure S3). In addition, no haplotype was shared between West and East Atlantic specimens and the two groups were also recognized by an ABGD analysis. Given the context of dramatic G. altavela population size declines and the recent reclassification of the species from vulnerable to threatened by the IUCN, further and deeper taxonomic investigations carried out using standard taxonomic analyses combined with nuclear marker analyses on samples from the entire geographic distribution of the species are required to infer species boundaries more securely within the genus.
Another finding highly relevant towards conservation efforts is that G. natalensis, to date, considered an endemic Southern Africa species (Figure 1), is phylogenetically closely related to G. altavela from the Mediterranean Sea. The possible synonymy of G. natalensis and G. altavela has already been suggested based on previous morphological analyses (Yokota et al., Reference Yokota, White, de Carvalho, Last, White, de Carvalho, Séret, Stehmann and Naylor2016) and, if future data support this hypothesis, G. natalensis should be reclassified as a junior synonym of G. altavela, consequently extending its distribution from the Mediterranean to Mozambique and Madagascar.
Linnaeus (Reference Linnaeus1758) first described the spiny butterfly ray from the Mediterranean Sea as Raja altavela. The species was later transferred to the Gymnura genus. The rays currently identified as G. altavela in the South-west Atlantic were originally identified as Gymnura binotata, described by Lunel (1879, as Pteroplatea binotata) from juvenile samples from Rio de Janeiro, the same area where our samples were collected. Interestingly, the main difference used to discriminate between G. altavela and G. binotata was ‘two white blotches on dorsal surface’. In Rio de Janeiro, we observed two G. altavela morphotypes, one classic (Figure S4A) and the other very rare (Figure S4B and Video S2), with only one specimen captured during 15 years of sampling efforts, with two white spots on the dorsum, which possibly gave rise to the description of G. binotata as a distinct species. More recently, we video recorded one individual with spots at Armação dos Búzios, in the state of Rio de Janeiro (22°46′06″S 41°47′26″W) (Video S2). Theses blotches are probably simply a polymorphism within the species, similar to the one also found in G. bimaculata (Norman, 1925). This species was differentiated from G. japonica based on the presence of a pair of white ocelli at the posterior part of the spiracule. Half a century later, Isouchi (Reference Isouchi1977) studied white ocelli pattern variations in G. bimaculata and described that one female with ocelli gave birth to pups with no ocelli, comprising strong evidence that the ocelli are a mere polymorphism. Alternatively, these blotches could be the sign of interspecific hybridization, and, thus, very rare and not detectable by a mitochondrial gene, or the result of reproductive secondary contacts overcoming species boundaries, which has been already recorded in batoids and sharks (Portnoy et al., Reference Portnoy, McDowell, Heist, Musick and Graves2010).
Recently, two new Gymnura species were distinguished (Yokota & de Carvalho, Reference Yokota and de Carvalho2017) within specimens formerly assigned as G. micrura. Alongside the evidence presented here, this underlines the high taxonomic instability of spiny butterfly rays. Therefore, morphological and molecular taxonomic studies are urgently required to help elucidate Gymnura species boundaries worldwide.
Consequences for conservation
Although the findings presented herein should be viewed with caution due to the use of a single mitochondrial gene, they are significant, as they indicate that a vulnerable and overfished species may be even more threatened than previously noted, since global population size estimates may be inflated by aggregating data from two different lineages that could, in fact, comprise distinct species. If confirmed by further morphological and genetic studies, the implications of this study will likely result in the uplisting of Gymnura altavela and, furthermore, indicate that the south-western Atlantic taxon is under even greater threat. These preliminary results should result in a greater engagement of the scientific community interested in fish conservation, which may be crucial for conservation efforts towards this important species.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S002531542200056X.
Data
Sequences used here can be retrieved from GenBank. Accession numbers are listed in Table S2.
Acknowledgements
The authors thank Rebeca A. Marques for aiding in laboratory procedures and two anonymous reviewers for valuable manuscript contributions. The authors thank the National Council for Scientific and Technological Development (CNPq) and Seed grants from IMAM-AquaRio. Thanks are also due to Eduardo Lukezic and the ‘Azul Profundo’ dive operator for ceding the video reported herein.
Author contributions
A.V.: contributed to data acquisition and analysis, interpreting the findings and writing the manuscript. F.L.: contributed to data acquisition and writing the manuscript. A.M.S.C.: contributed providing financial support, helping in data analysis, interpreting the findings and writing the manuscript. M.V.: contributed to formulating research questions, providing financial support, interpreting the findings and writing the manuscript.
Financial support
National Council for Scientific and Technological Development (CNPq) and Seed grants from IMAM-AquaRio.
Conflict of interest
The authors declare no conflict of interest.
