To be, or not to be: that is the hamletic question of cryptic evolution in the Eastern Atlantic and Mediterranean Raja miraletus species complex

Citation: Ferrari, A.; Crobe, V.; Cannas, R.; Leslie, R.W.; Serena, F.; Stagioni, M.; Costa, F.O.; Golani, D.; Hemida, F.; Zaera-Perez, D.; et al. To Be, or Not to Be: That Is the Hamletic Question of Cryptic Evolution in the Eastern Atlantic and Mediterranean Raja miraletus Species Complex. Animals 2023,13, 2139. https:// doi.org/10.3390/ani13132139 Academic Editor: James Albert Received: 31 May 2023 Revised: 25 June 2023 Accepted: 25 June 2023 Published: 28 June 2023 Copyright: © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). animals Article To Be, or Not to Be: That Is the Hamletic Question of Cryptic Evolution in the Eastern Atlantic and Mediterranean Raja miraletus Species Complex Alice Ferrari 1,† , Valentina Crobe 1,† , Rita Cannas 2, Rob W. Leslie 3, Fabrizio Serena 4, Marco Stagioni 5, Filipe O. Costa 6, Daniel Golani 7, Farid Hemida 8, Diana Zaera-Perez 9, Letizia Sion 10 , Pierluigi Carbonara 11 , Fabio Fiorentino 4,12 , Fausto Tinti 1,* and Alessia Cariani 1 1Department of Biological, Geological and Environmental Sciences, University of Bologna, 40126 Bologna, Italy; [email protected] (A.F.); valentina.cr[email protected] (V.C.); [email protected] (A.C.) 2 Department of Life and Environmental Sciences, University of Cagliari, 09126 Cagliari, Italy; r[email protected] 3Branch Fisheries Management, Department Agriculture, Forestry and Fisheries, Cape Town 8018, South Africa; r[email protected] 4 Institute for Biological Resources and Marine Biotechnology, National Research Council, 91026 Trapani, Italy; [email protected].it (F.S.); [email protected].it (F.F.) 5 Laboratory of Marine Biology and Fisheries, Department Biological, Geological and Environmental Sciences, University of Bologna, 61032 Fano, Italy; [email protected] 6Centre of Molecular and Environmental Biology (CBMA) and ARNET-Aquatic Research Network, Department of Biology, University of Minho, Campus de Gualtar, 4710-057 Braga, Portugal; [email protected] 7Department of Evolution, Systematics and Ecology, The Hebrew University of Jerusalem, Jerusalem 9190401, Israel; [email protected] 8 Ecole Nationale Supérieure des Sciences de la Mer et de l’Aménagement du Littoral, Campus Universitaire de Dely Ibrahim, Algiers 16320, Algeria; [email protected] 9Institute of Marine Research, 5817 Bergen, Norway; diana.zaera-per[email protected] 10 Department of Biosciences, Biotechnologies and Environment, University of Bari Aldo Moro, 70125 Bari, Italy; [email protected] 11 COISPA Technology and Research, 70126 Bari, Italy; [email protected] 12 Stazione Zoologica Anton Dohrn, 90149 Palermo, Italy *Correspondence: [email protected] † These authors contributed equally to this work. Simple Summary: The Raja miraletus species complex exhibits high levels of morphological and ecological stasis along with the antipodean distribution in the Eastern Atlantic and Indian Oceans. We investigated genetic variability and differentiation between taxa and geographical populations by integrating mitochondrial and nuclear DNA markers. The extraordinary occurrence of at least five different sibling taxa in the Northeastern Atlantic Ocean and Mediterranean Sea is documented, supporting cryptic speciation and stabilising selection. Abstract: Despite a high species diversity, skates (Rajiformes) exhibit remarkably conservative morphology and ecology. Limited trait variations occur within and between species, and cryptic species have been reported among sister and non-sister taxa, suggesting that species complexes may be subject to stabilising selection. Three sibling species are currently recognised in the Raja miraletus complex: (i) R. miraletus occurring along the Portuguese and Mediterranean coasts, (ii) R. parva in the Central-Eastern Atlantic off West Africa and (iii) R. ocellifera in the Western Indian Ocean off South Africa. In the present study, the genetic variation at mitochondrial and nuclear markers was estimated in the species complex by analysing 323 individuals sampled across most of its geographical distribution area to test the hypothesis that restricted gene flow and genetic divergence within species reflect known climate and bio-oceanographic discontinuities. Our results support previous morphological studies and confirm the known taxonomic boundaries of the three recognised species. In addition, we identified multiple weakly differentiated clades in the Northeastern Atlantic Ocean and Mediterranean, at least two additional cryptic taxa off Senegal and Angola, a pronounced Animals 2023,13, 2139. https://doi.org/10.3390/ani13132139 https://www.mdpi.com/journal/animals Animals 2023,13, 2139 2 of 19 differentiation of ancient South African clades. The hidden genetic structure presented here may represent a valuable support to species’ conservation action plans. Keywords: cartilaginous fish; brown skate; conservation biology; population genetics; mtDNA; microsatellite loci 1. Introduction The evolutionary debate on the nature of species boundaries [ 1 , 2 ] is based on paradigms such as Mayr’s discontinuous variation and reproductive isolation of species and Darwin’s continuity between varieties, geographical populations and species [ 2 ]. Natural hybrid zones and secondary contacts with gene introgression unequivocally show that species boundaries have a semipermeable nature [ 3 ] and that intrinsic barriers to gene flow (i.e., preand postzygotic barriers) are in some cases incomplete. Species life cycles, ecological features and adaptive phenotypes are particular key points influencing species distribution and dispersal within the marine environment [ 4 , 5 ], in which permanent and intermittent breaks (e.g., landmasses, unsuitable habitats, upwelling areas and oceanographic fronts) may isolate populations, enabling ecological differentiation [ 6 , 7 ], genetic divergence [ 8 ], reproductive isolation and speciation [ 9 – 11 ]. Molecular systematics and phylogenetics have greatly contributed to the assessment of relationships among taxa and effectively contributed to delineate the hierarchy of evolutionary frames in which recently diverged taxa exhibit, on average, lower divergence than taxa in the later stages of the speciation process [ 6 ]. Over the last two decades, increased evidence has emerged for speciation governed by entirely different mechanisms, leading to so-called sibling or cryptic species (sensu Bickford [ 12 ]). The idea that species can evolve into similar morphologies is well established [ 13 ], but the use of molecular delimitation methods has now brought cryptic species to the forefront in many research arenas [ 6 , 14 , 15 ]. Bickford et al. [ 12 ] identified at least two recurrent themes in animals wherein morphological distinctiveness and reproductive isolation are unpaired: in groups using non-visual mate-recognition signals (e.g., chemical, olfactory, acoustic, electric-field senses) and in groups living under environmental conditions that promote the stabilising selection of morphological traits (e.g., extreme habitats, specialised host–parasite relationships, deep-sea environments, fishery pressure). Among elasmobranchs, skates and rays exhibit highly effective modulation of electro-sensory signals depending on behaviour ([ 16 ] and references therein). At the same time, they display a marked conservation of ecological and morphological traits [ 17 – 23 ], especially between recently diverged species [ 12 ]; a strong evolutionary success in terms of resilience at the evolutionary scale [ 24 ] and a high degree of endemism [ 25 , 26 ] and species richness [ 27 – 29 ]. Investigations into the role of biogeographical barriers on the speciation of marine organisms have increasingly concentrated across several taxa [ 30 – 32 ]. Prior to 2016, the brown skate Raja miraletus Linnaeus, 1758 was thought to be distributed throughout the Mediterranean Sea and from northern Portugal, along the western and south-eastern coasts of Africa [ 33 , 34 ]. This distributional range is much wider than expected for a small-sized rajid, given the limited potential for dispersal in a species with a relatively sedentary behaviour of adults and juveniles [ 35 – 39 ] and the lack of egg dispersal [ 40 ]. Nominal R. miraletus exhibits a pronounced benthic ecology, with most records from 10 m to 150 m on sandy and hard bottoms [ 25 , 33 ] and a generalist feeding behaviour [ 41 , 42 ]. Due to its high and stable abundance over its distribution, small body size and early maturation (age at maturity estimated at 2.7 years; [ 43 ]), it was considered highly resilient to exploitation and was assessed globally as Least Concern in the International Union for Conservation of Nature Red List [ 44 , 45 ]. Since then, Raja ocellifera Regan 1906 has been resurrected for the southern African population [ 28 , 46 ] and was assessed as Endangered in 2020 [ 46 ]. The newly described R. parva Last and Séret 2016 from west Africa has not been assessed [28]. Animals 2023,13, 2139 3 of 19 The R. miraletus complex exhibits a high level of stasis of the general external morphology over its range; all populations exhibit a distinctive bright tricolored (blue, black and yellow) eyespot on the upper ochre-brownish surface at the base of each pectoral fin [ 34 , 47 ]. After the first identification of three parapatric or allopatric groups (Mediterranean, West Africa and South Africa) based on the variation of morphometric and meristic characters [ 48 ], preliminary evidence of cryptic speciation in R. miraletus was observed by integrating results from mitochondrial DNA analysis, morphology and host–parasite relationships from specimens collected in Central-Southern Africa [ 49 , 50 ]. Raja miraletus was then recognised as a species complex of at least four valid taxa based on combined data for COI and NADH2 [ 51 ]: (1) the northernmost R. miraletus, occurring in the Mediterranean Sea and adjacent North-Eastern Atlantic waters, (2) the southernmost and resurrected R. ocellifera (associated with mtDNA data from Naylor et al. [ 49 ]; Raja cf. miraletus 1, NCBI Accession Number JQ518895, [ 49 ]) occurring, in the Western Indian Ocean, off South Africa, and in the Indian Ocean, from False Bay to Durban, (3) the central African R. parva, distributed from Senegal to Angola (associated with mtDNA data from Naylor et al. [ 49 ]; Raja cf. miraletus 2, NCBI Accession Number JQ518890, [ 49 ]) and (4) a still undescribed taxon (Raja miraletus, NCBI Accession Number JQ518891, [ 49 ]), occurring from Mauritania to Senegal where it is therefore sympatric with R. parva. The advent of high-throughput DNA sequencing technologies and the launch of global DNA-based biodiversity assessments (e.g., DNA barcoding; [ 52 ]) has provided raw data, enabling the determination of taxonomic, ecological and evolutionary aspects of cryptic and sibling species, where the term “sibling” denotes a cryptic species with a recent common ancestor, implying a sister species relationship [ 53 ] and even more challenging conservation issues [ 12 ]. Moreover, molecular methods coupling markers obtained from mitochondrial DNA (mtDNA) and nuclear DNA (nuDNA) have improved the resolution of species boundaries and revealed gene introgression/hybridisation phenomena in marine fish including elasmobranchs [ 10 , 54 – 56 ]. Understanding the liaison between species life-history traits, ecology and the adaptive phenotypes leading to hidden population divergence and reproductive isolation is of utmost importance for skates, whose conservation is often hampered by the lack of species-specific data [57]. This study uses the integrative support of the mtDNA cytochrome oxidase subunit I barcode region (COI) and eight nuDNA EST-linked polymorphic microsatellite loci (ESTSSRs; [ 58 ]) to estimate the genetic variation among 323 specimens collected across almost the full geographic range of the R. miraletus species complex [ 40 , 47 ] and exhibiting the typical tricoloured eyespots. We tested hypotheses of the relationship between restricted gene flow and genetic divergence within the species complex, specifically in relation to climatic and oceanographic discontinuities. Additionally, we sought to establish parallel patterns between our findings and variations in morphology and parasite prevalence, which were independently assessed [ 48 , 50 ]. As compared to previous knowledge, our findings contributed to describe a richer scenario concerning the taxonomic units and zoogeographic boundaries characterising the R. miraletus complex. 2. Materials and Methods 2.1. Sampling A total of 323 specimens of from the R. miraletus species complex were collected between 2000 and 2014 (Tables 1and S1) during scientific research programs (South Africa, Angola and Mediterranean Sea) by contracted commercial fishermen (Senegal, Levantine Sea and Israel) or at local fish markets (Algeria). Scientific trawl surveys were carried out in South Africa (2006 and 2011, Africana cruises), Angola (2006, Nansen cruises), the whole Mediterranean Sea (2000–2014 Mediterranean International Trawl scientific Surveys, MedITS; [ 59 ]) and national scientific trawl surveys (2000–2010 Italian Gruppo Nazionale Demersali surveys, GruND; [ 60 ]; the 2007 Portuguese scientific surveys of the Instituto Português de Investigação do Mar) allowed for a comprehensive sampling, covering most of the wide geographical distribution of the R. miraletus complex (Figure 1). All individuals Animals 2023,13, 2139 4 of 19 were easily assigned to the R. miraletus complex based on their very distinctive morphotype and species-specific diagnostic characters [ 25 , 40 ]. Fin clips and muscle tissues were cut from each individual using sterile tweezers and clippers, transferred to a clean tube filled with 96% ethanol and stored at −20 ◦C for subsequent DNA analyses. Table 1. Sampling data and locations. “N (tot)”, refers to the total number of individuals sampled in this study, according to the methods indicated in the Source column. Of these “N (COI)” and “N (EST-SSRs)” refer to individuals COI-sequenced and genotyped in this study. N (tot) = 0 when available COI sequences were retrieved from public databases and integrated into the mtDNA dataset, as specified in the “Source” column. The last row refers to geographical samples previously described in McEachran et al. 1989 [ 48 ]. 1—Mediterranean group. 2—Mauritania and Senegal group. 3—Gulf of Guinea–equatorial African group. 4—Angolan sample. 5—South African sample. Sampling Area Macro Area Code N (tot) N (COI)N (EST-SSRs) Years Source (Trawl Survey Program) McEachran et al., 1989 [48] Central-Southern Africa (CSA) South Africa—South Coast SAF 8 6 8 2006 ST (Africana) 5 South Africa—South Coast 0 5 0 2007, 2012 GB/BOLD 5 South Africa—South Coast 32 30 31 2011 ST (Africana) 5 Namibia NAM 0 3 0 2009, 2010 GB/BOLD 5 Angola ANG 27 27 26 2006 ST (Nansen) 4 Senegal SEN 5 5 5 2007 CF 2 Northeastern Atlantic–Mediterranean Sea (NEAM) Portugal POR 3 3 3 2007 ST (IPIMAR) n.a. Portugal 0 7 0 2005, 2007 GB/BOLD n.a. Balearic Islands BAL 19 19 16 2006 ST (MedITS) 1 Balearic Islands 0 3 0 2013 ST (MedITS) 1 Sardinia SAR 11 11 8 2002, 2005 ST (MedITS) 1 Tuscany TUS 26 22 21 2005, 2006 ST (MedITS) 1 Tuscany 16 6 13 2008, 2010 ST (MedITS) 1 Algeria ALG 8 8 5 2002, 2003 FM (Algiers) 1 Algeria 10 9 8 2009, 2010 FM (Algiers) 1 Strait of Sicily—Adventura Bank SIC 22 22 22 2014 ST (MedITS) 1 Strait of Sicily—Maltese Bank 16 12 8 2000, 2002 ST (MedITS) 1 Strait of Sicily—Maltese Bank 0 11 0 2007 GB/BOLD 1 Adriatic Sea—Northern Italian coast ADR 39 31 20 2006, 2007 ST (MedITS; GruND) 1 Adriatic Sea—Croatian coast 24 24 8 2002, 2004 ST (MedITS; GruND) 1 Adriatic Sea—Southern Italian coast 19 16 19 2004 ST (MedITS; GruND) 1 Adriatic Sea—Albanian coast 19 13 17 2004 ST (MedITS; GruND) 1 Ionian Sea 4 3 4 2004 ST (MedITS; GruND) 1 Greece—Aegean coast GRE 0 3 0 2014 GB/BOLD 1 Israel ISR 8 7 7 2009 CF 1 Israel 0 3 0 2012 BOLD 1 Israel 0 4 0 2014 GB/BOLD 1 Levantine Sea LEV 0 2 0 2016 GB/BOLD 1 Levantine Sea 7 7 7 2009 CF 1 n.a.: not available. ST: Scientific Trawl survey. CF: Contracted Fishermen. FM: Fishery Market. GB: GenBank Database. BOLD: Barcoding of Life Database. Animals 2023,13, 2139 5 of 19 Animals 2023, 13, x FOR PEER REVIEW 5 of 20 Figure 1. Sampling sites of Raja miraletus species complex collected in the Northeastern Atlantic– Mediterranean and the Central–Southern Africa regions. A simplified representation of the main oceanic currents of the Eastern Atlantic is indicated. Acronyms of the Macro Area Codes are given as in Table 1. Colours in agreement with the haplotype network legend of Figure 2. Different dimensions of circles are related to sample size. Figure 1. Sampling sites of Raja miraletus species complex collected in the Northeastern Atlantic– Mediterranean and the Central–Southern Africa regions. A simplified representation of the main oceanic currents of the Eastern Atlantic is indicated. Acronyms of the Macro Area Codes are given as in Table 1. Colours in agreement with the haplotype network legend of Figure 2. Different dimensions of circles are related to sample size. Animals 2023,13, 2139 6 of 19 Animals 2023, 13, x FOR PEER REVIEW 6 of 20 Figure 2. TCS network of cytochrome oxidase subunit I (COI) haplotypes shown by Raja miraletus across most of its distribution area. CSA, Central–Southern Africa; NEAM, Northeastern Atlantic– Mediterranean Sea. Circles are proportional to haplotype frequencies. Black dots between branch nodes indicate substitutions. Black circles at network nodes represent unsampled haplotypes. Acronyms of the Macro Area Codes are given in Table 1. 2.2. Genetic Data Analysis Detailed protocols used for DNA extraction, PCR amplification, DNA sequencing and genotyping of mitochondrial and nuclear markers [56,58,61,62] are described in the Supplementary Material Text S1. 2.2.1. Genetic Diversity A total of 281 COI newly generated sequence electropherograms were manually edited and aligned by CLUSTAL W software [63] implemented in MEGA v.11 [64]. The presence of stop codons and sequencing error was verified through amino acidic translation [65]. Individual COI sequences were first compared with sequences deposited in public repositories in order to confirm their phylogenetic identity and rule out any error due to mishandling of samples on board or during the laboratory activities, namely GenBank (http://www.ncbi.nlm.nih.gov/genbank/, accessed on 21 May 2023) through the BLAST algorithm (http://blast.ncbi.nlm.nih.gov/Blast.cgi, accessed on 21 May 2023) and the Barcode of Life Data System (BOLD), using the BOLD identification engine ([66]; http://www.boldsystems.org, accessed on 21 May 2023). A total of 41 additional homologous COI sequences of R. miraletus complex were retrieved from both databases selecting records from different geographical locations (South Africa, Namibia, Strait of Sicily, Aegean Sea and Israel) when metadata giving the collection area were accessible Figure 2. TCS network of cytochrome oxidase subunit I (COI) haplotypes shown by Raja miraletus across most of its distribution area. CSA, Central–Southern Africa; NEAM, Northeastern Atlantic– Mediterranean Sea. Circles are proportional to haplotype frequencies. Black dots between branch nodes indicate substitutions. Black circles at network nodes represent unsampled haplotypes. Acronyms of the Macro Area Codes are given in Table 1. 2.2. Genetic Data Analysis Detailed protocols used for DNA extraction, PCR amplification, DNA sequencing and genotyping of mitochondrial and nuclear markers [ 56 , 58 , 61 , 62 ] are described in the Supplementary Material Text S1. 2.2.1. Genetic Diversity A total of 281 COI newly generated sequence electropherograms were manually edited and aligned by CLUSTAL W software [ 63 ] implemented in MEGA v.11 [ 64 ]. The presence of stop codons and sequencing error was verified through amino acidic translation [ 65 ]. Individual COI sequences were first compared with sequences deposited in public repositories in order to confirm their phylogenetic identity and rule out any error due to mishandling of samples on board or during the laboratory activities, namely GenBank (http://www.ncbi.nlm.nih.gov/genbank/, accessed on 21 May 2023) through the BLAST algorithm (http://blast.ncbi.nlm.nih.gov/Blast.cgi, accessed on 21 May 2023) and the Barcode of Life Data System (BOLD), using the BOLD identification engine ([ 66 ]; http://www.boldsystems.org, accessed on 21 May 2023). A total of 41 additional homologous COI sequences of R. miraletus complex were retrieved from both databases selecting records from different geographical locations (South Africa, Namibia, Strait of Sicily, Animals 2023,13, 2139 7 of 19 Aegean Sea and Israel) when metadata giving the collection area were accessible ([ 26 , 67 – 79 ]; Table 1; see Table S1 for more detail). The retrieved sequences were aligned with those newly generated and a final dataset of 322 homologous COI sequences was obtained. The number of polymorphic sites (S), the number of haplotypes (H), the haplotype diversity (Hd), the nucleotide diversity ( π ; [ 80 ]) and their standard deviations were calculated using DNASP v.6 [ 81 ]. The haplotype frequencies were estimated using ARLEQUIN v.3.5.2.2. [82]. The average genetic distances observed within and between the two identified Central– Southern African and the Northeastern Atlantic–Mediterranean clades of R. miraletus complex were calculated using the Tamura-Nei (1993) model implemented in MEGA as the best evolutionary substitution model following the corrected Akaike Information Criterion (AICc; [ 83 ]). Genetic distances were then compared with the range of those estimated among other congeneric species. For this, we retrieved homologous COI sequences of the following species public databases (NCBI and BOLD): Raja straeleni Poll 1951, Raja microocellata Montagu 1818, Raja asterias Delaroche 1809, Raja brachyura Lafont 1873, Raja clavata Linnaeus 1758, Raja montagui Fowler 1910, Raja polystigma Regan 1923, Raja radula Delaroche 1809 and Raja undulata Lacepède 1802 ([26,84]; Table S2). A total of 256 chromatograms for each of the eight EST-SSR loci were obtained and manually inspected using GENEMAPPER v.5 (Applied Biosystems, Waltham, MA, USA). Allele calling and binning were performed with GENEMAPPER. The presence of null alleles, stuttering and allele drop-out was tested using MICRO-CHECKER v.2.2.3 [ 85 ] with 1000 randomisations on Bonferroni correction. The multilocus EST-SSR genotypes were analysed using GENETIX v.4.05 [ 86 ] to estimate observed (H O ) and expected heterozygosity (H E ) and the number of alleles (A). Jack-knifing over loci was performed to assess the single-locus effects on Weir and Cockerham’s F-statistics estimators. The deviation from the Hardy–Weinberg equilibrium (HWE) and Linkage Disequilibrium (LD) was investigated using GENEPOP on the web v.4.2 [ 87 ]. Allelic richness (Ar) and inbreeding coefficient (F IS ) were estimated using FSTAT v.2.9.3.2 [88]. 2.2.2. Population Connectivity and Phylogenetic Signals The phylogenetic relationships among individual haplotypes were inferred by the TCS method implemented in the software POPART [ 89 ]. The graphical representation of the resulting network has been modified with Adobe Photoshop. Population connectivity within the R. miraletus species complex was investigated by estimates of Φ st and Fst values using ARLEQUIN with 10,000 permutations, p< 0.05. The Tamura-Nei (1993) substitution model was applied to the mtDNA dataset to estimate Φ st values. Genetic heterogeneity among the geographical samples was also assessed by a hierarchical analysis of molecular variance (AMOVA, [ 90 ]). Significance was assessed using a null distribution of the test statistic generated by 10,000 random permutations of the individuals in the samples. The significance threshold of the pairwise comparisons (p< 0.05) was adjusted with the sequential Bonferroni correction for multiple simultaneous comparisons [91] implemented in the R package “sgof” [92]. To unravel the individual-based genetic clustering, the COI and the EST-SSR datasets were analysed using the Bayesian algorithm implemented in BAPS v.6.0 [ 93 ] and STRUCTURE v.2.3.4 [ 94 ], respectively. The latter analysis on SSRs loci was carried out assuming an admixture ancestry model with the geographical origin of samples as prior information (LOCPRIOR models), associated with a correlated allele frequencies model. For each simulation of K (1–20), five independent replicates were run, setting a burn-in of 200,000 iterations and 500,000 iterations for the Markov Chain Monte Carlo (MCMC) simulation. Cluster matching and permutation were performed using CLUMPAK [ 95 ], while the most likely value for K was estimated from the mean log probability of the data using four alternative statistics (medmedk, medmeak, maxmedk and maxmeak) carried out using STRUCTURESELECTOR [96]. Animals 2023,13, 2139 8 of 19 Discriminant analysis of principal components (DAPC) using the R package Adegenet v.2.0.1 [97] was implemented in R v.4.0.5 (R Core Team [98]) using sampling locations as a priori groups (K = 5). Then, the optimal number of PCs to use in the DAPC was determined using the optim.a.score() command. Phylogenetic relationships between and within the Central–Southern African and Northeastern Atlantic–Mediterranean COI lineages were estimated using a Bayesian coalescent approach, implemented in BEAST v.1.10.4 [ 99 ]. Sequences of R. undulata (BOLD record ELAME177-09, NCBI accession numbers KT307412, KT307413, KT307414), the closest related species to the R. miraletus complex, were used as outgroups. The Bayesian reconstruction was obtained using the Hasegawa, Kishino and Yano (HKY + G) model of evolution [ 100 ], as the most appropriated model inferred by MEGA software, a strict molecular clock model, the Yule Process as species tree prior, the Piecewise linear and constant root as population size prior. To ensure convergence of the posterior distributions, an MCMC run of 60,000,000 generations sampled every 1000 generations with the first 25% of the sampled points removed as burn-in was performed. We analysed the log file using TRACER v.1.7.2 [ 101 ] to calculate the robustness of the posterior distributions for all parameters and recover average divergence time and 95% confidence intervals. The plausible trees obtained with BEAST were summarised using the program TREEANNOTATOR and the resulting phylogenetic relationships among population samples and the posterior probabilities at nodes were visualised with FigTree v.1.4.4 (available at http://tree.bio.ed.ac.uk/software/figtree/, accessed on 21 May 2023) and edited with the iTol v.6.7.5 online tool [102]. Cryptic species were also delimited by using two different methods: the distancebased method “Automatic Barcode Gap Discovery” (ABGD; [ 103 ]) computed on the online web application (http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html, accessed on 23 May 2023), using default values, and the phylogenetic-based method Bayesian Poisson Tree Process (bPTP, [ 104 ]), conducted on the web server (available at http://species.h-its. org/ptp/, accessed on 23 May 2023) with 100,000 MCMC generations, a thinning interval of 100% and 10% of burn-in. 3. Results 3.1. Genetic Diversity The COI dataset was a total of 322 sequences, while the EST-SSR dataset was made up of a total of 256 individuals overall distributed in 14 geographic samples (Tables 1and S1) . The final COI alignment consisted of 529 nucleotide positions and included 76 variable sites (14.3%) and 64 parsimony informative sites (12.1%). On average, COI polymorphism showed low estimates of nucleotide diversity ( π ) and very high haplotype diversity (Hd), with ANG being the most polymorphic sample (Hd = 0.858 ± 0.041 SD, π= 0.02543 ±0.00380 SD , K = 13.453; Table S3). Thirty-nine haplotypes were found and none were shared between samples from the Northeastern Atlantic–Mediterranean and Central–Southern African Regions (Figure 2and Table S4 in Supplementary Material). The average Tamura-Nei genetic distances (DTN) among geographical samples of the Northeastern Atlantic–Mediterranean were extremely low (DTN = 0.0025 ± 0.0011 SE; Table 2), while those observed among geographical samples of the Central–Southern Africa were an order of magnitude higher (mean DTN = 0.0183 ± 0.0029 SE). The DTN between Northeastern Atlantic–Mediterranean and Central–Southern Africa samples were much higher (mean DTN = 0.0733 ± 0.0117 SE) and is comparable between species distances recorded among species in the genus Raja (Table 2). Summary statistics of the eight polymorphic microsatellite loci per geographical sample and over all the loci considered are shown in Table S5 in the Supplementary Material. Mean allelic richness (Ar mean ) ranged from 1.242 (POR) to 1.724 (SEN). After Bonferroni correction, significant LD was not detected between any pair of loci, and the average mean observed and expected heterozygosity (H O /H E ) for the eight loci was 0.259/0.392. After applying the Bonferroni correction, significant HWE departures were found over all loci in Animals 2023,13, 2139 9 of 19 several geographical samples, apart from SEN, POR, SIC and ISR. The Portuguese sample was monomorphic at five loci (LERI 26, LERI 34, LERI 63, LERI 40 and LERI 44). Overall, MICRO-CHECKER results detected the presence of scoring errors such as stuttering and null alleles in all loci (Table S5), regardless of the geographical sample. Nevertheless, we did not exclude any of them since Jack-knife analysis did not reveal outliers outside of the confidence interval (Table S6). Table 2. Mean genetic distances within (in grey, in diagonal) and between congeneric Raja species. Standard error values for distances between species are reported above the diagonal. CSA, Central– Southern Africa; NEAM, Northeastern Atlantic–Mediterranean Sea. The mean genetic distance between NEAM and CSA is indicated in bold. Raja asterias Raja brachyura Raja clavata Raja microocellata Raja montagui Raja polystigma Raja radula Raja straeleni Raja undulata Raja miraletus (CSA) Raja miraletus (NEAM) Raja asterias 0.0025 ± 0.0012 0.0139 0.0106 0.0132 0.0135 0.0121 0.0095 0.0111 0.0121 0.0155 0.0189 Raja brachyura 0.0863 0.0030 ± 0.0016 0.0107 0.0094 0.0119 0.0114 0.0121 0.0125 0.0136 0.0141 0.0172 Raja clavata 0.0591 0.0514 0.0038 ± 0.0000 0.0115 0.0103 0.0113 0.0070 0.0059 0.0141 0.0133 0.0169 Raja microocellata 0.0877 0.0465 0.0598 0.0000 ± 0.0000 0.0125 0.0107 0.0114 0.0115 0.0127 0.0140 0.0172 Raja montagui 0.0829 0.0606 0.0511 0.0638 0.0000 ± 0.0000 0.0066 0.0126 0.0116 0.0128 0.0132 0.0172 Raja polystigma 0.0748 0.0526 0.0514 0.0596 0.0229 0.0000 ± 0.0000 0.0108 0.0112 0.0123 0.0119 0.0164 Raja radula 0.0487 0.0685 0.0280 0.0648 0.0646 0.0564 0.0010 ± 0.0012 0.0071 0.0119 0.0137 0.0173 Raja straeleni 0.0593 0.0593 0.0150 0.0601 0.0590 0.0558 0.0282 0.0019 ± 0.0010 0.0122 0.0131 0.0174 Raja undulata 0.0767 0.0736 0.0733 0.0799 0.0790 0.0711 0.0742 0.0735 0.0009 ± 0.0010 0.0122 0.0170 Raja miraletus (CSA) 0.1053 0.0857 0.0861 0.0947 0.0830 0.0778 0.0908 0.0907 0.0770 0.0183 ± 0.0029 0.0117 Raja miraletus (NEAM) 0.1072 0.1020 0.0989 0.1007 0.0921 0.0897 0.0974 0.1033 0.0959 0.0733 0.0025 ± 0.0010 3.2. Population Connectivity and Phylogenetic Signals Caution should be applied when interpreting the results obtained here due to the small sample size for some localities and the subsequent decrease in the discriminatory power of the analyses. The TCS network of the COI haplotypes (Figure 2) identified two main haplogroups, differentiated by at least 30 mutations and corresponding to the Central–Southern African and the Northeastern Atlantic–Mediterranean samples. The former haplogroup included 23 haplotypes that grouped into four largely differentiated geographic clusters located off Senegal, Angola/Namibia/South Africa and two off Angola. The Senegalese cluster formed only by the SEN sample (N = 5) showed three slightly differentiated private haplotypes. In contrast, the Angolan sample (ANG, N = 27) showed strongly differentiated haplotypes grouped into two endemic Angolan subclusters together and a third cluster shared with the South African (SAF, N = 40) and Namibian (NAM = 3) samples. The Northeastern Atlantic–Mediterranean haplogroup included 16 weakly divergent haplotypes (Figure 2 and Table S4). Four of them were shared by several samples and areas: (i) the haplotype Hap_24 was shared by Portuguese, Algerian and Strait of Sicily samples; (ii) the most frequent Hap_25 was shared by samples from Algeria, Balearic Islands, Sardinia, Strait of Sicily, Tuscany and Adriatic Sea; (iii) the Hap_28 was shared by samples from Algeria, Strait of Sicily and Adriatic Sea; iv) the Hap_32 was shared by Adriatic and Greek samples. In contrast, the Eastern Mediterranean samples from the Israeli coast (Hap_37 and Hap_38) and Levantine Sea (Hap_39) yielded only three endemic haplotypes. Most of the pairwise Φ st values among 14 geographical samples based on COI data were significant even after the Bonferroni correction was applied (Table S7). High levels of differentiation were observed between the African and Northeastern Atlantic– Mediterranean samples, as well as between the Western and Eastern Mediterranean. The Animals 2023,13, 2139 16 of 19 17. McEachran, J.D.; Dunn, K.A. Phylogenetic analysis of skates, a morphologically conservative clade of Elasmobranchs (Chondrichthyes: Rajidae). Copeia 1998,1998, 271–290. [CrossRef] 18. Tinti, F.; Ungaro, N.; Pasolini, P.; De Panfilis, M.; Garoia, F.; Guarniero, I.; Sabelli, B.; Marano, G.; Piccinetti, C. Development of molecular and morphological markers to improve species-specific monitoring and systematics of Northeast Atlantic and Mediterranean skates (Rajiformes). J. Exp. Mar. Biol. Ecol. 2003,288, 149–165. [CrossRef] 19. Iglésias, S.P.; Toulhoat, L.; Sellos, D.Y. Taxonomic confusion and market mislabelling of threatened skates: Important consequences for their conservation status. Aquat. Conserv. Mar. Freshw. Ecosyst. 2010,20, 319–333. [CrossRef] 20. Carugati, L.; Melis, R.; Cariani, A.; Cau, A.; Crobe, V.; Ferrari, A.; Follesa, M.C.; Geraci, M.L.; Iglésias, S.P.; Pesci, P.; et al. Combined COI barcode-based methods to avoid mislabelling of threatened species of deep-sea skates. Anim. Conserv. 2022 ,25, 38–52. [CrossRef] 21. Griffiths, A.M.; Sims, D.W.; Cotterell, S.P.; El Nagar, A.; Ellis, J.R.; Lynghammar, A.; McHugh, M.; Neat, F.C.; Pade, N.G.; Queiroz, N.; et al. Molecular markers reveal spatially segregated cryptic species in a critically endangered fish, the common skate (Dipturus batis). Proc. R. Soc. B Biol. Sci. 2010,277, 1497–1503. [CrossRef] 22. Cannas, R.; Follesa, M.C.; Cabiddu, S.; Porcu, C.; Salvadori, S.; Iglésias, S.P.; Deiana, A.M.; Cau, A. Molecular and morphological evidence of the occurrence of the norwegian skate Dipturus nidarosiensis (Storm, 1881) in the Mediterranean Sea. Mar. Biol. Res. 2010,6, 341–350. [CrossRef] 23. Carbonara, P.; Bellodi, A.; Zupa, W.; Donnaloia, M.; Gaudio, P.; Neglia, C.; Follesa, M.C. Morphological traits and capture depth of the Norwegian skate (Dipturus nidarosiensis (Storm, 1881)) from two Mediterranean populations. J. Mar. Sci. Eng. 2021 ,9, 1462. [CrossRef] 24. Domingues, R.R.; Hilsdorf, A.W.S.; Gadig, O.B.F. The importance of considering genetic diversity in shark and ray conservation policies. Conserv. Genet. 2018,19, 501–525. [CrossRef] 25. Serena, F.; Mancusi, C.; Barone, M. Field identification guide to the skates (Rajidae) of the Mediterranean Sea. Guidelines for data collection and analysis. Biol. Mar. Mediterr. 2010,17, 204. 26. Cariani, A.; Messinetti, S.; Ferrari, A.; Arculeo, M.; Bonello, J.J.; Bonnici, L.; Cannas, R.; Carbonara, P.; Cau, A.; Charilaou, C.; et al. Improving the conservation of Mediterranean Chondrichthyans: The ELASMOMED DNA barcode reference library. PLoS ONE 2017,12, e0170244. [CrossRef] 27. Ebert, D.A.; Compagno, L.J. Biodiversity and systematics of skates (Chondrichthyes: Rajiformes: Rajoidei). Biol. Skates 2007 ,27, 5–18. 28. Last, P.R.; White, W.T.; de Carvalho, M.R.; Séret, B.; Stehmann, M.F.W.; Naylor, G.J.P. (Eds.) Rays of the World; CSIRO Publishing: Clayton, Australia; Comstock Publishing Associates: Sacramento, CA, USA; pp. i–ix + 1–790. ISBN 978-0-643-10914-8. 29. Serena, F.; Abella, A.J.; Bargnesi, F.; Barone, M.; Colloca, F.; Ferretti, F.; Fiorentino, F.; Jenrette, J.; Moro, S. Species diversity, taxonomy and distribution of Chondrichthyes in the Mediterranean and Black Sea. Eur. Zool. J. 2020,87, 497–536. [CrossRef] 30. Grant, V. The systematic and geographical distribution of hawkmoth flowers in the temperate North American flora. Bot. Gaz. 1983,144, 439–449. [CrossRef] 31. Quintero, I.; Keil, P.; Jetz, W.; Crawford, F.W. Historical biogeography using species geographical ranges. Syst. Biol. 2015 ,64, 1059–1073. [CrossRef] 32. Henriques, S.; Guilhaumon, F.; Villéger, S.; Amoroso, S.; França, S.; Pasquaud, S.; Cabral, H.N.; Vasconcelos, R.P. Biogeographical region and environmental conditions drive functional traits of estuarine fish assemblages worldwide. Fish Fish. 2017 ,18, 752–771. [CrossRef] 33. Neat, F.; Pinto, C.; Burrett, I.; Cowie, L.; Travis, J.; Thorburn, J.; Gibb, F.; Wright, P.J. Site fidelity, survival and conservation options for the threatened flapper skate (Dipturus cf. intermedia). Aquat. Conserv. Mar. Freshw. Ecosyst. 2015,25, 6–20. [CrossRef] 34. Frisk, M.G.; Jordaan, A.; Miller, T.J. Moving beyond the current paradigm in marine population connectivity: Are adults the missing link? Fish Fish. 2014,15, 242–254. [CrossRef] 35. Wearmouth, V.J.; Sims, D.W. Movement and behaviour patterns of the critically endangered common skate Dipturus batis revealed by electronic tagging. J. Exp. Mar. Biol. Ecol. 2009,380, 77–87. [CrossRef] 36. Hunter, E.; Buckley, A.A.; Stewart, C.; Metcalfe, J.D. Repeated seasonal migration by a thornback ray in the southern North Sea. J. Mar. Biol. Assoc. UK 2005,85, 1199–1200. [CrossRef] 37. Hunter, E.; Buckley, A.A.; Stewart, C.; Metcalfe, J.D. Migratory behaviour of the thornback ray, Raja clavata, in the southern North Sea. J. Mar. Biol. Assoc. UK 2005,85, 1095–1105. [CrossRef] 38. Musick, J.A.; Ellis, J.K.; Hamlett, W. Reproductive evolution of chondrichthyans. In HAMLETT, WC, Reproductive Biology and Phylogeny of Chondrichthyes, Sharks, Batoids and Chimaeras; CRC Press: Boca Raton, FL, USA, 2005; pp. 45–71. 39. Compagno, L.J.V.; Ebert, D.A. Southern African Skate Biodiversity and Distribution. In Biology of Skates; Ebert, D.A., Sulikowski, J.A., Eds.; Developments in Environmental Biology of Fishes 27; Springer: Dordrecht, The Netherlands, 2007; pp. 19–39, ISBN 978-1-4020-9703-4. 40. Stehmann, M.F.W.; Bürkel, D.L. Rajidae. In Fishes of the North-Eastern Atlantic and the Mediterranean; Whitehead, P.J.P., Bauchot, M., Hureau, J., Nielsen, J., Eds.; Unesco: Paris, France, 1984; Volume 1, ISBN 9789230022150. 41. Kadri, H.; Marouani, S.; Bradai, M.N.; Bouaïn, A. Food habits of the brown ray Raja miraletus (Chondrichthyes: Rajidae) from the Gulf of Gabès (Tunisia). Mar. Biol. Res. 2014,10, 426–434. [CrossRef] Animals 2023,13, 2139 17 of 19 42. Šanti´c, M.; Radja, B.; Pallaoro, A. Feeding habits of brown ray (Raja miraletus Linnaeus, 1758) from the eastern central Adriatic Sea. Mar. Biol. Res. 2013,9, 301–308. [CrossRef] 43. Tsikliras, A.C.; Stergiou, K.I. Age at maturity of Mediterranean marine fishes. Mediterr. Mar. Sci. 2015,16, 5–20. [CrossRef] 44. Dulvy, N.K. Raja miraletus. IUCN Red List of Threatened Species. Available online: https://doi.org/10.2305/IUCN.UK.2019-3 .RLTS.T124569516A124512700.en (accessed on 17 May 2023). 45. Dulvy, N.K.; Walls, R.H.L.; Abella, A.; Serena, F.; Bradai, M.N. Raja miraletus (Mediterranean Assessment). The IUCN Red List of Threatened Species. Available online: https://doi.org/10.2305/IUCN.UK.2020-3.RLTS.T124569516A176535719.en (accessed on 17 May 2023). 46. Ebert, D.A.; Wintner, S.P.; Kynes, P.M. An annotated checklist of the chondrichthyans of South Africa. Zootaxa 2021 ,4947, 1–127. [CrossRef] 47. Compagno, L.J.V.; Ebert, D.A.; Smale, M.J. Guide to the Sharks and Rays of Southern Africa; Struik: Cape Town, South Africa, 1989. 48. McEachran, J.D.; Séret, B.; Miyake, T. Morphological variation within Raja miraletus and status of R. ocellifera (Chondrichthyes, Rajoidei). Copeia 1989,1989, 629–641. [CrossRef] 49. Naylor, G.; Caira, J.; Jensen, K.; Rosana, K.; Straube, N.; Lakner, C. Elasmobranch phylogeny: A mitochondrial estimate based on 595 species. In Biology of Sharks and Their Relatives, 2nd ed.; CRC Press: Boca Raton, FL, USA, 2012; pp. 31–56, ISBN 978-1-4398-3924-9. 50. Caira, J.N.; Rodriguez, N.; Pickering, M. New African species of Echinobothrium (Cestoda: Diphyllidea) and implications for the identities of their skate hosts. J. Parasitol. 2013,99, 781–788. [CrossRef] 51. Last, P.R.; Séret, B. A New eastern central Atlantic skate Raja parva sp. nov. (Rajoidei: Rajidae) belonging to the Raja miraletus species complex. Zootaxa 2016,4147, 477–489. [CrossRef] 52. Hebert, P.D.N.; Cywinska, A.; Ball, S.L.; de Waard, J.R. Biological identifications through DNA barcodes. Proc. R. Soc. Lond. B Biol. Sci. 2003,270, 313–321. [CrossRef] 53. Knowlton, N. Cryptic and sibling species among the decapod Crustacea. J. Crustac. Biol. 1986,6, 356–363. [CrossRef] 54. Morgan, J.A.T.; Harry, A.V.; Welch, D.J.; Street, R.; White, J.; Geraghty, P.T.; Macbeth, W.G.; Tobin, A.; Simpfendorfer, C.A.; Ovenden, J.R. Detection of interspecies hybridisation in Chondrichthyes: Hybrids and hybrid offspring between Australian (Carcharhinus tilstoni) and common (C. limbatus) blacktip shark found in an Australian fishery. Conserv. Genet. 2012 ,13, 455–463. [CrossRef] 55. Arlyza, I.S.; Shen, K.-N.; Solihin, D.D.; Soedharma, D.; Berrebi, P.; Borsa, P. Species boundaries in the Himantura uarnak species complex (Myliobatiformes: Dasyatidae). Mol. Phylogenet. Evol. 2013,66, 429–435. [CrossRef] 56. Frodella, N.; Cannas, R.; Velonà, A.; Carbonara, P.; Farrell, E.; Fiorentino, F.; Follesa, M.; Garofalo, G.; Hemida, F.; Mancusi, C.; et al. Population connectivity and phylogeography of the mediterranean endemic skate Raja polystigma and evidence of its hybridization with the parapatric sibling R. montagui.Mar. Ecol. Prog. Ser. 2016,554, 99–113. [CrossRef] 57. Siskey, M.R.; Shipley, O.N.; Frisk, M.G. Skating on thin ice: Identifying the need for species-specific data and defined migration ecology of Rajidae Spp. Fish Fish. 2019,20, 286–302. [CrossRef] 58. El Nagar, A.; McHugh, M.; Rapp, T.; Sims, D.W.; Genner, M.J. Characterisation of polymorphic microsatellite markers for skates (Elasmobranchii: Rajidae) from expressed sequence tags. Conserv. Genet. 2010,11, 1203–1206. [CrossRef] 59. Spedicato, M.T.; Massutí, E.; Mérigot, B.; Tserpes, G.; Jadaud, A.; Relini, G. The MEDITS trawl survey specifications in an ecosystem approach to fishery management. Sci. Mar. 2019,83, 9–20. [CrossRef] 60. Relini, G. Demersal trawl surveys in Italian Seas: A short review. Actes Colloq. IFREMER 2000 ,26, 76–93. Available online: https://archimer.ifremer.fr/doc/00716/82814/87636.pdf (accessed on 20 March 2023). 61. Ivanova, N.V.; Zemlak, T.S.; Hanner, R.H.; Hebert, P.D.N. Universal primer cocktails for fish DNA barcoding. Mol. Ecol. Notes 2007,7, 544–548. [CrossRef] 62. Catalano, G.; Crobe, V.; Ferrari, A.; Baino, R.; Massi, D.; Titone, A.; Mancusi, C.; Serena, F.; Cannas, R.; Carugati, L.; et al. Strongly structured populations and reproductive habitat fragmentation increase the vulnerability of the Mediterranean starry ray Raja asterias (Elasmobranchii, Rajidae). Aquat. Conserv. Mar. Freshw. Ecosyst. 2022,32, 66–84. [CrossRef] 63. Thompson, J.D.; Higgins, D.G.; Gibson, T.J. CLUSTAL W: Improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 1994 ,22, 4673–4680. [CrossRef] 64. Tamura, K.; Stecher, G.; Kumar, S. MEGA11: Molecular Evolutionary Genetics Analysis version 11. Mol. Biol. Evol. 2021 ,38, 3022–3027. [CrossRef] 65. Moulton, M.J.; Song, H.; Whiting, M.F. Assessing the effects of primer specificity on eliminating numt coamplification in DNA barcoding: A case study from Orthoptera (Arthropoda: Insecta). Mol. Ecol. Resour. 2010,10, 615–627. [CrossRef] 66. Ratnasingham, S.; Hebert, P.D.N. Bold: The Barcode of Life Data System (http://www.barcodinglife.org). Mol. Ecol. Notes 2007 ,7, 355–364. [CrossRef] 67. Crobe, V.; Ferrari, A.; Hanner, R.; Leslie, R.W.; Steinke, D.; Tinti, F.; Cariani, A. Molecular taxonomy and diversification of Atlantic skates (Chondrichthyes, Rajiformes): Adding more pieces to the puzzle of their evolutionary history. Life 2021 ,11, 596. [CrossRef] 68. Steinke, D.; Connell, A.D.; Hebert, P.D.N. Linking adults and immatures of South African marine fishes. Genome 2016 ,59, 959–967. [CrossRef] Animals 2023,13, 2139 18 of 19 69. Van Der Bank, H. DNA barcoding results for some Southern African elephantfish, guitarfish, rattails, rays, sharks and skates. Int. J. Oceanogr. Aquac. 2019,3, 000163. [CrossRef] 70. Serra-Pereira, B.; Moura, T.; Griffiths, A.M.; Serrano Gordo, L.; Figueiredo, I. Molecular barcoding of skates (Chondrichthyes: Rajidae) from the Southern Northeast Atlantic. Zool. Scr. 2011,40, 76–84. [CrossRef] 71. Costa, F.O.; Landi, M.; Martins, R.; Costa, M.H.; Costa, M.E.; Carneiro, M.; Alves, M.J.; Steinke, D.; Carvalho, G.R. A ranking system for reference libraries of DNA barcodes: Application to marine fish species from Portugal. PLoS ONE 2012 ,7, e35858. [CrossRef] 72. Ramírez-Amaro, S.; Ordines, F.; Picornell, A.; Castro, J.A.; Ramon, C.; Massutí, E.; Terrasa, B. The evolutionary history of mediterranean batoidea (Chondrichthyes: Neoselachii). Zool. Scr. 2018,47, 686–698. [CrossRef] 73. Ferrari, A.; Tinti, F.; Maresca, V.B.; Velonà, A.; Cannas, R.; Thasitis, I.; Costa, F.O.; Follesa, M.C.; Golani, D.; Hemida, F.; et al. Natural history and molecular evolution of demersal Mediterranean sharks and skates inferred by comparative phylogeographic and demographic analyses. PeerJ 2018,6, e5560. [CrossRef] 74. Landi, M.; Dimech, M.; Arculeo, M.; Biondo, G.; Martins, R.; Carneiro, M.; Carvalho, G.R.; Brutto, S.L.; Costa, F.O. DNA barcoding for species assignment: The case of Mediterranean marine fishes. PLoS ONE 2014,9, e106135. [CrossRef] 75. Vella, A.; Vella, N.; Schembri, S. A molecular approach towards taxonomic identification of elasmobranch species from Maltese fisheries landings. Mar. Genomics 2017,36, 17–23. [CrossRef] 76. Gkafas, G.A.; Megalofonou, P.; Batzakas, G.; Apostolidis, A.P.; Exadactylos, A. Molecular phylogenetic convergence within Elasmobranchii revealed by Cytochrome Oxidase Subunits. Biochem. Syst. Ecol. 2015,61, 510–515. [CrossRef] 77. Zambounis, A.G.; Ekonomou, G.; Megalofonou, P.; Batzakas, G.; Malandrakis, E.; Martsicalis, P.; Panagiotaki, P.; Neofitou, C.; Exadactylos, A. Molecular Phylogenetic Interrelations between Species of the Elasmobranchii subclass. 2010; unpublished; submitted to the EMBL/GenBank/DDBJ databases. 78. Shirak, A.; Dor, L.; Seroussi, E.; Ron, M.; Hulata, G.; Golani, D. DNA barcoding of fish species from the Mediterranean coast of Israel. Mediterr. Mar. Sci. 2016,17, 459–466. [CrossRef] 79. Yokes, M.B. DNA Barcoding of Marine Fish Species from Turkish Coastline. 2016; unpublished; submitted to the EMBL/GenBank/DDBJ databases. 80. Nei, M. Molecular Evolutionary Genetics; Columbia University Press: New York, NY, USA, 1987. [CrossRef] 81. Rozas, J.; Ferrer-Mata, A.; Sánchez-DelBarrio, J.C.; Guirao-Rico, S.; Librado, P.; Ramos-Onsins, S.E.; Sánchez-Gracia, A. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Mol. Biol. Evol. 2017,34, 3299–3302. [CrossRef] 82. Excoffier, L.; Lischer, H.E.L. Arlequin Suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 2010,10, 564–567. [CrossRef] 83. Akaike, H. A new look at the statistical model identification. Curr. Contents Eng. Technol. Appl. Sci. 1981,12, 42. 84. Knebelsberger, T.; Landi, M.; Neumann, H.; Kloppmann, M.; Sell, A.F.; Campbell, P.D.; Laakmann, S.; Raupach, M.J.; Carvalho, G.R.; Costa, F.O. A reliable DNA barcode reference library for the identification of the North European shelf fish fauna. Mol. Ecol. Resour. 2014,14, 1060–1071. [CrossRef] 85. Van Oosterhout, C.; Hutchinson, W.F.; Wills, D.P.M.; Shipley, P. MICRO-CHECKER: Software for identifying and correcting genotyping errors in microsatellite data. Mol. Ecol. Notes 2004,4, 535–538. [CrossRef] 86. Belkhir, K.; Borsa, P.; Chikhi, L.; Raufaste, N.; Bonhomme, F. GENETIX 4.05, Logiciel Sous Windows TM Pour la Génétique des Populations. 1996. Available online: https://www.scienceopen.com/document?_vid=7cfcd230-1958-4cfc-a571-ce0ba003e63f (accessed on 1 April 2022). 87. Rousset, F. GENEPOP’007: A complete re-implementation of the GENEPOP software for Windows and Linux. Mol. Ecol. Resour. 2008,8, 103–106. [CrossRef] 88. Goudet, J. FSTAT 2.9.3, a Program to Estimate and Test Gene Diversities and Fixation Indices. 2001. Available online: http: //www2.unil.ch/popgen/softwares/fstat.htm (accessed on 15 March 2023). 89. Leigh, J.W.; Bryant, D. POPART: Full-Feature software for haplotype network construction. Methods Ecol. Evol. 2015 ,6, 1110–1116. [CrossRef] 90. Excoffier, L.; Smouse, P.E.; Quattro, J.M. Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction Data. Genetics 1992,131, 479–491. [CrossRef] 91. Rice, W.R. Analyzing tables of statistical tests. Evolution 1989,43, 223. [CrossRef] 92. Castro-Conde, I.; de Uña-Álvarez, J. Sgof: An R package for multiple testing problems. R J. 2014,6, 96–113. [CrossRef] 93. Cheng, L.; Connor, T.R.; Sirén, J.; Aanensen, D.M.; Corander, J. Hierarchical and spatially explicit clustering of DNA sequences with BAPS software. Mol. Biol. Evol. 2013,30, 1224–1228. [CrossRef] 94. Hubisz, M.J.; Falush, D.; Stephens, M.; Pritchard, J.K. Inferring weak population structure with the assistance of sample group information. Mol. Ecol. Resour. 2009,9, 1322–1332. [CrossRef] 95. Kopelman, N.M.; Mayzel, J.; Jakobsson, M.; Rosenberg, N.A.; Mayrose, I. Clumpak: A program for identifying clustering modes and packaging population structure inferences across K. Mol. Ecol. Resour. 2015,15, 1179–1191. [CrossRef] 96. Li, Y.-L.; Liu, J.-X. StructureSelector: A web-based software to select and visualize the optimal number of clusters using multiple methods. Mol. Ecol. Resour. 2018,18, 176–177. [CrossRef] 97. Jombart, T.; Devillard, S.; Balloux, F. Discriminant analysis of principal components: A new method for the analysis of genetically structured populations. BMC Genet. 2010,11, 94. [CrossRef] Animals 2023,13, 2139 19 of 19 98. R Core Team R: A language and environment for statistical computing 2021. 99. Suchard, M.A.; Lemey, P.; Baele, G.; Ayres, D.L.; Drummond, A.J.; Rambaut, A. Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10. Virus Evol. 2018,4, vey016. [CrossRef] 100. Hasegawa, M.; Kishino, H.; aki Yano, T. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J. Mol. Evol. 1985,22, 160–174. [CrossRef] 101. Rambaut, A.; Drummond, A.J.; Xie, D.; Baele, G.; Suchard, M.A. Posterior summarization in Bayesian phylogenetics using tracer 1.7. Syst. Biol. 2018,67, 901–904. [CrossRef] 102. Letunic, I.; Bork, P. Interactive Tree Of Life (ITOL): An online tool for phylogenetic tree display and annotation. Bioinformatics 2007,23, 127–128. [CrossRef] 103. Puillandre, N.; Lambert, A.; Brouillet, S.; Achaz, G. ABGD, Automatic Barcode Gap Discovery for primary species delimitation. Mol. Ecol. 2012,21, 1864–1877. [CrossRef] 104. Zhang, J.; Kapli, P.; Pavlidis, P.; Stamatakis, A. A general species delimitation method with applications to phylogenetic placements. Bioinforma. Oxf. Engl. 2013,29, 2869–2876. [CrossRef] 105. Wallace, J.H. The batoid fishes of the east coast of southern Africa. Part III: Skates and electric rays. S. Afr. Assoc. Mar. Biol. Res. Investig. Rep. 1967,17, 1–62. 106. Hirschfeld, M.; Dudgeon, C.; Sheaves, M.; Barnett, A. Barriers in a sea of elasmobranchs: From fishing for populations to testing hypotheses in population genetics. Glob. Ecol. Biogeogr. 2021,30, 2147–2163. [CrossRef] 107. Sandoval-Castillo, J.; Beheregaray, L.B. Oceanographic heterogeneity influences an ecological radiation in elasmobranchs. J. Biogeogr. 2020,47, 1599–1611. [CrossRef] 108. Henriques, R.; Potts, W.M.; Sauer, W.H.; Santos, C.V.; Kruger, J.; Thomas, J.A.; Shaw, P.W. Molecular genetic, life-history and morphological variation in a coastal warm-temperate sciaenid fish: Evidence for an upwelling-driven speciation event. J. Biogeogr. 2016,43, 1820–1831. [CrossRef] 109. Henriques, R.; Potts, W.M.; Santos, C.V.; Sauer, W.H.; Shaw, P.W. Population connectivity and phylogeography of a coastal fish, Atractoscion aequidens (Sciaenidae), across the Benguela current region: Evidence of an ancient vicariant event. PLoS ONE 2014 ,9, e87907. [CrossRef] 110. Chevolot, M.; Hoarau, G.; Rijnsdorp, A.D.; Stam, W.T.; Olsen, J.L. Phylogeography and population structure of thornback ray (Raja clavata L., Rajidae). Mol. Ecol. 2006,15, 3693–3705. [CrossRef] 111. Valsecchi, E.; Pasolini, P.; Bertozzi, M.; Garoia, F.; Ungaro, N.; Vacchi, M.; Sabelli, B.; Tinti, F. Rapid Miocene-Pliocene dispersal and evolution of mediterranean rajid fauna as inferred by mitochondrial gene variation. J. Evol. Biol. 2005 ,18, 436–446. [CrossRef] 112. Patarnello, T.; Volckaert, F.A.M.J.; Castilho, R. Pillars of Hercules: Is the Atlantic–Mediterranean transition a phylogeographical break? Mol. Ecol. 2007,16, 4426–4444. [CrossRef] 113. Melis, R.; Vacca, L.; Cariani, A.; Carugati, L.; Charilaou, C.; Di Crescenzo, S.; Ferrari, A.; Follesa, M.C.; Mancusi, C.; Pinna, V.; et al. Baseline genetic distinctiveness supports structured populations of thornback ray in the Mediterranean Sea. Aquat. Conserv. Mar. Freshw. Ecosyst. 2023,33, 458–471. [CrossRef] 114. Serena, F. Field Identification Guide to the Sharks and Rays of the Mediterranean and Black Sea; Food & Agriculture Organisation of The United Nations: Rome, Italy, 2005; 97 pp. + 11 colour plates + egg cases. Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.