Phylogenomics and species delimitation for effective conservation of manta and devil rays

Practical biodiversity conservation relies on delineation of biologically meaningful units. Manta and devil rays (Mobulidae) are threatened worldwide, yet morphological similarities and a succession of recent taxonomic changes impede the development of an effective conservation strategy. Here, we generate genome‐wide single nucleotide polymorphism (SNP) data from a geographically and taxonomically representative set of manta and devil ray samples to reconstruct phylogenetic relationships and evaluate species boundaries under the general lineage concept. We show that nominal species units supported by alternative data sources constitute independently evolving lineages, and find robust evidence for a putative new species of manta ray in the Gulf of Mexico. Additionally, we uncover substantial incomplete lineage sorting indicating that rapid speciation together with standing variation in ancestral populations has driven phylogenetic uncertainty within Mobulidae. Finally, we detect cryptic diversity in geographically distinct populations, demonstrating that management below the species level may be warranted in certain species. Overall, our study provides a framework for molecular genetic species delimitation that is relevant to wide‐ranging taxa of conservation concern, and highlights the potential for genomic data to support effective management, conservation and law enforcement strategies.


| INTRODUCTION
Biodiversity conservation relies on the delimitation of species units, particularly with respect to global conventions and regulatory frameworks. Advances in sequencing technologies now make it possible to uncover fine-scale patterns of genetic variation with unprecedented resolution (Helyar et al., 2012;Morin et al., 2004). However, there is concern this could lead to over-splitting of species, and therefore constrain management options while artificially increasing extinction events (Frankham et al., 2012;Isaac et al., 2004;Mace, 2004;Zachos, 2013;Zachos et al., 2013). On the other hand, under-splitting species has the potential to leave novel evolutionary genetic lineages unrecognized (Gippoliti et al., 2018;Morrison et al., 2009).
When using genomic data, it has therefore been recommended that a robust discovery and validation framework is applied, employing multiple methods (Carstens et al., 2013) and taking additional lines of evidence and expert opinion into account (Coates et al., 2018;Dayrat, 2005;Padial et al., 2010;Stanton et al., 2019). Such an approach provides an objective yet conservative evidence base to differentiate within and among species of conservation concern.
Alongside recent developments in genome sequencing, analytical methods in phylogenetics have also advanced (Fujita et al., 2012).
Nevertheless, multispecies coalescent approaches have a tendency to delimit population structure over speciation (Sukumaran & Knowles, 2017) and when used in isolation, a distinction between these two scenarios is often not possible (Chambers & Hillis, 2020).
This runs the risk of assigning species to the smallest diagnosable unit and is particularly problematic both for taxa where gene flow is restricted by environmental or geographical barriers (Sukumaran & Knowles, 2017), and for taxa with high levels of gene flow but for which limited geographical sampling has taken place (Chambers & Hillis, 2020;Pante et al., 2015). Under the general lineage concept, species represent independently evolving lineages that remain intact when in contact with close relatives (de Queiroz, 1998(de Queiroz, , 2007. Accordingly, when it is possible to sample from regions where individuals come into geographical contact, a distinction between independently evolving lineages and population-level variation can be made (Chambers & Hillis, 2020;Leaché et al., 2019). Analysing samples from regions such as these, together with consideration of similarities and a succession of recent taxonomic changes impede the development of an effective conservation strategy. Here, we generate genome-wide single nucleotide polymorphism (SNP) data from a geographically and taxonomically representative set of manta and devil ray samples to reconstruct phylogenetic relationships and evaluate species boundaries under the general lineage concept. We show that nominal species units supported by alternative data sources constitute independently evolving lineages, and find robust evidence for a putative new species of manta ray in the Gulf of Mexico. Additionally, we uncover substantial incomplete lineage sorting indicating that rapid speciation together with standing variation in ancestral populations has driven phylogenetic uncertainty within Mobulidae. Finally, we detect cryptic diversity in geographically distinct populations, demonstrating that management below the species level may be warranted in certain species. Overall, our study provides a framework for molecular genetic species delimitation that is relevant to wide-ranging taxa of conservation concern, and highlights the potential for genomic data to support effective management, conservation and law enforcement strategies.
The importance of judiciously defined species and conservation units is particularly apparent in fisheries management, where overexploitation has led to population declines worldwide (Agnew et al., 2009;Burgess et al., 2013;Hutchings & Reynolds, 2004;Worm et al., 2009). One group of heavily targeted fishes are the manta and devil rays (Mobula spp.; family Mobulidae; collectively, mobulids). Despite substantial economic value through tourism (O'Malley et al., 2013), these highly mobile, circumglobally distributed megafauna are threatened by targeted and bycatch fishing often driven by a demand for their gill plates (Couturier et al., 2012;O'Malley et al., 2017). Mobulids have among the most conservative life history traits of all fishes, including low fecundity, long gestation periods and late maturation (Couturier et al., 2012;Dulvy et al., 2014), and as a result, current levels of exploitation are considered unsustainable (Croll et al., 2016 Adnet, et al., 2020;Stewart et al., 2018). To develop effective conservation measures, there is an urgent need for both species and intraspecies boundaries within mobulids to be resolved.
The first major taxonomic revision of the genus Mobula was undertaken in 1987, where a total of nine extant species were described (Notarbartolo di Sciara, 1987). Over two decades later, the genus Manta was revised, where a second species, Manta alfredi, was resurrected based on morphological and meristic data (Marshall et al., 2009). In 2018, a reclassification of the family Mobulidae was carried out using sequence data from mitogenomes and over 600 nuclear exons (White et al., 2018). White et al., (2018) proposed the genus Manta (manta rays; previously Manta alfredi and Manta birostris) to be subsumed into Mobula (devil rays), and for the total number of species to be reduced to eight, by recognizing M. japanica, M.
rochebrunei and M. eregoodootenkee as junior synonyms of M. mobular, M. hypostoma and M. kuhlii, respectively. Both the synonymization of M. japanica and M. mobular, and the monophyly of the genus Mobula are supported by additional studies (Adnet et al., 2012;Naylor et al., 2012;Paig-Tran et al., 2013), including an independent phylogenetic analysis using mitogenome data and two nuclear markers (Poortvliet et al., 2015). Nevertheless, as previous work has been geographically restricted and reliant on few samples, the relationships between other species groups remain problematic and patterns of intraspe- Furthermore, for over a decade, a third putative species of manta ray has been hypothesized to occur in the Atlantic and Caribbean (Hinojosa-Alvarez et al., 2016;Marshall et al., 2009). However, molecular support has remained inconclusive due to reliance on relatively few genetic markers displaying conflicting evolutionary signals (Hinojosa-Alvarez et al., 2016). Similar patterns of discordance have been observed between M. alfredi and M. birostris, which despite showing evidence of separation at nuclear markers, are indistinguishable based on mitochondrial DNA White et al., 2018). Through using phylogeographical analysis and coalescent models across multiple populations of manta rays, Kashiwagi et al. (2012) attributed these observations to recent speciation and post-divergence gene flow. To explore the extent of such phenomena across the entire genus, there is a need to apply similar coalescent-based approaches to additional species groups. Moreover, by employing genome-wide markers in a comprehensive set of samples, it should be possible to resolve problematic species boundaries and determine classifications that encapsulate the full extent of diversity across devil rays.
Here, we use double-digest restriction site-associated DNA (ddRAD) sequencing to generate a dataset of genome-wide markers in an extensive geographically and taxonomically representative set of mobulid individuals. The resulting data were used to: (a) reconstruct phylogenetic relationships within Mobulidae, (b) evaluate support for nominal species boundaries using coalescent-based approaches and (c) explore patterns of intra-specific variation across taxa. We examine our results in light of recent morphological and ecological evidence, and in doing so, demonstrate utility in delimiting units for conservation, detecting cryptic diversity, and improving our understanding of associated evolutionary processes in a global radiation of socio-economically important marine megafauna.

| Sampling and DNA extraction
Tissue samples were obtained from existing collections and sampling initiatives of researchers and organizations worldwide, yielding a total of 116 samples from all nominal species (Figure 1 and Table S1). These kuhlii cf. eregoodootenkee (Notarbartolo di Sciara, Adnet, et al., 2020). An additional five samples from Rhinoptera bonasus were obtained as an outgroup. Samples were identified to species based on characteristics described by Stevens et al. (2018), using species names valid at the time of collection. To distinguish between original species names and those currently recognized, we refer to the recently synonymized species M. japanica and M. rochebrunei using the abbreviated expression "cf." throughout the paper. Genomic DNA was extracted using the Qiagen DNeasy Blood and Tissue Kit and DNA yield was measured using a Qubit 3.0 Broad Range Assay. Extracts were quality assessed on 1% agarose gels stained with SafeView. The single sample of M. hypostoma cf. rochebrunei was from a museum specimen stored in formalin and yielded no detectable DNA.  (Asahida et al., 1993;Chang et al., 1995), together with our required sequencing coverage. Unique P1 and P2 barcode combinations were ligated to resulting DNA fragments, which were then size-selected between 400 and 700 bp using gel electrophoresis and then PCR (polymerase chain reaction)-amplified. A pilot ddRAD library was sequenced on an Illumina MiSeq at the Institute of Aquaculture, University of Stirling. Subsequent ddRAD libraries were sequenced by Edinburgh Genomics on an Illumina HiSeq High Output version 4 using the 2 × 125PE read module.

| Data quality control and filtering
Data was quality assessed using faStqc and processed in StackS version 1.46 (Catchen et al., 2011). The prOcESS _ raDtaGS.pl module in StackS was used to demultiplex the data, filter for adaptor sequences (allowing two mismatches), remove low-quality sequence reads (99% probability) and discard reads with any uncalled bases. To minimize linkage disequilibrium in the single nucleotide polymorphism (SNP) data, only forward reads were retained for subsequent analyses.
Short fragments not removed through size-selection were filtered out using a custom bash script (8.5% of reads).
The DEnOvOmap.pl program in StackS was then used to assemble loci and call SNPs. We selected the three main assembly parameter values that generated the largest number of new polymorphic loci shared across 80% of individuals, following Paris et al. (2017).
Four identical reads were required to build a stack (-m), stacks differing by up to four nucleotides were merged into putative loci (-M) and putative loci across individuals differing by up to five nucleotides were written to the catalogue (-n), giving an average locus coverage of 105X. To examine the effect of missing data on phylogenomic inference we used the pOpulatiOnS.pl program in StackS to generate two datasets with varying levels of missing data. The first contained SNPs present in at least 10 individuals (Dataset A), and the second contained SNPs present in at least 90 individuals (Dataset B). The former reflects a more relaxed filter resulting in a higher genotyping rate, while the latter represents a stricter threshold. To remove paralogous loci and mitigate for allelic dropout (Arnold et al., 2015;Gautier et al., 2013), loci sequenced at greater than twice or less than one-third the standard deviation of coverage, respectively, were identified and excluded from both datasets using vcftOOlS (Danecek et al., 2011). The remaining loci were then assessed for excess heterozygosity using vcftOOlS, and those exhibiting a significant probability of heterozygote excess were excluded. Next, because StackS ignores indels, SNPs in the last five nucleotide positions were assumed erroneous and excluded. Finally, the remaining loci in each dataset were written to a whitelist and filtered for a single random SNP per locus to minimize linkage using the pOpulatiOnS.pl program in StackS. This resulted in 7,926 SNPs with 47.1% missing data in Dataset A, and FIGURE 1 Geographic locations of samples used in this study. Samples were identified to species based on characteristics described by Stevens et al. (2018), using species names valid at the time of collection. Species are represented by coloured circles, scaled for number of samples. Total numbers of samples for each species is provided in the legend. Further details are provided in Table S1. A small amount of variation was applied to the location of each sampling site to avoid over-plotting and improve interpretation. [Colour figure can be viewed at wileyonlinelibrary.com] 1,762 SNPs with 14% missing data in Dataset B (Table S2). The sequences in which SNPs were observed were concatenated for each individual and used as input for phylogenetic analysis.

| Phylogenetic reconstruction and clustering
To infer relationships among individuals and examine the extent of monophyly, we first carried out maximum likelihood phylogenetic analysis using concatenated ddRAD Datasets A and B. A maximum likelihood phylogeny was estimated using raxml version 8.2.11 (Stamatakis, 2014) under the GTRGAMMA model of rate heterogeneity. This model was selected following an assessment of best fit using both Akaike's and Bayesian information criterion (AIC and BIC) in jmODEltESt2 (Darriba et al., 2012). Branch support was assessed with 1,000 bootstrap replicates. raxml identified four highly supported clades separated by long branches. We next calculated F ST values among these inferred clusters using the pOpulatiOnS.pl program in StackS. To further assess how individuals clustered within groups and investigate possible cases of population structure, we also performed a principal components analysis (PCA) using the R package aDEGEnEt (Jombart, 2008) on Datasets A and B. To improve visualization, this was carried out separately for each clade (Table S3). After assessment of 10 PCs, three were retained in all cases.

| Coalescent-based species delimitation
To objectively compare alternative species delimitation models against the taxonomy proposed by White et al. (2018), we used Bayes factor delimitation (BFD*, Leaché et al., 2014). This method allows for direct comparison of marginal likelihood estimates (MLEs) for alternative species delimitation hypotheses, or models, under the multispecies coalescent. BFD* was conducted using the modified version of Snapp (Bryant et al., 2012), implemented as a plug-in to bEaSt version 2.4.8 (Bouckaert et al., 2014). Due to computational constraints, we ran this analysis separately for each clade, as previously described, using Dataset B only. However, to evaluate support for interaction from higher phylogenetic levels, we included four random individuals from a sister species. Our null model matched the taxonomy proposed by White et al. (2018) and alternative models were informed by the literature and analyses herein.
Full descriptions of our species delimitation models are described in Tables S4-S7. For each clade, we also assessed a model in which individuals were randomly assigned to two or three species.
For all models, path sampling involved 10 steps (1,000,000 MCMC [Monte Carlo Markov chain] iterations, 20% burn-in), implementing the log-likelihood correction. Because MLEs are affected by improper prior distributions, a gamma distribution was implemented on the lambda (tree height) parameter. To assess the effect of priors on the ranking order of models, models were also assessed retaining the default 1/X distribution on lambda, implementing upper and lower bounds (10,000 and 0.00001 respectively), for a proper prior.
Bayes factors (2log e BF) were calculated from the MLE for each model for comparison (Kass & Raftery, 1995;Leaché et al., 2014), as follows: Positive 2log e BF values indicate support for the null model (>10 is decisive; Leaché et al., 2014) and negative values favour the tested model. The best supported species units were then used to build consensus sequences and generate maximum likelihood consensus trees for Datasets A and B. For this, we provided a population map to the pOpulatiOnS.pl program in StackS where each individual was assigned a species based on the results from our BFD* analysis and individual-based phylogeny. Phylogenetic analysis was then performed on the resulting consensus sequences using raxml (Stamatakis, 2014) as described above.

| Species tree analysis
To test tree topology and evaluate uncertainty due to incomplete lineage sorting, we estimated species trees using Snapp based on the species units best supported by both BFD* and our individualbased phylogeny (Bryant et al., 2012). This method allows each SNP to have its own history under the multispecies coalescent, while bypassing the need to sample individual gene trees. Due to the computational capacity required to run Snapp, three individuals per species were randomly selected from Dataset B while maximizing geographical coverage within species. Random sampling of individuals with replacement was repeated a further three times, resulting in four subsampled alignments (Table S8). MCMC chains consisted of 5,000,000 iterations, sampling every 1,000 and retaining default priors on lambda and theta for each independent analysis. Convergence to stationary distributions was observed after 20% burn-in in tracEr (Rambaut et al., 2018). The distribution of trees was visualized in DEnSitrEE version 2.2.6 (Bouckaert, 2010) and maximum clade credibility (MCC) trees were drawn using trEEannOtatOr version 2.4.7 (Bouckaert et al., 2014). Alternative prior combinations produced highly concordant results.
Multispecies coalescent-based approaches assume that any discordance of topologies among loci results from incomplete lineage sorting, and do not consider introgression as a source. We therefore used trEEmix on Dataset A to evaluate evidence for significant introgression events within the Mobulidae. This method examines the extent to which variation between user-defined groups is explained by a single bifurcating tree. Given the topological uncertainty identified using Snapp, specifically regarding the placement of M. mobular, the three-population test (Reich et al., 2009) was additionally used to test for "treeness" between clades. Like trEEmix, the three-population test estimates covariance of allele frequencies among groups 2log e BF=2 * MLE null − MLE test but is simpler and less parameterized and therefore potentially more powerful for identifying introgression. For this analysis, M. mobular, M. alfredi and M. thurstoni were randomly selected to represent their respective clades.

| Phylogenetic reconstruction and clustering
To infer relationships among individuals and examine the extent of monophyly and clustering, we carried out maximum likelihood phylogenetic analysis and PCA. Phylogenetic trees based on ddRAD Datasets A and B were highly congruent, with species groups forming well-supported clades separated by long branches indicating that missing data had little effect on our results (Figure 2; Figure S1).
PCAs based on Datasets A and B were also highly similar, and revealed clustering within each clade that mirrored patterns in both phylogenetic trees (Figure 3; Figure S2).

| Coalescent-based species delimitation
To compare alternative species delimitation models against the for further details). We therefore considered the single species models to be the most reasonable, in line with the results of our maximum likelihood phylogenetic analysis. In all cases, models assessing interaction from higher phylogenetic levels and random individual assignments performed comparatively poorly (Tables S4-S7). Furthermore, marginal likelihood estimates were unaffected by lambda priors, with no change in the rank order of models (Tables S4-S7).

| Species tree analysis
To examine relationships among species, we estimated species trees based on the best supported species groups from our maximum likelihood and BFD* analysis. Consensus species trees estimated under the multispecies coalescent exhibited broadly consistent topologies and theta estimates across independent runs, indicating that subsampling did not significantly bias the tree topology inferred with Snapp ( Figure S5). Therefore, the species tree for individual subsample one is presented here (Figure 4; see Figure S5 for all subsets). Relationships between sister species within clades remained particularly stable across alternative topologies. However, some uncertainty with respect to the placement of the putative new mobulid species was observed, although a sister relationship with M. birostris was the most common. This is in line with the maximum likelihood consensus trees produced using raxml, where the putative new species is sister to M. birostris, with 100% bootstrap support Figures S1 and S6). Theta estimates inferred from the Snapp analysis varied across lineages, where they tended to be larger for deeper branches of the tree ( Figure S7). trEEmix inferred an admixture graph similar to trees produced with raxml ( Figure S8), explaining 99.86% of the variance, indicating mobulid species placement is unaffected by admixture. Furthermore, we found no evidence of introgression between clades containing M. alfredi, M. mobular and M. thurstoni, through three-population tests (Table S9).

| DISCUSSION
Conservation management relies on the classification of diversity into species units. Phylogenetic inference and species delimitation are therefore of fundamental importance to biodiversity conservation in assigning strategies to the appropriate level of biological integrity. In this study, we combined ddRAD sequencing together with coalescent-based species delimitation to evaluate species boundaries and explore intraspecific patterns of variation among the globally threatened Mobulidae. Our study has broad-reaching implications for practical conservation and provides a framework for carrying out molecular genetic species delimitation in wide-ranging taxa of conservation concern.
The phylogenetic relationships and species boundaries that we infer in this study largely support those proposed in a recent revision of Mobulidae by White et al. (2018). However, we find compelling evidence to refute the synonymization of M. eregoodoo and M. kuhlii. Our findings coincide with morphological, ecological and behavioural evidence reported in a recent study by Notarbartolo di Sciara, Adnet, et al., (2020). Here, distinct differences in coloration, cephalic fin length, and tooth and branchial filter plate FIGURE 4 Tree clouds of sampled species trees from Snapp analysis to visualise the range of alternative topologies. The tree presented here is based on 30 individuals from subsample one and 1242 SNPs from Dataset B. Individuals were assigned to the ten species units supported by our analyses herein. See Figure S5 for tree clouds based on individual subsamples twofour and Table S8 for further details on individual subsets. Illustrations © Marc Dando, Wild Nature Press, and are not to scale. The illustration of the putative new mobulid species is based on images of hundreds of individuals off the Yucatan Peninsula, Gulf of Mexico [Colour figure can be viewed at wileyonlinelibrary.com] morphology were characterized using specimens collected globally. Our results, together with this morphological evidence, therefore strongly suggest these taxa should be considered distinct species under an integrative taxonomic framework. Mobula eregoodoo and M. kuhlii partially share a geographical range across a region with intense fishing pressure (Notarbartolo di Sciara, Fernando, Adnet, Cappetta, & Jabado, 2017). It is therefore critical that they are monitored and managed as separate units, to avoid both overstating abundance estimates and underestimating perceived extinction risk (Coates et al., 2018;Isaac et al., 2004;Mace, 2004).
In contrast to this species split between M. eregoodoo and M. kuhlii, our results support the synonymization of M. mobular and M. mobular cf. japanica as described by White et al. (2018). Such an assertion is also in line with a previous study of mitogenome and nuclear data of this species (Poortvliet et al., 2015), as well as morphological similarities that have been documented in these taxa for several decades (Notarbartolo di Sciara, 1987; Notarbartolo di Sciara, . However, through analysing a larger number of geographically representative samples and genome-wide molecular markers, our work addresses the constraints of previous studies, and therefore provides the most robust genetic evidence to date that M. mobular should be considered a single species unit for conservation. Unfortunately, we were unable to assess the taxonomic status of M. hypostoma cf. rochebrunei, as the museum sample we had access to yielded no detectable DNA. Given the extent of illegal, unreported and unregulated fishing pressure in the West African region (Agnew et al., 2009) combined with the scarcity of M. hypostoma cf. rochebrunei sightings, there is a risk that the synonymization outlined by White et al. (2018) could lead to an underestimate of extinction vulnerability. We therefore recommend that high priority be given to a comprehensive evaluation of mobulid diversity in West Africa (see Stewart et al., 2018).
Another important outcome of our study is support for a putative new species of manta ray in the Gulf of Mexico based on multilocus nuclear SNP markers. Importantly, samples from both M. birostris and the putative new species were collected from within sites in the Gulf of Mexico, suggesting that these independent lineages co-occur in sympatry. This work adds to the growing body of morphological, ecological and mitochondrial DNA evidence that this unique genetic lineage may constitute an undescribed species worthy of independent management. Nevertheless, although we did not formally test for hybridization among manta ray species, we identified one individual that could be considered genetically intermediate between M. birostris and the putative new species. This indicates that, as in M. birostris and M. alfredi (Walter et al., 2014), hybridization may be occurring. Management of these taxa as independent units will therefore be challenging, potentially requiring blanket protection in regions where sympatry and/or hybridization occur. To develop suitable protective measures, it is critical that a formal description of the putative new species is carried out. This will allow researchers to confidently determine its full range, which may extend into international waters or span areas with high fishing pressure.
Recent studies have demonstrated that caution must be applied when using the multispecies coalescent as it has a tendency to over-split taxa into the lowest possible denomination (Chambers & Hillis, 2020;Leaché et al., 2019;Sukumaran & Knowles, 2017). This is particularly true when factors including population structure, isolation by distance and sampling design are not taken into account (Chambers & Hillis, 2020). Consequently, the approach can make it difficult to disentangle species boundaries from population-level processes, presenting a genuine risk to conservation, where over-splitting can negatively impact assignment of management units and estimates of extinction vulnerability (Frankham et al., 2012;Isaac et al., 2004;Mace, 2004;Zachos, 2013;Zachos et al., 2013). Nevertheless, under the general lineage concept, there is little controversy in assigning species status to genetically distinct populations that occur in sympatry (Leaché et al., 2019;Rannala et al., 2020), such as the case here for M. birostris and the putative new mobulid species, and for M. kuhlii and M. eregoodoo. For the former, samples from both species were collected from the same site in the Gulf of Mexico, while for the latter, samples from both species were collected from within the same ~ 120km stretch of South African coastline. The patterns we uncover are therefore unlikely to be attributable to population structure, and we can be confident that these taxa reflect genetically distinct and separately evolving lineages.
For M. alfredi, M. kuhlii and M. mobular, models containing splits below the species level received the highest level of support in our comparison of species delimitation hypotheses. However, because these groupings reflected geographical separation in sample origin, attributing these patterns to species boundaries is less clear cut. In particular, it has been shown that limited geographical sampling in widely distributed species can give rise to the appearance of distinct lineages (Chambers & Hillis, 2020). This pattern can occur regardless of whether a species is made up of fragmented populations with limited gene flow or exhibits isolation by distance (Barley et al., 2018;Hedin et al., 2015). Although we analysed samples from multiple populations of M. alfredi, M. kuhlii and M. mobular, large areas of their geographic range remain under-represented. Together with the fact that support for models recognizing geographically separated species only marginally out-performed those that did not, we interpret such patterns as evidence for population-level processes as opposed to speciation events.
Nevertheless, these BFD* results, alongside the geographical clusters characterized in the PCA, indicate that isolated populations represent unique genetic diversity and that management below the species level may be warranted (Coates et al., 2018;Stanton et al., 2019). In line with this, the geographic signal we have uncovered lays the foundation for developing tools to support fisheries traceability and population monitoring (Ogden, 2008(Ogden, , 2011. For example, diagnostic markers with the power to discriminate between populations could be employed for determining the geographic origin of devil ray samples. Such traceability tools could be used to understand the impact of fisheries on individual populations or to identify geographic regions that may be routinely supplying the gill plate trade (Fields et al., 2020;O'Malley et al., 2017). The extent of population structuring described within M. alfredi, M. kuhlii and M. mobular indicates that future traceability studies will be feasible in these species. Interestingly, for M. thurstoni and M. birostris, both of which are circumglobally distributed, no evidence for population structure was found, despite samples originating from across ocean basins. This was in some respects surprising given the reproductive biology of mobulids and localized movement patterns that have been reported in certain species (Braun et al., 2015;Deakos et al., 2011;Jaine et al., 2014;Stewart et al., 2016). However, both M. thurstoni and M. birostris occupy more pelagic and offshore habitats than their coastal counterparts, and are hypothesized to exhibit highly migratory behaviour, presenting opportunities for gene flow across ocean basins. To understand the biological significance of these patterns in more detail, and contribute further information to support management, a more comprehensive assessment of population structure within Mobulidae using high-density markers is required.
We used the outcomes of our species delimitation and maximum likelihood phylogenetic analysis to provide insights into the evolutionary history of devil rays. A common problem when inferring species relationships is the impact of conflicting phylogenetic signal across the genome on gene-tree discordance (Degnan & Rosenberg, 2009;Schrempf & Szöllösi, 2020). Fortunately, the multispecies coalescent model provides a framework with which to incorporate topological uncertainty in phylogenetic reconstruction (Rannala et al., 2020). Through employing an extension of the multispecies coalescent developed for SNP datasets (Bryant et al., 2012), we find evidence for discordance among tree topologies in manta and devil rays. In particular, we find substantial uncertainty in the placement of M. mobular, with this species shifting position relative to the manta rays and the other species. This is in contrast to both our maximum likelihood phylogeny and consensus species trees, where M. mobular was sister to the manta rays with ≥95% bootstrap support. Nevertheless, when taken together, our results are in line with those of White et al. (2018) suggesting the previously recognized genus Manta is nested within Mobula, and therefore provide further justification for the associated change in nomenclature. We also un- The two primary biological causes of gene-tree discordance are introgression and incomplete lineage sorting (Hahn, 2018;Mallet et al., 2016). In the case of M. mobular, we find no evidence for admixture and therefore attribute the topological uncertainty to incomplete lineage sorting. This phenomenon is more likely to occur when the time between speciation events is small, because there is reduced opportunity for lineages to coalesce (Maddison, 1997;Pamilo & Nei, 1988). Mobulidae have undergone a recent and rapid burst of speciation (Poortvliet et al., 2015) that is likely to have increased the level of incomplete lineage sorting. Similarly, levels of lineage sorting are proportional to ancestral population size, given that the rate of coalescence is lower when effective population size (N e ) is large. Interestingly, we find that N e , as inferred from theta estimates, is larger for deeper branches of the mobulid tree. Thus, it is also possible that standing genetic variation in ancestral populations may have led to lower levels of lineage sorting in mobulid rays. Similar patterns have been identified across a wide range of taxa (Brawand et al., 2014;Pollard et al., 2006;Suh et al., 2015;Thompson et al., 2019;Wang et al., 2018) and our results highlight how the use of multilocus approaches within a coalescent framework can improve our understanding of the evolutionary processes shaping genomic variation.

| CONCLUSIONS
In summary, we combine genome-wide data together with multiple methods to characterize lineage diversity and examine species units in globally threatened manta and devil rays. We find that two recently synonymized species constitute independently evolving lineages and uncover robust support for a putative new species of manta ray in the Gulf of Mexico. We further characterized differences between geographically separate populations, highlighting the need for management below the species level. Our findings have important implications for practical conservation and are relevant to the enforcement of CITES and CMS regulations by laying the groundwork for species identification and regional traceability of products in trade. Furthermore, we demonstrate the power of genomic data to resolve and identify diversity within organismal radiations and improve our understanding of the evolutionary processes generating biodiversity. Consequently, we provide a framework for molecular genetic species delimitation of relevance to wide-ranging taxa of conservation concern, and endorse the value of genomic research to inform conservation, management and law enforcement when integrated with complementary biological data.