Puzzling over spurdogs: molecular taxonomy assessment of the Squalus species in the Strait of Sicily

The actual occurrence of Squalus megalops in the Mediterranean Sea has recently been questioned. Several research works which sought to assess available morphological and meristic features that differentiate S. megalops from other Squalus species in the Mediterranean Sea, revealed poor discriminatory power and high variability of the assessed characters, especially when comparing S. megalops and S. blainville . The application of molecular tools does not support the presence of S. megalops . In the present study, we screened spurdog species from the Strait of Sicily using a molecular taxonomy approach based on two mitochondrial DNA markers and we report the occurrence of two Squalus lineages characterizing specimens collected from the stretch of sea between Tunisia, southern Sicily, Malta and Libya. The results support the hypothesis that a common species, S. blainville , currently inhabits the Mediterranean Sea, while a second and rare species is probably an occasional visitor with high morphological similarity to the S. megalops and S. blainville but is genetically distinct from both. Within this perspective, the occurrence of S. megalops in the Mediterranean Sea is not confirmed and our study highlights the taxonomic uncertainties in relation to the occurrence and distribution of Squalus species in this region. We encourage the establishment of a coordinated international effort to implement a comprehensive and integrated taxonomic assessment on this genus which represents an irreplaceable component of the biodiversity of the area.


Introduction
The intrinsic low variation of morphological characters specific to elasmobranchs hinders their taxonomic identification at the species level and consequently undermines their conservation at different geographical scales (McEachran & Dunn 1998;Quattro et al. 2006;Aschliman et al. 2012). The lack of well-preserved holotypes for many shark species (e.g. Centrophorus spp.), misidentifications in databases and in the literature, and challenges in retrieving representative series of specimens for comparison are top-down impediments to the proper taxonomic identification and the potential revision of genera (Veríssimo et al. 2014).
In particular, the genus Squalus Linnaeus, 1758 is distributed worldwide (Ebert & Stehmann 2013) and includes about 26 different species of longlived sharks (Viana et al. 2016) inhabiting the waters of the continental shelf and upper slope, between 300-700 m of depth (Whitehead et al. 1984;Cannizzaro et al. 1995;Serena et al. 2009), as well as some seamounts and the waters around oceanic islands (Compagno 1984;Veríssimo et al. 2017). The genus is highly represented in bycatch and sev-eral studies have focused on the mitigation of the fishery impact on this group (Mandelman & Farrington 2007;Stoner & Kaimmer 2008;Tallack & Mandelman 2009). In addition, most spurdog shark species are included in the IUCN Red List of Threatened Species and are currently classified in categories from Data Deficient to Endangered (i.e. the Greeneye spurdog S. chloroculus Last et al. 2007;IUCN 2019). For this reason, highlighting the need to increase our knowledge about taxonomy, biology and ecology of these species appears now essential. Difficulties due to the problems of distinguishing between morphologically similar species and to the lack of effective and user-friendly identification field guides have been commented upon by many authors (Garrick 1960;Compagno 1984;Munoz-Chapuli & Ramos 1989;Veríssimo et al. 2017). With the progressive expansion of molecular taxonomy (i.e. DNA barcoding) and its integration with morphological approaches, unravelling the taxonomy of confounding groups of sharks has become a priority for their conservation (Ebert et al. 2010;Veríssimo et al. 2014;Pfleger et al. 2018). In more recent times, the well-established DNA barcoding technique based on the mitochondrial DNA (mtDNA) cytochrome c oxidase subunit I (COI) has been coupled with analysis of faster-evolving mitochondrial genes or long ribosomal DNA to improve the resolution and support of inferred phylogenetic relationship, up to the description of the evolutionary history of species (Avise 2004;Naylor et al. 2012;Krehenwinkel et al. 2019).
Recently, Veríssimo et al. (2017) reported uncertainties in the identification of Squalus specimens caught along the eastern Atlantic Ocean and in the Mediterranean Sea. In the latter area, Squalus acanthias Linnaeus, 1758 and Squalus blainville (Risso, 1827) are considered as resident species, while the presence of S. megalops (MacLeay, 1881) has been reported from Tunisia (Marouani et al. 2012). From the latter region, several specimens were analysed with a multidisciplinary approach and significant differences were reported for the two identified groups defined respectively as S. blainville and S. megalops (Marouani et al. 2012). Sequence data from both identified species were included later in the comprehensive assessment of the Squalus genus carried out in Veríssimo et al. (2017). Using two mtDNA markers, three main groups and four main Clades corresponding to S. acanthias (Clade A), S. blainville/S. megalops (Clade B), S. megalops (Clade C) and Squalus mitsukurii Jordan & Snyder, 1903 (Clade D) were identified in the Eastern Atlantic and the Mediterranean Sea (Veríssimo et al. 2017). Individuals classified under Clade C, which is highly divergent from both S. megalops from Australian waters and Clade B, originates from tropical West African coasts, with the exception of the single specimen from Tunisia, previously identified as S. blainville by Marouani et al. (2012). The sequence data associated with the individual from Tunisian waters described as S. megalops in Marouani et al. (2012) fitted in Clade B. Given the above described results and since specimens classified as potential S. megalops are included in both Clade B and Clade C, these findings further support current inconsistencies in species identification within the genus Squalus and the need for an accurate redescription of Squalus species, especially in the Mediterranean Sea, to stabilize the systematic and facilitate specimens identification.
Extensive studies conducted in other areas of the Mediterranean Sea considering both morphological (chondrocranium and other body traits) and genetic (COI sequences) analyses, revealed the presence of only one spurdog species, S. blainville, in the Ionian, Lybian, Aegean Seas and in the Sardinian waters (Kousteni et al. 2016;Bellodi et al. 2018). These findings spotlighted the stretch of sea between Tunisia, southern Sicily, Malta and Libya, known as the Strait of Sicily, as a more interesting area for spurdog species. Differently from Marouani et al. (2012), Bonello et al. (2016) did not identify diagnostic features (e.g. dermal denticles) which could effectively discriminate between S. blainville and S. megalops. The presence of S. blainville in the Maltese waters was assessed through the use of the DNA barcoding approach (Bonello et al. 2016). In the same region, Vella et al. (2017) collected and analysed individuals belonging to the nominal S. blainville and genetically clustering within Clade B (sensu Veríssimo et al. 2017), while three individuals were classified as Squalus sp. by the authors as clustering in the genetic Clade C (sensu Veríssimo et al. 2017).
In the present study, we screened all spurdogs caught in the Strait of Sicily between 2016 and 2017 to investigate the species composition in this stretch of sea. Comparing new and publicly available mtDNA gene sequences (COI and NADH dehydrogenase subunit 2; NADH2) we investigated the possible co-occurrence of two or more species of spurdog within the study area. Furthermore, since previous studies revealed that nominal species characterized by a wide geographical distribution share mitochondrial lineages, we evaluated the relationship between the genetic lineages retrieved in the Strait of Sicily and other Clades previously identified within the genus Squalus.

182
A. Ferrari et al.  Bonello et al. 2016) and ten specimens collected off the Tuscany coasts (GFCM-GSA 9) in 2017 were included in this study. All specimens were preserved at −20°C until the collection of main biological data and tissues samples for genetic analyses. In particular, individual muscle tissue or fin clips were preserved in 96% ethanol for laboratory analyses (see Table S1 for details on the field species assignment and geographical origin of specimens). Total genomic DNA (gDNA) was extracted from approximately 20 mg of tissue using the Wizard® SV Genomic DNA Purification System by Promega, according to the manufacturer's instructions. The quality of the extracted gDNA was assessed on a 0.8% agarose gel electrophoresis.
Similarly, a fragment of the NADH2 (~1100 bp) gene was amplified for all individuals using the Met-F and Trp-R primers (Vella et al. 2017). PCR conditions consisted of an initial denaturation of 2 min at 95°C, followed by 28 cycles of 45 s at 95°C, 45 s at 54°C and 60 s at 72°C and a final extension step for 7 min at 72°C. PCR outcomes were evaluated on a 2% agarose gel and preserved at 4°C until purification. All PCR products were purified using the ExoSAP-IT™ Express PCR Product Cleanup Reagent (ThermoFisher Scientific) following the manufacturer's protocol. All amplicons were then sequenced with the same primers used for the amplification by the external provider Macrogen Europe (Amsterdam, the Netherlands).

Data analysis
The mtDNA COI and NADH2 sequences electropherograms were imported in MEGA v.X (10.1; Kumar et al. 2018), and carefully assessed and edited. All the sequences for each marker were aligned with the CLUSTAL W algorithm (Thompson et al. 1994) as implemented in MEGA and the correct amino acidic translation was verified to exclude nuclear mitochondrial pseudogenes (Song et al. 2008).
The software DnaSP v.6 (Rozas et al. 2017) has been used to compute the number of haplotypes, the number of polymorphic and parsimony-informative sites, the haplotype (Hd) and nucleotide diversity (Pi) were estimated for each of the newly obtained COI and NADH2 datasets.
Available COI and NADH2 sequences were retrieved for S. acanthias, Squalus albifrons Last et al. 2007 Table S2). In both datasets, when possible, retrieved data had different geographic origin to properly assess intraspecific variability.
For each mtDNA marker, a Neighbour-joining (NJ) tree (Saitou & Nei 1987) was computed with MEGA using p-distance (Collins & Cruickshank 2013) and, although no nucleotide gap was detected across both datasets, the 'pairwise deletion' option for the treatment of gaps and missing data was selected. A bootstrap test (BS) with 10,000 replicates (Felsenstein 1985) was performed to evaluate the robustness of reconstructions. The average intra and interspecific genetic distances among the clades observed were calculated with MEGA.
In order to obtain more robust and comparable results in relation to previous studies (Marouani et al. 2012;Bonello et al. 2016;Vella et al. 2017;Veríssimo et al. 2017), a concatenated dataset was created merging the newly obtained COI and NADH2 individual sequences with those published by Veríssimo et al. (2017) and Vella et al. (2017). The two mtDNA genes were concatenated in a congruent dataset (see Table S2), using Phyutility v.2.2 (Smith & Dunn 2008). A second NJ tree was computed on the concatenated dataset as reported above.
The software IQ-TREE v.1.6.12 (Nguyen et al. 2015) was firstly used to identify the best substitution model to be applied to the concatenated dataset and then to compute the Maximum-Likelihood (ML) reconstruction using the TN+F+ G4 model and an ultrafast BS (Hoang et al. 2018) with 10,000 replicates.
The Bayesian Inference (BI) reconstruction was unravelled with MrBayes v.3.2.7 (Ronquist et al. 2012) using two independent runs of 1,000,000 generations and a 25% burn-in cut-off. Run convergence was assessed considering a mean standard deviation of split frequencies of <0.01 between runs.
For all the topology reconstruction methods, Cirrhigaleus asper (Merrett, 1973) (COI and NADH2 Accession numbers MN982926 and JQ518974, respectively) and Cirrhigaleus australis White, Last & Stevens, 2007 (NC_024059) were chosen as outgroups for the analyses. The acquired trees were summarised in one congruent topology and edited using FigTree v.1.4.2 (Rambaut & Drummond 2012). BS values for NJ and ML trees were reported near nodes, as well as the posterior probability (P) obtained for the BI reconstruction.
The phylogeographic relationship of species haplotypes was inferred with the Median Joining Tree clustering algorithm (Bandelt et al. 1999) implemented in the software PopART (Leigh & Bryant 2015). The network was obtained considering all the Squalus spp. nominal species and lineages detected across the con-

Results
DNA was successfully extracted from all tissue samples, although PCR amplification and sequencing were successful for 121 and 107 individuals for COI and NADH2, respectively. These newly obtained sequences of Squalus specimens collected in the Strait of Sicily and Tuscany coasts, have been deposited in GenBank under the Accession Numbers (COI Accession Numbers MW537886-MW537998; NADH2 Accession Numbers MW557187-MW557293).
After retrieving additional sequences from public data repositories, the final dataset for COI included 735 sequences (499 bp long) representing 18 nominal species, while the NADH2 dataset included 352 sequences (523 bp long) representing 20 nominal species. The concatenated dataset consisted of 222 sequences (1022 bp long).
Considering the congruence among the NJ, ML and BI tree topologies obtained using the concatenated dataset (Table S2), all reconstructions were summarised in one topology (Figure 1), which assigned most of the individuals collected off the southern coast of Sicily, Malta and Tuscany to a bigger cluster including S. blainville from the Mediterranean Sea and Squalus sp. Clade B (sensu Veríssimo et al. 2017) from South Africa, Angola, Namibia, Morocco, Portugal and the Mediterranean Sea (BS = 100% for both NJ and ML and posterior probability, P = 1 for BI). Conversely, two specimens of Squalus (S7 and S9) are included within the bigger group composed by individuals of Squalus sp. collected in Gabon and Guinea, representing the genetic Clade C (sensu Veríssimo et al. 2017; BS = 100% for both NJ and ML and P = 1 for BI). Furthermore, the three Squalus sp. collected around Malta by Vella et al. (2017) were grouped within Clade C (Figure 1).
The NJ tree topologies reconstructed considering a large number of nominal species characterised by a wide geographic distribution (e.g. Australia, China, Indonesia, New Zealand, Japan, Taiwan, UK, USA) have consistently highlighted the existence of three distinct and well-supported groups ( Figure S1, for COI; Figure S2, for NADH2). A first group (group I; Table S2) 184 A. Ferrari et al. included the two main Clades of S. acanthias and S. suckleyi with high BS values in both COI and NADH2 reconstructions (100% and 99%, respectively), a second one (group II; Table S2) included the nominal S. blainville, S. brevirostris, S. megalops and S. raoulensis (BS = 100% in COI tree topology and BS = 99% in NADH2), and a last one (group III in Table S2) included all the other species considered with BS = 99% in COI. Considering the COI mitochondrial gene, the genetic distance within groups was generally low, ranging from 0% (Squalus albifrons, S. grahami, S. griffini, Squalus cf. megalops, Squalus cf. mitsukurii from Hawaii) to 1.47% (S. nasutus) and 1.49% (S. megalops; Table S3). This was also confirmed by the NADH2 data, for which the genetic distance measured in S. mitsukurii complex sp. D and S. raoulensis was also 0%, while in S. brevirostris the genetic distance reached 3.02% (Table S4). COI genetic distances between groups ranged from 0.10% (S. clarkae vs S. cubensis) to 7.90% (S. acanthias vs S. megalops or S. acanthias vs S. raoulensis; Table S3, for COI), while they ranged from 0% (S. mitsukurii complex sp. D vs Squalus sp. Clade D) and 8.78% (S. edmundsi vs S. suckleyi) for NADH2 (Table S4). COI genetic distances increased when comparing Squalus individuals caught in the Strait of Sicily with S. megalops from Australia (1.40%), with Squalus cf. megalops from Mauritius (6.50%) and with Squalus sp. Clade C (sensu Veríssimo et al. 2017), which included individuals S7 and S9 collected in the Strait (6.60%; Table S3). NADH2 genetic distances showed similar results. In particular, distances increased when comparing Squalus sp. caught in the Strait of Sicily with S. blainville (0.52%), with Squalus sp. Clade B (sensu Veríssimo et al. 2017; 0.56%), with Squalus cf. megalops from Mauritius (4.95%) and with Squalus sp. Clade C (sensu Veríssimo et al. 2017), which included S7 and S9 individuals (5%; Table S4).
The haplotype network in Figure S3 showed the presence of two most frequent haplotypes shared by almost all the samples ascribable to S. blainville. Clusters of fewer individuals sharing the same haplotype, down to one individual displaying only a different variant, are separated from one to five mutational steps. The S. blainville haplogroups were separated from Squalus samples collected in Sicily and Malta and shared haplotypes with Squalus sp. Clade C from Gabon and Guinea, by at least 68 mutations. Furthermore, 59 mutational steps separated the S. blainville haplogroups from Squalus cf. megalops (Mauritius), while the latter was separated from Squalus sp. Clade C by 24 mutations.

Discussion
The use of molecular taxonomy is largely considered a powerful tool for the correct assessment of specimens identification, the discovery of new species and, in some cases, also the identification of cryptic Corrigan et al. 2008;Liu et al. 2013) and/or intricate complexes of species (Ward et al. 2008;Arlyza et al. 2013). Indeed, among sharks and batoids, new species have been described (De Astarloa et al. 2008;White & Iglesias 2011;Daly-Engel et al. 2018;Pfleger et al. 2018), old species resurrected (Ebert et al. 2010;Viana & de Carvalho 2018), and identification issues have been resolved with DNA barcoding when the morphological methods gave misleading results (Bonello et al. 2016;Cariani et al. 2017).
The lack of robust data from original descriptions compromises the correct identification of specimens with implications for species conservation and management purposes. This is particularly true for many genera, Squalus included (Veríssimo et al. 2017). The integration of morphological and molecular taxonomy techniques has been proved helpful (Henderson et al. 2016). In the last years, few studies in the context of integrative taxonomy were successfully performed on the genus Squalus aiming at the integration of new molecular taxonomy techniques to more classical morphological analyses with the aim to clarify taxonomic ambiguities between some of the species (Marouani et al. 2012;Bonello et al. 2016).
The evidence obtained in the present study, which followed in the footsteps of previous case studies in terms of genetic methodologies and analytical approach (Ward et al. 2007;Naylor et al. 2012;Vella et al. 2017;Veríssimo et al. 2017;Bellodi et al. 2018), confirmed the branching of the genus Squalus into three main lineages as referring to Squalus acanthias/S. suckleyi (group I), S. blainville/ S. megalops/S. raoulensis/S. brevirostris (group II) and a heterogeneous group of species mainly represented by the S. mitsukurii species complex (group III). At a finer resolution, the assignment of public data of Squalus individuals to Clades A-D (sensu Veríssimo et al. 2017) was confirmed by robust tree topologies. Almost all the individuals sampled in the Strait of Sicily and classified as Squalus sp. clustered in Clade B (sensu Veríssimo et al. 2017) and thus we suggest that they could be assigned to the nominal species S. blainville. On the other hand, two individuals sampled in the Strait of Sicily and classified as Squalus sp. fell in Clade C, together with four more specimens collected in the adjacent Tunisia (BOLD: FOAI329-09) and Malta waters (data from Vella et al. 2017) and showed a genetic distance from the former Clade B as high as the one existing between well-distinguished nominal species (e.g. S. mitsukurii and S. raoulensis from group III and II, respectively). Similar distances between genetic lineages were described in Veríssimo et al. (2017).
Therefore, here we confirm the hypothesis that specimens of Clade C are not related to the Australian S. megalops, since these lineages are separated by a distance of 6.20%. In fact, individuals of Squalus sp. from Clade C are closer to S. megalops from Mauritius (2.30% according to the COI gene). This was confirmed by the species-haplotype network, where at least 68 mutational steps separated Squalus sp. clustering in Clade B (ascribable to S. blainville) from Squalus sp. clustering in Clade C from tropical Western Africa. Similar haplotype differentiation between these same lineages was observed by Vella et al. (2017), who reported 83 mutations across a concatenated dataset 1600 bp long.
Hence, here we support previous studies that suggest that S. megalops does not occur in the eastern Atlantic and Mediterranean waters and that individuals composing Clade C should be considered a new species that needs formal description and proper taxonomic assessment (Last & Stevens 2009;Veríssimo et al. 2017). As a matter of fact, the occurrence of S. megalops appears more likely limited to the Australian waters, since the tree topologies obtained herein showed Squalus blainville from the Strait of Sicily (this study) clustering with longnose spurdog S. blainville from Sardinian coasts (Bellodi et al. 2018), Malta (Bonello et al. 2016;186 A. Ferrari et al. Vella et al. 2017), Tuscany (this study), Spain, Libyan and Ionian Seas, Greece (Kousteni et al. 2016) along with Squalus sp. Clade B (Veríssimo et al. 2017). Unfortunately, the sequence data from Marouani et al. (2012) are not publicly available and thus a direct comparison was not possible. The results discussed here, reinforced by both the genetic distances measured between groups and clades and the haplotype network, support the hypothesis that a common species, S. blainville, is currently inhabiting the Mediterranean Sea, while a second and extremely rare one, which has not yet been extensively described, is probably an occasional visitor of the Strait. This second species shows a strikingly high morphological similarity to S. blainville, but it is genetically distinct from both S. blainville and S. megalops. Within this perspective, the absence of records of S. megalops in the Atlantic Ocean and Mediterranean Sea (Straube et al. 2013;Veríssimo et al. 2017;Viana & de Carvalho 2018) can be confirmed. As previously mentioned, cryptic speciation among elasmobranchs is very common (Borsa et al. 2016(Borsa et al. , 2018 and the number of new descriptions, re-descriptions and resurrections of species is growing with the increasing application of molecular tools and integrated taxonomic methodologies. Starting from the accurate morphological data registered in Marouani et al. (2012) a dedicated effort is needed to identify and assess diagnostic features that characterize the individuals ascribed to Clade C exhibiting high morphological similarity with both species (S. megalops and S. blainville). What is more, a comparative approach involving both mtDNA markers and highly polymorphic nuclear DNA loci, besides morphology, is recommended to enhance the power of analyses aiming at unravelling the gene flow and migratory patterns of such a rare and geographically limited species.
To conclude, current inconsistencies in species identification within the genus Squalus need to be resolved and an accurate redescription of Squalus species is advised, especially in the Mediterranean where, along with commercial fishery, the bycatch of sharks may lead to the erosion of local biodiversity including genetic diversity. For all these reasons, the growing concern about the vulnerability of sharks to fishing pressure (Dulvy et al. 2014) and overexploitation (Simpfendorfer & Kyne 2009) now requires strong measures for species protection and management, especially in this exploited stretch of Sea.
To continue the acquisition of new information and the resolution of old problems, the establishment and strengthening of an international network of collaborations and the maximisation of sampling effort would go a long way towards filling the gaps in our knowledge of these shark species, which represent an irreplaceable component of biodiversity, in terms of both species and genetic richness, to be protected and conserved before they are "Lost Before Found" (White et al. 2019).

Funding
This research was funded by RFO and Canziani grants given to FT and AC, and by the FFO grant of the University of Bologna for the research fellowships of AF. JJB, LB and PJS were supported by funds from the University of Malta Research Committee. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript

Disclosure statement
No potential conflict of interest was reported by the authors.

Supplementary material
Supplemental data for this article can be accessed here.