Introduction

Speciation, or the origin of species, was famously referred to as that “mystery of mysteries” (Darwin 1859; Cannon 1961) and understanding its process remains one of the greatest challenges in evolutionary biology. The biological species concept (BSC) regards species as “groups of interbreeding natural populations which are reproductively isolated from other such groups” (Mayr 1963, 1995). Geographic isolation is considered to play a significant role in the process of speciation, which in this case can be considered a dynamic continuum throughout which a number of stages of divergence may be observed (Mayr 1954). Where populations are only recently diverged they can appear morphologically indistinguishable. This is more likely to occur in taxa that demonstrate a high degree of morphological conservatism and do not rely on visual cues for reproductive interactions (Bickford et al. 2007).

Skates (Order Rajiformes, Class: Chondrichthyes) are characterised by a highly conserved morphology. The identification of species boundaries, fundamental for effective conservation management, is therefore challenging. This issue is reflected in market mis-labelling of commercial landings (Iglésias et al. 2010; Griffiths et al. 2013). Complementary molecular approaches can assist with defining these boundaries and verifying species identifications. In particular, the mitochondrial ‘barcoding’ gene, cytochrome c oxidase subunit I (COI), can reveal differentiation at fine taxonomic levels, yet usefully diverges little within species (Hebert et al. 2003). COI barcoding has recently proven a successful tool for validating skate identification in a number of fisheries (Spies et al. 2006; Serra-Pereira et al. 2011; Coulson et al. 2011; Mabragaña et al. 2011; Lynghammar et al. 2014).

Like many elasmobranchs worldwide, skate populations in Europe have experienced significant declines, and in some cases even local extirpation, owing to historic and current high levels of fishing pressure (Brander 1981; Dulvy et al. 2000). The reproductive strategy of skates is characterised by late maturity and low fecundity, which coupled with susceptibility to capture, makes them arguably some of the most vulnerable of all exploited marine fishes (Dulvy et al. 2014). Their commercial importance is often low compared to other demersal fisheries, resulting in historically poor management of this group (Dulvy et al. 2000). Despite this, catches remain significant, with the total weight of all landings of the most common European skate species amounting to 16,600 tonnes (FAO 2009). Reporting species-specific landings became European law in 2009 (EC No 43/2009) however, such protective legislation requires an unambiguous definition of a species, which formally encompasses all natural variation.

Raja maderensis (Lowe, 1838) is a poorly known endemic species, reportedly restricted to waters close to the island of Madeira, however a similar form has also been recorded from the Azores (Stehmann 2009). It is classified as ‘Data Deficient’ by the IUCN red list of threatened species (Stehmann 2009), which considers factors such as geographic extent of range and abundance to determine extinction risk. Its vulnerability is therefore currently unknown and more data concerning its biology and ecology is required to establish responsible management. In particular, the taxonomic status of skates identified as R. maderensis in the Azores needs to be resolved, especially as recent molecular evidence suggested little difference between R. clavata and a specimen identified as R. maderensis landed in a Portuguese fish market (Serra-Pereira et al. 2011). The distinct dorsal pattern (Fig. S1) of R. maderensis reportedly distinguishes it from the morphologically highly similar thornback ray (R. clavata). However, the thornback ray is noticeably variable with regard to its general morphology when compared to many other European rajiids, with polychromatism recorded in the Mediterranean (Mnasri et al. 2009) and a number of recognised patterns relatively widely known in the UK (Pritchard 1977), as well as documented cases of albinism (Ball et al. 2013). It also occupies a much greater geographical range, particularly in terms of latitude, than many other skates. The thornback ray is distributed from Norway and Iceland to Northwest Africa, including the Mediterranean and Black Seas, and until recently was considered to inhabit waters as far south as southern Africa (Froese and Pauly 2011). The latter records most likely represented mis-identifications of the recently diverged and morphologically similar biscuit skate, Raja straeleni, distributed in southern Africa (Pasolini et al. 2011). Similarly, the rough ray (Raja radula), restricted to the Mediterranean (Morey et al. 2009) is regarded as a sister species to R. clavata. This species is considered to have split from the main lineage (Valsecchi et al. 2005), following re-colonisation of the western Mediterranean basin by R. clavata after the Late Miocene Messinian salinity crisis, around 5 Myr ago (Rogl 1998; Duggen et al. 2003). The question of whether R. maderensis represents a distinct species, or merely another morphological variant of the polytypic thornback ray, therefore remains unanswered. Clearly, there is a need to sample R. clavata and R. maderensis across their geographical range and assess their genetic relatedness in comparison to other closely related species, such as R. straeleni and R. radula.

DNA barcoding using COI has proven successful in distinguishing commonly captured species of European skate, but showed little difference between R. clavata and a single specimen of R. maderensis (Serra-Pereira et al. 2011). This suggests the need for further investigation using additional samples and a more variable genetic marker that affords better resolution. The mitochondrial control region (CR) has been shown to have a higher degree of nucleotide polymorphism in skate species (Valesecchi et al. 2005) than COI. In the present study both markers were used to determine the genetic distinctiveness of specimens identified as R. maderensis and R. clavata over a considerable portion of their geographic range. Their taxonomic and evolutionary status (based on the COI gene) was also examined in the context of 67 other Atlantic species of the family Rajidae.

Materails and methods

Sample collection

The majority of tissue samples were collected between June 2006 and December 2013 from a range of scientific research cruises throughout Europe and the Mediterranean (Fig. 1 and summarised in Table S1, S2). A minority of tissue samples, most notably all those from Portugal and South Africa, were collected at landing or at fish markets.

Fig. 1
figure 1

Map of sampling locations for R. clavata, R. maderensis and R. radula (samples of R. straeleni collected from a fish market in Cape Town, South Africa are not shown)

Full size image

Photographs of specimens were not collected with tissue samples. However, sampling was undertaken by experienced researchers following detailed biological descriptions. Samples of R. maderensis collected by the Universidade dos Açores and by Instituto Português do Mar e da Atmosfera were identified following Stehmann and Burkel (1984). Photographs of reference specimens are included in Fig. S1. In particular, the lack of dark crossbars on the tail and the smoothness of the skin were considered distinguishing features of R. maderensis, relative to the common R. clavata caught off the Azores (G. Menezes, pers. comm).

DNA sequencing

Genomic DNA was isolated using the ‘Wizard’ protocol (Promega Madison, WI). Approximately 1000 bp of the mitochondrial control region was amplified using primers RayDloopFor (5′-CATTAATCGACTRTCAACTATTTCATT-3′) and Ray12s (5′-TACTGAGGCTAGGACCAAAC-3′), following Griffiths et al. (2010). A subset of individuals, with particular focus on R. maderensis and R. clavata from around the Azores and Northern Africa, were also amplified for COI (the ‘DNA barcode’) using the primers FishF2 (5′-TCGACTAATCATAAAGATATCGGCAC-3′) and FishR2 (5′-ACTTCAGGGTGACCGAAGAATCAGAA-3′) (Ward et al. 2005). Reaction conditions followed Serra-Pereira et al. (2011). PCR products were sequenced by Macrogen Inc (Korea) and Source Bioscience (UK), and results checked by eye in BIOEDIT version 7.1.11 (Hall 1999) for errors and ambiguous base calls.

Control region analysis: divergence between species

Sequences of the CR from R. maderensis (n = 37), R. clavata (n = 196), R. straeleni (n = 6) and R. radula (n = 3) were aligned in MEGA version 5 (Tamura et al. 2011). Nucleotide diversity indices: number of segregating sites (S), number of haplotypes (Nh), haplotype diversity (Hd) and nucleotide diversity (p) were calculated using DnaSP version 5.10 (Rozas et al. 2003). The evolutionary divergence between pairs of species was estimated in Mega 5 using the K80 + I substitution model (Kimura 1980), determined as appropriate by JModelTest (ΔAICc = 3.86, Darriba et al. 2012). Standard error estimates were obtained by bootstrapping (1000 replicates). Haplotype genealogy was reconstructed from trees using a maximum likelihood (ML) approach with the above substitution model implemented in PHYML 3.0 (Guindon and Gascuel 2003) and visualised using Haploviewer (Salzburger et al. 2011). To assist with interpretation, samples were labelled on the haplotype network according to species. In the case of R. clavata, they were additionally grouped by broad geographic sampling region based on previous findings (Chevolot et al. 2006): Northern north east Atlantic; Portuguese coast; Mediterranean and Azores.

Control region analysis: differentiation among sampling regions

To compare genetic differentiation between geographically distinct sample collections of R. clavata (grouped as above) with that observed between species (R. maderensis and R. clavata), the pairwise fixation index ΦST (which takes into account nucleotide diversity) was computed using the Kimura 2-P distance model in Arlequin 3.5 (Excoffier et al. 2005). In addition, a hierarchical analysis of molecular variance (AMOVA) was computed using the same model in Arlequin 3.5 (Excoffier et al. 2005). Using this analysis, total variance is divided into variance components representing: (1) differences among groups (ΦCT), (2) differences among populations within groups (ΦSC) and (3) differences within populations (ΦST). The best geographic/species subdivisions are assumed to have maximum among group variance (ΦCT values), with a low within-group component (ΦSC), which significantly differs from a random distribution. A number of hypothesised groupings of samples were therefore considered in the AMOVA (summarised in Table 1) with samples identified as R. maderensis along with the a priori geographic sampling regions of R. clavata used as a basis for comparison.

Table 1 Analysis of molecular variance (AMOVA) for five hypothesised groupings comprising Raja maderensis and four a priori geographic sampling regions of Raja clavata
Full size table

Cytochrome oxidase I analysis

Sequences produced from tissue samples of R. maderensis (n = 37), R. straeleni (n = 6) and a subset of 30 R. clavata were aligned alongside those from the closely related R. radula downloaded from BOLD database. A haplotype genealogy was constructed from trees using a ML approach. A HKY (Hasegawa et al. 1985) nucleotide substitution model determined by JModelTest (ΔAICc = 1.42, Darriba et al. 2012) was implemented in PHYML 3.0, and genealogy was visualised using Haploviewer (Salzburger et al. 2011).

To build a comparative phylogeny, sequences of 65 additional Atlantic skate species were downloaded from the Barcode of Life (BOLD) and Genbank databases. Only sequences originating from the following specialist publications were utilised to minimise inaccuracy with species identification: Spies et al. (2006), Serra-Pereira et al. (2011), Coulson et al. (2011), Mabragaña et al. (2011) and Lynghammar et al. (2014). Additional sequences from the starry ray (R. asterias, a Mediterranean endemic species) were also downloaded, as this species is considered closely related to R. clavata (summarised in Table S3). Evolutionary divergence between pairs of species and genera was estimated using the TrN + I + Γ substitution model (Tamura and Nei 1993). Standard error estimates were obtained by bootstrapping (1000 replicates) in Mega 5. The number of sequences was subsequently collapsed so only unique haplotypes were represented using online fasta toolbox software FaBox (Villesen 2007). Phylogenetic relationships were then inferred using maximum-likelihood (ML) tree methods with the HKY + I + Γ model (Hasegawa et al. 1985) This was determined as the most appropriate by JModeltest (ΔAICc = 35.11) and implemented in PhyML 3.0 with 1000 bootstrap replications to establish the robustness of internal tree branches. The tree was rooted using sequences of ancestral chimaera species: Chimaera monstrosa, Chimaera opalescens and Chimaera fulva.

Results

Control region analysis: divergence between species

A total of 242 samples of R. clavata, R. maderensis, R. straeleni and R. radula were sequenced for the CR (Table S1). According to the diversity measures calculated, R. clavata and R. maderensis showed comparable haplotype diversity (Table 2). Of the four haplotypes identified for R. maderensis, three were shared with R. clavata and the fourth, which was unique, differed only by a single base. The estimated evolutionary distance between R. clavata and R. maderensis was very low (0.0076 ± 0.0029). In contrast, all haplotypes for R. straeleni and R. radula were species-specific. The calculated genetic distance between R. clavata and R. straeleni was 0.018 (±0.0046) and between R. clavata and R. radula, 0.023 (±0.0052, Table 3).

Table 2 Summary of genetic parameters for the mtDNA control region by species
Full size table
Table 3 Estimates of evolutionary divergence between species based on mtDNA control region sequences
Full size table

The CR haplotype network revealed five main haplogroups, with the vast majority of samples found in two of these haplogroups (Fig. 2). The first haplogroup (group A) contained the most common R. clavata haplotype (68 % of individuals), which was found in every region sampled, except the Azores. The second haplogroup (group B) contained the second most common R. clavata haplotype (6 % of individuals), found only in the Azores and the west coast of Ireland (one individual). This haplotype also represented the majority of individuals identified as R. maderensis (57 %) and all haplotypes of R. maderensis were contained within this haplogroup. Haplogroups A and B were only separated by two stepwise mutations. An additional haplogroup of R. clavata was also identified that included nine individuals from three different sampling locations. The individuals in this group, which was separated by nine stepwise mutations from the next nearest R. clavata haplotype (group E), originated from the southern North Sea, the west coast of Scotland and the Portuguese coast. The remaining two haplogroups represented specimens of R. straeleni (group C) and R. radula (group D), separated by nine and 14 mutational steps, respectively, from the nearest R. clavata haplotype.

Fig. 2
figure 2

Haplotype genealogy based on mtDNA control region sequences for R. clavata, R. maderensis, R. straeleni and R. radula. Samples are patterned according to species and those from R. clavata are also labelled according to geographic region. Haplotype size is proportional to the number of individuals represented. Solid black dots represent mutational steps, dashed boxes represent separate haplogroups

Full size image

Control region analysis: differentiation among sampling regions

Little genetic divergence was evident between the Northern north east Atlantic and the Portuguese coast (pairwise ΦST = 0.061, P = 0.026), or the Mediterranean (ΦST = 0.025, P = 0.055, Table 4). The Portuguese coast and the Mediterranean were more strongly differentiated (ΦST = 0.18, P = < 0.001). However, pairwise comparisons of Azores R. clavata samples with other sample sites demonstrated a highly significant degree of separation (ΦST = 0.64–0.95, P < 0.0001). More importantly, samples identified as R. maderensis were not so strongly differentiated from proximate Azorean R. clavata, (ΦST = 0.21, P = 0.02, Fig. 3) despite being similarly differentiated from Atlantic shelf and Mediterranean groups. Overall the AMOVA results best supported the existence of only two groups, with among group variation constituting 73.9 % of the total (P = 0.02). The first group was represented by R. clavata distributed across the northern north east Atlantic, the Portuguese coast and the Mediterranean; and the second group was represented by R. clavata collected from the Azores and individuals identified as R. maderensis (Table 1).

Table 4 Phi (Φ)ST for four sampling regions of R. clavata (NNEA Northern North East Atlantic, PC Portuguese coast, MED Mediterranean Sea), contrasted with R. maderensis (MAD)
Full size table
Fig. 3
figure 3

Pairwise ΦST (Phi) plotted against ranked ΦST for population and species comparisons of R. clavata and R. maderensis established from sequence variation at the mtDNA control region. Cluster A = contrasts: MAD ~ AZ, MED ~ PC, MED ~ NNEA, NNEA ~ PC, Cluster B = contrasts: MAD ~ PC, MAD ~ NNEA, MAD ~ MED, AZ ~ PC, AZ ~ NNEA, AZ ~ MED where MAD = individuals identified as Raja maderensis and for geographic sampling regions of R. clavata: NNEA Northern north east Atlantic, PC Portuguese coast, MED Mediterranean Sea, AZ Azores Islands

Full size image

Cytochrome oxidase I analysis

Sequences of the COI gene for R. clavata (n = 30) and R. maderensis (n = 38) revealed substantial overlap, with haplotypes shared between species. Only one specimen identified as R. maderensis had a unique haplotype, separated by a single mutation from R. clavata (Fig. 4). Specimens of R. straeleni and R. radula were clearly separated by six and 15 stepwise mutations, respectively from the nearest R. clavata haplotype.

Fig. 4
figure 4

Haplotype genealogy based on mtDNA COI region sequences for R. clavata, R. maderensis, R. straeleni and R. radula. Haplotype size is proportional to the number of individuals represented. Solid black dots represent mutational steps

Full size image

Considering sequences of 69 species from the Atlantic, the overall mean sequence divergence between species was 0.157 ± 0.035, but within the genus Raja it was only 0.1045 ± 0.0642. Samples from closely related species R. clavata, R. maderensis and R. straeleni showed high levels of similarity (Table 5). In particular, among specimens of R. maderensis and R. clavata, divergence ranged from 0.00 to 0.016. Average divergence between species within genera ranged considerably, with the highest levels seen in the Leucoraja (0.1094 ± 0.0165, Table S4).

Table 5 Estimates of evolutionary divergence between species based on mtDNA COI sequences calculated by Tamura-Nei nucleotide substitution model
Full size table

The phylogenetic analysis produced a tree topology supported by high bootstrap values. In particular the main evolutionary lineages (families Rajidae, Blainville 1816 and Arhynchobatidae, Fowler 1934) were clearly identifiable (Ebert and Compagno 2007). Within Rajidae, clades with strong support were the Amblyrajini, Rajini, and a third containing Raja eglanteria and Rostroraja alba of the Rostrorajini (Chiquillo et al. 2014), as well as genera Neoraja and Malacoraja of the Gurgesiellini (McEachran and Dunn 1998). Within the Rajini, the genus Raja was not monophyletic, although there was strong support for reciprocal monophyly of most species. In particular, R. radula, R. straeleni and R. clavata were clearly separated, while R. maderensis and R. clavata had overlapping haplotypes (Fig. 5).

Fig. 5
figure 5

Maximum likelihood tree using a HKY substitution model based on a common 556 bp COI sequence from 69 skate species of the Atlantic. Only bootstrap values (based on 1000 replicates) >70 % are shown. Branches where species formed monophyletic groups were collapsed and the shape of the triangle represents the number of haplotypes (vertically) and their diversity (horizontally). Individuals highlighted as mis-identifications in the original studies were relabelled accordingly. If identities were not resolved, they were removed from the analysis. Star denotes area of the tree detailing R. clavata and R. maderensis

Full size image

Discussion

One of the key issues in evolutionary biology lies in identifying the patterns and causes of reproductive isolation and speciation among closely-related taxa (Mayr 1963; Coyne and Orr 2004). Skates represent a particularly speciose group constituting nearly a quarter of all known chondrichthyan fishes (Ebert and Compagno 2007). Perhaps somewhat paradoxically, they are extremely morphologically conserved, maintaining a dorso-ventrally flattened body form typical of all batoids and brown-grey coloration. Despite evidence of speciation in the form of accumulated genetic differences (Pasolini et al. 2011), species identification remains challenging, potentially due to the selection of a well-adapted phenotype (Futuyma 1979). In this study it was considered whether the endemic R. maderensis truly denotes an individual species or whether it represents a morphotype of the polytypic and widely distributed R. clavata. Using molecular approaches, there is little evidence to suggest they are distinct species, but R. maderensis may constitute part of a diverging population.

Control region analysis

Based on sequence variation of the CR, there were no fixed differences between R. clavata and R. maderensis. Only a single R. maderensis specimen possessed a unique haplotype, and this differed from a R. clavata haplotype by a single base. This sharply contrasts with R. straeleni and R. radula, where species boundaries were clearly defined in the haplotype network and numerous mutational steps separated these groups. In fact, the separation of clades for R. clavata and R. maderensis appears to more closely follow geographic origin of samples, rather than species identity. This meant individuals identified as R. maderensis actually grouped with geographically proximate R. clavata captured around the Azores. Conversely, comparisons of population sub-division between R. clavata from the Azores or R. maderensis to continental shelf or Mediterranean groups, demonstrated concordant high levels of divergence. If R. maderensis has been reproductively isolated from R. clavata for an extended period, it would generally be expected that levels of divergence between the species would greatly exceed those observed between populations from within the same species (reviewed in Mallet 1995). This is clearly not the case and data suggest that individuals inhabiting areas around islands and seamounts in the eastern Atlantic (Macaronesia) potentially represent a distinct demographic unit, which is isolated from populations inhabiting the continental shelf of Europe. The low nuclear allelic diversity in R. clavata from the Azores region and its high genetic differentiation is certainly consistent with a strong bottleneck and physical isolation (Chevolot et al. 2006). Similarly, Naylor et al. (2012) also observed genetic differences between R. clavata from the Azores and the continental shelf of mainland Europe. A number of studies considering the population genetic structure of European fish species illustrate the distinctiveness of this region, such as the white sea bream, Diplodus sargus, (González-Wangüemert et al. 2010) and the shanny, Lipophrys pholis (Stefanni et al. 2006). This is perhaps not surprising when considering the Azores islands represent the most remote archipelago of the North Atlantic, some 1300 km from Portugal suggesting both depth and oceanic distance are significant dispersal barriers for a variety of fish species. This is likely to be even more of an obstacle in demersal species such as skate that lack larval dispersal and have a restricted depth distribution.

Significantly, the haplotype most commonly associated with R. maderensis was also isolated in a small number of R. clavata from Ireland (with other highly similar haplotypes from within the same clade also identified in R. clavata from around the North Sea, Scotland, Ireland and Portugal). Previously, it has been suggested that the Azores represented a refugium for R. clavata during the last glacial maxima (Chevolot et al. 2006) and that waters around the British Isles were directly re-colonised by fish from this region. Our results demonstrate that these areas share haplotypes, yet the exact route of re-colonisation is unclear, given that very similar haplotypes are represented along the African seamounts and Portuguese coast.

A highly divergent group of CR haplotypes was also identified in nine of the specimens of R. clavata from across the continental shelf, specifically the southern North Sea; western Scotland; the Irish Sea and the Portuguese coast. Given the large number of stepwise mutations that separate these from other haplotypes, the possibility of mis-identification was considered. This seems unlikely as the closest known species to R. clavata, R. straeleni and R. radula are also represented in the network. Other species that could have been mistaken for R. clavata in British waters include Raja brachyura and Raja montagui. These species, along with representatives of all the commonly captured species of skate around the UK, have been previously sequenced at the same CR region and R. clavata remains the closest match in terms of sequence homology (Griffiths et al. 2010). Hybridisation can also be ruled out for the same reason, as the haplotypes do not appear to originate from any other contemporary sympatric species of skate. No photographs accompany these tissue samples, but cryptic speciation is also doubtful as five of the individuals were also genotyped at a suite of microsatellite loci and all had alleles falling within the range identified for R. clavata (Griffiths, unpublished data). Interestingly, the majority of the mutations accounting for this group are located within a 50 bp region (~340–390 bp), suggesting that these changes could have arisen through a single gross mutational event. Alternatively, this group of haplotypes could have arisen within an additional, unidentified refugium, where R. clavata persisted during the last glacial maxima. There is growing evidence from a number of marine species that many organisms survived in isolated populations much closer to the ice sheets than previously thought (Provan et al. 2005; Hoarau et al. 2007; Olsen et al. 2010; Finnegan et al. 2013). This has been suggested in previous studies of R. clavata to explain the diversity of haplotypes present in northern Europe (Chevolot et al. 2006), putatively suggesting Hurd Deep (i.e. the region of the western English Channel) as the refugium. However, the sequences identified in the current study remain highly divergent (separated from the most common R. clavata haplotypes by a similar genetic distance as R. radula or R. straeleni), they are extremely rare and do not appear to be exclusive to any specific region of the continental shelf that was sampled.

Cytochrome oxidase I analysis

The numbers of COI haplotypes identified in R. maderensis and R. clavata were low (2 and 6 respectively), with a high frequency shared haplotype common to both species in the Azores and surrounding seamounts. This essentially supports the results of the CR sequencing, that individuals identified as R. maderensis and R. clavata from around Macaronesia represent a single demographic unit, which is distinct from R. clavata inhabiting areas of the European continental shelf and Mediterranean.

The substantial phylogeny reconstructed here reveals that R. clavata and R. maderensis do not form monophyletic groups, yet the closely related R. straeleni and R. radula both branch with strong bootstrap support (81 and 86 % respectively), with R. straeleni demonstrating the closest relationship with the R. clavata and R. maderensis. This analysis extends the previous work of Chiquillo et al. (2014) by including genera of both hard-nosed (upper clade) and soft-nosed (lower clade) skates and supports the polyphyletic nature of the genus Raja. In particular, Raja eglanteria certainly requires re-arrangement as it aligns closely with Rostroraja alba, forming part of a newly proposed tribe: Rostrorajini (Chiquillo et al. 2014). These taxonomic issues highlight how understudied the phylogenetics of batoids remains and the tree presented here represents one of the most comprehensive phylogenetic reconstructions of skates to date.

The COI gene has proven to be a useful species diagnostic tool owing to its diversity between, but conservation within, species (Hebert et al. 2003). As a result, numerous attempts have been made at establishing boundaries which usefully define species. An extensive analysis of chordates demonstrated that 98 % of species pairs were separated by a COI gene sequence divergence greater than 2 %, although this work also highlighted that the appropriate cut-off should be calibrated with known species of the taxon under study (Hebert et al. 2003). A detailed analysis of fish suggested that 3.5 % divergence was appropriate for the group (Ward et al. 2009), however rates of nucleotide substitution have been shown to be lower in elasmobranchs (Martin et al. 1992). Using a sample set of 69 skate species from the Atlantic, the average genetic distance between species was 15.7 % and within the genus Raja, 10.4 % (or 8.6 % if R. eglanteria was not included). This sharply contrasts with that calculated between R. maderensis and R. clavata (0.4 %). Meier et al. (2008) considered that the lowest distance to the nearest neighbour might be a more appropriate cut-off. In this case, comparison with the levels of divergence between R. clavata and its close relatives R. straeleni and R. radula, may be informative. We calculated genetic distances of 1.1 and 2.2 % between these species. Therefore, the distance between R. maderensis and R. clavata is considerably below any of the species cut-offs considered, even the most conservative and taxonomically appropriate. Nevertheless, such cut-offs remain contentious as taxonomic groups exist on an evolutionary continuum, and are unlikely to be fixed (Trewick 2008).

Distinctiveness of Raja maderensis

The combined analysis of both CR and COI sequences failed to support the genetic distinctiveness of R. maderensis. New animal species most commonly evolve through allopatric speciation, which results from geographic isolation, beginning with highly structured populations (Coyne and Orr 2004; Gavrilets 2004). This had been suggested as the predominant mode of evolution in European skates (Pasolini et al. 2011). Raja clavata and R. straeleni are thought to have separated during the early Pleistocene as a result of glacial activity in equatorial Africa around 1.3 Ma. Using the highly variable CR, we calculated a mean sequence divergence of 1.8 % between these two species, which equals that demonstrated by Pasolini et al. (2011). Similarly, the sequence divergence we calculated between R. clavata and the endemic Mediterranean skate R. radula (2.3 %) compares well to previous estimates used to date the evolution of R. radula to less than 5 million years ago (2.9 %, estimated from data from Valsecchi et al. 2005). This result supports the possibility that the independent demographic unit comprising specimens identified as R. maderensis and geographically proximate Azorean R. clavata represents an earlier stage of such a process of divergence.

The lack of genetic distinctiveness of R. maderensis is difficult to reconcile with the morphological differences ascribed to it (Stehmann and Burkel 1984; Ebert and Stehmann 2013). Perhaps this distinct morphology corresponds to examples of local adaptation or phenotypic plasticity in response to varying environmental conditions, indeed R. clavata is known to exhibit considerable variation in coloration and skin texture across its distribution in the UK and Mediterranean Sea (Pritchard 1977; Mnasri et al. 2009). It is acknowledged that difficulties distinguishing between R. clavata and R. maderensis may have potentially led to misidentification, which would contribute to apparent overlap of shared haplotypes between species. Despite this, wholescale misidentification seems very unlikely given the number of samples analysed and the experience of the researchers. Morphological identification remains challenging, although appropriate taxonomic guides were utilised (Stehmann and Burkel 1984; Ebert and Stehmann 2013), suggesting the conclusions are robust. Further biological and ecological information concerning R. maderensis is required to assess the taxonomic validity of this group. However, this work is the result of many years of focused sample collection, and represents a substantial effort to collect data from remote and inaccessible areas of Europe and Africa. This difficulty is also reflected in the fact that no samples were collected from the coastal waters of Madeira, a key part of the described range of R. maderensis. The region remains difficult to trawl due to its rough topography, but samples were collected in near-by waters, approximately 100 nautical miles north east of the archipelago. Indeed, further sampling along the coast of north Africa is not only important for assessing the genetic distinctiveness of R. maderensis and establishing the distribution of R. clavata and R. straeleni, but also the elasmobranch fauna of the whole region generally remains rather poorly described (Bonfil 1994; Compagno 1994; Kroese and Sauer 1998).

A recent review by Dulvy et al. (2014) highlights that most chondrichthyan catches worldwide are unregulated, misidentified, unrecorded, aggregated or discarded at sea, resulting in a lack of species specific landing information. Skates are particularly challenging to identify and a lack of species composition monitoring has masked severe population declines (Brander 1981; Walker and Hislop 1998; Iglésias et al. 2010). They are slow growing, late maturing fishes with a relatively low fecundity; a set of life history traits which renders them highly susceptible to over-exploitation (Dulvy and Reynolds 2002). Endemic species with restricted ranges are of particular concern and there is a growing need to gather more information pertaining to their biology and ecology to assess their risk of extinction, used as a factor in determining IUCN classification. Therefore, determining the taxonomic status of R. maderensis has important implications for conservation and management. The application of genetic markers to this question seems especially appropriate, given that this approach has been highly successful in unambiguously distinguishing skates across Europe and elucidating taxonomic issues in the group (Spies et al. 2006; Griffiths et al. 2010; Iglésias et al. 2010; Serra-Pereira et al. 2011; Coulson et al. 2011; Mabragaña et al. 2011; Lynghammar et al. 2014). In terms of IUCN assessment this study highlights polytypic nature of R. clavata, further complicating its assessment for red-listing. While these results do not appear to support the occurrence of R. maderensis as a distinct species, our extensive sampling suggests that if it does exist it must have an extremely restricted range and small population size, rendering it extremely vulnerable to extinction.

Regardless of their morphological differences, the individuals collected from across the Azores and seamounts proximal to North Africa appear to belong to a single demographic unit, distinct from other populations on the European continental shelf and Mediterranean. However, our analysis of the CR revealed these groups were not entirely reciprocally monophyletic, which remains an important criterion in defining Evolutionary Significant Units (reviewed in Moritz 1994). Despite this, previous analysis with nuclear markers has also revealed the genetic distinctiveness of Azorean R. clavata (Chevolot et al. 2006). Recognising the evolutionary significance and conservation value of this group, would certainly serve to protect the unique morphology associated with R. maderensis, regardless of its underlying origin.