During scientific surveys on fishing vessels operating various fishing gears from 2016 to 2020, muscle tissue samples were opportunistically collected from 300 individuals (86 Mm and 214 Mp), representing populations from the Strait of Sicily (SIC; GSA 16) and the Adriatic Sea (ADRI; GSA 17 and 18). Within ADRI, samples were taken from the Italian Northern-Central coasts (N-ADRI; GSA 17) and the Montenegrin coasts (S-ADRI; GSA 18), with S-ADRI samples collected exclusively in 2020 ( Figure 1 ; Supplementary Figures S1 , S2 ). Sampled individuals were identified using the latest diagnostic morphological features from Marino et al. (2018) such as black spot presence, distance of the nostrils and mouth shape, and, when possible, dermal denticle observation under stereoscopes. Total length (TL; in centimeters), sex, and haul coordinates were recorded for individuals: Mm ranged from 40- cm to 150-cm TL (24 females and 14 males) and Mp from 30-cm to 120-cm TL (58 females and 68 males) in N-ADRI. In SIC, Mm ranged from 30- cm to 150-cm TL (21 females and 17 males) and Mp from 30- cm to 90-cm TL (30 females and 21 males). In S-ADRI, only female Mp were sampled, ranging from 92-cm to 138-cm TL ( Figure 1 ); only three Mm were sampled in S-ADRI and were not included in the analyses. From each individual, a muscle sample was collected and stored in pure grade ethanol at 4 degree Celsius (°C) for further analyses.

Genomic DNA was extracted from 30-mg to 40-mg tissue samples by salting-out protocol (Patwary et al., 1994), and the extract quality was checked by 1% agarose gel in TBE buffer (1×) electrophoresis, with GelRed stain (0.025 μL/mL; Biotium). Extracted DNA was conserved at −20°C.

A total number of 19 loci were used ( Supplementary Table S1.1 ), starting from previously tested microsatellites (Boomer and Stow, 2010; Chabot and Nigenda, 2011; Giresi et al., 2012; Marino et al., 2015a; Maduna et al., 2017). Among these, four loci were amplified in only one of the two species, according to previous successful amplification and testing (Marino et al., 2015a; Maduna et al., 2017), namely, Mmu9 and Gg22 in Mm and MaND5 and Mh29 in Mp. Thus, for each species, 17 microsatellites for each species were included throughout the analysis. The microsatellites were divided into three groups and amplified by a Multiplex PCR kit (QIAGEN) according to the published amplification profiles (Marino et al., 2015a; Hull et al., 2019). For each microsatellite, fluorophores were chosen according to length to avoid overlapping ( Supplementary Table S1.1 ). A separate single locus PCR was carried out for the Mmu11 locus, and PCR products were assembled before sending to the genotyping service (see Supplementary Materials for Mmu11 amplification profile). This was necessary only in Mmu11 because the allelic peaks did not match when comparing single and multiplex amplifications. After checking the successful amplification by electrophoresis in 1.8% agarose gel, PCR products were sent to the BMR Genomics (Padua, Italy) for genotyping service where an ABI Prism 3100 Genetic Analyzer and LIZ500 as size standard were used. Scoring to examine the allelic profiles was performed for each sample by the software PEAK SCANNER v1.0 (Applied Biosystems). The binning was done by FLEXIBIN (Amos et al., 2007).

For the CR amplification, the primer pair MaCYB/MaDLP was used (Hull et al., 2019), following the protocols reported therein. After confirming the amplification by electrophoresis, amplicons were purified by EuroSAP - PCR Enzymatic Clean-up (Euroclone, Italy) and then sent for Sanger sequencing at the Eurofins Genomics (Colone, Germany). All sequences were manually checked for quality by FinchTV (Geospiza Inc.) and then trimmed, edited, and forward primer–removed. All sequences were visualized in Mega 6 (Tamura et al., 2013) and aligned by Muscle algorithm (Edgar, 2004) with default parameters.

The sequenced fragment encompassed the last part of cytochrome b, two Transfer Ribonucleic Acid (tRNA), and the first part of the CR (orientation 5′-3′) when compared to the full mitochondrial genome of Mm, available on GenBank (accession number MH559351.1) by Blastn (Altschul et al., 1990) and Clustal Muscle alignment (Edgar, 2004). The initial sequenced fragment was composed of 713 bp and the first 230 bp encompassed the cytochrome b and tRNAs. The remaining 483 bp included the beginning of the CR region, which started at base position 15638 of the Mm complete genome ( Supplementary Figure S3 ). All the downstream analyses were performed considering only the CR fragment (483 bp).

Microsatellite genotypes were first checked for null alleles using MicroChecker v2.2.3 (Van Oosterhout et al., 2004), and their presence was evaluated on the estimation of genetic structure by FreeNA (Chapuis and Estoup, 2007). GENEPOP ON THE WEB v4.2 (Raymond and Rousset, 1995; Rousset, 2008) was used for testing linkage disequilibrium (LD) and deviations from Hardy–Weinberg equilibrium (HWE). For statistical significance of multiple tests, Benjamini–Hochberg correction (B-H; Benjamini and Hochberg, 1995) was applied with the function p.adjust on RStudio 1.3.1093-1 (RStudio Team, 2020). Bayesian structure analysis by STRUCTURE (Pritchard et al., 2000) was run to obtain the most reliable species identification for each specimen, assuming an admixture ancestry model with independent allelic frequencies and without a prior on sample origin (see Supplementary Materials for more details). Using the species identification based on STRUCTURE, for each species, the nDNA genetic diversity was calculated using the R package DiveRsity (Keenan et al., 2013), as the number of alleles observed per locus per population sample (A), the allelic richness (A_R), the observed and expected heterozygosity (H_O and H_E), and the inbreeding coefficient (F_IS). For the mtDNA, the CR diversity was calculated by DnaSP (Rozas et al., 2017) on the basis of the same species identification. Total number of haplotypes (H), haplotype diversity (h), and nucleotide diversity (π) were obtained. The haplotype network was created using TCS (Clement et al., 2000) and edited with PopART (Leigh and Bryant, 2015). To avoid that population structure was influenced by family structure, Colony v2.0.6.7 (Jones and Wang, 2010) was used for each species in separate runs selecting for female and male polygamy with inbreeding and cloning, by full-likelihood method at 95% and no prior. Only one sample for each full-sibling pair was kept for further analysis (see Supplementary Materials ).

To assess the level of genetic structure among the sampling sites in the two species, a Bayesian clustering analysis by STRUCTURE was performed with correlated allelic frequencies together with the abovementioned settings and software to account for uneven sample size.

For both nDNA and mtDNA, genetic differentiation, pairwise F_ST, and Φ_ST respectively, was determined by ARLEQUIN (Excoffier and Lischer, 2010) with 10⁴ permutations across sampling sites. Only for the Mp, the comparisons were performed between the ADRI and SIC (pooling together N-ADRI and S-ADRI), and among N-ADRI, S-ADRI, and SIC. When appropriate, B-H correction for multiple comparisons was performed with experiment-wide significance at p-values < 0.01. To assess whether unbalanced sample size in Mp may have affected the genetic differentiation, a random subsampling was performed, taking 50 random individuals from N-ADRI (N = 37) and S-ADRI (N = 13) while also balancing the sex ratio; the metrics of genetic differentiation were computed as described above.

The discriminant analysis of principal components (DAPC) (Jombart et al., 2010) was used to investigate population genetic structure. By summarizing between-group variation and ignoring within-group variation, DAPC effectively identified genetic structure within the samples. This multivariate approach assigned individuals to distinct groups and assessed inter-population differentiation without relying on specific population genetics models (e.g., HWE and LD).

A low genetic diversity emerged in this study, particularly evident in the nDNA of Mp with an expected heterozygosity of about 0.25 and only about two alleles expected per locus. For Mm, the genetic diversity was slightly higher, with about 35% expected heterozygotes and 3.5 alleles per locus, and, to a small extent, greater than what previously found in Mm with only nine microsatellites (Hull et al., 2019). The low nuclear diversity found in our study is in line with the one found in other benthic elasmobranchs in the Mediterranean populations such as the Galeus melastomus in the southern Adriatic Sea and in Sicily (A_R = 2.7–2.8 and H_E = 0.36–0.4; Di Crescenzo et al., 2022) and the Raja asterias in the northern Adriatic Sea (A_R = 3.35 and H_E = 0.56; Catalano et al., 2022). On the other hand, a higher nuclear diversity (A_R = 5.2–6.5, H_E = 0.56–0.76) was reported in the northern Adriatic Sea for other benthic elasmobranch species, such as the Schyliorhinus canicula (Gubili et al., 2014), the R. clavata (Melis et al., 2023), and the epibenthic Squalus acanthias (Gračan et al., 2020). The mtDNA diversity found in our study mirrors the pattern observed with nDNA, showing a very low variability in Mp (h = 0.16–0.31 and π = 0.0003–0.001) and higher values in Mm (h = 0.52–0.54 and π = 0.001); in Mm, despite the shorter sequenced fragment in our study, CR showed similar diversity values to those previously found (Hull et al., 2019). CR diversity, however, is known to vary between closely related species and within the same species (Subramanian and Lambert, 2011). Accordingly, several studies in other elasmobranch species reported very different patterns with a moderate to high CR genetic diversity in some species [reviewed by Hull et al. (2019); h = 0.78–0.99 and π = 0.004–0.35], including S. canicula from the Adriatic Sea (h = 0.71 and π = 0.003; Gubili et al., 2014).

For the aforementioned benthic elasmobranch species, whose genetic diversity was investigated in many Mediterranean subareas (Gubili et al., 2014; Catalano et al., 2022; Di Crescenzo et al., 2022; Melis et al., 2023), including the two smooth-hound species of the present study, their steep decline found support in literature (Colloca et al., 2017; Ferretti et al., 2013; Dulvy et al., 2021). This calls for further research to deepen the understanding of fishery-induced effect of the genetic diversity of over-exploited elasmobranch species (Domingues et al., 2018)

The evidence of genetic structure found in our study for Mm confirms and extends what was previously known at the global level for this species (Hull et al., 2019). In fact, although, as expected, our extended microsatellite panel confirms the existence of genetic differentiation between the Adriatic Sea and the Strait of Sicily, our study adds new support for the occurrence of a genetic structure between these two Mediterranean areas that, thanks to our higher sample size, were detected, for the first time, using also mtDNA. In Mm, the occurrence of one high-frequency shared haplotype in SIC and ADRI could be due to a recent common ancestor, whereas the other shared haplotypes could be the result of recent immigration. However, the existence of the unique haplotype (N10) points out to the occurrence of some level of isolation reflecting either localized adaptation or bottlenecks in N-ADRI. In Mp, although one prevailing ancestral haplotype was found, multiple unique haplotypes were detected in both N-ADRI, S-ADRI, and SIC, highlighting a higher degree of isolation and genetic drift. Through sequencing a longer CR fragment, focused research could discern whether mtDNA divergence could also result from a sex-biased dispersal (Hirschfeld et al., 2021; Phillips et al., 2021). However, the haplotype network and the value of Φ _ST and F_ST could potentially be explained by alternative hypothesis such as female-mediate gene flow, incomplete lineage sorting, or strong population genetic declines, as also discussed at global level in Hull et al. (2019). This result of genetic differentiation is in line with what reported for this species at the global scale, Mm displayed a significant level of isolation shaped by the effect of oceanic currents or other biogeographical barriers, despite the capability to undertake rare and long migrations (Mann and Bullen, 2009; da Silva et al., 2013; Maduna et al., 2016). In South African water, the combination of tagging and genetic methods in Mm contributed to find a contrasted dynamic between a strong site fidelity in a relatively short period and a significant gene flow at evolutionary scale (Klein et al., 2022).

With regard to Mp, our study provided the first evidence for a significant structure between SIC and ADRI, with a pattern that mirrors the Mm structure. It is worth noting that seasonal and permanent currents occur in this area; in particular, the ADRI circulation consisted in wind-driven currents, one from south to north along the eastern coast and another one from north to the south in the western coast, and three seasonal gyres (Russo and Artegiani, 1996). The ADRI hydrodynamics did not appear to influence the genetic differentiation between N-ADRI and S-ADRI in Mp while being responsible for the genetic structure of many marine species with larval dispersal (Papetti et al., 2013; Matić-Skoko et al., 2018). Similarly, other populations of Mustelus species appeared to be genetically structured, including the gummy shark (M. antarcticus) in Australian waters despite a wide movement range (238–900 km) (Braccini et al., 2017) and the brown smooth-hound (M. henlei) in the Gulf of California (Sandoval-Castillo and Beheregaray, 2015). Our analyses of re-assignment using DAPC further supported the concept of two distinct and yet slightly connected populations, highlighting the presence of a semi-permeable barrier between the Adriatic Sea and the Mediterranean Sea.

Despite being opportunistic, our sampling represents the first coordinated subregional effort, yielding the largest genetic dataset for these exploited species. Combining microsatellite genotyping with CR sequencing, we detected genetic diversity and a statistically significant population structure between the Adriatic Sea and the Strait of Sicily, confirming two distinct genetic stocks for Mm (SIC) and Mp (ADRI) for the first time. Our results were in line with previous analyses on Mediterranean and benthic elasmobranch species that showed significant genetic structuring between Mediterranean areas (Gubili et al., 2014; Catalano et al., 2022; Di Crescenzo et al., 2022; Melis et al., 2023). The diversity of habitats present between ADRI and SIC, encompassed by deep water and the presence of the Strait of Otranto, seemed to act as an only semi-permeable dispersal barrier for many populations of elasmobranch species, highlighting the importance of identifying genetic stock. Behaviors like philopatry and sex-biased dispersal may play a role in shaping the genetic structure and contribute to the connectivity in elasmobranch (Chapman et al., 2015; Phillips et al., 2021), but, due to the opportunistic sampling, this study was not able to address this goal and needs to be further investigated. Considering the stock definition (Ihssen et al., 1981), fishery management of these two stocks could be improved tailoring stock specific strategies that consider different growth rate, reproduction and nursery areas, and harvesting dynamic (Cadrin et al., 2014). The commercial relevance and the conservation issues of Mm and Mp in the two Mediterranean populations require the delineation of management units. Albeit a holistic approach is needed, our study provides solid evidence in relation to the genetic stocks of these exploited and threatened shark species.