Barcode Reveals Hidden Diversity and Cryptic Speciation among Butterfly Rays Distributed in the Americas

: The taxonomic status of butterfly rays within the genus Gymnura remains a subject of ongoing debate among researchers. Some authors recognize up to five valid species for the Americas, while others considered several to be synonyms, which has posed a persistent challenge. We aimed to shed light on this complexity by employing molecular operational taxonomic units (MOTUs) based on the mitochondrial gene cytochrome oxidase I (COI). Genetic sequences were obtained from fresh muscle tissue collected in the marine ecoregions corresponding to the type locality from all the nominal butterfly ray species distributed along the Eastern Tropical Pacific (ETP). Our results unveiled compelling findings; all the species delimitation models used consistently identified seven MOTUs for the American continent and an extra G. altavela MOTU restricted to Africa. In addition, our results and models exceeded the worldwide accepted interspecific threshold of 2.0%. Remarkably, our results support the taxonomic reinstatement of Gymnura afuerae (Hildebrand, 1946) as a valid species, with a range expanding into the ETP in the Southern Hemisphere. Similarly, our data support the recent suggestion of resurrecting Gymnura valenciennii (Dum é ril, 1865) as a valid species in the western Atlantic. These findings urge a reassessment of the conservation status and a comprehensive taxonomic revision of American butterfly rays.


Introduction
The family Gymnuridae encompasses a single genus, Gymnura van Hasselt 1823, commonly known as butterfly rays [1].These distinctive rays are characterized by their rhomboidal discs, typically more than 1.5 times wider than it is long, and they possess a slender, whip-like tail, which may or may not bear serrated spines [2,3].With a maximum disc width (Dw) of up to 2.6 m, these benthopelagic batomorphs predominantly inhabit the tropical continental shelves of the Atlantic, Indian, and Pacific oceans [1,2,4,5].Despite being part of a monogeneric family, the precise number of valid Gymnura species remains debated among taxonomists, accounting for up to 16 butterfly rays worldwide.For example, the authors of [6] recognized 16 species, the authors of [1] listed 14, the authors of [7] included 12, the author of [4] identified 11, and the authors of [3] acknowledged only 10 distinct species.Recently, the authors of [8], based on molecular markers, asserted that almost half of the Gymnura species worldwide are unknown, indicating a subestimation of its current richness.The taxonomic challenge of butterfly rays is further complicated by their remarkable phenotypic similarity, often accompanied by substantial intraspecific variation and sexual dimorphism [3,9,10], adding complexity to their accurate taxonomic Taxonomy 2024, 4 identification and resulting in cryptic speciation [8].Historically, 10 nominal species of butterfly rays have been described from the Americas [11], but most were synonymized.Using new mitochondrial cytochrome oxidase subunit I (COI) sequence data, a widely used gene also known as the barcode gene for species delimitation, we aim to explore its potential for understanding the alpha diversity and evolutionary history of the Gymnura in the Neotropical region.
In the Americas, traditionally, two butterfly ray species, Gymnura altavela (Linnaeus, 1758) [12] and Gymnura micrura (Bloch & Schneider, 1801) [13], have been reported along the Atlantic coast [11,14,15].Meanwhile, at least three species, namely Gymnura marmorata (Cooper, 1863) [16], Gymnura crebripunctata (Peters, 1869) [17], and Gymnura afuerae (Hildebrand, 1946) [18], have been mentioned to occur across the eastern Pacific coast by several authors [2,[19][20][21][22][23].However, other authors currently consider the latter species to be a synonym of G. crebripunctata [3,4,11].Interestingly, several of the authors mentioned above noted and quoted the presence of butterfly ray complex species.Notably, recent efforts spanning the last 15 years have been dedicated to unravelling these complex species within the Neotropical region.For instance, on the Atlantic coast, the authors of [10] identified two new species, Gymunura lessae Yokota & Carvalho, 2017 [10] and Gymnura sereti Yokota & Carvalho, 2017 [10], within the G. micrura complex.Consequently, the current consensus recognizes G. lessae, G. micrura, and G. altavela distributed in the western Atlantic, with G. sereti restricted to the coast of Africa [10].In a separate study, the authors of [9] utilized a combination of morphometrics and molecular techniques, specifically the mitochondrial cytochrome b gene, to distinguish the California butterfly ray G. marmorata and the Longsnout butterfly ray G. crebripunctata as separate species in the northeastern Pacific.Notably, in this work, the authors failed to determine the taxonomic status of G. afuerae due to an insufficient sample size of the Peruvian butterfly ray entity described in [18] in 1946.
Questions regarding the Neotropical butterfly rays' richness and taxonomic validity persist (despite the recent effort mentioned above), particularly for those across the Eastern Tropical Pacific (ETP).For example, G. afuerae had a geographic distribution restricted to the Southern Hemisphere of the ETP, particularly in Ecuador and Peru, which is the type locality of the latter (Isla de Lobos; see [18,24,25]).In addition, this butterfly ray entity stands as one of the most contentious taxa within this genus on this continent even today, where its unresolved taxonomic status directly impacts the overall richness and conservation assessment of this group across their geographic range.In this sense, various studies have presented conflicting perspectives on G. afuerae, with some authors considering it a senior synonym [6, 9,22,[24][25][26] while others suggest it as a junior synonym of G. marmorata [19,27] or G. crebripunctata [3,4,11].The limited representation or absence of G. afuerae when focusing on this ray group underscores the imperative need for a comprehensive reevaluation of its taxonomic status by leveraging newly available techniques.
A similar case has been hypothesized for the Spiny butterfly ray G. altavela, where its presence and distribution on both Atlantic coasts have been questioned [4,8,10].The considerable morphological and phenotypic resemblance among the American butterfly rays, along with the population-level data gaps and the ongoing taxonomic and systematic controversy, may have led to mistakes and misidentification of these batomorphs in the published sources reverberating through the full extent of their ranges and geographic distributions.
Molecular analyses based on the mitochondrial gene COI have proven instrumental in resolving taxonomic uncertainties among Neotropical American endemic batomorphs or with a limited number of preserved individuals [28,29], similar to the challenges addressed in our study.Analysis using the COI has facilitated the differentiation of phenotypically similar [30,31] or taxonomically understudied species [29,[32][33][34], shedding light on discrepancies in sequences found in online databases.Furthermore, in the Americas, barcode analysis based on COI data has played a pivotal role in confirming synonymization of two or more elasmobranch species [35,36] or cryptic speciation of morphologically similar species [8,31,[36][37][38][39][40].Understanding the biological diversity and richness within the Ameri-can continent is crucial for formulating effective conservation and management strategies at both the continental and national levels, which is a primary research priority concerning elasmobranchs in Latin America [41].In this sense, our study aims to elucidate the taxonomic status of Neotropical butterfly rays, particularly in the ETP region, by leveraging barcode species delimitation data and models.

Materials and Methods
The muscle tissue samples were collected from eight individuals representing all nominal species distributed along the ETP region (Figure 1).These individuals included four G. afuerae (accession codes: PP756688-PP756691), two G. crebripunctata (accession codes: PP756686 and PP756687), and two G. marmorata (accession codes: PP756684 and PP756685).All specimens were obtained from artisanal fishing activities.Collection locations were chosen to be geographically close to the type locality of each nominal species (i.e., Isla de Lobos, Peru; Mazatlan, Mexico; and San Diego, USA, respectively) and within the same marine bioregion, as described in [42].These specimens were deposited into the ichthyological collections at the Centro Interdisciplinario de Ciencias Marinas del Instituto Politécnico Nacional in México (http://coleccion.cicimar.ipn.mx/,accessed on 12 February 2021) and the Instituto del Mar Perú sede Tumbes in Peru.
synonymization of two or more elasmobranch species [35,36] or cryptic speciation of morphologically similar species [8,31,[36][37][38][39][40].Understanding the biological diversity and richness within the American continent is crucial for formulating effective conservation and management strategies at both the continental and national levels, which is a primary research priority concerning elasmobranchs in Latin America [41].In this sense, our study aims to elucidate the taxonomic status of Neotropical butterfly rays, particularly in the ETP region, by leveraging barcode species delimitation data and models.

Materials and Methods
The muscle tissue samples were collected from eight individuals representing all nominal species distributed along the ETP region (Figure 1).These individuals included four G. afuerae (accession codes: PP756688-PP756691), two G. crebripunctata (accession codes: PP756686 and PP756687), and two G. marmorata (accession codes: PP756684 and PP756685).All specimens were obtained from artisanal fishing activities.Collection locations were chosen to be geographically close to the type locality of each nominal species (i.e., Isla de Lobos, Peru; Mazatlan, Mexico; and San Diego, USA, respectively) and within the same marine bioregion, as described in [42].These specimens were deposited into the ichthyological collections at the Centro Interdisciplinario de Ciencias Marinas del Instituto Politécnico Nacional in México (http://coleccion.cicimar.ipn.mx/,accessed on 12 February 2021) and the Instituto del Mar Perú sede Tumbes in Peru.The total DNA was extracted using a DNeasy Blood and Tissue kit (QIAGEN, Hilden, Germany) ® following the manufacturer's instructions.A segment of the COI gene was amplified using the universal FISH-F1 and FISH-R1 primers [43].The PCR mixture was prepared in a final volume of 25 µL, including 5 µL BSA (0.22%), 5 µL Taq buffer (5x), 0.5 µL dNTPs, 1.2 µL (10 mM) of each primer, 0.125 µL Taq polymerase (5 U/µL), 8.3 µL H2O, and 1.5 µL of the DNA sample.The PCR cycling conditions were as follows: initial denaturation at 94 • C for 120 s, followed by 35 cycles of denaturation at 94 • C for 30 s, annealing at 54 • C for 30 s, and extension at 72 • C for 60 s.A final extension step was carried out at 72 • C for 10 min.The amplified products were sequenced in both directions in South Korea (Macrogen, Inc.Seoul, South Korea).
Sequence Analysis: The obtained sequences were visualized and edited using Chro-masPro 1.3.3(Technelysium Pty Ltd., South Brisbane, Queensland, Australia).Our analysis included sequences retrieved from the Boldsystem (https://www.boldsystems.org/index.php accessed on 5 March 2024) and GenBank (https://www.ncbi.nlm.nih.gov/genbank/accessed on 8 March 2024) databases.These sequences were sourced from butterfly ray species geographically distributed along the coasts of the American continent.Specifically, we reviewed 37 Gymnura sequences from the Neotropical continent (30 from the Western Atlantic and 7 from the Eastern Pacific; Supplemental Table S1).Four sequences from G. altavela collected in the Eastern Atlantic were also included since this species has been mentioned to have a distribution on both Atlantic coasts (Supplemental Table S1).In addition, the Malaysian sequence of Gymnura poecilura (Shaw, 1804) [44] (MW498634) was included as the farthest congeneric species.Finally, the Urotrygon rogersi (Jordan & Starks, 1895) [45] OL604516 sequence was used as an outgroup.
Phylogenetic tree estimation: All sequences were aligned using MEGA X v10 software [46] through the MUSCLE algorithm and were subsequently checked for quality in terms of length and consistency, with sequences from other individuals included in the analysis.To determine our dataset's most suitable nucleotide substitution model, we utilized the Akaike criterion and Bayesian information, estimated within MEGA X v10 software.The HKG+I nucleotide substitution model was the best fit for the data.The phylogenetic tree topologies were constructed using the maximum likelihood and neighbor-joining methods [47] with the Kimura two-parameter model [48].We conducted 1000 bootstrap replicates, using MEGA X v10 software to assess the tree's robustness.A genetic distance matrix was computed among the assessed sequences using the parameters described above for each phylogenetic tree.
For the Bayesian inference tree, we employed MrBayes 3.2.0software [49].The parameters were configured with a Markov chain Monte Carlo (MCMC) chain length of 10,000,000 generations.We conducted this analysis with an optimized relaxed clock and HKG+I substitution model, where tracelog and treelog sampling was established every 1000 generations.The quality of the posterior estimates related to the effective sample size (after a 20% burn-in configuration) of the resulting MCMC chain was examined, and values exceeding the rule of thumb of 200 were confirmed for the parameters with MCMC Tracer Analysis Tool Version 1.7.2 software [50].Subsequently, we employed the standard configuration (10% burn-in, posterior probability limit = 0.0) within TreeAnnotator 2.7.1 software to generate the maximum clade credibility tree with node mean heights.The final tree was observed using the software FigTree version 1.4.4 [51].
Species delimitation: This was achieved by identifying the molecular operational taxonomic units (MOTUs) using six distinct models.These models were (1) automatic barcode gap discovery (ABGD; https://bioinfo.mnhn.fr/abi/public/abgd/abgdweb.htmlaccessed on 23 April 2024), (2) assemble species by automatic partitioning (ASAP; https: //bioinfo.mnhn.fr/abi/public/asap/accessed on 25 April 2024), which is similar to ABGD and was recently proposed by the authors of [52], (3) the poison tree process (PTP) [53] (https://species.h-its.org/ptp/accessed on 28 April 2024), (4) the maximum likelihood poison tree process (mPTP), and (4) the Bayesian inference poison tree process (bPTP).The first two models relied on genetic distances between the haplotypes to identify the threshold between interspecific (speciation) and intraspecific (populational) processes within the dataset [52,54].We used the following pre-established parameter values for these models: model = K2P, Pmin = 0.001, Pmax = 0.1, steps = 10, and number of bins = 20.We used the Kimura two-parameter model for the ASAP model, and the first two ASAP models with the lowest score were used as suggested by the authors.
As the poison tree process speciation delimitation model is based on a probabilistic framework, a phylogenetic tree was required as input to simulate the speciation and coalescent processes.This identified the nodes where the tree structure changed due to speciation and coalescent events (called the "transition points") in the phylogenetic tree.This model employed an MCMC approach to estimate the probability of these transition points (number of MOTUs within the tree).Finally, the optimal clustering of each individual or sequence was assigned based on maximizing the likelihood of the observed data given in the model [53,55].The bPTP is similar to the mPTP but relies on a Bayesian inference approach for speciation delimitation.For our butterfly ray case, a maximum likelihood tree was constructed using the best nucleotide substitution model (HKG+I), chosen based on the Akaike information criterion (AIC) in MEGA X v10.This was followed by adjustment of the parameters for this model, including the MCMC of 100,000 generations, to ensure convergence.Finally, the reliability of Bayesian support was confirmed by assessing the convergence equilibrium distribution [55].

Results
In our comprehensive analysis, which encompassed 51 COI sequences, we consistently observed a congruent and similar phylogenetic tree generated through two different methods: the Kimura two parameters and Bayesian approaches (Figure 2).The tree's robustness, along with the support for each node, was evident through high bootstrap and posterior probability values (Figure 2).
Our findings revealed an intriguing and compelling latitudinal gradient within various western Atlantic Gymnura species.For example, within G. micrura, two distinct sister clades emerged-one comprising the Mexican sequences and the other encompassing Brazilian sequences-revealing a pronounced geographical pattern (Figure 2).A similar latitudinal effect was observed for G. altavela (Figure 2).Interestingly, the American G. altavela sequences formed a cohesive clade, with the African G. altavela sequences as their sister clade (Figure 2).These findings underscore the significance of integrating both geographical and genetic factors in the study of butterfly rays.
In the context of the eastern Pacific species, our results contributed valuable insights for confidently assigning online sequences labeled as Gymnura sp. from Peru.Specifically, these sequences clustered within the same clade as our sequences of G. afuerae (Figure 2), affirming their accurate identification at the species level.Additionally, we identified a previously misidentified sequence, Gymnura marmorata PHANT702-08, which distinctly grouped with our sequences of G. crebripunctata (Figure 2).

Genetic Distances
The average net genetic distances among the Neotropical MOTUs assessed in this study from the trees described above ranged from 0.00% to 23.6% (Table 1).A lack of genetic distance was observed between the three groups, namely (1) Gymnura sp. and our G. afuerae, (2) G. marmorata (MX)1 and our group G. crebripunctata, and (3) among the three G. altavela species from the eastern Atlantic region, including Mauritania, Tunisia, and the Mediterranean Sea (Table 1).Notably, the MOTU used as an outgroup, G. poecilura, consistently exhibited the highest genetic distances, ranging from 17.2% (G.altavela ME) to 23.6% (our G. crebripunctata) (Table 1).

Genetic Distances
The average net genetic distances among the Neotropical MOTUs assessed in this study from the trees described above ranged from 0.00% to 23.6% (Table 1).A lack of genetic distance was observed between the three groups, namely (1) Gymnura sp. and our G. afuerae, (2) G. marmorata (MX)1 and our group G. crebripunctata, and (3) among the three G. altavela species from the eastern Atlantic region, including Mauritania, Tunisia, and the Mediterranean Sea (Table 1).Notably, the MOTU used as an outgroup, G. poecilura, consistently exhibited the highest genetic distances, ranging from 17.2% (G.altavela ME) to 23.6% (our G. crebripunctata) (Table 1).
The models based on genetic distance were similar and congruent, estimating an interspecific threshold value ranging between 1.5% and 2.2%.For instance, the ABGD species delimitation model consistently showed an initial partition of nine groups with a barcode gap distance of 0.015.However, the nine recursive partitions estimated between 14 and 6 groups, with the prior maximal distance ranging from 0.001 to 0.060, respectively.Interestingly, partition numbers four, five, and six consistently found nine groups, showing the highest group mode among all nine partitions calculated by this model (Figure 2).Congruently, the first two ASAP model partitions showed a threshold ranging between 0.014 (ASAP score = 3.0, p-value = 7.49 × 10 −3 , and W = 6.46 × 10 −4 ), resulting in nine subsets (MOTUs) (Figure 2).Consistently, for the ASAP2 model (ASAP score = 5.50, p-value = 1.16 × 10 −1 , and W = 4.58 × 10 −4 ), the number of subsets estimated was nine (Figure 2), with a threshold distance value of 0.022.Table 1.Net evolutionary divergence estimates and standard errors.This table presents net evolutionary divergence estimates (below the diagonal) and their corresponding standard errors (above the diagonal) for the distinct groups assessed, categorized by country of origin of the sampled butterfly ray collections.This analysis was performed using 50 nucleotide sequences and employed the Kimura two-parameter model.Country codes used are as follows: BR = Brazil, MT = Mauritania, ME = Mediterranean, MX = Mexico, PE = Peru, TU = Tunisia, and US = the United States of America.The group labeled G. marmorata (MX)1 represents an online sequence, while G. marmorata (MX)2 constitutes the new sequences generated in this study.

Species Delimitation Models and Grouping
Regardless of their algorithms, all species delimitation models consistently revealed a shared outcome, identifying nine distinct MOTUs among the butterfly rays assessed (Figure 2).Intriguingly, seven of these MOTUs formed clades for the American continent.At the same time, one cluster was specific to the African continent sequences of G. altavela.The remaining MOTU represented the outgroup sequences (Figure 2).Surprisingly, all species delimitation models indicated cryptic speciation for both coasts of the Neotropic region, revealing four and three MOTUs for the eastern and western coasts, respectively.For instance, in the western Atlantic region, two closely related yet distinct sister clades were identified for the eastern USA coast and the Gulf of Mexico sequences of the current known species G. lessae (Figure 2).Additionally, G. micrura from Brazil formed a distinct clade (Figure 2).Remarkably, all western Atlantic sequences of the spiny butterfly ray G. altavela, irrespective of their collection site (USA or Brazil), constituted a single clade distinct from the G. altavela from the eastern Atlantic region (Figure 2).Along the eastern Pacific coast, three distinct MOTUs were identified and consistently supported by all species delimitation models.Among these, G. crebripunctata and G. afuerae were closely related, while the California butterfly ray G. marmorata sequences from the USA and Mexico formed a single clade (Figure 2).

Discussion
Our results emphasize the substantial sequence dataset (n = 37) generated for the barcode gene available for Neotropical butterfly rays.However, it is noteworthy that these sequences exclusively originate from individuals collected in four countries across the American continent (i.e., USA, Mexico, Brazil, and Peru).It is notable that available American butterfly ray sequences generally represent the species' most extreme geographic distribution, creating genetic gaps and limiting insights along their known distribution.Remarkably, the number of sequences previously available for the butterfly ray entities distributed along the TEP was seven, and they were duplicated here, increasing the information available for G. afuerae and G. crebripunctata by 200% and that for G. marmorata by 50%.Our findings support the existence of cryptic speciation in Neotropical Gymnura, aligning with previously published research based on different genetic markers [8,9,56].As a result, the knowledge regarding accepting five Neotropical Gymnura as valid species was underestimated, and it needs to be updated.For instance, across the western Atlantic region, it seems to inhabit four Gymnura MOTUs, including one which still awaits a formal description, and for the eastern Pacific region, three butterfly ray entities are supported, including the taxonomic resurrection of a historically controversial Gymnura species.However, numerous biological and geographic aspects remain unknown for these butterfly rays, particularly for the newly identified MOTUs on both American coasts.Therefore, as butterfly rays are exploited as a marine resource along their distribution, we strongly encourage the need for additional studies to address their current population dynamics and conservation status.
Genetic distance, often crucial for species delimitation, has conventionally relied on a fixed threshold, commonly set around 2.0% for interspecific differences in barcode-based analyses [43,57].However, this approach may be unsuitable as the threshold varies across taxa.For example, elasmobranchs exhibit a slower COI gene mutation rate than bony fish (or actinopterygian fish) and other high-level taxonomic groups [58].Consequently, studies on sharks and rays, including Carcharhinus, Mobula, and Squalus, reported barcode genetic differences of around 1.0% between distinct and valid species [35,40,[58][59][60][61]. Overall, our COI gene genetic distances estimated among the assessed butterfly ray groups exceed both the minimal threshold values established for elasmobranch and the worldwide standardized 2.0% for other higher taxa [43,57].These findings support that the MOTUs estimated here are confident for the species-level taxon.As a corollary, our results suggest potential new taxonomic entities congruent with the previous results [8,56].Furthermore, taxonomic misidentification of COI sequences available in public databases has been mentioned for several American batomorph groups, such as Aetobatidae [33], Myliobatidae [31], as well as Urotrygonidae [28,29,40], which were also found for the butterfly rays in the present study (PHANT702-08), highlighting the need to be cautious when taxonomic determination is exclusively based on genetic sequences.This emphasize the importance of using up-to-date information on the target study group, including detailed taxonomic determination based on the external and morphometric traits of the individuals being studied.In addition, storing specimens collected during research in ichthyological collections for future consultation and review is always an excellent practice.
Four valid butterfly ray species are recognized in the western Atlantic region [8,10,62].However, our findings suggest that the Florida land extension and Greater Antilles islands could serve as a physical barrier, restricting the genetic flow between the northeastern USA and the Gulf of Mexico populations of the recently described G. lessae to the extent that they represent distinct and independent lineages.Considering G. lessae was described from a North Carolina specimen by the authors of [10] in 2017, our results raise the possibility that the Gulf of Mexico lineage warrants a formal description as a new species.A similar situation occurred with G. altavela, initially described from a Mediterranean Sea specimen by Carolus Linnaeus in 1758 [12].Numerous authors have discussed the currently accepted presence of G. altavela on both Atlantic coasts [56,63], aligning with our results, which support the uncertainties surrounding the taxonomy of G. altavela and echo the hypothesis in [8] suggesting that Gymnura valenciennii (Duméril, 1865) [64] should be taxonomically resurrected as the only butterfly ray with a serrated spine in the tail distributed along the western Atlantic coast while G. altavela is confined to the eastern Atlantic coast.
On the eastern Pacific coast, all species delimitation models identified the same three distinctive MOTUs (Figure 2).The first MOTU identified was G. marmorata, the only species exhibiting a serrated tail spine, which is consistent with previous studies [3,9].The distribution of this species encompasses the area from California to Peru [65], although its presence needs to be confirmed in several countries along its distribution.The second clade belongs to G. crebripunctata from the Northern Hemisphere, which formed a distinct and independent MOTU from its Southern Hemisphere sister clade, including both public sequences labeled Gymnura sp. and our genetic data generated for the nominal species G. afuerae.Our results affirm the previous hypothesis regarding the close relationship between these two entities or morphotypes [3][4][5]9].It seems that the oceanographic phenomena occurring in the equatorial area, such as gyres, current and countercurrent systems, including the equatorial undercurrent, the primary southern subsurface countercurrent, and the farther secondary southern subsurface countercurrent [66,67], may act as a genetic flow barrier to these two closely related butterfly ray species.In addition, the lack of genetic variation, as well as the collection location (Peru) of both online sequences labeled as Gymnura sp.(MH194465 and MH194463) and our G. afuerae barcode data forming into a single and independent MOTU from G. crebripunctata, support the taxonomic resurrection of G. afuerae as a senior synonym.Remarkably, after 77 years, our results corroborate Hildebrand's caution when describing G. afuerae, stating "It seems advisable at present to consider the Peruvian butterfly ray as specifically distinct from the Mexican" [18] (p.73).
Our barcode results shed light on the debated richness of butterfly rays in the Neotropics, a longstanding subject discussed by various authors [3,4,9,10], particularly in the ETP.Beyond its impact on species richness, our findings have broader implications for the conservation, distribution, and management of Neotropical butterfly rays, which rank high on the priority list of elasmobranch conservation in Latin American countries [41].Furthermore, our study underscores the necessity of a morphological taxonomic review within this genus in this area, particularly along the ETP.Identifying useful criteria for morphological species differentiation becomes crucial, especially without relying on characters linked to sexual dimorphism [40] and particularly in species which sympatrically coexist along the coasts of the Americas.Finally, considering our results, we strongly advocate for a reassessment of the conservation status of American butterfly rays.These insights contribute to the region's scientific understanding and practical conservation efforts.

Conclusions
Barcode gene and species delimitation models help understand a cryptic speciation.We suggest the need for the taxonomic resurrection of Gymnura afuerae (Hildebrand, 1946) [18] and Gymnura valenciennii (Duméril, 1865) [64].Our findings emphasize the need for a morphological taxonomic review within the American Gymnura species.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/taxonomy4030027/s1,Table S1: Identification code and accession link of the butterfly rays (Gymnura) Cytochrome Oxidase subunit I (COI) sequences available in the two online databases.

Figure 1 .Figure 1 .
Figure 1.Geopolitical map of a section of the American continent and the marine ecoregions of the Eastern Pacific according to [42], indicating the collection locations of the sequences by species (different forms) obtained in this study.The grey-filled figures indicate the sequences obtained in the present investigation, the black-filled figures indicate the sequences available from online databases, Figure 1.Geopolitical map of a section of the American continent and the marine ecoregions of the Eastern Pacific according to [42], indicating the collection locations of the sequences by species (different forms) obtained in this study.The grey-filled figures indicate the sequences obtained in the present investigation, the black-filled figures indicate the sequences available from online databases, and the blank figures indicate the type localities of the entities evaluated in this study.The same geometric figure may include more than one sample or nearby areas.Finally, this figure also shows the lack of genetic information throughout the known distribution of butterfly rays in the Neotropics (e.g., south from the Cortezian to the Guayaquil marine ecoregion and the Caribbean Sea in the Western Atlantic).

Figure 2 .
Figure 2. Phylogenetic analysis of the Gymnura species documented throughout the American continent using the cytochrome oxidase subunit I (COI) mitochondrial gene.Bootstrap support values (>95%) and Bayesian posterior probabilities are indicated above and below each branch, respectively, providing confidence in the tree topology.Vertical black bars on the right illustrate the results from the species delimitation models evaluated, highlighting distinct molecular operational taxonomic units (MOTUs).The new sequences generated in this study are denoted in a black font.Abbreviations: BR = Brazil, E.P. = eastern Pacific, E.A. = eastern Atlantic, I.O.= Indian Ocean, MA = Malaysia, MT = Mauritania, ME = Mediterranean, MX = Mexico, PE = Peru, US = the United States of America, TU = Turkey, and W.A. = western Atlantic.Clusters of butterfly ray MOTUs with tailserrated spines are indicated with asterisks (*).

Figure 2 .
Figure 2. Phylogenetic analysis of the Gymnura species documented throughout the American continent using the cytochrome oxidase subunit I (COI) mitochondrial gene.Bootstrap support values (>95%) and Bayesian posterior probabilities are indicated above and below each branch, respectively, providing confidence in the tree topology.Vertical black bars on the right illustrate the results from the species delimitation models evaluated, highlighting distinct molecular operational taxonomic units (MOTUs).The new sequences generated in this study are denoted in a black font.Abbreviations: BR = Brazil, E.P. = eastern Pacific, E.A. = eastern Atlantic, I.O.= Indian Ocean, MA = Malaysia, MT = Mauritania, ME = Mediterranean, MX = Mexico, PE = Peru, US = the United States of America, TU = Turkey, and W.A. = western Atlantic.Clusters of butterfly ray MOTUs with tail-serrated spines are indicated with asterisks (*).