Species Delimitation of Southeast Paciﬁc Angel Sharks ( Squatina spp.) Reveals Hidden Diversity through DNA Barcoding

: Angel sharks are distributed worldwide in tropical to subtropical waters. Across the Eastern Paciﬁc Ocean (EPO), two valid species are reported: The Paciﬁc angelshark Squatina californica and the Chilean angelshark Squatina armata ; however, there is still uncertainty about their geographic distribution, mainly along the northern Peru coast where the species have been reported to be sympatric. The aim of this study is to describe the genetic differences between the genus Squatina from the EPO, including samples from northern Peru, and using DNA barcoding and three species delimitation models: Poisson tree processes (PTP) model, Bayesian implementation of the PTP (bPTP) model and the general mixed Yule coalescent (GMYC) model. The three approaches summarized 19 nominal Squatina species in 23 consensus Molecular Taxonomic Units (MOTU). Only 16 of them were in accordance with taxonomic identiﬁcations. From the EPO, four Squatina MOTUs were identiﬁed, one from North America ( S. californica USA/Mexico) and three sampled in northern Peru, S. californica Peru, S. armata and Squatina sp. (a potential new species). This study contributes to the management and conservation policies of angel sharks in Peru, suggesting the presence of an undescribed species inhabiting the northern Peruvian coast. The use of molecular approaches, such as DNA barcoding, has the potential to quickly ﬂag undescribed species in poorly studied regions, including the Southeast Paciﬁc, within groups of ecologically and economically important groups like angel sharks.


Introduction
One of the most diverse groups of marine predators is the subclass Elasmobranchii (i.e., sharks, skates, and rays). Among them, there is a small but highly distinctive group of bizarrely-shaped benthic sharks, commonly called angel sharks. This group of ray-like sharks belongs to the monophyletic family Squatinidae (Bonaparte, 1838) [1][2][3] encompassing a unique genus, Squatina (Dumeril, 1805), with 22 extant species described based on morphological characters or molecular information [4][5][6][7]. Although, there are two additional species from the Gulf of Mexico described [8], Squatina mexicana Castro-Aguirre, Espinosa-Pérez and Huidobro-Campos, 2007, and Squatina heteroptera Castro-Aguirre, Espinosa-is considered as an operational definition that groups individual DNA sequences (i.e., cluster of sequences) based on an explicit algorithm and is used to estimate diversity at the species level. A powerful tool that uses standardized gene regions (e.g., cytochrome c oxidase-COI) to delimit and identify species is DNA barcoding. This molecular tool couples a comprehensive dataset of COI DNA sequences together with validated identified voucher specimens to support taxonomic studies [35,36]. Among their benefits, DNA barcoding can assist in defining phylogenetic relationships and species geographic boundaries. Furthermore, DNA barcoding together with single locus species delimitation methods have recently shown to be an effective tool for validating elasmobranch identification and describing new species in a number of fisheries [37][38][39][40][41][42][43].
Stelbrink et al. [4] employed COI and 16S rRNA markers to provide a comprehensive global phylogeny of 17 Squatina species, including the two species found along the EPO, S. armata and S. californica. In that study, S. armata is the first species to branch off the clade that includes the North and South American species, indicating the existence of two different species. However, in that study, only samples from the ends of both distributions, in California and Chile, were used. Thus, the detailed distribution of both species along the EPO remains uncertain, especially in the areas of North of South America where both species have been reported [14,15,18,44]. Additionally, angel sharks are target species and are captured by small-scale fisheries [28]. However, reports of landings are rather general and likely include several species of angel sharks, thus molecular tools could serve well for the accurate identification of the species landed, their distributions, and the fisheries with which they interact [28].
In this regard, the aim of the present study is to describe the genetic differences within the genus Squatina from the EPO, including samples of angel sharks from an area not previously covered (i.e., northern Peru), and integrate them with mtDNA COI sequences data from Barcode of Life Data (BOLD) System to evaluate species boundaries using species delimitation methods and determine MOTUs.

Morphological Identification and Sample Collection
As part of the sampling campaign of the consortium for DNA barcoding of Peruvian marine species (PeMar), a total of eight specimens of angel sharks were sampled between 2017 and 2018 from fish markets and landing sites in northern Peru (Figures 1 and 2). All specimens collected were identified using external morphological characteristics following a literature review [1,10,23]. Each specimen was photographed, and one specimen was fixed and deposited in the fish collection of the San Marcos-Natural History Museum (UNMSM) for further analyses. Muscle tissue samples were extracted from all specimens and preserved in 96% ethanol at room temperature (17-20 • C). Sequences, sample records, and voucher numbers can be viewed in Table 1.

DNA Extraction, Amplification and Sequencing
Total genomic DNA was isolated using Tissue Kit (Thermo Scientific) following the manufacturer's instructions. A partial fragment (~655 base-pair) of the mitochondrial Cytochrome Oxidase subunit I (COI) gene was amplified through Polymerase Chain Reaction (PCR) using primers FishF1-FishR1 or FishF2-FishR2 [35], that amplify an overlapping region from the 5' region of the COI gene. The PCR was performed with a final volume of 25 µL containing 16.35 µL distilled water, 2.5 µL dNTP (8 nM), 0.6 µL of each primer (5 µM), using just one pair of primers (i.e., F1/R1 or F2/R2) and 0.6 µL of Taq polymerase (5 µ/µL). PCR conditions consisted of an initial denaturation at 95 • C for 2 min, followed by 30 cycles including denaturation at 95 • C for 45 s, annealing at 52 • C for 45 s, and extension at 72 • C for 60 s, and a final extension at 72 • C for 5 min. Amplified products were checked on 1% agarose gel and both strands per amplicon were sent to Macrogen (Rockville, MD, USA) for Sanger sequencing. The sequencing was carried out using the same set of primers that was used in the PCR, however a greater number of samples were amplified and sequenced using the Fish F1 and Fish R1 primers, since these had better efficiency for our samples. Sequences were cleaned and contigs were assembled using the software CodonCode Aligner (CodonCode Corporation, Dedlham, MA, USA). Multiple alignments were done using a ClustalW algorithm [45], implemented in the software MEGA 7 [46] and were checked manually for misalignments and trimmed to the shortest common sequence length.
We performed and compared three molecular species delimitation methods using the pipeline SPdel (https://github.com/jolobito/SPdel, accessed on 8 March 2021) added to a delimitation species method based on classical taxonomy. The first method implemented was the general mixed Yule coalescent (GMYC) model [57] with a single threshold. The GMYC model identifies the threshold value at the transition of branching patterns that are characteristics of the speciation process [58] versus coalescence process [59], and identifies significant changes in the branching rates in a time-calibrated ultrametric tree. The other methods used were the Poisson tree processes (PTP), and the Bayesian implementation of the PTP model (bPTP) [60]. These two methods directly use the number of substitutions as opposed to the GMYC method that uses time. For all models, an ultrametric tree was used as an input file which was generated in BEAUti v2.1 [61], with a normal relaxed clock and a birth and death model, and a HKY+G substitution model suggested by the Bayesian Information Criterion in jModeltest 2 [62]. The analysis was implemented with 60,000,000 million MCMC generations and with a burn-in of 10%. To run the analysis, we used BEAST v2.5 [61] implemented in the Cyber Infrastructure for phylogenetic Research (CIPRES; https://www.phylo.org, accessed on 8 March 2021) [63]. Convergence and adequate sample size (greater than 130) were evaluated in Tracer v1.6.0. The pipeline SPdel used for our analysis performs a comparison between the four chosen methods and ultimately generates a consensus species delimitation considering three molecular approaches. Only MOTUs supported by at least two of the applied delimitation methods were considered consensus MOTUs. Furthermore, we quantified the degree of genetic divergences for nominal species and for each one of the species delimitation methods used, calculating the values of intragroup and intergroup Kimura 2-parameter (K2-P) genetic distances in SPdel.

Taxonomic Identification
Most of our sampled specimens were not retained because they were part of fishers' daily catch. For this reason, identification was done mainly in the field (i.e., with fresh specimens) and through photographs allowing for diagnostic taxonomic characters to be assessed. The four specimens collected at the ECOMPHISA fish market and close to the Bayovar Port were identified as S. armata (Table 1, Figure 1). For the identification we used an illustrated guide [10] and a taxonomic key [15], distinguishing the following combination of characteristics: reddish-brown to blackish dorsal surface (Figure 2a), white ventral surface with dark brown edged of pectoral and pelvic fins (Figure 2b), narrow and simple barbels (Figure 2c), anterior nasal flaps fringed (Figure 2d), and thorns on snout, between eyes and spiracles (Figure 2e). One angel shark pup was collected off the coast of Punta Pico (Table 1, Figure 1). It was identified as S. californica following the aforementioned illustrated guide, distinguishing the following characteristics: reddishbrown dorsal surface with scattered light spots (Figure 2k), white edged pectoral and pelvic fins (Figure 2l), presence of thorns (Figure 2m,n), and pale dorsal fins (Figure 2o). Finally, three specimens collected in Mancora (Table 1, Figure 1) were identified as S. californica in the field, however after the molecular analysis, they were allocated to the taxon Squatina sp. (Figure 2f-j).

MOTU Delimitation Analyses
We obtained eight sequences (610 base-pair) from two nominal species, S. armata (n = 4), S. californica (n = 1), and Squatina sp. (n = 3), from northern Peru with 30 parsimony informative sites. The final alignment of mtDNA COI sequences resulted in 591 bp (shortest common sequence length) comprised of a total of 19 nominal species (Result S1).
The species delimitation analyses showed that PTP and the bPTP method delimited the same 23 Squatina MOTUs (Figure 3), with a maximum intra-MOTU divergence of 0.99% (for the MOTU of S. squatina) and a minimum inter-MOTU divergence of 0.49% (between the MOTUs of S. californica collected in Peru and USA/Mex). On the other hand, the GMYC method delimited Squatina 25 MOTUs with a maximum intra-MOTU divergence of 0.99% (for the MOTU of S. squatina) and a minimum inter-MOTU divergence of 0.16% (between the MOTUs of S. squatina collected in Ireland and Turkey). Single-locus species delimitation results from PTP, bPTP, and GMYC approaches were summarized by using the pipeline SPdel, in 25 consensus MOTUs. The maximum intra-MOTU distance was 0.99% (for the MOTU of S. squatina) ( Table 3) and the minimum inter-MOTU distance was 0.49% (between the MOTUs of S. californica collected in Peru and USA/Mexico) ( Table 3). In contrast, if considering species delimited through traditional taxonomy, the maximum intraspecific distance was 2.51% (between specimens morphologically identified as S. africana) ( Table 3) and the minimum interspecific distance was 0 between specimens of S. formosa (collected in Thailand) and S. nebulosa (collected in China and Southern China Sea) ( Table 3).    Table 3). The S. californica samples formed a polyphyletic group, with two MOTUs closely related, one from the Northeast Pacific (samples collected within the Gulf of California, Mexico, and in California, USA), and the other composed by samples collected in northern Peru (Figure 3). The minimum intra-MOTU divergence calculated for S. californica from Northeast Pacific was 0% and the maximum intra-MOTU value was 0.83%. Furthermore, the minimum and the maximum inter-MOTU divergences between the MOTU of S. californica from the northern hemisphere and the southern hemisphere were 0.5 and 0.83%, respectively. Samples identified initially as S. californica (Pemar_V0209, Pemar_V0210, and Pemar_V0211) from northern Peru, were grouped into a clade, separately from S. californica and S. armata and for this reason they were renamed as Squatina sp. The minimum and maximum intra-MOTU divergences for Squatina sp. were 0 and 0.33%, respectively.

Discussion
Only two valid species had been previously described for the EPO: Squatina californica and Squatina armata [7,10,28], and both species have been reported as sympatric for the northern Peruvian coast [15]. Nonetheless, there is controversy about the validity of these species due to the lack of taxonomic studies confirming their presence across the whole range of their geographic extents [10]. Our results show a new scenario, reporting the existence of four MOTUs along the EPO (Figure 3), revealing a hidden diversity that may include at least one new species for the genus Squatina in the Southeast Pacific.
The Bayesian phylogenetic tree obtained shows a similar topology compared with the comprehensive phylogenetic analysis carried out by Stelbrink et al. [4], and all species delimitation analyses performed (i.e., PTP, bPTP, and GMYC) yielded mostly the same result. The four MOTUs found in the EPO were grouped within the clade of North and South American species described by Stelbrink et al. [4]. However, the four MOTUs are not phylogenetically closely related with the exception of the MOTU of S. californica from North America and the MOTU of S. californica from northern Peru that are sister MOTUs (Figure 3).
All the algorithms used split the cluster of S. californica into two MOTUs. One group includes specimens collected along the Pacific coast of the United States and specimens from the Gulf of California, the other group includes sequences from Punta Pico located in the region of Tumbes, Peru ( Figure 3). In some previous studies, two different population lineages originated in the Northern Pacific and in the Gulf of California were found [4,13,64], but this division was not supported by the coalescent species delimitation methods herein used, instead of reciprocal monophyly. Nonetheless, the minimum K2-P genetic distances within the MOTUs of S. californica from the northern hemisphere (0.5%) are comparable to the lowest genetic distances found between Squatina nominal species (0.5% for S. guggenheim and S. occulta). In regard to S. californica from the southern hemisphere, it was discriminated as a different MOTU from the clade of the northern hemisphere, even though the maximum K2-P genetic distances between both clades is 0.5%, similar to the values observed between the population from North America and other Squatina nominal species. This significant heterogeneity between S. californica populations may be promoted by their ecological behavior (e.g., limited ability for sustained swimming due to their morphological and anatomical characteristics) [13], their reproductive behavior (e.g., philopatric behavior) [28] but also due to geographic barriers (e.g., deep marine basins as barriers to dispersal of these populations) [13]. To elucidate if the degree of separation of these MOTUs of S. californica corresponds to a recent divergent species or a strong population structure, more studies including specimens of different ontogenetic stages collected along the Eastern Pacific, their morphological, anatomical, and ecological traits and nuclear markers to evaluate the degree of separation of these MOTUs, assessing introgression and gene flow, are necessary.
A third EPO MOTU corresponds to samples matching the original S. armata description (e.g., the presence of narrow and simple barbels and anterior nasal flaps fringed, and thorns on the head). This MOTU clade is the first to branch off from the American clade and includes samples from Chile and Peru. The last EPO MOTU corresponds to Squatina sp. from northern Peru, a new lineage that represents a potential undescribed new species. The MOTU was related to the clade formed by S. dumeril and the two MOTUs of S. californica, but with a low posterior probability support and their exact evolutionary relationships remain unknown ( Figure 3). Additionally, the nearest species to Squatina sp. using K2-P distance was S. guggenheim, with a higher value (1.5%) compared to other nominal Squatina species. An integrative taxonomic study to evaluate if this MOTU corresponds to a new species is imperative. In addition, since there are no taxonomic studies providing a direct comparison of both EPO species, besides the two original descriptions, a revision is necessary of S. armata and S. californica to clarify characteristics and dichotomous identification keys that can help with field identification along the distribution.
The exact distributions of S. californica and S. armata or even of the four MOTUs remain uncertain. Published studies reporting the presence of angel sharks off the coasts of South America are based on literature reviews [16,17] or governmental reports that do not clearly differentiate between both species, reporting angel sharks as Squatina spp. [65]. Similarly, previous distributions reports could be erroneous due to the presence of the cryptic species herein described or by the difficulties during species assignment due to the lack of clear morphological taxonomic characters to distinguish these species.
The limited movement, site fidelity, and preference for coastal areas, as well as other characteristics, such as their large body size and low reproductive output, make angel sharks susceptible to overexploitation by fisheries [28]. In Peru, both S. californica and S. armata have been reported as caught or landed by small scale fisheries by the Institute of Peruvian Sea (Instituto del Mar del Peru in Spanish; IMARPE) [19,65,66]. Before 1996, landings of angel sharks were reported under their common name in Spanish "angelotes" [65]. Nonetheless, some fishing expedition reports, catalogues and identification guides mentioned only the presence of S. armata [19,66,67].
From 1996 to 2010, all fishery landings of angel sharks along the coast were reported by IMARPE as S. californica [65]. These reports showed a marked decline over the period 1996 to 2010 [65]. Currently, the fishery continues and there is still a lack of detailed landing reports of angel sharks in Peru [68]. Managing angel sharks in groups (e.g., genus level) can mask population declines and can represent an impediment to fisheries research but also it may hamper the national enforcement regulation for conservation and management [69]. Landings information along the Peruvian coast (2010-2019) is still reported as Squatina sp. by the Ministry of Production (Ministerio de Producción in Spanish; PRODUCE, Nro Registro 00090925-2020) indicating, at least for the northern region of Peru (i.e., Tumbes, Piura, Lambayeque, and La Libertad regions) a decline of caught specimens reported in tons. Due to the complicated taxonomy of angel sharks added to a potential new, undescribed species, landings reports might be underestimating the population decline of all the species distributed in Peruvian waters. As well as in Peru, in the United States and Mexico, a decline has also been observed [28,[70][71][72], for instance, reported landings from US fisheries showed a large decline influenced by fisheries regulation (e.g., minimum landing size and ban on gillnets and trammel nets) [28].
In our species delimitation analysis, that included 19 from the 22 nominal Squatina species, some differences were observed between MOTUs and nominal data that deserve further discussion. All species delimitations split S. africana into four MOTUs. A recent study which included sequences from South Africa and Indian Ocean specimens, found a high number of haplotypes within S. africana [73] using DNA fragments of COI gene, and with a high K2-P genetic distances (up to 2.5%), when compared with values reported in this study between other pairs of nominal Squatina species (e.g., 0.99% between specimens of S. albipunctata and S. pseudocellata). Similar to S. californica, the split of S. africana into 4 MOTUs may be explained due to the life-traits of angel sharks and the geographical barriers. Along the South African coast several barriers to gene flow have been detected, which often coincide with patterns of ocean currents [74]. For example, the South-East and East coasts include three marine ecoregions, the Agulhas Ecoregion, Natal Ecoregion, and part of the Delagoa Ecoregion. Biogeographic forcing agents defining these ecoregions include bathymetric or coastal complexity and currents [75]. This is evident along the east side of South Africa, where a very narrow extension of the continental shelf and wide, deep areas are observed between the northern part of the Agulhas Ecoregion and the Southern part of the Natal Ecoregion. The South-East and East coasts appear to have several population genetic boundaries for several rocky shore and estuarine species (e.g., the fishes Clinus cottoides, Clinus superciliosus, and Muraenoclinus dorsalis) [74,76]. Although previous studies showed a genetic break in these areas, the effect of biogeographic barriers on demersal shark species is still largely unknown and needs further investigation [77]. Finally, the MOTU S. africana 1 was collected in the Indian Ocean close to the Maldives, far from the sampling area of MOTUs S. africana 3 and S. africana 4 which is characterized by different oceanographic and geomorphologic features.
Samples of S. nebulosa (sampling location in Japan and Thailand) and S. formosa (sampling location in China and Sea of South China) from BOLD System were clustered in the same MOTU sharing the same haplotype. Diagnosing both species is particularly challenging, their subtle external morphological differences may cause an overlap among many characteristics used [31]. Due to the difficulties in the correct identification of these species, the source from where we retrieved the sequences (i.e., GenBank which lacks physical voucher information) and the fact the taxonomy for these individuals were not revised, our result needs to be tested in an integrative taxonomic study including morphology and genetic analyses.
Finally, our results have far-reaching implications for the management and conservation policies of angel sharks in EPO, suggesting the presence of an undescribed species inhabiting the northern Peruvian coast, living in sympatry with the two species already reported. Many species of angel sharks are affected by artisanal and industrial fisheries and their populations have declined around the world [28]. All efforts to quantify the impact of fisheries in South American angel sharks will be unsuccessful without a precise identification of the species. The use of molecular approaches, such as DNA barcoding, has a potential to quickly identify undescribed species in poorly studied regions like the Southeast Pacific [78]; and in ecologically and economically important groups like Elasmobranchii that hold a high level of taxonomic uncertainties [78].
Supplementary Materials: The following are available online at https://www.mdpi.com/1424-2 818/13/5/177/s1, Table S1: Complete list of COI sequences of angel sharks retrieved from BOLD System and used in this study, including collection location and access numbers Bold System and GenBank Result S1: Final matrix containing all the DNA sequences worked on this paper.