A tale of two seas: contrasting patterns of population structure in the small-spotted catshark across Europe

Elasmobranchs represent important components of marine ecosystems, but they can be vulnerable to overexploitation. This has driven investigations into the population genetic structure of large-bodied pelagic sharks, but relatively little is known of population structure in smaller demersal taxa, which are perhaps more representative of the biodiversity of the group. This study explores spatial population genetic structure of the small-spotted catshark (Scyliorhinus canicula), across European seas. The results show significant genetic differences among most of the Mediterranean sample collections, but no significant structure among Atlantic shelf areas. The data suggest the Mediterranean populations are likely to have persisted in a stable and structured environment during Pleistocene sea-level changes. Conversely, the Northeast Atlantic populations would have experienced major changes in habitat availability during glacial cycles, driving patterns of population reduction and expansion. The data also provide evidence of male-biased dispersal and female philopatry over large spatial scales, implying complex sex-determined differences in the behaviour of elasmobranchs. On the basis of this evidence, we suggest that patterns of connectivity are determined by trends of past habitat stability that provides opportunity for local adaptation in species exhibiting philopatric behaviour, implying that resilience of populations to fisheries and other stressors may differ across the range of species.


Introduction
Molecular genetic markers have had a profound impact in conservation and management [1][2][3]. They allow inferences to be made about the scale over which genetic connectivity occurs [4], the behavioural mechanisms leading to gene flow [5], the estimation of effective population sizes [6] and how historical processes have impacted on populations [7]. This information has importance in terms of potential for local adaptation [8,9], but is also valuable for sustainable management and conservation of these natural resources [10][11][12][13]. In fisheries management, genetic tools can help to identify discrete populations which represent demographically independent stocks that require individual management to ensure fisheries sustainability [10,14].
The effects of marine overexploitation have been clearly shown by severe declines in many shark species [15]. Understandably, much of the work investigating genetic diversity and population structure of elasmobranchs has been focused on the large pelagic sharks that are considered particularly vulnerable to exploitation in high seas fisheries [16,17]. However, these species represent a small fraction of the biodiversity of the group. The small-spotted catshark (Scyliorhinus canicula L., 1758) is a relatively small demersal species belonging to one of the largest families of sharks, the Scyliorhinidae. It is generally considered to be the most abundant catshark in European shelf seas [18] and occurs from Norway and the British Isles, south to Senegal, including the Mediterranean Sea [19]. It is an oviparous species that breeds most of the year and is relatively fecund for an elasmobranch. It has remarkable variation in reproductive parameters, but in British waters it has been shown to lay between 29 and 62 eggs from November to July each year [18,20]. In the Atlantic, it is often caught as by-catch in demersal fisheries, but its commercial importance is growing, particularly through its use as whelk bait [21], and it is also significant for recreational fishing in some regions [22]. In the Mediterranean, catsharks have been fished since ancient times, as documented by mosaics from the Roman age [23], and S. canicula is still targeted for consumption today [24]. Recent studies have shown very dramatic localized reductions in abundance [25]. For example, in the Adriatic Sea it has been estimated that the species has declined in abundance by up to 90% since the 1940s [26].
Investigations of elasmobranch population structure focusing on wide-ranging pelagic sharks have often revealed genetic differentiation over broad inter-or intra-oceanic scale [27,28]. By contrast, work on coastal and demersal species suggests they can have more highly divided population structure, which has implications for management and conservation [29,30]. Scyliorhinus canicula has a range of traits associated with a low dispersal potential, including internal fertilization and deposition of demersal eggs. In addition, mark-recapture studies suggest adults do not generally make long migrations [31]. These factors could potentially lead to population genetic structure in this species, a concept that has some support from apparent differences among populations in growth rates, habitat/depth preference and reproductive biology that could have arisen from local adaptation [20,25]. Indeed, populations within the Mediterranean show such marked changes from those in the Atlantic that they have historically been suggested as a different subspecies [32,33]. Similarly, a more recent study of sexual dimorphism in S. canicula noted significant morphological differences in dentition between west African, Mediterranean and west European populations, indicating the west African group could represent a distinct taxon [34].
Documented sex-biased dispersal and philopatry have also been shown to have serious effects on elasmobranch population structure. Female sharks generally make far greater investment in reproduction than males, potentially leading to discrepancies between optimal male-and female-fitness strategies, and generating sex-specific differences in behaviour [27,35]. Evidence from a growing number of studies suggests female philopatry in sharks may be widespread [36]. This is also significant as sexbiased differences in movement behaviour could potentially lead to sexual segregation in sharks, and in turn affect the sustainability of marine harvesting [37].
Molecular data have the power to infer historic processes, such as population expansions or contractions, locations of refugia and patterns of recolonization [38]. Inferences of this type are particularly enlightening within Northeast Atlantic marine ecosystems, as organisms in this region have been significantly impacted by the Pleistocene glacial cycles [39]. During the last glacial maximum, about 20 000 years ago, ice sheets dominated the majority of the United Kingdom and Ireland, while permanent sea ice may have extended as far south as the Bay of Biscay [40,41]. Therefore, the distributions of marine organisms might have been forced southwards into refugia, including the Mediterranean, north African coast and the Iberian Peninsula [42][43][44]. Along the Atlantic coast, there is also evidence for refugia further north and much closer to the ice sheets [43,45,46]. Subsequently, as the ice retreated, organisms were able to recolonize the more northerly regions that were previously glaciated. Phylogeographic investigations of a variety of marine taxa have shown a division between the Mediterranean Sea and the Atlantic, although the degree and geographical scale of the biogeographic separation varies among even closely related species [47]. Therefore, it is plausible that these locations also acted as refugia for the small-spotted catshark.
Here, we test for population genetic structure among populations of S. canicula collected across European seas. A particular focus is made of the Atlantic-Mediterranean transition, as it is often considered to be an important phylogeographic break. We also use these data to test for sex-biased differences in dispersal and philopatry. We discuss the results in the light of published work on the behavioural ecology of the species, and highlight the conservation and management implications. 3. Material and methods 3. 1

. Sample collection and DNA extraction
Tissue samples were collected between 2007 and 2011 (with the exception of the Western English Channel collection site, see below), primarily from research cruises throughout Europe and the Mediterranean (figure 1; electronic supplementary material, S1). Temporally replicated samples were collected from the same approximate region of the Western Channel (2003, n = 45; 2008, n = 26; 2010, n = 39). A minority of samples were also collected at landing or at fish markets, most notably all those from Portugal, Sardinia and Crete. Genomic DNA was isolated from S. canicula using the Wizard technique (Promega Madison, WI, USA).

Mitochondrial DNA sequencing
An approximately 900 base pair (bp) section of the mitochondrial DNA (mtDNA) control region was amplified using the newly designed primers ScyD1pF (ATGACATGGCCCACATATCC) and Scan2R (TTCTCTTCTCAAGACCGGGTA), using PCR conditions described in Griffiths et al. [50]. PCR products were cleaned and sequenced by Macrogen (Korea), using the forward primer ScyD1pF, which yielded a rsos.royalsocietypublishing.org R. Soc. open  Sample collections with inappropriately small sample sizes for estimating allelic frequencies from microsatellite loci (i.e. from Tenerife n = 1 and Norway n = 4; figure 1) were removed. Estimation of pairwise F ST values was conducted in ARLEQUIN using 10 000 randomizations, and 95% CI were estimated across loci in GDA 1.0 [58] with 10 000 bootstraps. Sequential Bonferroni corrections were used to minimize type I errors [59]. Within GENALEX 6.5 [60], Nei's pairwise genetic distances were calculated using the default settings and visualized by principal coordinate analysis (PCoA). Mantel tests were used to test for significant correlations between genetic distances and geographical distances, also in GENALEX.
A hierarchical analysis of molecular variance (AMOVA) was performed in ARLEQUIN, to test for significance between groups of sample collections. Data were grouped according to geographical location, with the following hierarchy: Northeast Atlantic (Scotland, North Sea, Bristol Channel, Western Channel, Ireland and Portugal) and the Mediterranean (Mallorca, Sardinia, the Adriatic and Crete).

Microsatellite data: clustering analysis
The complete microsatellite dataset was analysed in the software STRUCTURE v. 2.3.4 [61]. A 'hierarchical' approach [62] with multiple rounds of analysis was employed in order to capture the major structure within the data. Run-lengths included a 100 000 burn-in and 1 000 000 total length, with five iterations. The model incorporating sampling locations as prior information was employed and the numbers of clusters, K, varied between 1 and 12 in each run. The K method of Evanno et al. [63] was applied to judge the most likely values of K. In the analysis of each cluster, plots of the absolute values of ln Pr(X|K) and K were generated by STRUCTURE HARVESTER [64]. 3. 6. Mitochondrial DNA: summary statistics and population differentiation Summary statistics were primarily calculated using ARLEQUIN. Rarefied haplotype richness, standardized to a sample size of 22, was calculated using CONTRIB 1.02 [2]. Genetic differences among localities were estimated in ARLEQUIN, using both the genetic distance-based Φ ST and the frequencybased mtF ST , with significance of differences estimated with 10 000 permutations. A hierarchical AMOVA was performed in ARLEQUIN, using the same scenarios as implemented with the nuclear markers. Tamura and Nei's [65] genetic distances were calculated between all sample collections (excepting the lone African/Canaries sample) using MEGA v. 5.1 [66] and visualized using PCoA in GENALEX. Mantel tests were used to look for significant correlations between genetic distances and geographical distances, also in GENALEX.

Mitochondrial DNA: haplotype network and demographic analysis
A median joining network was used to investigate genealogical relationships between mtDNA haplotypes with NETWORK 4. 6 [69] tests were calculated in ARLEQUIN [53], with significant negative values indicative of recent population expansion. The demographic history was also evaluated by mismatch distribution analysis. Typically, a population of constant size is characterized by a multimodal distribution; alternatively, one that has experienced expansion usually shows a unimodal distribution [70]. Bayesian skyline plots (BSPs) were generated in the software package BEAST v. 1.7 [71], and plotted using the upper 95% highest posterior density. Generally the programme defaults were used, except the HKY + Γ + I mutation model was selected and the Markov chain Monte Carlo (MCMC) was set between 50 and 200 million iterations, depending on length required for convergence. A fixed clock was set using a divergence rate of 0.361% estimated from homologous control sequences from S. canicula and Scyliorhinus stellaris (AM Griffiths 2010, unpublished data) and the divergence date (22 Ma) from Sorenson et al. [72]. This is the closest calibration point available, but it has resulted in a comparatively slow estimate rate of divergence, albeit one broadly similar to those calculated for mtDNA in other sharks (0.8% [73], 0.67% [74]). It is important to highlight the uncertainty in estimation of the clock rates, which could have very significant effects on the phylogeographic reconstructions in the BSP.

Sex-biased dispersal
Initially, to assess differences between male and female dispersal a test comparing pairwise genetic distances between populations for microsatellite and mitochondrial data was employed. Specifically, with male-biased dispersal, lower genetic distances were expected to be present in biparentally inherited (microsatellite) markers than in maternally inherited (mitochondrial) markers. Following Daly-Engel et al. [35], paired t-tests in R 2. 15.1 were used to compare mtF ST calculated using mtDNA haplotype frequencies and F ST calculated using microsatellite allele frequencies.
Tests for sex-biased dispersal based on the microsatellite data alone were conducted in FSTAT. Samples from Sardinia were removed from analyses as no information on sex was available. Five methods based on differences on the inbreeding coefficient (F IS ), fixation index (F ST ), degree of relatedness, mean assignment indices (mAIc) and variance of the assignment indices (vAIc) between the philopatric and dispersing groups were used [75]. In principle, unequal levels of gene flow between males and females would lead to a Wahlund effect, and a heterozygote deficit resulting in a higher F IS in the most dispersing sex, while also leading to correspondingly lower F ST and relatedness values. The assignment index statistic indicates the probability of a genotype occurring in a population, and unequal levels of gene flow between males and females could lead to negative values of the mAIc in the most dispersing sex (as the distribution is centred on zero), while also increasing the corresponding vAIc.

Microsatellite data: summary statistics and population differentiation
No significant differences were detected between temporal samples from the Western Channel (pairwise F ST ranged from less than 0.000 to 0.006 and none were significant at the 99% CI) and all individuals were grouped together into a single sample collection. The mean observed heterozygosity across all populations was 0.632, varying between 0.571 (Mallorca) and 0.659 (North Sea). The mean expected heterozygosity across populations was 0.615, varying between 0.562 (Mallorca) and 0.661 (Portugal;  table 1; electronic supplementary material, S2). There was no evidence of deviation from HWE after sequential Bonferroni correction (electronic supplementary material, S2). Furthermore, there was no consistent evidence of scoring issues or null alleles. There were 37 (of 858, 4.3%) significant tests of LD across sample sites at the 95% CI, none of which remained significant after sequential Bonferroni correction. Mean allelic richness across loci varied between 5.320 (Mallorca) and 6.200 (Portugal; table 1). Across all 12 loci, significant genetic differentiation was observed (global F ST = 0.039, p < 0.001). Pairwise F ST values ranged from −0.005 to 0.070 and were significant for all combinations that included the Mediterranean sampling sites, while genetic differentiation was absent among all Northeast Atlantic populations, after Bonferroni correction (table 2). The PCoA plot of pairwise genetic distance between sample collections clearly separated the Atlantic and the Mediterranean, but showed a greater level of division among the Mediterranean than the Atlantic sampling sites (figure 2a). There was a significant association between genetic and geographical distance across all samples (R 2 = 0.357, p = 0.002, figure 3a). This was also observed within the Atlantic (R 2 = 0.333, p = 0.047), but not within the Mediterranean (R 2 = 0.087, p = 0.383). Hierarchical AMOVA (table 3) showed significant variation between the Atlantic and the Mediterranean groups (2.09%), and significant variation among populations within these regions (1.78%), although within population genetic variation was greatest (96.13%).

Microsatellite data: clustering analysis
Analysis of the complete microsatellite dataset in STRUCTURE initially divided the individuals into two clusters (electronic supplementary material, S3). The first cluster included the samples from the Atlantic, with all individuals demonstrating admixture coefficients showing more than 80% membership to this cluster, reflecting clear common ancestry. Similarly, individuals from Sardinia, Crete and the Adriatic also demonstrated admixture coefficients showing more than 80% membership to a second cluster, suggesting clear division between these two clusters. The Mallorca samples demonstrated some degree

Mitochondrial DNA data analysis: haplotype network and demographic analysis
The haplotype network demonstrated that two common haplotypes predominate in the Atlantic, but are also present in the Mediterranean (figure 4). There was little evidence of population structure in the Northeast Atlantic, as haplotypes did not appear to assort by sample location. However, a number of closely related haplotypes were unique to the eastern Mediterranean, supporting the distinctiveness of the sample collections from this region, particularly those from the eastern basin (Crete and Cyprus

Sex-biased dispersal
The mtF ST was significantly greater than F ST across the sampling area (paired t-test; t = 6.487, p < 0.001; tables 2 and 5) consistent with male dispersal and female philopatry. In direct comparisons of male and female microsatellite data, males showed a higher F IS , lower F ST , lower relatedness values, lower mAIc and higher vAIC, all consistent with male-biased dispersal. However, only F ST and mAIc showed a significant difference between the sexes (table 6).

Discussion
This study represents the most detailed analysis of genetic variation from a species in the largest family of sharks, the Scyliorhinidae, incorporating samples from a wide trans-European and Mediterranean area. The results showed high gene flow in the Northeast Atlantic Ocean, a region of connectivity between the Atlantic and the Mediterranean populations in the western part of the Mediterranean Sea, and genetic  differences among populations across the Mediterranean basin. Indications of male-biased dispersal and female philopatry were also identified, supporting growing evidence of sex-based differences in dispersal and behaviour in elasmobranchs [76,77].
A lack of temporal variation in the Western English Channel over an 8-year sampling period supported the stability in population structure in this region over time, which remains a key factor in interpreting genetic data. By contrast, the null hypothesis of panmixia in S. canicula was rejected by analyses for both mtDNA control region sequence and nuclear microsatellite markers when considering all the sample collections across the Northeast Atlantic and the Mediterranean. This indicates that traits associated with limited dispersal potential may have played an important role in limiting gene flow, a finding that is becoming increasingly common in coastal and demersal sharks [13,30]. A significant pattern of isolation by distance (IBD) was also identified by both sets of markers, primarily related to significant genetic differences between the Atlantic and the Mediterranean, especially eastern Mediterranean, populations. Such regional variation could also be linked to ontogenetic differences of S. canicula, with length and age at sexual maturity attained earlier in the Mediterranean than in the Atlantic [20,25,33], so that genetic and morphological differences appear to coincide.
The mtDNA data identified two main groups of haplotypes; the first included the highest frequency haplotypes, central to the network that were predominantly associated with samples from the Northeast Atlantic and the Balearic Islands. The second includes haplotypes unique to specimens found in central and eastern parts of the Mediterranean ( figure 4). Furthermore, pairwise Φ ST revealed highly significant differences between biogeographically distinct regions of the Mediterranean: the western Mediterranean, eastern Mediterranean and Adriatic Sea. Such a trend runs counter to previous studies on elasmobranchs in these waters, where a lack of divergence has been found [44,78,79]. However, similar patterns of population subdivision have been described by the only other investigation of S. canicula [23], where variation within the cytochrome oxidase I gene was analysed from predominantly Mediterranean individuals. That study did not find evidence of population structure at such a fine geographical scale, perhaps due to the relatively conservative nature of the protein coding gene region used [80]. Additionally, Barbieri et al. [23] grouped distant areas of the Mediterranean owing to sample size restrictions, potentially reducing the resolution of their analysis. Application of microsatellite loci in the current study also supported evidence of genetic sub-division between S. canicula in the eastern and western Mediterranean. Indeed, the results suggest differentiation at an even smaller scale, between sample collections separated by less than 500 km in the western Mediterranean (Balearic Islands and Sardinia).
The pattern of highly divided population structure across the Mediterranean contrasts very sharply with results from the Northeast Atlantic. Regardless of the markers used, or the analytical approach, there was no evidence of significant genetic differences between sample collections originating from this region. Individuals from Norway and Africa demonstrated microsatellite alleles and haplotypes that were generally common across the Atlantic, suggesting little evidence of population structure, even at the most extreme latitudinal ranges considered. This lack of genetic differences across the Atlantic waters does not conform to expectations from typical life-history characteristics of small, demersal elasmobranchs that suggest low dispersal potential [81]. Nevertheless, similar patterns of little genetic evidence of population structure have also been reported for the starry ray (Amblyraja radiata) in the north Atlantic over comparable scales [81,82]. Perhaps both of these species correspond to a more typical pattern of marine species population structure with large effective population sizes and high gene flow within the Atlantic that minimize the effects of genetic drift and lead to low levels of population structure that are difficult to detect [83,84].
While nuclear and mtDNA data generally exhibited highly concordant results, differing signals of population structure were identified around Mallorca, with mtDNA suggesting a close similarity with the Atlantic group Balearic Islands representing an important region of secondary contact. There is strong evidence that many marine species have invaded the western Mediterranean through the Strait of Gibraltar since the last glacial maximum, bringing previously allopatric lineages into contact [47,[85][86][87]. The differences in the modes of inheritance between the two marker types may also explain the differing patterns of population structure identified, with the smaller effective population size of mtDNA potentially leading to more rapid genetic drift and the fixation of Atlantic haplotypes within the populations of the western Mediterranean.

Phylogeography
There have been long-standing hypotheses suggesting that the Strait of Gibraltar (the so-called 'Pillars of Hercules') or the Almeria-Oran Front represents phylogeographic barriers that shape the biogeographical patterns of marine Atlantic-Mediterranean organisms (reviewed in [47]). A number of molecular studies corroborate these scenarios with genetic discontinuity between the Atlantic and the Mediterranean populations occurring at the Strait of Gibraltar [88][89][90][91][92]. However, since the opening of the Strait of Gibraltar occurred at the end of the Messinian salinity crisis, its status as a current barrier to gene flow has been questioned. This study does not support a genetic discontinuity across the Strait, as an important zone of secondary contact between these populations is actually present in the Balearic Sea (followed by increasing genetic distinctiveness in the more isolated and semi-enclosed regions of the eastern Mediterranean, perhaps corresponding to relictual populations). Similar zones of secondary contact between the Mediterranean and the Northeast Atlantic stocks are also observed in other fish, such as the Atlantic bonito and the swordfish [93,94], suggesting this may be a common phylogeographic scenario.
The Mediterranean and the Northeast Atlantic share a history of inter-connectedness and have a large number of species in common. However, some evidence of contrasting demographic signatures was detected between the Atlantic and the Mediterranean populations of S. canicula. The Atlantic populations typically showed patterns of population expansion, while the Mediterranean appeared much more stable. This could be explained by the long-term stability and the consistent existence of suitable habitat in the Mediterranean during the Pleistocene glacial and interglacial cycles, as has been observed for other species [95,96]. Many of the Mediterranean basins are relatively deep, potentially providing suitable habitat for S. canicula during climatically driven sea-level fluctuations, allowing its persistence [97]. Recolonization of North Atlantic shelf habitats may have been rapid over the current interglacial cycle, meaning that it would not have promoted genetic discontinuity among regions. Interestingly, no evidence of northern refugia was supported by our data, in contrast to other marine taxa [43][44][45]98,99]. This lack of spatial association is also shown in the haplotype distribution where the two main haplotypes are shared by individuals caught along the Atlantic coast (figure 4). Together these data are suggestive of a relatively sudden population expansion. The apparent absence of a northward decrease in genetic diversity (table 1) is, however, not generally supportive of the 'leading edge hypothesis' [100], where latitudinal genetic variation is reduced in recently colonized populations owing to stochastic processes.

Sex-biased dispersal
Despite obvious congruence in patterns of population structure identified between the mtDNA and microsatellite data (tables 2 and 5), global testing across all populations demonstrated significant differences that were consistent with male-biased dispersal and female philopatry in S. canicula. This supports growing evidence of sex-based differences in dispersal and behaviour in elasmobranchs more widely [76,77]. However, there are limitations with the tests employed in this study. In the comparison of mtDNA and nuclear DNA markers, the reduced effective population size of mtDNA leads to expectations of greater F ST values, regardless of differences in sex-biased dispersal. Therefore, the recent approach of Daly-Engel et al. [35] was used; they suggested this comparison is robust with the use of numerous polymorphic nuclear markers that provide strong statistic power to detect population structure [101]. If population sizes are stable, it will also decrease the chance that differences in mtDNA and nuclear DNA are driven purely by variation in marker effective population size [35].
In order to overcome the limitations with comparisons of mtDNA and nuclear DNA, an additional suite of tests for sex-biased dispersal that focus on the microsatellite markers alone was conducted. These tests have a number of assumptions, including non-overlapping generations, equal sample size, juvenile dispersal and equal sex ratios [75] that are not necessarily satisfied in this case. However, the methods do provide a valuable way of interrogating the data and have been used in similar studies [102][103][104].
The apparent male-biased dispersal supports a long history of work demonstrating strong sex-based differences in the behaviour of the small-spotted catshark. It was the first elasmobranch species for which systematic analyses of unequal sex ratios in trawl catches provided clear evidence for unisexual aggregations [105,106], with recent work suggesting that this is the result of sexual harassment by males [107,108]. Finally, our results also support the commonly made association between the reproductive strategies of sharks and marine mammals that much greater investment in reproduction by females than males is at the evolutionary root of these differences in behaviour [27,35].

Conservation implications
The assessment and management of shark stocks is not well established, in part, because many characteristics such as dispersal or migratory behaviour are not fully understood [37]. Additional factors such as female philopatry, sexual segregation and sex-biased dispersal should be better considered in any management regimes, as spatially focused fisheries could result in the differential exploitation of sexes. The results of this study clearly show the potential for S. canicula to form multiple stocks within its distributional range. This has important implications for sustainable management, as effective conservation measures may need to be implemented at the level of the demographic unit to ensure longterm stock viability in the face of exploitation [109]. This is especially relevant in the case of the northern Adriatic Sea stock, which appears to have undergone dramatic declines in abundance [26], but for which recovery may not be as simple as immigration from proximate regions.