Genetic differentiation and phylogeography of Mediterranean-North Eastern Atlantic blue shark (Prionace glauca, L. 1758) using mitochondrial DNA: panmixia or complex stock structure?

Background The blue shark (Prionace glauca, Linnaeus 1758) is one of the most abundant epipelagic shark inhabiting all the oceans except the poles, including the Mediterranean Sea, but its genetic structure has not been confirmed at basin and interoceanic distances. Past tagging programs in the Atlantic Ocean failed to find evidence of migration of blue sharks between the Mediterranean and the adjacent Atlantic, despite the extreme vagility of the species. Although the high rate of by-catch in the Mediterranean basin, to date no genetic study on Mediterranean blue shark was carried out, which constitutes a significant knowledge gap, considering that this population is classified as “Critically Endangered”, unlike its open-ocean counterpart. Methods Blue shark phylogeography and demography in the Mediterranean Sea and North-Eastern Atlantic Ocean were inferred using two mitochondrial genes (Cytb and control region) amplified from 207 and 170 individuals respectively, collected from six localities across the Mediterranean and two from the North-Eastern Atlantic. Results Although no obvious pattern of geographical differentiation was apparent from the haplotype network, Φst analyses indicated significant genetic structure among four geographical groups. Demographic analyses suggest that these populations have experienced a constant population expansion in the last 0.4–0.1 million of years. Discussion The weak, but significant, differences in Mediterranean and adjacent North-eastern Atlantic blue sharks revealed a complex phylogeographic structure, which appears to reject the assumption of panmixia across the study area, but also supports a certain degree of population connectivity across the Strait of Gibraltar, despite the lack of evidence of migratory movements observed by tagging data. Analyses of spatial genetic structure in relation to sex-ratio and size could indicate some level of sex/stage biased migratory behaviour.

control region) amplified from 207 and 170 individuals respectively, collected from six localities across the Mediterranean and two from the North-Eastern Atlantic. Results. Although no obvious pattern of geographical differentiation was apparent from the haplotype network, st analyses indicated significant genetic structure among four geographical groups. Demographic analyses suggest that these populations have experienced a constant population expansion in the last 0.4-0.1 million of years. Discussion. The weak, but significant, differences in Mediterranean and adjacent North-eastern Atlantic blue sharks revealed a complex phylogeographic structure, which appears to reject the assumption of panmixia across the study area, but also supports a certain degree of population connectivity across the Strait of Gibraltar, despite the lack of evidence of migratory movements observed by tagging data. Analyses of spatial genetic structure in relation to sex-ratio and size could indicate some level of sex/stage biased migratory behaviour.

INTRODUCTION
The blue shark (Prionace glauca, Linnaeus 1758; BS henceforth) is one of the most abundant epipelagic sharks that is found in all oceans from 60 • N to 50 • S (Compagno, 1984). Blue sharks are rarely targeted by commercial fishing, but feature prominently as by-catch of fisheries targeting large pelagic fish, especially swordfish and tuna longlines (Fowler et al., 2005). BS populations trend data are available only for a part of the geographic range and stock assessments are highly uncertain (Dulvy et al., 2014;Coelho et al., 2017); due to the huge amount of by-caught BS (approx. 20 million per annum, Stevens, 2009), the species has being categorized worldwide as ''Near Threatened'' in the IUCN Red List (Stevens, 2009). Based on recent assessment (ICCAT, 2015), the North Atlantic stock is unlikely to be currently overfished. The Mediterranean BS, on the other hand, is estimated to have undergone a 90% decline over three generations, primarily due to overfishing (Ferretti et al., 2008), and is now categorized as ''Critically Endangered'' (Sims et al., 2016). Given the vast amount of poorly reported by-catch, the increasing commercial value of the species (Megalofonou, Damalas & Yannopolous, 2005) and the persistent issue of the global trade in shark fin products, of which BS is the main component (Clarke et al., 2006), a more explicit management is needed for this species, which should be underpinned by robust knowledge of its population structure.
In the Atlantic, BS is distributed from Canada to Argentina, on the western side, and from Norway to South Africa on the eastern side, including the Mediterranean Sea (Compagno, 1984). The population structure and dynamics of Atlantic BS is still poorly known, despite several long-term tagging studies, which revealed extensive movements of BS tagged in the western side of the North Atlantic (henceforth NA), with well documented eastward trans-Atlantic migrations (Kohler, Casey & Turner, 1998;Kohler et al., 2002;Kohler & Turner, 2008;Vandeperre et al., 2014). Sexual segregation was also evident, with a concentration of mature females in more temperate waters of the northernmost NA, and immature males predominant in the southernmost NA (Sampaio da Costa, 2013). Mature BS of both sexes seemed to be distributed in the southern part of NA, while immature individuals of both sexes and sub-adult females are usually distributed in the northern areas (Kohler et al., 2002). Conversely, a prevalent occurrence of immature juveniles is reported in the Mediterranean Sea (Megalofonou, Damalas & De Metrio, 2009;Kohler et al., 2002). A significant genetic heterogeneity among potential BS nurseries from the Atlantic Ocean (Portugal and Azores) and those from South Africa was detected by Sampaio da Costa (2013) from mitochondrial and nuclear marker variation. Their finding indicated a deeper separation between the northern and the southern NA nurseries and supported a male philopatry behaviour to mating areas exclusively contributing to a single nursery ground. Contradictorily, a recent genetic survey (Veríssimo et al., 2017) carried out on the same dataset (i.e., young-of-year and <2 years juveniles) collected from the same nurseries, enriched with more samples from different areas (i.e., coasts of Brazil), and using the same type of markers, showed a lack of spatio-temporal genetic differentiation, suggesting the presence of a panmictic population in the whole Atlantic.
To date, no genetic data are available for the Mediterranean BS population and population structure and dynamics of BS in the Mediterranean are presently inferred only by Atlantic-Mediterranean integrated tagging studies and fishing data assessments (Kohler, Casey & Turner, 1998;Kohler et al., 2002;Ferretti et al., 2008;Kohler & Turner, 2008;Megalofonou, Damalas & De Metrio, 2009).
Irrespective of the small recapture rate (out of the 91,450 BS specimens tagged in the north western Atlantic, only 5.9% were recaptured), extensive tag-recapture surveys carried out from 1962 to 2000, indicated that North Atlantic BS form a single stock and that trans-Atlantic migratory movements were quite frequent, likely favoured by the oceanic current system (Kohler et al., 2002). Focusing on the Atlantic-Mediterranean connectivity, the reproductive migratory movements of Atlantic BS towards Mediterranean and the degree of population connectivity between the two areas are still unknown, because only one adult BS male tagged in the north-western Atlantic and one sub adult female tagged in the North-Eastern Atlantic were recaptured in the Mediterranean (Kohler et al., 2002). The large majority of BS tagged in the Mediterranean Sea were immature and remained in the tagging area, with the only exception of a subadult female that moved a short distance to the adjacent north-eastern Atlantic area. Most of the BS caught in the Mediterranean (99% and 98% for males and females, respectively) are immature, indicating that the Mediterranean BS stock consists primarily of small immature BS of both sexes, with a sex-ratio skewed toward females or males, depending on different geographical areas (Kohler et al., 2002;Megalofonou, Damalas & De Metrio, 2009). A high number of pregnant females was observed in the Adriatic, North Ionian Sea and Ligurian Sea, suggesting potential nursery grounds for BS (Megalofonou, Damalas & De Metrio, 2009;F Garibaldi, pers. comm., 2017). On the other hand, the adjacent South-Eastern North Atlantic BS was prevalently composed by primarily mature individuals of both sexes with male-based sex ratio.
The primary aims of this study is to test the null hypothesis of panmixia between North Atlantic and Mediterranean BS, by comparing the mtDNA genetic variation of two gene regions, the control region (CR) and the Cytochrome b (Cytb) among four population samples collected in the North-Eastern and South-Eastern North Atlantic and in the Western and Eastern Mediterranean. Given the female philopatry observed in other carcharhiniformes (Mourier & Planes, 2013;Tillett et al., 2012), mtDNA markers are likely to be useful to spot localised groups due to site-fidelity. Accordingly, this work aims to provide further and needed data on matrilineal genetic structure, female philopatry and demography of Mediterranean BS. These, previously lacking, data will contribute to a better understanding and inclusion of the Mediterranean BS dynamics in the wider North Atlantic population model, to improve assessment and management of BS stocks in the area.  (Table S1). The BS individuals were grouped according the Total Length (TL) in three size categories (Pratt, 1979;Vandeperre et al., 2014): juveniles (J, TL ≤ 120 cm), young (Y, TL = 120-180 cm) and large (L, TL ≥ 180 cm).

Blue shark sampling
A unique and transparent sampling documentation tool was developed within the project, in order to render data public. This tool can be used by everyone as an interactive map visiting the website: https://fishreg.jrc.ec.europa.eu/web/medbluesgen/sampling-data.

Molecular methods
Individual fin clips or skeletal muscle tissue samples were collected and preserved in 96% ethanol and kept at −20 • C until laboratory analyses. DNA extraction was carried out using the Invisorb R Spin Tissue Kit, Invitek (STRATEC Molecular, Birkenfeld, Germany) and the Wizard R Genomic DNA Purification Kit (Promega, Madison, WI, USA) following the manufacturers' protocols.
Species-specific primer pairs for the amplification of the mitochondrial control region (CR) and cytochrome b (Cytb) genes were designed. Homologous complete CR and Cytb sequences of Prionace glauca available in GenBank were retrieved and aligned using ClusterW algorithm implemented in MEGA ver.7.0 (Tamura et al., 2013). Primer pairs were designed using the online software PRIMER3 (ver.0.4.0) (Untergasser et al., 2012), minimizing the propensity of oligos to form hairpins or dimers or to hybridize or prime from unintended The designed primer pairs (control region: CR-Blues-F 5 AAACACATCAGGGGAAGGA G3 , CR-Blues-R 5 CATCTTAGCATCTTCAGTGCC3 ; Cytochrome-b: Cytb-Blues-F 5 TCCTCACAGGACTCTTCCTAGC3 , Cytb-Blues-R 5 GTCGAAAGATGGTGCTTCGT3 ) were tested using a temperature gradient to identify the most suitable melting temperatures (T m = from 50 • C to 60 • C) according to PCR cycling conditions described by Ovenden et al. (2009).
Once the optimal melting temperature was identified, the PCR thermal profile was adjusted and the PCR reactions were performed for both markers in a final volume of 50 µL containing 31.75 µL of distilled sterile H 2 O, 8 µL of Buffer 10× (Tris-HCl; final 1×), 3 µL of MgCl 2 (25 mM; final 1.5 mM), 2 µL of dNTPs (10 mM; final 0.37 mM), 2.5 µL (10 µM; final 0.46 µM) of each primer, 0.25 µL (5U/µL; final 1.5U) of Taq polymerase and 2 µL of template DNA(10-20 ng). The temperature profile included an initial denaturation at 94 • C for 2 min, followed by 35 cycles of denaturation at 94 • C for 30 s, annealing at 60 • C for 30 s, elongation at 72 • C for 30 s and a final elongation step at 72 • C for 5 min. PCR amplicons were sequenced using the external service provider MACROGEN R Europe.

Data analysis
The CR and Cytb nucleotide sequences obtained were validated with the homologous gene sequences deposited in the GenBank with the BLASTn search implemented in the NCBI website (Altschul et al., 1990), and aligned using the ClusterW algorithm implemented in MEGA ver.7.0 (Tamura et al., 2013). When aligned to the complete BS mitochondrial genome, Cytb sequences mapped from nucleotide position 14,530 to 15,291 and CR from 15,651 to 16,397.
Given the high potential of geographical dispersal of the species, sequence data were grouped according to the four geographical areas: EMED, WMED, SNEATL and NNEATL (Fig. 1). The software DNAsp v.5.10.01 (Librado & Rozas, 2009) was used to assess the genetic diversity parameters at both markers: the number of haplotypes (Nh), the number of polymorphic sites (S), the haplotype (h) and nucleotide diversity (π ) with associated standard deviation (stdev).
Demographic history was investigated using the mismatch distribution as implemented in the DNAsp software (Librado & Rozas, 2009).
Furthermore, historical demographic trend of the four groups was investigated using Bayesian Skyline Plot (BSP) analysis implemented in the software BEAST v.1.8.2 (Drummond et al., 2005;Drummond et al., 2012), using the best evolutionary models for both Cytb and CR markers inferred using JModelTest 2.1.1 (Darriba et al., 2012), and the average mutation rate for sharks, 0.62% and 0.31% for CR and Cytb respectively (Martin & Palumbi, 1993;Galván-Tirado et al., 2013). The same software and parameters, with associate software TreeAnnotator and FigTree, were used to define the phylogeny of the Mediterranean and Eastern Atlantic BS populations.
The Cytb sequence dataset exhibited 16 polymorphic segregating sites while CR dataset showed 27 polymorphic segregating sites. The Cytb haplotype diversity ranged from 0.784 to 0.835, and that of the CR from 0.894 to 1.000. The Cytb nucleotide diversity ranged from 0.001 to 0.002, and that of the CR from 0.004 to 0.008. Detailed genetic diversity of BS samples collected from the four macro areas and all sampling locations is presented in Table 1 and Table S2, respectively.
The Cytb and CR haplotype networks highlighted the distribution of haplotypes irrespective of the geographical origin of BS samples, indicating the lack of phylogeographical structure in the Mediterranean and adjacent North Atlantic BS (see Fig. 2, Fig. S1). In the Cytb network, the four main frequent haplotypes were shared by BS from all the four geographical areas, except for the most frequent haplotype which was shared by BS from the three geographical areas, SNEATL, WMED and EMED. In the CR network, six most frequent haplotypes (No. individuals ≥ 10) were observed. Although these six haplotypes were shared by all geographical areas, three of them were shared by Mediterranean and SNEATL, one by Mediterranean and NNEATL, and two within the Mediterranean. In both networks, most of the NNEATL haplotypes were singletons (Fig. 2). The AMOVA (Table 2) revealed a significant overall ST among population samples for both markers. Significant partition of molecular variance among areas was observed when BS sampling locations were grouped according to the four geographical areas in both markers (AMOVA4), according to three areas (NEATL (NNEATL+SNEATL) vs WMED vs EMED; AMOVA3), and according to two areas (NEATL (NNEATL+SNEATL) vs MED (WMED+EMED), for both dataset. However, the grouping that best described the partitioning of genetic variance is when the different sampling locations are subdivided into four areas showing the lowest partition of molecular variance among populations within group.

Figure 2 Cytochrome-b (A) and Control Region (B) Maximum Likelihood Haplotype Network of Mediterranean/North East Atlantic Blue
With the Cytb sequence data, all pairwise ST values among the four geographical areas were significant except that between the two Atlantic groups ( ST = 0.1152; p = 0.019) that became non-significant after the Bonferroni correction for multiple tests (Martin & Douglas, 1995) (Table 3). Unlike the CR dataset, only the pairwise ST values between SNEATL and the two Mediterranean areas and between WMED and EMED remained significant after the Bonferroni correction for multiple tests (Table 3).
AMOVA and pairwise ST analyses were performed on a reduced dataset, selecting only juvenile and immature specimens from each sampling site. Despite the reduced sample sizes and the complete absence of data from the site NNEATL, the results obtained are in agreement with the values observed with the complete dataset (Tables S3 and S4). The Cytb distribution of sequence mismatch pairwise differences showed a skewed unimodal distribution in all four BS macro areas suggesting a recent bottleneck or sudden population expansion (Fig. S2). A unimodal mismatch distribution was obtained with CR dataset in the NNEATL BS. The CR mismatch distribution of EMED, SNEATL and NNEATL BS resulted to a slightly ragged pattern (Fig. S2) that could suggest a more constant population size of the Mediterranean BS over generations.

Table 2 Analysis of molecular variance (AMOVA) of Cytochrome b (Cytb) and Control Region (CR) of the Mediterranean and North-eastern Atlantic Blue Sharks (Prionace glauca).
Both BSP analyses suggested a constant population size increase of Mediterranean and North-eastern Atlantic BS, starting more recently in the Mediterranean than in the North-eastern Atlantic (∼0.02-0.15 Mya vs 0.15-0.4 Mya; Fig. 3). Divergence time analysis based on both markers (Fig. S3) highlights a similar pattern of separation between two main groups, composed by BS from all regions, without any evidence of separation between defined geographic areas. The separation between the two clades, which is strongly supported of Posterior Probability (PP = 1.0) in both markers, is dated back to 1.24 Mya and 0.94 Mya using Cytb and control region, respectively.

DISCUSSION AND CONCLUSIONS
The BS is probably the most mobile shark species in the world (Stevens, 1990) and past research works, using both mitochondrial and nuclear markers, have struggled to find genetic structure at interoceanic scale (Sampaio da Costa, 2013;King et al., 2015;Li et al., 2016;Veríssimo et al., 2017). This high level of gene flow make it difficult to define clear BS population units. In the Pacific Ocean, the lack of structure may be the result of the combination of high potential of migration and the lack of effective barriers to gene flow (Veríssimo et al., 2017).
Experimental data have indicated that no significant genetic structure is detected in spatially distant BS samples (King et al., 2015;Li et al., 2016;Veríssimo et al., 2017). Our results revealed significant signals of geographical structuring for Mediterranean and adjacent Atlantic BS, with several frequent mtDNA haplotypes that are exclusive of the Mediterranean BS and other that are shared with the Atlantic population samples.
While both haplotype networks failed to evidence a clear geographical structure, either between Mediterranean and North Atlantic BS or within the Mediterranean, the results of AMOVA revealed a significant partition of molecular variance among all population samples and when they were grouped according to the four geographical areas with both mitochondrial markers (8.87% for Cytb and 7.93% for CR). Previous studies carrying out AMOVA on the Atlantic BS using the control region variation, showed a significance variance among groups formed by the North Atlantic BS collected from Portugal and Azores and by the South African (See Table 7 of Sampaio da Costa, 2013) or Brazilian BS (Veríssimo et al., 2017). On the contrary, the global population genetics carried out by Fitzpatrick (2012), using concatenated fragments from: 16S, tRNA, COII, ATPase and control region genes, showed no significance variation among oceans, based upon comparisons between North Atlantic and all sampling locations combined (See Table  5.7 of Fitzpatrick, 2012). Although BS exhibits high potential of dispersal and migration, our results seem to reject an absence of geographical structure in the Mediterranean and adjacent North-eastern Atlantic BS. The pairwise st analysis revealed a geographical structuring between the two Mediterranean groups and Southern North-eastern Atlantic BS, with a closer genetic similarity of the Southern North-eastern Atlantic with the Western Mediterranean BS rather than with the Eastern Mediterranean BS. This pattern of differentiation seems to suggest that reproductive movements, such as female philopatry, may occur between the Western Mediterranean and the Southern North-eastern Atlantic BS. In addition, pairwise st values highlighted that the EMED BS are the more divergent from the NATL BS Given that SNNEATL specimens are from a previously identified nursery site (Veríssimo et al., 2017), the pairwise ST values could suggest that specimens from WMED can be reproductively related to the SNNEATL, while EMED could represent a nursery site in itself (Megalofonou, Damalas & De Metrio, 2009).
Our sampling work has also preliminarily revealed significant differences between North-eastern Atlantic and Mediterranean BS by sex-ratio and size. This pattern could be the result of a sex-biased reproductive migratory behaviour that could contribute to explain the significant phylogeographical structure. Similarly, size differences were observed between WMED and EMED BS, with the large and sexually mature individuals abundant in the easternmost Mediterranean sampling location (Aegean Sea) while the sub-adult and juvenile BS frequent in the Adriatic and Ionian Seas. The great abundance of juvenile BS in the Adriatic Sea seemed to confirm the nursery role of this area for BS (Megalofonou, Damalas & De Metrio, 2009). The biological data reported in Megalofonou and colleagues (2009) describe a larger amount of big female in the easternmost Mediterranean (e.g., Aegen Sea) which is in agreement with the pattern inferred from our dataset. Conversely, using data on size and maturity stages, Kohler and colleagues (2002) observe that the majority of sharks from the Mediterranean Sea were juvenile and immature (99% of males and 98% of females; mean = 65 cm of fork length). The difference may be related to the different sampling design and fishing gear used in the studies. In fact, the majority of the data collected by Kohler and colleagues (2002) came from volunteer recreational fishermen, while the individuals from Megalofonou, Damalas & De Metrio (2009) and from this work, originated principally as by-catch from commercial fisheries, such tuna and swordfish longline.
Overall, Mediterranean and adjacent North-eastern Atlantic BS displayed a complex geographical structure in which weak but significant differences proved that a certain degree of population connectivity across the Strait of Gibraltar occurred. These results are in contrast with those obtained by tagging data in the past (Kohler, Casey & Turner, 1998;Kohler et al., 2002;Kohler & Turner, 2008;Poisson et al., 2015). Similar findings of genetic differences were observed in other shark species, more related to a benthic environment, such the small-spotted catshark, Scyliorhinus canicula, and the velvet belly lanternshark, Etmopterus spinax (Gubili et al., 2014;Gubili et al., 2016;Kousteni et al., 2015). The reported evidence of genetic structure in the blue shark analyzed in this study are associated with geographical differences in sex-ratio and size. Our results suggest BS in the NE Atlantic and the Mediterranean are not panmictic. There is still no direct observations of mating events take place in the Eastern Mediterranean, but the biological data analysis results support the Eastern Mediterranean as an important nursery area for this species (Megalofonou, Damalas & De Metrio, 2009). Such microevolutionary pattern of differentiation of Mediterranean and North-eastern Atlantic BS prompt the need for a deeper population genetic analysis on the same population samples with more powerful markers for investigating potential subtle structure of BS populations (e.g., microsatellites or SNPs) to provide robust data on BS population structure that are of priority for the BS stock management. High genetic diversity values are usually related to large population size (Frankham, 1996), and the high genetic diversity showed by both Mediterranean and North-eastern Atlantic BS at the two mitochondrial makers advocates in favour of a large size of these populations. Mediterranean and North-eastern Atlantic BS showed higher Cytb gene polymorphism than Pacific BS (Mediterranean and North-eastern Atlantic: h = 0.777-0.814; π = 0.002-0.004; Pacific: h = 0.517-0.768; π = 0.0007-0.0011, Li et al., 2016).
Based on nuclear markers, similar values of observed heterozygosity were detected between Pacific and North Atlantic BS (Sampaio da Costa, 2013;King et al., 2015;Veríssimo et al., 2017). High genetic diversity in abundant species is likely due to a combination of demographic factors, such as local population sizes, fast generation times and high rates of gene flow with other populations (Hague & Routman, 2016). The high genetic diversity shown by Mediterranean and North-eastern Atlantic BS could be a consequence of the short time elapsed, in proportion to the relative generation time, since the population started to suffer overexploitation. In fact, the abundance of the Mediterranean BS has declined by ∼78-90% over the past 30 years (Ferretti et al., 2008), approximately corresponding to three generations; the BS generation time was estimated at 8.2 and 9.8 years for South African and North Atlantic populations, respectively (Cortés et al., 2015). Furthermore, biological characters such as the large size of litters, the low nucleotide substitution rate compared to other vertebrates (Martin, Naylor & Palumbi, 1992), the high potential of migration and the high gene flow between geographical distant populations, may have affected the relationship between genetic diversity and population size, masking the sudden potential population bottleneck of the last three decades, without genetic erosion.
Otherwise, the mismatch distributions of the different macro areas appear to be slightly skewed unimodal, related to a recent bottleneck or a sudden population expansion (Fig. S2), and given the Bayesian skyline plots (Fig. 3), there is overwhelming evidence that the Mediterranean and North East Atlantic populations have undergone a constant population expansion during the last 400-200 Kya, especially within the Mediterranean samples.
The data we show here represent a novelty for the knowledge of Mediterranean blue shark, and our findings highlight the importance of the Mediterranean Sea as nursery area for this species, with direct implication to specific conservation measures for the species.
This work sheds new light on the understudied BS of the Mediterranean Sea, and emphasizes the need of conducting further population genetic surveys on this population. With ongoing efforts, (i.e., https://fishreg.jrc.ec.europa.eu/web/medbluesgen/) a greater understanding of the genetic diversity, spatial population structure and gene flow in this species will be achieved, which will enable us to devise more effective strategies for the management of this increasingly exploited ocean predator.