Commercial sharks under scrutiny: Baseline genetic distinctiveness supports structured populations of small-spotted catsharks in the Mediterranean Sea

The present study, based on microsatellite markers, describes a population genetic analysis of the small-spotted catshark Scyliorhinus canicula (Linnaeus, 1758), representing one of the most abundant and commonly caught cartilaginous fishes in the Mediterranean Sea and adjacent areas. The analyses were performed to unravel the genetic features (variability, connectivity, sex-biased dispersal) of their relative geographic populations, both at the small (around the coast of Sardinia, Western Mediterranean Sea) and at a larger spatial scale (pan-Mediterranean level and between the Atlantic Ocean and the Mediterranean Sea). Individual clustering, multivariate and variance analyses rejected the hypothesis of genetic homogeneity, with significant genetic differences mainly within the Mediterranean between the Western and Eastern basins, as well as between the Mediterranean and the NE Atlantic Ocean. In detail, our results seem to confirm that the Strait of Gibraltar could not represent a complete barrier to the exchange of individuals of small-spotted catshark between the Atlantic Ocean and the Mediterranean Sea. In the latter area, a complex genetic structuring for S. canicula was found. Apart from differences among the Western, Eastern and Adriatic sites, within the Western basin the small-spotted catsharks around Sardinian waters are strongly differentiated from all others (both from the eastern Tyrrhenian Sea and southernmost part of the Algerian basin) and are demographically stable. Several possible mechanisms, both biological and abiotic (e.g., migratory behavior, waterfronts, and oceanographic discontinuities), are discussed here to explain their peculiar characteristics. Overall, the genetic data presented, both at the local and regional level, could represent a baseline information, useful for the temporal monitoring of populations, and to assess the effects of present or future fishing/management/conservation measures.


Introduction
Elasmobranchs are by far the most endangered group of fishes in the Mediterranean Sea, with over 40% of species classified by the IUCN as threatened: Vulnerable (8.1%), Endangered (12.8%) and Critically Endangered (19.8%) . The remarkable increase in fishing effort and the inherent low resilience of this taxon to harvesting have led to the decrease of Mediterranean shark and batoid populations during the last 50 years (Serena et al., 2020). Characteristic life-history traits such as slow growth, late maturity, low fecundity and productivity, long gestation periods, and a long-life span make elasmobranchs vulnerable and highly impacted by fisheries. Other factors, such as habitat loss and environmental degradation, have contributed to abundance and diversity decline in the Mediterranean area (Peristeraki et al., 2020) and highlight the urgent need for a conservation action plan for chondrichthyan populations of the Basin (Dulvy et al., 2014;Dulvy et al., 2017; SPA/RAC-UN Environment/MAP, Jorgensen et al., 2022). However, specific management measures to improve their status are still limited due to insufficient information on many species (Elliott et al., 2020;Dulvy et al., 2021). In the Mediterranean Sea, nearly 30% of species are still classified as Data Deficient (18.6%) or Not Evaluated (10.5%) .
Intense harvesting and related population declines may induce genetic effects such as the loss of genetic variation, inbreeding, genetic drift, and selective genetic changes, reducing the resilience of overfished species (Domingues et al., 2018;Pacoureau et al., 2021). Therefore, genetic studies, describing the intraspecies diversity of sharks and rays, are urgently needed.
It has been suggested that, apart from relying only on studies on the basic biology, life history and population ecology, shark and ray fisheries management and conservation policies should also consider genetic diversity metrics. Because the long-term survival of a species is strongly dependent on the levels of genetic diversity within and between populations, knowing the effective population size, observed heterozygosity, and allelic richness is of priority interest (Ovenden et al., 2015;Domingues et al., 2018). Thus, urgency should be placed on genetic studies of shark and ray species that are currently heavily exploited, conducting not only analyses essential to determining management priorities, but a genetic monitoring program producing temporal replicates to assess genetic variation over time. These data series would enable the assessment of temporal stability of the population structure, the loss of genetic diversity over time, and/or changes in population size (Domingues et al., 2018;Jorgensen et al., 2022).
In this context, the present paper deals with population genetic analysis of an exploited elasmobranch which is widespread and abundant in the Mediterranean and the adjacent eastern Atlantic Ocean: the small-spotted catshark Scyliorhinus canicula (Linnaeus, 1758) (Cashion et al., 2019;Follesa et al., 2019). S. canicula inhabits the continental shelf and uppermost slope waters throughout the Mediterranean Sea as well as the continental shelves of the Northeastern Atlantic Ocean, from Norway and the British Isles to Senegal and possibly Benin, from nearshore to a depth of 800 m (Finucci et al., 2021). Recent fishery data indicate it is a commonly caught species in the western Mediterranean and encountered less regularly in the eastern basin (Follesa et al., 2019). The small-spotted catshark is frequently captured accidentally but is often retained as a valuable by-catch of fisheries that focus on more productive teleost fish species (Follesa et al., 2019). However, there is no evidence of significant population decline in S. canicula; on the contrary, stable or even increasing population trends are reported in some parts of their range, justifying the inclusion of the species in the global IUCN Red List in the lower-risk categories as Least Concern (Finucci et al., 2021). The acquisition of genetic data, both at the local and regional level, could represent useful baseline information for the temporal monitoring of S. canicula populations, which constitute a valuable tool to assess the effects of present or future fishing/management/ conservation measures.
Previous microsatellite-based studies on S. canicula have consistently identified genetic structuring within the Mediterranean Sea but not among Atlantic sites (Gubili et al., 2014;Kousteni et al., 2015b;Ramıŕez-Amaro et al., 2018). In Kousteni et al. (2015b) data were exclusively limited to the Mediterranean Sea, with a denser sampling in the Eastern Mediterranean Sea. In Gubili et al. (2014) the samples were more abundant in the Atlantic Ocean, while Ramıŕez-Amaro et al. (2018) concentrated the sampling effort in the Alboran area, the Balearic Islands and northern Mediterranean coasts of Spain.
Recently, using genome wide nuclear SNPs, Manuzzi et al. (2019) revealed fine-scale genetic differentiation across the northeastern Atlantic, between the British Isles and southern Iberia.
The present study deals with the analysis of genetic diversity, population structure and connectivity of small-spotted catsharks within the Mediterranean Sea, and between the Mediterranean Sea and the Atlantic Ocean. The main aim of the study is to complement and compare the genetic data obtained in previous studies with a very similar set of microsatellite loci (Gubili et al., 2014;Kousteni et al., 2015b;Ramıŕez-Amaro et al., 2018), but implementing quite a different sampling scheme, with a denser sampling on the Central-Western Mediterranean Sea and a focus around the Sardinia Island. A tested panel of nuclear microsatellite markers has been used to evaluate the dispersal capabilities and the occurrence of potential physical barriers to gene flow and hence to improve the knowledge of S. canicula populations for appropriate management and protection of the species. The outcomes are discussed under an overview of the population structure of the species. The new genetic data represent the main strength of the study, as they guarantee a reference point for future monitoring.
2 Material and methods

Sampling
A total of 229 S. canicula specimens were collected between 2008 and 2012 during scientific bottom trawling campaigns (e.g., MEDITS; Spedicato et al., 2020) or by commercial fishermen, from 11 sites located in the northeastern Atlantic Ocean (Nor) and throughout the Mediterranean Sea ( Figure 1 and Tables 1, S1).
Laboratory procedures for tissue sampling and storage are described in the Supplementary Information (SI).
2.2 DNA extraction, loci amplification, and microsatellite data variability DNA extraction protocols, and PCR conditions with primer pairs are described in SI and Table S2, respectively. Individuals of S. canicula were genotyped at 12 microsatellite loci (Griffiths et al., 2011). Locus Scan7 failed to amplify in Adriatic sample (AdrN) and hence it was excluded from the final dataset. The microsatellites' allele size was scored using GENEMARKER v.1.8 (SoftGenetics Inc.) after setting the proper panels for binning analysis. Micro-Checker v.2.2.3 was used to check the data for genotyping errors and null alleles (Van Oosterhout et al., 2004) whose frequency for each locus and population was inspected with the program FreeNA (Chapuis and Estoup, 2007), using the algorithm by Dempster et al. (1977). Mean number of alleles (Na), observed (Ho) and unbiased expected heterozygosity (UHe) were calculated using GenAlEx v.6.5 (Peakall and Smouse, 2012). Total number of alleles per sample and Allelic richness (Ar) were estimated using the R package PopGenReport v.3.0.4 (Adamack and Gruber, 2014). All loci and samples were tested for linkage disequilibrium (LD) and deviations from Hardy-Weinberg expectations (henceforth HWE, null hypothesis H1 = heterozygote deficiency) using exact tests in Genepop v.4.2 (Rousset, 2008) with default settings. The inbreeding coefficient F IS was computed over all loci for each sample and for single locus over all samples with GenoDive v. 3.04 (Meirmans, 2020).

Genetic differentiation and population structure
The software POWSIM v.4.1 (Ryman and Palm, 2006) was used to assess the statistical power to detect genetic population differences with the sets of microsatellite loci (SI).

Marttinen
Genetic differentiation among population samples was investigated based on pairwise fixation indices, whose significance was calculated with 1,000 permutations using GenoDive. Pairwise F ST values were also calculated using the ENA (Excluding Null Alleles) method implemented in FreeNA to account for the possible effect of null alleles. Given that F ST tends to underestimate differentiation when heterozygosity is high (as in our case), we also calculated an allelic differentiation measure (Jost's D EST ;Jost, 2008) with GenoDive. The two indices, F ST and D EST , are reported to quantify different but complementary aspects of population structure (Jost et al., 2018). All statistical tests with multiple pairwise comparisons were corrected, calculating the false discovery rate (FDR) according to Benjamini and Hochberg (1995) as implemented in Myriads (Carvajal-Rodrıǵuez, 2017). A priori grouping of sampling sites were also tested with the Analysis of Molecular Variance (AMOVA) in Arlequin (Excoffier and Lischer, 2010); further details on the groups tested are provided in SI. To determine whether genetic differentiation was driven by geographical distance creating a pattern of isolation by distance (IBD), linearized pairwise F ST estimates (F ST /1-F ST ) were correlated against log-transformed geographical distances between samples (Rousset, 1997) in GENODIVE using Mantel tests (standard and stratified; Meirmans, 2020) with all the sites and the Mediterranean locations only. Given that null alleles can have important consequences on inferences, we used the pairwise F ST ENA values as recommended by Sere et al. (2017). Geographical distances were estimated as the minimum linear distance between pairs of locations by sea on Google Earth.

Demography, population size and connectivity
Microsatellite data were checked for the occurrence of recent demographic changes, namely genetic bottlenecks and population growth using the software BOTTLENECK v.1.2.02 (Piry et al., 1999).
NEESTIMATOR v.2.1 (Do et al., 2014) was used to calculate effective population size (Ne) and the number of effective breeders (Neb) for each sampling site. The LD method and the molecular coancestry method of Nomura (2008) were used to evaluate Ne and Neb, respectively.
To reconstruct source-sink population dynamics and the evolutionary processes leading to the present genetic diversity distribution, relative migration rates among all populations were estimated based on allele frequency data according to the method of Sundqvist et al. (2016) using the divMigrate function implemented in the R package diveRsity (Keenan et al., 2013), using G ST (Nei, 1973) as a differentiation measure based on allele fixation.
Detection of first-generation-migrants was performed using GENECLASS 2.0 (Piry et al., 2004). For each individual, the likelihood of being a resident (i.e., born in the sampling location) or a migrant from another reference population was estimated using the Bayesian method (Rannala and Mountain, 1997). Monte Carlo resampling with 10,000 simulated individuals following a simulation algorithm .

Sex-biased dispersal
Tests for sex-biased dispersal were conducted following Goudet et al. (2002). To test for possible effects of sex-biased dispersal on partitioning of genetic variation, a corrected assignment index (AIc) (Paetkau et al., 1995) was computed in GenAlEx. The difference in AIc values between males and females was tested using Wilcoxon's rank-sum test (Wilcoxon, 1945). Pairwise F ST and F IS values were calculated for both females and males to detect differences in their gene flow and inbreeding. Furthermore, relatedness (Re), the mean assignment index (mAIc) and variance of the assignment index (vAIc) were tested for significance for each sex using 10,000 permutations in FSTAT v2.9.4 (Goudet, 2003).
Details of each of the aforementioned analyses are provided in the S1.

Comparison of gene diversity among areas
Considering that in recent years S. canicula populations from different areas have been analyzed using 9 microsatellite loci used also in the present study, gene diversity estimates from previous papers (Gubili et al., 2014;Kousteni et al., 2015b;Ramıŕez-Amaro et al., 2018) were recalculated and compared with our data. After the elimination of areas with low sampling size or reduced representativity (Central Mediterranean and Adriatic), samplings were grouped in macro-areas and differences were tested with a nonparametric Kruskal-Wallis test, while pairwise comparisons between macro-areas were analyzed for differences between means in XLSTAT v2019.1.1 following the Dunn's procedure with a Bonferroni correction.

Genetic diversity
The final dataset was composed of 229 individuals genotyped at 11 loci. After false discovery rate correction, none of the loci was found to be in linkage disequilibrium. The statistics of the microsatellite data for both population samples and loci are described in Tables 1, S3. The microsatellite loci showed a moderate polymorphism varying from six alleles at locus Scan5 to 14 alleles at Scan17. The number of alleles per sample varied from 39 in the Ionian (Ion) to 74 in Algeria (Alg). Allelic richness ranged from 2.956 in the Ionian (Ion) to 3.807 in Morocco (Mor). The level of observed heterozygosity was slightly lower than the unbiased expected heterozygosity (uHe) in most populations. Significant deviations from probability for the multilocus Hardy-Weinberg test occurred at nine of the 11 loci and in all populations except for SicC and Ion sampling locations (Table 1). However, many of the population-by-locus combinations tests showed no deviation of Hardy-Weinberg equilibrium. The inbreeding coefficient F IS ranged from 0.012 to 0.180; it was significantly higher than zero only in four sites (SarSW, Mor, SicC, AdrN). The proportion of null alleles (Fnu) ranged from 0 to 17.50%, with an average of 3.77%, and the estimation of overall F ST gave comparable results (0.048 vs. 0.051 with vs. without the ENA correction, with overlapping 95% CI), consequently, as above, all the loci were considered as adequate for the following analyses.

Genetic differentiation and population structure
The simulation obtained with POWSIM highlighted the highresolution power of the dataset. Both Chi-square and Fisher's tests revealed genetic differentiation in most runs and were able to detect F ST values of as low as 0.0075 in > 99% of the runs (t 45 : c 2 = 0.99, F = 1.000, P value < 0.05).
Pairwise F ST and D EST values were mostly statistically significant, with F ST ranging from 0 to 0.133 and D EST from 0 to 0.297 (Table 2).
Pairwise F ST and D EST values for genetic differentiation were mostly statistically significant, with F ST ranging from 0 to 0.133 and D EST from 0 to 0.297 ( Table 2).
The Bayesian clustering analysis in STRUCTURE revealed the highest support for the presence of two or three genetic clusters (Figures 2, S1A, B). Delta K statistics, and Likelihood peaked at K = 2: one cluster coinciding with Sardinian samples, and the second one containing the remaining locations. However, K = 3 was also estimated as a possible number of genetic clusters, according to MaxMedK and MasMeanK Puechmaille's estimators (Figures S1C, D); along with the cluster including all the Sardinian samples, the other two clusters grouped the samples from the Atlantic Ocean and the North African coast (Nor+Mor+Alg), and the third one with all the remaining sites east of the Sicilian Strait (SicC, Ion, AdrN, Lev). No substantial differences were found when the combination of the alternative ancestry priors and an uncorrelated frequencies model was tested as suggested by Wang (2017) to account for the unbalanced sampling ( Figure S2). Similar to STRUCTURE, BAPS results indicated the best scenario at K = 2 (log(ml) = -7635.434), followed by K = 3 (log(ml) = -7638.007), while further substructures had lower support ( Figure 2). DAPC yielded overall congruent results with STRUCTURE K = 3, with Sardinian population differentiated from all other populations along the first axis, while Atlantic and North African samples (Nor+Mor+Alg) differed from all the remaining sites east of the Sicilian Strait (SicC, Ion, AdrN, Lev), separated along the second axis ( Figure 3). When the DAPC analysis was restricted to the Mediterranean locations, the eastern-central Mediterranean (Lev +AdrN+Ion+Sic) appeared distinct from the southern-western (Mor+Alg) and Sardinian samples ( Figure S3). The 'K-clustering means' results indicated two possible outcomes: the best clustering at K = 2 (SarN+SarSE+SarNW+SarSW/all other locations) according to pseudo-F, or, according to BIC, at K = 8 (Nor/SarN +SarSE+SarNW+SarSW/Mor/Alg/SicC/Ion/AdrN/Lev).
Hierarchical AMOVA showed a global F ST = 0.040 (p = 0.000), rejecting the hypothesis of genetic homogeneity (Table 3). However, less than 3.96% of the total variance was explained by differences among populations and 96.04% was within populations. AMOVA arranged by the two clusters identified by previous analyses (STRUCTURE, BAPS, K-clustering) confirmed the structuring (F CT = 0.031). However, when additional structures were tested, some of them indicated the occurrence of separation among clusters, especially the three groups structure (F CT = 0.035), and the eight groups, the last one characterized by a differentiation among groups (F CT = 0.0444) and the lowest differentiation within groups (F SC = 0.007) ( Table 3).
The Mantel test, performed to find out whether the population differentiation observed with the F ST pairwise analyses was due to an isolation by distance, showed that the variation observed could be attributable to isolation by distance both for the whole dataset (r 2 = 0.354, p = 0.012) and the Mediterranean sites only (r 2 = 0.218, p = 0.036). However, due to the limitations of the Mantel test (Meirmans, 2012;Meirmans, 2015), we also performed stratified Mantel tests within each of the identified clusters (see above), finding no significant isolation by distance in all cases.

Demography, population size and connectivity
Bottleneck tests failed to show statistical evidence that S. canicula populations from different origins had undergone a recent reduction in population size, except for Nor, even when locations were grouped in clusters as identified in previous analyses (Table S4).
High disparity of contemporaneous effective population size (Ne) in respect to location was recorded, with the lowest values in   Table 1 for geographical location codes. Abbreviations are in accordance with Figure 1.   (Table 4). The lowest estimate for the number of effective breeders within a cohort was recorded in AdrN and Alg (Table 4). Directional relative migration networks, based on the proportion of genetic diversity that resides among populations (G ST ) of S. canicula populations constructed with divMigrate, show that gene flow exists among the studied localities at different levels ( Figure 4A). Migration networks show that the predominance of gene flow was in general symmetrical, since directional migration was found statistically significant in a limited number of cases ( Figure 4B). When applying arbitrary filter thresholds to retain the most informative connectivity values, the highest gene flow was detected among the neighboring localities in the Western Mediterranean basin (i.e., among Sardinian sites, and Alg/Mor), whereas migration to the most distant localities seems to be limited or very scarce (Figures 4C, D).
Results from the first-generation-migrants detection analysis (Table 5) showed that about 80.70% (60.00-96.30%) of the individuals were residents (i.e., born in the sampling locations). Individuals detected as migrants were not necessarily from neighboring locations. The highest percentage of migrants (40%) were found in SicC and Ion, while the lowest values were recorded in SarNW and AdrN. However, the very small sample size of some location (i.e., SicC, Ion, Nor) call for a cautious interpretation of results.

Sex-biased dispersal
Analyses of sex-biased dispersal were conducted on a reduced dataset consisting only of individuals for which sex data and genotypes without missing data were available: S. canicula (86 females and 96 males, most of them from Sardinia, 61 females and 59 males). In the current study, patterns of molecular variation across sexes trended toward signals of female-biased dispersal in S. canicula. Indeed, the inbreeding coefficient, the fixation index, the relatedness, and the mean assignment indices indicate signals of female-biased dispersal, but in all cases the results are not statistically supported (Table 6).

Comparison of gene diversity among areas
The highest values of gene diversity were found in the Atlantic and Western Mediterranean (North African coasts) and the lowest in Alboran, Balearic basin, and Sardinia (Table S5). Differences in gene diversity between the macro-areas: (Atlantic = ATL; ANA = Alboran +NAfrica; BSA = Balearic + Sardinia; EMED = East Mediterranean) were statistically significant (Kruskal-Wallis P = 0.001; Table S6). In particular, the Atlantic and the most western Mediterranean samples (North African+Alboran) were found to be significantly different from all the other groups (Table S7).

Discussions
In this study, the population genetic structure of the smallspotted catshark was assessed at different spatial scales (within the Mediterranean Sea, and between the Mediterranean and the Atlantic Ocean, respectively) using polymorphic nuclear microsatellite loci. The new data acquired in this study complement the information obtained in previous studies (Gubili et al., 2014;Kousteni et al., 2015b;Ramıŕez-Amaro et al., 2018) with a very similar set of loci but a different sampling coverage, especially focused on the Western Mediterranean basin, around the Sardinia Island. Primarily, we used these data to test the hypothesis of genetic heterogeneity among populations. In fact, the life history traits (slow growth rate, late maturity, low fecundity, and no pelagic larval phase) suggest limited dispersal capabilities of the species, and the complex geomorphology of the Mediterranean Sea indicates the occurrence of potential physical barriers to gene flow. Results confirmed the occurrence of significant genetic differentiation in the investigated areas.
A weak to moderate but significant population differentiation was detected in S. canicula. Most of the pairwise F ST values were statistically significant except for few comparisons involving the Sardinian and Sicilian Channel samples (SicC-Alg and SicC-AdrN). However, the results for the Sicilian Channel, as well as Ionian Sea, and North Sea samples should be interpreted cautiously, due to their low sample size (<10 individuals).
The number of genetic clusters detected (from two to eight clusters) depended on analytical approaches. Individual-based analyses (STRUCTURE and BAPS) revealed the highest support for two clusters corresponding to the Sardinian samples in relation to all the remaining locations. Discriminant Analysis of Principal Components indicated that a subdivision in three clusters was also possible (Sardinian populations, as well as Atlantic and North African samples differentiated from other populations east of the Sicilian channel). In the same way, AMOVA rejected the hypothesis of genetic homogeneity and confirmed the occurrence of structuring at K = 2, K = 3, and K = 8 (each location different from others, apart from Sardinian populations grouped together).
Genetic differentiation of the Sardinian population has already been reported, unfortunately, no information on the precise origin of Sardinian small-spotted catshark was available in Gubili et al. (2014), so a direct comparison with our study is difficult. Genetic distinctiveness of Sardinian small-spotted catsharks has been reported in both close species (i.e., longnosed skate Dipturus oxyrinchus (Linnaeus, 1758) Melis et al., 2020), and in very distant taxa (i.e., mollusks Gopel et al., 2022;and corals Cannas et al., 2015;Cannas et al., 2016). Our data clearly indicate that S. canicula from all over Sardinia, both from the east and west side located in different basins, are plausibly connected by movements of adults in the coastal waters of the island, while the deepest waters, separating Sardinia from the Balearic archipelago and the North African coasts, cannot be easily crossed, potentially acting as a barrier to dispersal. Indeed, the genetic distinctiveness between Sardinian samples and Algerian ones could be ascribed to the role of the Channel of Sardinia-Tunisia, a stretch of deep (∼1900 m) and wide sea (∼70 km at 500 m depth) (Bouzinac et al., 1999;Testor and Gascard, 2005), that is the passage between the Algerian basin, on the west side, and the Tyrrhenian basin and the Sicily Channel on the east side (Bouzinac et al., 1999). Similarly, significant divergence between the sample collections from the Balearic Islands and the Algerian Basin was ascribed by Kousteni et al. (2015b) to the deep waters of the Balearic Abyssal Plain, a strong natural barrier to dispersal. Indeed, S. canicula is an obligate bottom-dweller, inhabiting the continental shelves; it is known for its non oceanodromous behaviour, avoiding open seas or abyssal plains (Kousteni et al., 2015a). Hirschfeld et al. (2021) performed a global overview of population genetic studies to investigate how marine barriers and dispersal ecology shape genetic connectivity in Elasmobranchs. They found that the three most common barriers are related to ocean bathymetry, with the maximum depth of occurrence, maximum body size and habitat of each species recognized as predictors of genetic connectivity, along with environmental tolerance and reproductive behavior (Hirschfeld et al., 2021). In particular, depth was described as the most common barrier generating genetic structure at intermediate down to surprisingly small geographical scales in the Mediterranean Sea were the depths of sub-basins (i.e., Algerian, Balearic, Sicilian Strait; Hirschfeld et al., 2021) form barriers at distances of less than 500 km in benthopelagic sharks and skates (maximum depth of occurrence of 800 and 630 m, respectively) that maintain connectivity along continuous shelf habitat (Gubili et al., 2014;Kousteni et al., 2015b;Frodella et al., 2016). Similarly, in Galeus melastomus Rafinesque, 1810, biotic factors such as bathymetry and water circulations might be contributing to the differentiation of samples (Di Crescenzo et al., 2022). Furthermore, Catalano et al. (2022) reported the genetic structuring of populations of the starry ray, Raja asterias Delaroche 1809, to be related to the extended deep marine waters in the western and eastern Mediterranean, acting as dispersal barriers for the species.
Apart from bathymetry, other abiotic drivers such as water circulations and temperature were indicated as affecting the Individuals on the diagonal are identified as residents. Sample location codes as in Table S1 and Figure 1. dispersion of benthic sharks (Di Crescenzo et al., 2022 and reference therein). In the present study, the known complex hydrography of the area, with the formation of gyres and fronts in the water southwest of Sardinia and between the Baleares and Sardinia (Bouzinac et al., 1999;Millot and Taupier-Letage, 2005;Testor and Gascard, 2005) could represent additional obstacles to dispersal, even if the role of fronts as factors affecting connectivity was found not strongly relevant for benthic-vagile species such as the species under study as it was for pelagic species (Pascual et al., 2017). The Sardinian area is peculiar also for its substantial demographic stability in the recent past and relatively high Ne values, where the recorded variance in Ne estimates across a range of P crit values suggests a history of gene flow and/or the presence of first-generation dispersers (Waples and England, 2011). Our bottleneck tests failed to show statistical evidence that S. canicula populations from different origins had undergone a recent reduction in population size. These results are in line with a previous study in which sample collections from Sardinia and Crete showed demographic stability (Gubili et al., 2014), but in contrast with the findings of other studies in different areas as highlighted by Ramıŕez-Amaro et al. (2018) who observed signals of past population bottleneck events in both the close Alboran and Balearic basins through analyses based on microsatellite data. This population event could have been caused by climate change during the Pleistocene glacial/interglacial periods, when much of the outer continental shelf in the area was exposed because of the sea-level drop, severely restricting the habitat available to S. canicula and causing a decrease in their populations (Ramıŕez-Amaro et al., 2018). Similarly, a weak decline has been described in the Aegean Sea related to the restricted habitat availability during the Pleistocene (Kousteni et al., 2015b). On the other hand, mitochondrial data indicated S. canicula experienced a slow and constant increase in population size over the last 350,000 years all over its distributional range (Ferrari et al., 2018), as well as a slight increase of the population size in the Ionian Sea (Kousteni et al., 2015b). However, the reliability of genetic bottleneck tests for detecting recent population declines has been questioned in long lived species such as S. canicula (Lippéet al., 2006;Peery et al., 2012;Bradke et al., 2021), and hence these results are to be interpreted cautiously.
The dispersal of species could be influenced by sea bottom topography and current patterns, morphological traits, and even by their reproductive behaviour. S. canicula is oviparous, producing egg cases anchored to macroalgae or various objects in the substrate, with behavioural strategies characterized by site fidelity and philopatry, as documented by biological, tagging, and several studies (Gubili et al., 2014;Kousteni et al., 2015b and references therein). According to Goudet et al. (2002) in general members of the most dispersing sex should display the following combination: a higher F IS and vAIC values, lower F ST , relatedness and mAIc values than the resident sex. Genetic data has provided evidence of malebiased dispersal and female philopatry in S. canicula populations (Gubili et al., 2014), but regrettably the sample from Sardinia was removed from their analyses as no information on sex was available. In the current study, patterns of molecular variation across sexes trended toward signals of female-biased dispersal in S. canicula, both in Sardinia and all over the Mediterranean, however, these were not statistically supported. Therefore, we should hypothesize that, among our samples, patterns of gene flow were symmetric across females and males, as estimated in about 48% of the elasmobranch species investigated so far (Phillips et al., 2021). To date, possible female biased dispersal (FBD) has been documented only in a single Carcharhiniformes species but only in a single population of its entire geographical range: the Galapagos Islands population of Carcharhinus galapagensis (Snodgrass and Heller, 1905) represents an uncommon case of documented female biased dispersal in shark that appears a rare condition in elasmobranchs (Phillips et al., 2021).
Alternatively, the remarkable genetic diversity could be explained by the fact that the Sardinian populations correspond to or are close to some putative glacial refugia for the species in the past (at the Last Glacial Maximum), as already proposed for red corals , with an opportunity for divergence due to local adaptation. Similarly, in the NE Atlantic Ocean, the three genetic clusters have been proposed to derive from three putative ancestral populations of S. canicula (Manuzzi et al., 2019). As a matter of fact, in species with limited dispersal capability, such as the small-spotted catshark, due to the limited ability of sustained swimming related to the tail morphology (Scacco et al., 2010), gene flow may be slow in breaking down the genetic differences accumulated during previous periods (Manuzzi et al., 2019).
It is largely demonstrated that geographic distance between sampling locations may be not sufficient to explain the presence or absence of connectivity in marine species. Nevertheless, in the present study, a significant pattern of isolation by distance (IBD) was identified, contrary to previous studies where no pattern of isolation by distance was detected (Barbieri et al., 2014) or it was found significant only when Atlantic and Mediterranean S. canicula were analysed together but not within the Atlantic or within the Mediterranean separately (Gubili et al., 2014). Although we detected a pattern of isolation by distance with standard Mantel tests only (when all locations were analyzed together), stratified Mantel tests (within each of the identified clusters) were not significant and suggested that other factors, such as oceanographic barriers impeding the gene flow among certain areas, could be playing a more prominent role than isolation by distance.
Overall, our results are in line with the previous genetic studies carried out on S. canicula, based both on mitochondrial and microsatellite markers, used separated or combined (Barbieri et al., 2014;Gubili et al., 2014;Kousteni et al., 2015b). In our study, despite the small sample size, the Atlantic sample was found statistically different from all the Mediterranean samples in all pairwise comparisons. Several authors concordantly indicated the distinction between the Atlantic and Mediterranean S. canicula (Barbieri et al., 2014;Gubili et al., 2014;Kousteni et al., 2015b). Intraspecific variation between the Atlantic and Mediterranean has also been described for biological features of populations such as adult size, mean length for the egg deposition, first maturity age, and sexual dimorphism (Gubili et al., 2014 and references therein;, suggesting they could represent different subspecies, especially for geographically distant subpopulation as West Africa (Gubili et al., 2014 and references therein). Additional differences such as teeth morphology, and color pattern, were listed to suggest the existence of a second species in the Mediterranean Sea distinct from S. canicula, named Scyliorhinus duhameli (Garman, 1913), distributed in the Adriatic and Mediterranean seas, along the continental shelves of Croatia, Greece, Tunisia, and Algeria (Soares and Carvalho, 2019).
However, DAPC and AMOVA analyses grouped the North Sea (Nor) sample with those from the North African coasts (Morocco and Algeria); similarly, the results from the first-generationmigrants analysis suggest the possible gene flow brought by individuals entering from the Atlantic Ocean and connecting locations along the African coasts, even if the very small sample sizes call for a cautious interpretation of these results. Considering the data from this and previous studies (Gubili et al., 2014;Kousteni et al., 2015b;Ramıŕez-Amaro et al., 2018), the Atlantic and North African samples are similar also in respect to the gene diversity, being characterized by higher values than other Mediterranean areas (Balearic, Sardinia, East Mediterranean).
Contrary to the findings of previous studies both on bony (e.g., Lo Brutto and Arculeo, 2004;Sěgvic-Bubićet al., 2016) and cartilaginous fishes (e.g., Catarino et al., 2015;Gubili et al., 2016) that report a clear division between Atlantic and Mediterranean populations, our results seem to confirm that the Strait of Gibraltar could not represent a complete barrier to the exchange of individuals of this species between the Atlantic Ocean and the Mediterranean Sea as suggested by Ramıŕez-Amaro et al. (2018). According to Ramıŕez-Amaro et al. (2018), the Alboran Sea, the westernmost portion of the Mediterranean located east of the Strait of Gibraltar, could represent a transitional zone between the Atlantic and the Mediterranean S. canicula populations, being genetically different from the adjacent Mediterranean (Spanish coasts and Balearic Islands), but similar to the southern Portuguese Atlantic area (Ramıŕez-Amaro et al., 2018). In the past, the occurrence of at least four main groups of genetically distinct small-spotted catshark have been proposed: North Atlantic, Atlantic Transition (Portugal and Alboran), Western and Eastern Mediterranean Sea (Gubili et al., 2014;Kousteni et al., 2015b;Ramıŕez-Amaro et al., 2018).
Recently, in contrast with previous studies that found no evidence of genetic structure across the northeastern Atlantic using traditional molecular markers (Gubili et al., 2014), additional fine-scale genetic differentiation was detected using SNPs in the mentioned area, with marked genetic divergence described for the Celtic Sea and South Portugal collections from their closest neighbors (Manuzzi et al., 2019). Similarly in the Mediterranean Sea, our results indicate that further structuring exists. For instance, in the Western basin, the Sardinian samples (four sites, two located in the Algero-Provencal basin and two in the western Tyrrhenian Sea) were strongly differentiated from all others (both from the eastern Tyrrhenian Sea and southernmost part of the Algerian basin). Additionally, the Adriatic samples were genetically distinct from the populations both from the Eastern and Western Mediterranean. The distinctiveness of the Adriatic area has been previously described (Barbieri et al., 2014;Gubili et al., 2014), as well as the differentiation between the Balearic and Algerian basins (Kousteni et al., 2015b;Ramıŕez-Amaro et al., 2018), and between the islands of Mallorca and Sardinia (Gubili et al., 2014).
The role of the oceanographic discontinuities can be invoked for geographically isolating and hence genetically differentiating the Adriatic samples and the easternmost from the Levantine Sea. Indeed, the Otranto Channel has been involved in limiting gene flow between the Adriatic and adjacent seas (Pascual et al., 2017), while the Sicily Channel could represent a barrier for the genetic exchanges between the eastern and the western Mediterranean, similarly to the deep Hellenic Trench (4,000 m deep), isolating the S. canicula populations both from the Aegean and Ionian Seas from the central-western Mediterranean (Kousteni et al., 2015b). Similar patterns of genetic differentiation across the Siculo-Tunisian area were also detected in other chondrichthyan species such as Raja miraletus Linnaeus, 1758 (Ferrari et al., 2018) and Raja asterias Delaroche, 1809 (Cariani et al., 2017;Catalano et al., 2022) as well as in bony fishes (Mejri et al., 2009;Mejri et al., 2011).
Yet, although the evolutionary history of the species has more to be revealed, here we support that the genetic differentiation described for S. canicula is probably linked to bathymetry, water currents, past events, coupled with life-history traits and ecology.

Final conclusions and management implications
In this research, the sampling was designed to obtain geographical coverage at both small (around the coast of Sardinia, Western Mediterranean Sea) and larger spatial scales (at the pan-Mediterranean level, and between the Atlantic and the Mediterranean Sea). Additionally, experimental analyses were designed to primarily acquire useful information to support potential management, describing the genetic features such as variability, connectivity, sex-biased dispersal of small-spotted catshark populations.
In agreement with the results of the different analytical approaches carried out (individual clustering, multivariate and variance analyses), the hypothesis of genetic homogeneity was rejected, finding significant genetic differences between the Mediterranean and Atlantic sharks, as well as within the Mediterranean between the Western and Eastern basins.
Our results confirm that structuring exists in the Mediterranean Sea, especially evident in the Western Mediterranean where smallspotted catsharks from the Sardinian waters were strongly differentiated from all others (both from the eastern Tyrrhenian Sea and southernmost part of the Algerian basin). Special attention and close monitoring of these populations are strongly necessary and warmly suggested, also due to the fishing pressure regularly exerted on them, to fill the gaps in knowledge and acquire data useful for the effective management and conservation plan for this commercial species.
Future research actions should be directed to gathering information to investigate the occurrence of discrete groups of populations in areas not covered or underrepresented in previous studies to fill the gaps in our knowledge and acquire data useful for the effective management and conservation plan for this commercial species.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Ethics statement
Ethical review and approval was not required for the animal study because the following information was supplied relating to ethical approvals (i.e., approving body and any reference numbers): samples of shark 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. 1967No. /2006. This Regulation is de facto the unique authorisation needed to conduct this type of activity.