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

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.

Material and Methods

Blue shark sampling

Sampling was carried out between 2003 and 2016 and tissue specimens and individual data were collected by means of commercial fishermen and scientific surveys. Mediterranean BS were collected from multiple locations in the Eastern (Central Adriatic, CADR, 21; Ionian Sea, IONI, 15; Aegean Sea and Levantine Sea, AEGE, 20) and Western areas (South Tyrrhenian, TYRR, 10; Ligurian Sea, LIGU, 57; Balearic Islands, BALE, 42). North Atlantic BS were caught from the North Eastern Atlantic Ocean off the coasts of Portugal (SNEATL, 33) and Celtic Sea (NNEATL, 16) (Fig. 1). A total of 214 BS individuals were collected (N = 91 males, N = 101 females and N = 22 unsexed) (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).

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.

Results

Among sexed individuals (N = 192; Table S1), BS females significantly outnumbered males in the NA samples (sex-ratio 0.34, χ² test: 10.256 P_2tail = 0.001; P_1tail = 0, d.f. 1) while in the two Mediterranean BS groups a weak and not significant predominance of males was observed (WMED: 1.19, χ² test: 0.786 P_2tail = 0.375; P_1tail = 0.188, d.f. 1; EMED: 1.09, χ² test: 0.087 P_2tail = 0.768; P_1tail = 0.384, d.f. 1). Sized BS (N = 209) were composed by 63 juvenile, 82 young and 64 large individuals (Table S1). In the NA and WMED the young BS (TL = 120–180 cm; 48% and 50%, respectively) were predominant, while in the EMED a large predominance of juveniles was observed (TL ≤ 120 cm; 63%). Noticeably 67% of the BS sampled in the Ionian Sea and 95% of those sampled in the Adriatic Sea were juveniles. Large BS are similarly represented in the geographical groups with percentages varying from 25% (EMED) to 34% (NA), full details presented in Table S1.

A total of 207 and 170 BS individuals were sequenced for Cytb (762 bp) and CR (747 bp), respectively. Haplotype sequences (Cytb, N = 23 and CR, N = 55) were deposited in GenBank under the accession numbers MG515900–MG516106 and MG545732–MG545901 for Cytb and control region, respectively.

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.

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.

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.

Supplemental Information