Next Article in Journal
Effects of Individual and Pair Housing of Calves on Short-Term Health and Behaviour on a UK Commercial Dairy Farm
Next Article in Special Issue
Editorial: Sharks and Skates—Ecology, Distribution and Conservation
Previous Article in Journal
Intermittent and Mild Cold Stimulation Maintains Immune Function Stability through Increasing the Levels of Intestinal Barrier Genes of Broilers
Previous Article in Special Issue
Aggregative Behaviour of Spiny Butterfly Rays (Gymnura altavela, Linnaeus, 1758) in the Shallow Coastal Zones of Gran Canaria in the Eastern Central Atlantic
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Animals
Volume 13
Issue 13
10.3390/ani13132139
animals-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
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

by
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
1
Department of Biological, Geological and Environmental Sciences, University of Bologna, 40126 Bologna, Italy
2
Department of Life and Environmental Sciences, University of Cagliari, 09126 Cagliari, Italy
3
Branch Fisheries Management, Department Agriculture, Forestry and Fisheries, Cape Town 8018, South Africa
4
Institute for Biological Resources and Marine Biotechnology, National Research Council, 91026 Trapani, Italy
5
Laboratory of Marine Biology and Fisheries, Department Biological, Geological and Environmental Sciences, University of Bologna, 61032 Fano, Italy
6
Centre of Molecular and Environmental Biology (CBMA) and ARNET-Aquatic Research Network, Department of Biology, University of Minho, Campus de Gualtar, 4710-057 Braga, Portugal
7
Department of Evolution, Systematics and Ecology, The Hebrew University of Jerusalem, Jerusalem 9190401, Israel
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
9
Institute of Marine Research, 5817 Bergen, Norway
10
Department of Biosciences, Biotechnologies and Environment, University of Bari Aldo Moro, 70125 Bari, Italy
11
COISPA Technology and Research, 70126 Bari, Italy
12
Stazione Zoologica Anton Dohrn, 90149 Palermo, Italy
 
add Show full affiliation list
 
remove Hide full affiliation list
*
Author to whom correspondence should be addressed.
†
These authors contributed equally to this work.
Animals 2023, 13(13), 2139; https://doi.org/10.3390/ani13132139
Submission received: 31 May 2023 / Revised: 25 June 2023 / Accepted: 25 June 2023 / Published: 28 June 2023
(This article belongs to the Special Issue Sharks and Skates: Ecology, Distribution and Conservation)
Browse Figures
Review Reports Versions Notes

Abstract

:

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 differentiation of ancient South African clades. The hidden genetic structure presented here may represent a valuable support to species’ conservation action plans.

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., pre- and 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,10,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,18,19,20,21,22,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,28,29].
Investigations into the role of biogeographical barriers on the speciation of marine organisms have increasingly concentrated across several taxa [30,31,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,36,37,38,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].
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,55,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 (EST-SSRs; [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 (Table 1 and Table 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 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.

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 ([26,67,68,69,70,71,72,73,74,75,76,77,78,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 (HO) and expected heterozygosity (HE) 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 (FIS) 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].
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 distance-based 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 (Table 1 and Table 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 2 and 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 (Armean) 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 (HO/HE) for the eight loci was 0.259/0.392. After applying the Bonferroni correction, significant HWE departures were found over all loci in 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).

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 EST-SSR data showed a similar pattern of genetic differentiation in terms of Fst calculated over 12 macro areas, even after the Bonferroni correction was applied (Table S8).
The hierarchical AMOVA found the highest percentage of molecular variation among groups with four sample groupings (AMOVA 3: Southern Africa vs. Angola vs. Senegal vs. Northeastern Atlantic–Mediterranean Sea; Table S9) when using the COI dataset, and with five sample groups (40.06%, AMOVA 4: Southern Africa vs. Angola vs. Senegal vs. Portugal–West Mediterranean vs. Eastern Mediterranean) when using the AMOVA 5 (six groups) explained the genetic variation among samples with the proportion of the genetic variation among populations within very low groups (2.58% with mitochondrial data and 7.72% with the EST-SSR data).
The Bayesian clustering analysis based on mtDNA data (Figure 3a) revealed six genetic clusters: the first cluster (green) included individuals from SAF, NAM and ANG; the second (purple) and third clusters (light blue) ANG; the fourth cluster (red) was unique to SEN; the fifth (yellow) and sixth (blue) clusters were unique to the Northeastern Atlantic and Mediterranean with individuals from POR, BAL, SAR and TUS in the fifth cluster and those from GRE, ISR and LEV in the sixth cluster. Individuals from ALG, SIC and ADR were randomly associated with the fifth and the sixth clusters, suggesting the southern Mediterranean, especially the Siculo-Tunisian area, as a potential admixture zone of the latter clusters.
The outputs of the STRUCTURE analysis based on EST-SSR data analysed with STRUCTURESELECTOR did not provide clear-cut evidence of the most likely number of clusters using four alternative statistics (medmedk, medmeak, maxmedk and maxmeak), while the maximum value of ΔK was verified with K = 3 (Figures S1 and S2). Thus, the results from K = 2 to K = 7 were assessed with CLUMPAK (Figure S3). The barplot of the clustering K = 2 supported the separation of the Central–Southern Africa samples from Northeastern Atlantic–Mediterranean Sea and with an admixed genetic composition of the Senegalese individuals (Figure S3). The clustering K = 3 further discriminated between the Angolan sample from those of South Africa, as well as samples from the Western and Eastern Mediterranean Sea (Figure S3). At the same time, the Angolan sample displayed an intermediate genetic composition between the South African and Senegalese clusters. This trend was resolved by the clustering K = 4 (Figure 3b and Figure S3), corresponding to the best grouping revealed by AMOVA (Table S9). The plot showed a deep differentiation of the Northeastern Atlantic–Mediterranean samples east to the Strait of Sicily.
The DAPC computed on EST-SSR data with sampling locations used as a priori group identified 10 optimal numbers of PCs and separated five main clusters: South African, Angolan, Senegalese, Northeastern Atlantic–Western Mediterranean and Eastern Mediterranean clusters, with the South African and Angolan clusters partially overlapping (Figure S4).
Bayesian approach using MCMC simulation was used to test for a speciation signal (Yule process) between lineages from Central–Southern Africa and Northeastern Atlantic–Mediterranean (Figure 4, see Table S4 for haplotype distribution among samples). All effective sample size (ESS) values exceeded 200, indicating a solid evaluation of all parameters. The model based on the substitution rate estimated for mtDNA showed a clear separation between haplotypes from Central–Southern Africa (Hap_1 to Hap_23) and the Northeastern Atlantic–Mediterranean (Hap_24 to Hap_39; Figure 4). The phylogenetic relationships among lineages and haplotypes were congruent with the relationships obtained with the parsimony network results (Figure 2). Furthermore, within the main Central–Southern African lineage, four clusters of haplotypes were reconstructed with high posterior probability (PP = 1): the most basal Angolan haplotype H_20; a second cluster formed by six Angolan haplotypes (Hap_14 to Hap_19); a Senegalese cluster (Hap_21 to Hap_23); and a South African/Angolan/Namibian cluster formed by all the South African haplotypes (Hap_1–9), the Namibian haplotype Hap_10, the Angolan haplotypes (Hap_12 and Hap_13) and the Angolan/Namibian haplotype (Hap_11).
The two species delimitation approaches (the ABGD and bPTP methods) yielded the same result (Figure 4). Five groups have been identified: four are formed by samples from the Central–Southern Atlantic, in agreement with the phylogenetic reconstruction, and one formed by samples from the Northeastern Atlantic and Mediterranean.

4. Discussion

The R. miraletus species complex is widely distributed, occurring from the Mediterranean Sea down the west coast of Africa to South Africa. The geographically isolated population off the south coast of South Africa was originally described as R. ocellifera Regan, 1906, but was synonymised with R. miraletus [105]. An extensive morphological study by McEachran et al. [48] found a marked differentiation between specimens from South Africa and those from Mediterranean elements, whilst the West African samples, in particular those from Angola, showed intermediate meristic features. However, they considered the differences between the Mediterranean and South African populations to be clinal and concluded that R. miraletus is a polymorphic species with three partially separated populations. Strikingly stable gross morphology has always been misleading for taxonomists. Subsequent integrated molecular and morphometric studies have shown that the three partially separated populations are valid species [51]. Raja ocellifera has been resurrected [44] and Last and Séret [51] described a new species, R. parva from Liberia and Angola. Last and Séret [51] stated that R. parva differed from the Senegalese Raja cf. miraletus 1 (sensu Naylor et al. [49]) in total body length, a shorter snout and a smaller tail, and they suggest that two further putative species (or taxa) occur off Senegal, Guinea, Liberia, down to Angola. Among the West African samples, the R. parva were the most distinct, even though some characters were intermediate between the Mediterranean and South African specimens [48].
The recent designation of the R. miraletus species complex [51] increased interest in an evolutionary and phylogenetic investigation of the complex based on more extensive sampling and analysis of combined nuclear and mitochondrial genetic data. Ferrari et al. [73] inferred population structure within R. miraletus (sensu stricto) across the Mediterranean Sea based on an analysis of nucleotide variation in three mtDNA markers, while Crobe et al. [67] preliminarily recognised four divergent COI lineages from the Eastern Atlantic populations of this morphologically conserved taxon. It should be noted that our study (i) was based on an unprecedently high number of specimens unequivocally assigned to the R. miraletus species complex; (ii) included specimens collected from the localities within the distribution ranges of the four putative taxa occur and (iii) was based on an integrated analytical approach that combines sequence variation of both mtDNA (the universal COI barcode region) and nuDNA(allelic variation in eight polymorphic EST-linked microsatellite loci) markers. This methodology enables our study to advance the molecular characterisation of this Hamletic taxon and to increase the knowledge on its status in the Mediterranean Sea.
The mitochondrial and microsatellite data consistently agreed in genetically defining the taxonomic and geographical boundaries of R. miraletus (sensu stricto), which are distributed in the entire Mediterranean Sea and the adjacent Northeastern Atlantic Ocean, at least in the Portuguese coastal waters. The estimated mean Tamura-Nei genetic distance between the Northeastern Atlantic–Mediterranean clade and the Central–Southern African clade (DTN = 0.0735) was in the upper range of corresponding pairwise interspecific estimates among several congeneric Raja species supporting a specific level of differentiation.
Surprisingly, the host–parasite specificity established between different African population of the R. miraletus species complex, with species from the genus Echinobothrium (Cestoda: Diphyillidea) highlighted by Caira et al. [50], supports the clades identified by Naylor et al. [49] based on molecular data. Specimens of R. ocellifera (R. cf. miraletus 1 from South Africa from [49]) hosted E. yiae Caira, Rodriguez and Pickering, 2013, those of R. miraletus from Senegal hosted E. mercedesae Caira, Rodriguez and Pickering, 2013 and specimens of a second clade off Senegal likely belonging to R. parva (R. cf. miraletus 2 sensu Naylor et al. [49]) hosted two additional new species [50]. These findings were partially supported by our study, where the Senegalese sample (SEN) is a distinct subunit (Figure 2 and Figure 3; Hap_21, Hap_22 and Hap_23 in Figure 4; Φst in Table S7). Unfortunately, the very limited sampling along the Senegalese coast (N = 5) prevents any definitive conclusions.
Overall, our results emphasised the strong differentiation between South Atlantic-Indian R. ocellifera and the Northeastern Atlantic–Mediterranean R. miraletus and justifies the resurrection of the former taxon. The complex oceanographic conditions along the African coast with alternating cold and warm currents, from north to south cold Canary Current, warm Angola Current, cold Benguela Current and warm Agulhas Current, undoubtedly played a role in speciation along the African coast. Likewise, the complexity oceanographic and geological discontinuities characterising the Eastern Atlantic and the Mediterranean Basin may likely influence phylogeography, population structure and connectivity as well as evolution at multiple taxonomic levels [48,106]. Oceanographic heterogeneities, such as current systems, play a key role for ecologically divergent natural selection in elasmobranchs, such as the ecological radiation of the genus Pseudobatos Last, Séret and Naylor, 2016 in the Gulf of California, strongly influenced by habitat heterogeneity and the geological history of the region [107]. This condition seemed to be true not only for skates, but even more so for African coastal bony fish (i.e., genus Argynosomus De la Plyaie, 1835), whose evolutionary histories, including the dispersal phase, have been influenced by the Benguela Current [108,109]. The Benguela Current region ranges from Cape Agulhas to Cape Frio, where the north-flowing cold Benguela Current meets the south-flowing warm Angola current (see Hirschfeld et al. [106], Figure 2d,f) for current model map details). The cold waters of the Benguela system are likely to have strongly reduced gene flow between R. ocellifera and Raja cf. miraletus (sensu Naylor et al. [49]). This diversification was evident from the level of mean sequence divergence observed between the two geographical populations estimated at 7.3%, a value comparable to, or higher than, any divergence measured between R. undulata and other congeneric taxa. On the other hand, the intertropical Canary current inflowing from the northeast could have influenced the diversification of the Senegalese taxa, whose migration southwards would be hampered by the intermittent Cape Blanc upwelling area. This upper boundary could have influenced and limited gene flow in Northern Mauritania. Like other skates [10,110,111] and teleost species [112], no genetic differentiation was observed between Northeastern Atlantic and Mediterranean populations of R. miraletus. This suggests that the Strait of Gibraltar did not represent an effective barrier to gene flow, rather than an accession gate to ancient refugia. In contrast, Melis et al. [113] found moderate significant population differentiation between the Mediterranean and the Atlantic Ocean in the congeneric thornback skate R. clavata, suggesting an effective role of the Strait in limiting the dispersal of individuals.
The slight genetic population structure observed within the Mediterranean Sea represents a true novelty for this species complex. The unforeseen regional East–West Mediterranean structure highlighted by nuclear markers (Figure 3b and Figure S4) could be linked to bathymetry and hydrogeological fronts or discontinuities. In particular, the shallow bathymetry characterising the Southwest part of the Mediterranean, coupled with the species’ preferences for continental shelf habitats, may likely enhance the dispersal of brown skate. The area ranging from the easternmost part of Sicily and the adjacent geo-morphological depression of the Calabrian Arc (down to 3000 m) is dominated by cyclonic/anticyclonic inversions of water masses. The combination of these features could have limited gene flow between western and eastern populations and driven the differentiation of the Eastern Mediterranean samples (ADR, ISR and LEV; Figure 1). The specific habitat and depth preferences, the less pronounced migratory behaviour and the limited dispersal capability of taxa belonging to the R. miraletus species complex are rather common characteristics among skates [62,73,113,114], although other congeneric species did not show such evidence of deep differentiation at both nuDNA and mtDNA markers [10,56,72,110]. In detail, the Mediterranean starry ray R. asterias showed a strongly structured population with three geographical clades corresponding to the western, central–western and central–eastern Mediterranean areas [62]. R. clavata showed a weak but detectable phylogeographic structure in the Levantine Sea [73] and a finer structuring located off the Algerian coasts and Tyrrhenian basins, suggesting the occurrence of additional barriers to dispersal [113]. On the contrary, R. polystigma showed a slightly differentiated Adriatic haplotype but a near panmictic population in the central–western part of the basin [56]. These different patterns of population structure in such closely related species can be explained by their bathymetric range which drive different ecological features [72].

5. Conclusions

This study coupled a massive and extensive sampling effort covering the full distribution of the R. miraletus species complex with analysis of genetic variability in both mtDNA and nuDNA markers, individual clustering, phylogeography and variance at different population levels. The results were partially congruent with previous taxonomic and meristic analyses. The use of both nuclear and mitochondrial markers resulted in identifying signals of species differentiation and in supporting the existence of at least five cryptic taxa within the R. miraletus species complex, four of which have been previously suggested with scattered genetic data. In addition to the evolutionary meaning of this evidence, genetics is shown to aid conservation efforts by revealing hidden diversity that deserves special attention and in the monitoring of taxa that are important fishing bycatch species. The new insights highlighted in the present research paper suggest that the extraordinary intraspecific diversity observed across such a wide geographical scale should be carefully considered to update or set dedicated and effective measures to reduce the impact of skate bycatch during fishing activities and improve their conservation.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/ani13132139/s1, Text S1: Molecular methods; Table S1: List of samples with detail of individual code, location, geographical coordinates, biological data and indication of available sequences integrated to the original dataset (BOLD records and NCBI Accession Numbers); Table S2: List of public sequences used for the estimation of genetic distances among the genus Raja. Details of BOLD Record and NCBI Accession Numbers are also listed; Table S3: mtDNA polymorphism and its parameters. N, number of individuals; Nh, number of haplotypes; S, number of polymorphic sites; Hd, haplotype diversity; π, nucleotide diversity; k, average number of nucleotide differences and SD, standard deviation. The locations are indicated by codes as given in Table 1; Table S4: Distribution of the COI haplotypes in the Raja miraletus complex population samples; Table S5: Summary statistics of the EST-SSRs polymorphism per geographical sample and over all the loci considered. Table S6: Jack-knifing over all loci at 95% confidence, 1000 bootstraps; Table S7: Pairwise Φst values (below the diagonal) and associated significance (above the diagonal). Significant p-values are highlighted in bold (p < 0.05), * p significant after sequential Bonferroni; Table S8: Pairwise Fst values (below the diagonal) and associated significance (above the diagonal). Significant p-values are highlighted in bold (p < 0.05), * p significant after sequential Bonferroni; Table S9: Five-levels AMOVA performed from two and six sample groupings on both mtDNA and nuDNA; Figure S1: Processing of STRUCTURE results from the analysis conducted on R. miraletus individuals using medmedk, medmeak, maxmedk and maxmeak statistics obtained with STRUCTURESELECTOR; Figure S2: Absolute value of derived delta K (ΔK) for each cluster number (K) from one to 12; Figure S3: Results of the Bayesian clustering using STRUCTURE from K = 2 to K = 7 obtained with CLUMPAK; each vertical bar represents one individual, in which a different colour represents the estimated cluster membership. Acronyms of the samples are given in Table 1; Figure S4: Discriminant analysis of principal components (DAPCs) scatter plot conducted on 256 Raja miraletus complex individuals based on EST-SSRs loci. Clusters refer to SAF, South Africa; ANG, Angola; SEN, Senegal; NATL-WMED, Northeastern Atlantic Ocean and Western Mediterranean Sea; EMED, Eastern Mediterranean Sea. Colours refer to the genetic clusters observed in BAPS plot (Figure 2).

Author Contributions

Conceptualisation, F.T. and A.C.; data curation, A.F.; investigation, A.F. and V.C.; methodology, A.F.; resources, R.C., R.W.L., F.S., M.S., F.O.C., D.G., F.H., D.Z.-P., L.S., P.C. and F.F.; supervision, F.T. and A.C.; visualisation, A.F.; writing—original draft, A.F. and V.C.; writing—review and editing, V.C., R.C., R.W.L., F.S., M.S., F.O.C., D.G., F.H., D.Z.-P., L.S., P.C., F.F., F.T. and A.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no dedicated funding and was supported by institutional funding from University of Bologna assigned to FT and AC (RFO and Canziani grants).

Institutional Review Board Statement

The authors declare that the sample of brown skate individuals analysed in the present work were obtained from commercial and scientific fisheries. The activity was conducted with the observation of the Regulation of the European Parliament and the Council for fishing in the General Fisheries Commission for the Mediterranean (GFCM) Agreement area and amending Council Regulation (EC) No. 1967/2006. This Regulation is de facto the unique authorisation needed to conduct this type of activities. Dead fish vertebrates are out of scope of “the University of Bologna—Animal Care and Use Committee” since it authorises/not authorises experimental procedures on living vertebrates.

Informed Consent Statement

Not applicable.

Data Availability Statement

The newly produced sequence data that support the findings of this study are openly available in GenBank at https://www.ncbi.nlm.nih.gov/genbank (accessed on 27 June 2023) reference number OR193802–OR194004. The EST-SSR data that support the findings of this study are available as a genotype matrix from the corresponding author upon request.

Acknowledgments

We are deeply grateful to all the Mediterranean Trawl Survey, the South African Trawl Survey and the EAF-Nansen Programme partners for their collaboration in the R. miraletus complex sample collection. Moreover, we thank Samira Vinjau and Alessia Crosara for their contribution in the analysis of specimens.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Harrison, R.G. Linking evolutionary pattern and process. In Endless Forms; Oxford University Press: Oxford, UK, 1998; pp. 19–31. [Google Scholar]
  2. Mallet, J. Hybridization, ecological races and the nature of species: Empirical evidence for the ease of speciation. Philos. Trans. R. Soc. B Biol. Sci. 2008, 363, 2971–2986. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Harrison, R.G.; Larson, E.L. Hybridization, introgression, and the nature of species boundaries. J. Hered. 2014, 105, 795–809. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Norris, R.D. Pelagic species diversity, biogeography, and evolution. Paleobiology 2000, 26, 236–258. [Google Scholar] [CrossRef]
  5. Pinheiro, H.T.; Bernardi, G.; Simon, T.; Joyeux, J.-C.; Macieira, R.M.; Gasparini, J.L.; Rocha, C.; Rocha, L.A. Island biogeography of marine organisms. Nature 2017, 549, 82–85. [Google Scholar] [CrossRef] [PubMed]
  6. Palumbi, S.R.; Lessios, H.A. Evolutionary animation: How do molecular phylogenies compare to Mayr’s reconstruction of speciation patterns in the sea? Proc. Natl. Acad. Sci. USA 2005, 102, 6566–6572. [Google Scholar] [CrossRef] [Green Version]
  7. Lynghammar, A.; Christiansen, J.S.; Griffiths, A.M.; Fevolden, S.-E.; Hop, H.; Bakken, T. DNA barcoding of the Northern Northeast Atlantic skates (Chondrichthyes, Rajiformes), with remarks on the widely distributed starry ray. Zool. Scr. 2014, 43, 485–495. [Google Scholar] [CrossRef] [Green Version]
  8. Palumbi, S.R. Genetic divergence, reproductive isolation, and marine speciation. Annu. Rev. Ecol. Syst. 1994, 25, 547–572. [Google Scholar] [CrossRef]
  9. Avise, J.C.; Neigel, J.E.; Arnold, J. Demographic influences on mitochondrial DNA lineage survivorship in animal populations. J. Mol. Evol. 1984, 20, 99–105. [Google Scholar] [CrossRef]
  10. Pasolini, P.; Ragazzini, C.; Zaccaro, Z.; Cariani, A.; Ferrara, G.; Gonzalez, E.G.; Landi, M.; Milano, I.; Stagioni, M.; Guarniero, I.; et al. Quaternary geographical sibling speciation and population structuring in the eastern Atlantic skates (Suborder Rajoidea) Raja clavata and R. straeleni. Mar. Biol. 2011, 158, 2173–2186. [Google Scholar] [CrossRef] [Green Version]
  11. Milá, B.; Van Tassell, J.L.; Calderón, J.A.; Rüber, L.; Zardoya, R. Cryptic lineage divergence in marine environments: Genetic differentiation at multiple spatial and temporal scales in the widespread intertidal goby Gobiosoma bosc. Ecol. Evol. 2017, 7, 5514–5523. [Google Scholar] [CrossRef] [Green Version]
  12. Bickford, D.; Lohman, D.J.; Sodhi, N.S.; Ng, P.K.; Meier, R.; Winker, K.; Ingram, K.K.; Das, I. Cryptic species as a window on diversity and conservation. Trends Ecol. Evol. 2007, 22, 148–155. [Google Scholar] [CrossRef] [PubMed]
  13. Mayr, E. Systematics and the Origin of Species from the Viewpoint of a Zoologist; Columbia University Press: New York, NY, USA, 1942; p. 372. ISBN 9780674862500. [Google Scholar]
  14. Knowlton, N. Sibling species in the sea. Annu. Rev. Ecol. Syst. 1993, 24, 189–216. [Google Scholar] [CrossRef]
  15. Nygren, A. Cryptic polychaete diversity: A Review. Zool. Scr. 2014, 43, 172–183. [Google Scholar] [CrossRef]
  16. Bellono, N.W.; Leitch, D.B.; Julius, D. Molecular tuning of electroreception in sharks and skates. Nature 2018, 558, 122–126. [Google Scholar] [CrossRef] [PubMed]
  17. McEachran, J.D.; Dunn, K.A. Phylogenetic analysis of skates, a morphologically conservative clade of Elasmobranchs (Chondrichthyes: Rajidae). Copeia 1998, 1998, 271–290. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar]
  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. [Google Scholar] [CrossRef] [Green Version]
  27. Ebert, D.A.; Compagno, L.J. Biodiversity and systematics of skates (Chondrichthyes: Rajiformes: Rajoidei). Biol. Skates 2007, 27, 5–18. [Google Scholar]
  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. [Google Scholar] [CrossRef]
  30. Grant, V. The systematic and geographical distribution of hawkmoth flowers in the temperate North American flora. Bot. Gaz. 1983, 144, 439–449. [Google Scholar] [CrossRef]
  31. Quintero, I.; Keil, P.; Jetz, W.; Crawford, F.W. Historical biogeography using species geographical ranges. Syst. Biol. 2015, 64, 1059–1073. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar]
  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. [Google Scholar] [CrossRef]
  42. Šantić, 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. [Google Scholar] [CrossRef]
  43. Tsikliras, A.C.; Stergiou, K.I. Age at maturity of Mediterranean marine fishes. Mediterr. Mar. Sci. 2015, 16, 5–20. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar]
  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. [Google Scholar] [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. [Google Scholar]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  53. Knowlton, N. Cryptic and sibling species among the decapod Crustacea. J. Crustac. Biol. 1986, 6, 356–363. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  64. Tamura, K.; Stecher, G.; Kumar, S. MEGA11: Molecular Evolutionary Genetics Analysis version 11. Mol. Biol. Evol. 2021, 38, 3022–3027. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar]
  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. [Google Scholar] [CrossRef] [Green Version]
  79. Yokes, M.B. DNA Barcoding of Marine Fish Species from Turkish Coastline. 2016; unpublished; submitted to the EMBL/GenBank/DDBJ databases. [Google Scholar]
  80. Nei, M. Molecular Evolutionary Genetics; Columbia University Press: New York, NY, USA, 1987. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [CrossRef]
  83. Akaike, H. A new look at the statistical model identification. Curr. Contents Eng. Technol. Appl. Sci. 1981, 12, 42. [Google Scholar]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [CrossRef]
  91. Rice, W.R. Analyzing tables of statistical tests. Evolution 1989, 43, 223. [Google Scholar] [CrossRef]
  92. Castro-Conde, I.; de Uña-Álvarez, J. Sgof: An R package for multiple testing problems. R J. 2014, 6, 96–113. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  102. Letunic, I.; Bork, P. Interactive Tree Of Life (ITOL): An online tool for phylogenetic tree display and annotation. Bioinformatics 2007, 23, 127–128. [Google Scholar] [CrossRef] [Green Version]
  103. Puillandre, N.; Lambert, A.; Brouillet, S.; Achaz, G. ABGD, Automatic Barcode Gap Discovery for primary species delimitation. Mol. Ecol. 2012, 21, 1864–1877. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar]
  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. [Google Scholar] [CrossRef]
  107. Sandoval-Castillo, J.; Beheregaray, L.B. Oceanographic heterogeneity influences an ecological radiation in elasmobranchs. J. Biogeogr. 2020, 47, 1599–1611. [Google Scholar] [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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [CrossRef] [Green Version]
  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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar] [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. [Google Scholar]
Animals 13 02139 g001 550
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 13 02139 g001
Animals 13 02139 g002 550
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.
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.
Animals 13 02139 g002
Animals 13 02139 g003 550
Figure 3. Bayesian admixture analysis among individuals belonging to the Raja miraletus species complex across most of its geographical distribution area. (a) Distribution of the cytochrome oxidase subunit I clades in the population samples inferred with BAPS; each colour represents one distinct haplogroup (cluster), and each bar represents a different individual. (b) Results of the Bayesian individual clustering using STRUCTURE results for K = 4; each vertical bar represents one individual, in which a different colour represents the estimated cluster membership. Acronyms of the Macro Area Codes are given in Table 1.
Figure 3. Bayesian admixture analysis among individuals belonging to the Raja miraletus species complex across most of its geographical distribution area. (a) Distribution of the cytochrome oxidase subunit I clades in the population samples inferred with BAPS; each colour represents one distinct haplogroup (cluster), and each bar represents a different individual. (b) Results of the Bayesian individual clustering using STRUCTURE results for K = 4; each vertical bar represents one individual, in which a different colour represents the estimated cluster membership. Acronyms of the Macro Area Codes are given in Table 1.
Animals 13 02139 g003
Animals 13 02139 g004 550
Figure 4. Bayesian coalescent tree summarising phylogenetic relationships between the Central–Southern African (CSA) and Northeastern Atlantic–Mediterranean (NEAM) COI lineages. Raja undulata was used as an outgroup (BOLD record ELAME177-09, NCBI accession numbers KT307412, KT307413, KT307414). Posterior probability (≥0.5) values are reported on nodes. Coloured bars near haplotype nodes refer to the genetic clusters identified by BAPS results (Figure 3a). The outcomes of species delimitation analyses using ABGD and bPTP methods are shown as vertical bars on the right.
Figure 4. Bayesian coalescent tree summarising phylogenetic relationships between the Central–Southern African (CSA) and Northeastern Atlantic–Mediterranean (NEAM) COI lineages. Raja undulata was used as an outgroup (BOLD record ELAME177-09, NCBI accession numbers KT307412, KT307413, KT307414). Posterior probability (≥0.5) values are reported on nodes. Coloured bars near haplotype nodes refer to the genetic clusters identified by BAPS results (Figure 3a). The outcomes of species delimitation analyses using ABGD and bPTP methods are shown as vertical bars on the right.
Animals 13 02139 g004
Table
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.
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 AreaMacro Area CodeN (tot)N (COI)N (EST-SSRs)YearsSource (Trawl Survey Program)McEachran et al., 1989 [48]
Central-Southern Africa (CSA)
South Africa—South CoastSAF8682006ST (Africana)5
South Africa—South Coast0502007, 2012GB/BOLD5
South Africa—South Coast3230312011ST (Africana)5
NamibiaNAM0302009, 2010GB/BOLD5
AngolaANG2727262006ST (Nansen)4
SenegalSEN5552007CF2
Northeastern Atlantic–Mediterranean Sea (NEAM)
PortugalPOR3332007ST (IPIMAR)n.a.
Portugal0702005, 2007GB/BOLDn.a.
Balearic IslandsBAL1919162006ST (MedITS)1
Balearic Islands0302013ST (MedITS)1
SardiniaSAR111182002, 2005ST (MedITS)1
TuscanyTUS2622212005, 2006ST (MedITS)1
Tuscany166132008, 2010ST (MedITS)1
AlgeriaALG8852002, 2003FM (Algiers)1
Algeria10982009, 2010FM (Algiers)1
Strait of Sicily—Adventura BankSIC2222222014ST (MedITS)1
Strait of Sicily—Maltese Bank161282000, 2002ST (MedITS)1
Strait of Sicily—Maltese Bank01102007GB/BOLD1
Adriatic Sea—Northern Italian coastADR3931202006, 2007ST (MedITS; GruND)1
Adriatic Sea—Croatian coast242482002, 2004ST (MedITS; GruND)1
Adriatic Sea—Southern Italian coast1916192004ST (MedITS; GruND)1
Adriatic Sea—Albanian coast1913172004ST (MedITS; GruND)1
Ionian Sea4342004ST (MedITS; GruND)1
Greece—Aegean coastGRE0302014GB/BOLD1
IsraelISR8772009CF1
Israel0302012BOLD1
Israel0402014GB/BOLD1
Levantine SeaLEV0202016GB/BOLD1
Levantine Sea7772009CF1
n.a.: not available. ST: Scientific Trawl survey. CF: Contracted Fishermen. FM: Fishery Market. GB: GenBank Database. BOLD: Barcoding of Life Database.
Table
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.
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 asteriasRaja brachyuraRaja clavataRaja microocellataRaja montaguiRaja polystigmaRaja radulaRaja straeleniRaja undulataRaja miraletus (CSA)Raja miraletus (NEAM)
Raja asterias0.0025 ± 0.00120.01390.01060.01320.01350.01210.00950.01110.01210.01550.0189
Raja brachyura0.08630.0030 ± 0.00160.01070.00940.01190.01140.01210.01250.01360.01410.0172
Raja clavata0.05910.05140.0038 ± 0.00000.01150.01030.01130.00700.00590.01410.01330.0169
Raja microocellata0.08770.04650.05980.0000 ± 0.00000.01250.01070.01140.01150.01270.01400.0172
Raja montagui0.08290.06060.05110.06380.0000 ± 0.00000.00660.01260.01160.01280.01320.0172
Raja polystigma0.07480.05260.05140.05960.02290.0000 ± 0.00000.01080.01120.01230.01190.0164
Raja radula0.04870.06850.02800.06480.06460.05640.0010 ± 0.00120.00710.01190.01370.0173
Raja straeleni0.05930.05930.01500.06010.05900.05580.02820.0019 ± 0.00100.01220.01310.0174
Raja undulata0.07670.07360.07330.07990.07900.07110.07420.07350.0009 ± 0.00100.01220.0170
Raja miraletus (CSA)0.10530.08570.08610.09470.08300.07780.09080.09070.07700.0183 ± 0.00290.0117
Raja miraletus (NEAM)0.10720.10200.09890.10070.09210.08970.09740.10330.09590.07330.0025 ± 0.0010
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.

Share and Cite

MDPI and ACS Style

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

AMA Style

Ferrari A, Crobe V, Cannas R, Leslie RW, Serena F, Stagioni M, Costa FO, 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(13):2139. https://doi.org/10.3390/ani13132139

Chicago/Turabian Style

Ferrari, Alice, Valentina Crobe, Rita Cannas, Rob W. Leslie, Fabrizio Serena, Marco Stagioni, Filipe O. Costa, Daniel Golani, Farid Hemida, Diana Zaera-Perez, and et al. 2023. "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 13, no. 13: 2139. https://doi.org/10.3390/ani13132139

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop