Paleoclimate Shaped Bluefish Structure in the Northern Hemisphere

ABSTRACTBluefish (Pomatomus saltatrix), a highly migratory cosmopolitan predator, is the only extant representative of the family Pomatomidae. It has been the subject of many studies due to its commercial and recreational value, but much less research has been conducted on its global population structure. Here we investigate the population structure of this species and the effects of present and past oceanographic barriers to dispersal in its North Atlantic, Mediterranean, Marmara, and Black sea populations. We employed mitochondrial (cytochrome b and cytochrome oxidase subunit I genes) and nuclear (eight microsatellite loci) DNA as molecular markers. Three main genetic units of Bluefish were identified: American (West Atlantic waters), Spanish (East Atlantic–Western Mediterranean regions), and Turkish (Eastern Mediterranean, Marmara, and Black seas). Our results suggested that Bluefish is panmictic in the northwest Atlantic Ocean but not in the Mediterranean Sea. The common ancestor of the studied popula...


INTRODUCTION
Ocean currents and the apparent lack of physical barriers in the marine realm seem to facilitate extensive gene flow among marine fish populations (Palumbi 1994). Pelagic and demersal fishes are expected to exhibit little intraspecific genetic structuring even over large geographic distances (Ward et al. 1994), and marine environments are often seen as open habitats in which isolation by distance is the main mechanism that may promote speciation (Palumbi 1994). However, several studies have demonstrated the existence of marine physical barriers that produce intraspecific genetic fragmentation in marine systems. Population structuring of highly migratory marine fish can be promoted by currents (Machado-Schiaffino et al. 2010), salinity gradients (Nielsen et al. 2004), temperature boundaries (Crow et al. 2007), convergence of distinct water masses (Borrero-Perez et al. 2011), behavior (Campos-Telles 2011), or historical past events (Shen et al. 2011). Some barriers to gene flow are well defined by coastal shapes and features. Examples are the Gibraltar Strait for some Sparidae species (Bargelloni et al. 2003), the Siculo-Tunisian Strait for Sea Bass (Bahri-Sfar et al. 2000), the Florida Keys for gobies (Avise 1992), and the hydrographic isolation of the Aegean and Ionian and Adriatic seas for numerous species (Patarnello et al. 2007;Perez-Losada et al. 2007). On the other hand, behavioral traits such as homing can also account for population differentiation in some species (e.g., Castillo et al. 2005). Finally, though the genetic variability and population genetic structure of a species are shaped by both present and historical marine barriers, they are also affected by paleoecological history (Patarnello et al. 2007;Perez-Losada et al. 2007). Hence, to better understand speciation mechanisms in the marine environment, it is important to characterize not only population dynamics and structure but also life history strategies, environmental factors from past and present, and physical barriers to dispersal (Zardoya et al. 2004).
El paleoclima como determinante de la estructura poblacional del anjova en el hemisferio norte RESUMEN: la anjova (Pomatomus saltatrix), un depredador cosmopolita y altamente migratorio, es el único representante vivo de la familia Pomatomidae. Ha sido sujeto de numerosos estudios dado su valor comercial y recreativo, pero poco se sabe acerca de su estructura poblacional a nivel mundial. En este trabajo se investiga la estructura poblacional de esta especie y los efectos que tienen las barreras oceanográficas pasadas y presentes en la dispersión de  Bluefish (Pomatomus saltatrix) is a cosmopolitan, migratory, pelagic predator distributed over continental shelves and in estuaries of temperate waters of the Atlantic, Indian, and Pacific oceans and adjacent seas, including the Mediterranean, Aegean, and Black seas (Briggs 1960;Tortonese 1986;Pottern et al. 1989;Juanes et al. 1996). As a consequence of its broad distribution and the existence of potential oceanographic barriers, the species may be composed of multiple different populations, but there is limited information available (Goodbred and Graves 1996;Turan et al. 2006;Pardiñas et al. 2010). Bluefish is one of the most important recreational and commercial species along the east coast of the United States (Robillard et al. 2008) and is the target species of an important artisanal fishery in all Turkish seas: Marmara, Aegean, and Black (Ceyhan et al. 2007). Failure to detect and recognize population units can lead to local overfishing and ultimately to severe declines of fisheries (Hutchings 2000;Knutsen et al. 2009); thus, an accurate definition of Bluefish population structure is particularly important and necessary for fisheries management (Utter 1991;Wilson 2003).
The contemporary distribution of Bluefish is coincident with sea surface temperatures of 18-27°C (Juanes et al. 1996), and it has been suggested that shifts in its ranges and contacts between populations have resulted from historical changes in water temperature (Goodbred and Graves 1996). The sensitive behavioral response of Bluefish to temperature variations may provide new insights into the evolutionary consequences of the glacier-interglacier cycles and migrations into refuges for marine migratory species.
The objective of this study was to document the present population structure of Bluefish in its northern distribution across the North Atlantic ocean and Mediterranean, Marmara, and Black seas; identify possible ocean barriers to dispersal; and reconstruct the phylogeography of Bluefish to understand the role of climate for determining historical and present barriers to gene flow along the Atlantic and Mediterranean basins.

Sampling
A total of 120 samples of Pomatomus saltatrix collected from eight different locations ( Figure 1) between 2004 and 2009 were analyzed. These samples represent at most two consecutive generations because the maturity age is 2.4 in males and 1.9 in females (Dhieb et al. 2006). Four locations were in the northwest (NW) Atlantic Ocean, on the U.S. coast-New Jersey, Maryland, North Carolina, and Florida-and four locations were in the East Atlantic and in the Mediterranean basin-Cadiz (Atlantic Spanish coast) and Barcelona (West Mediterranean, Spanish coast) and Canakkale and Istanbul (Marmara and Black seas, Turkish coast).

DNA Extraction, Amplification, and Sequencing
DNA was extracted with Chelex (Bio-Rad, Berkeley, CA) following Estoup et al. (1996). Two mitochondrial genes and eight microsatellite loci were amplified. The mitochondrial cytochrome b (cyt b) gene sequences were obtained following the protocol described by Kocher et al. (1989), with the primers H151 and L148 described therein. The cytochrome oxidase subunit I (COI) gene was amplified with primers designed for Pomatomus saltatrix-COI-R-Pom: 5′-AAGAATGGGGTCTCCTCCAC-3′ and COI-F-Pom: 5′-TTGGTGCATGAGCTGGTATG-3′-with the software PRIMER3 (Rozen and Skaletsky 2000). Several potential primer sets were generated; we selected one set of primers that covered the maximum number of base pairs. The polymerase chain reactions (PCRs) to obtain the COI sequences were performed using the GeneAmp PCR system 2700 (Life Technologies, Waltham, MA). Total reaction volume was 40 μL and the reaction mix contained approximately 50 ng of DNA, 20 pmol of each primers 10 mM Tris-HCL pH 8.8, 250 μM of each dNTP, 5 U μL-1 of DNA Taq polymerase (Promega, Madison, WI), and 2.5 mM MgCl 2 . The PCR conditions were as follows: initial denaturing at 95°C for 5 min, 35 cycles of denaturing at 95°C for 30 s, annealing at 58°C for 30 s, an extension of 72°C for 30 s, and a final extension at 72°C for 20 min. PCR products were visualized in 2% agarose gels with 10 mg/mL of ethidium bromide. Stained bands were excised from the gel and DNA fragments were purified with an Eppendorf Perfect-Prep Gel CleanUp kit (Hamburg, Germany). Then, purified DNA was precipitated using standard 2-propanol precipitation and resuspended in formamide prior to sequencing. Sequencing was performed in an ABI PRISM 3100 Genetic Analyzer (Life Technologies), with a BigDye 3.1 Terminator system, in the Unit of Genetic Analysis of the Scientific-Technical Services of the University of Oviedo (Spain). Eight tetranucleotide microsatellites were assayed: elf 17, elf 19, elf 37, elf 39, elf 44, elf 46, elf 49, and elf 50 (Dos Santos et al. 2008). Reactions and conditions were modified from those described by the authors. PCRs consisted of the following: 95°C for 5 min, 35 cycles of 95°C for 30 s, the annealing temperature (Table 1) for 30s, 72°C for 30s, followed by 72°C for 20 min. PCRs were performed in a total volume of 20 μL in individual PCR reactions for each locus with the GeneAmp PCR system 2700 and 2720 Thermal Cycler (Applied Biosystems). Products were visualized in 50 mL 2% agarose gels with 2.5 μL of ethidium bromide (10 mg/mL) for verification and genotyped using an ABI PRISM 3100 Genetic Analyzer (Applied Biosystems), with GS500 LIZ 3130 size standard, in the Unit of Genetic Analysis of the Scientific-Technical Services of the University of Oviedo (Spain).

Genetic Diversity
Sequences of each mitochondrial gene (cyt b and COI) were aligned using ClustalW (Thompson et al. 1994) from the BioEdit Sequence Alignment Editor (Hall 1999) and were visually inspected to avoid base-calling errors. The incongruence length difference tests (Farris et al. 1995) were performed in PAUP* 4.0 (Swofford 1999) with 1,000 replicates and a p-value of 0.05. The two mtDNA gene sequences were concatenated and haplotypes defined with DNAsp v.4.50.3 (Rozas et al. 2003). MtDNA haplotype (H) and nucleotide diversity (pi) were calculated for each location using the software Arlequin version 3.0 (Excoffier et al. 2005).
Microsatellite allele sizes were estimated using the Gen-eMapper Software Version 4.0 (Applied Biosystems). Loci with scoring errors, large allele dropout, and null alleles were discarded employing the program MICROCHECKER (Van Oosterhout et al. 2004). Conformity to Hardy-Weinberg equilibrium was calculated using the exact probability test with GENEPOP software (Raymond and Rousset 1995). Microsatellite variation (number of alleles per locus, allelic richness, and observed and expected heterozygosity) was calculated with the programs GENETIX Version 4.03 (Belkhir et al. 2004) and FSTAT Version 2.9.3.2 (Goudet 2001).

Population Differentiation and Structure
A median-joining (Bandelt et al. 1999) haplotype network was constructed using the concatenated mtDNA genes to visualize the intraspecific relationships of the different haplotypes and their relative frequencies in the sampled populations with the program Network 4.5.1.6 (fluxus-engineering.com) with default settings. Network software reconstructed all possible, shortest, least complex, phylogenetic trees (maximum parsimony) from a data set under different algorithms.
Genetic divergence between populations was estimated using population pairwise F ST values calculated using nuclear and mtDNA data with the program Arlequin version 3.0 (Excoffier et al. 2005). The statistical significance of F ST values between samples was calculated with 1,000 permutations and 10,000 steps in a Markov chain. This software was also employed for analysis of molecular variance (locus by locus and standard AMOVAs) using 1,000 permutations.
Fu's F S test (Fu 1997) for selective neutrality was calculated in Arlequin v.3.0. For neutral markers, this test can be employed to detect changes in population size. Significant and negative F S values can be interpreted as signatures of population expansion (Dodson et al. 2007). Mismatch analysis was also used to explore the spatial and demographic evolution of the studied Bluefish populations with the raggedness index (Harpending 1994) and the sum of squared deviations (SSD; Schneider and Excoffier 1999). Both demographic and spatial mismatch analysis were calculated with Arlequin v.3.0 and were based on the null hypothesis of expansion; thus, nonsignificant values reveal population expansion.
Population structure across the study area was assessed using the program STRUCTURE 2.3.1 (Pritchard et al. 2000) by using microsatellite loci data. This software estimates the minimum number of population units with genetic identity in the data set under a Bayesian framework. This data set was analyzed under the "Admixture model," which assumes that individuals may have mixed ancestry. The parameter set consisted of a burn-in period of 30,000 steps followed by 300,000 Markov chain Monte Carlo (MCMC) iterations and seven runs for each K (number of genetic units estimated); the number of K was estimated following Evanno et al. (2005).
Mantel tests of association between genetic differentiation values (pairwise F ST values) and geographical Euclidean distances (linear distances between sampling locations) were calculated with Arlequin version 3.0 (Excoffier et al. 2005). Mantel tests were done to determine whether the considered populations follow an isolation-by-distance model. To determine possible geographical barriers to dispersal, we employed the software BARRIER v.2.2 (Manni et al. 2004), which can identify geographically continuous and discontinuous assemblages of samples from a spatial landscape. Geographical coordinates of each sampling area were mapped into a matrix connected by Delauney triangulation (Brassel and Reif 1979). Barriers in the triangulation were identified using mitochondrial and nuclear pairwise F ST distances. The analysis employed was based on Monmonier's maximum distance algorithm (Manni et al. 2004) to identify regions with sharp genetic change or discontinuity.

Phylogeny and Evolutionary History of Pomatomus saltatrix
Population divergence time estimations were done under a Bayesian MCMC framework using the two mtDNA genes with the software BEAST version 1.6.1 . Following a burn-in of 3 million cycles, rates were sampled once every 1,000 cycles from 30 million MCMC steps for an extended Bayesian skyline tree prior with a stepwise model for mitochondrial DNA and strict clock model. Bayesian intraspecific phylogenies are based on coalescent theory (Kingman 1982) and allow the inference of past population dynamics and parameters from contemporary gene sequences. The best evolutionary model of both sequences (COI and cyt b) and their priors (kappa, gamma-shape, proportion of invariant sites, etc.) were inferred by using jModelTest software version 0.11 (Posada 2008) and its implementation of the Akaike information criterion (Akaike 1974). The mutation rates employed were 2% per million years (MY) for cyt b (Brown et al. 1979) and 1.2% per MY for COI (Bermingham et al. 1997). Tracer version 1.5 ) was used to check that chains had converged to a stationary distribution. The analysis was repeated with longer runs (50 million MCMC steps) when data sets did not accomplish this condition.

Genetic Diversity
A total of 46 haplotypes for the concatenated COI-cyt b genes were detected among the studied samples. Most (68.75%) were observed solely among U.S. samples (Northwest Atlantic Ocean) and the rest were from the eastern area (Table 1; GenBank ID: JQ039400-JQ039435for COI and JQ039436-JQ039465 for cyt b haplotypes). High haplotype diversity and low nucleotide diversity were found in general for all regions, but differences between localities revealed a gradient from west to east across the studied area. U.S. samples exhibited higher diversities than European ones, and western Mediterranean samples were more diverse than those from the eastern Mediterranean (Table 1). The highest number of haplotypes and haplotype and nucleotide diversities corresponded to northwest Atlantic samples.
The eight microsatellite loci considered were amplified for the same 120 individuals across the sampling area. MICRO-CHECKER did not detect dropouts or scoring errors, but null alleles were found for one loci (elf 44) in different populations. Therefore, elf 44 was excluded from the data set. High genetic variability was found at the seven loci. The number of alleles per locus ranged from 18 (elf 17) to 29 (elf 37). All sampling areas were in Hardy-Weinberg equilibrium and there were no significant differences between expected and observed heterozygosities in any location (Table 1).

Population Differentiation and Population Structure
Differentiation between the two sides of the Atlantic Ocean was observed in the haplotype network (Figure 2), where the main separation between mitochondrial lineages corresponded clearly to the geographical differentiation between the two sides of the Atlantic Ocean, giving separate U.S. and European clades, with further internal division. The presence of unique haplotypes may be caused by recent differentiation, and geological or demographic factors and mutations would explain the haplotype structure (Verheyen et al. 2003).
The population structure of Bluefish based on both mitochondrial and nuclear DNA was consistent with three different genetic units in the sampling area: northwest Atlantic and west and east Mediterranean. F ST population pairwise comparisons revealed significant differences between northwest Atlantic and other locations and also between east and west Mediterranean samples (Table 2) for both types of markers. The northwest Atlantic group consisted of all of the North American localities: Florida, Maryland, North Carolina, and New Jersey. The west Mediterranean cluster included Cadiz (from the Atlantic Ocean) and Barcelona (west Mediterranean), and the east Mediterranean group was composed of the two Turkish localities, Istanbul and Canakkale.
Bayesian analysis confirmed that Bluefish from this study belonged to three different genetic units (K = 3), corresponding to the same three clusters detected with the genetic distances but with a moderate level of admixture between them ( Figure  3). There is also a secondary maximum at K = 8, supported by much lower likelihood than K = 3, suggesting a hierarchical island model (Evanno et al. 2005) of eight different populations clustered in three main sets. The AMOVA confirmed these three genetic units. Both AMOVAs, standard (for mtDNA) and locus by locus (for microsatellite loci; Table 3), showed significant differences among the three groups and within populations but not among populations within the three previously defined groups.
In the standard AMOVA, the highest percentage of variation was observed between groups (71%; P < 0.001), whereas the highest percentage of variation in the locus by locus AMOVA was within populations (93.72%; P < 0.001). These results were consistent with STRUCTURE and mitochondrial and nuclear F ST results.
Demographic analyses indicated that all of the populations except Canakkale were in expansion (Table 4). In the Mantel test, a high and significant correlation was detected between genetic and Euclidean distances between samples (r = 0.734; P < 0.001), suggesting a pattern of isolation by distance for Pomatomus saltatrix. The software BARRIER (Manni et al. 2004) identified two main boundaries coincident with the two longer geographical distances between samples. The two detected barriers were supported with both types of molecular markers ( Figure 4). However, all of the previous population structure results (STRUCTURE, mitochondrial, and nuclear F ST and AMOVAs) suggested some permeability across the three genetic units. The Gibraltar Strait was not a barrier to gene flow for Pomatomus saltatrix because the Spanish samples on the two sides of the strait (Cadiz and Barcelona) were not significantly different.

Phylogeny and Evolutionary History of Pomatomus saltatrix
The time to the most recent common ancestor (TMRCA) of the sampled Bluefish was estimated to be 480,000 years ago (95% highest posterior density [HPD] = 0.28-0.72 MY), and the TMRCA for the NW Atlantic Ocean was dated as 252,000 years (95% HPD = 0.13-0.39 MY) and 148,300 years (95% HPD = 0.05-0.27 MY) for the Mediterranean samples. The population growth rate estimated with the extended Bayesian skyline model for the U.S. samples reached its maximum approximately 40,100 years ago, whereas for the European clade the maximum occurred 23,420 years ago ( Figure 5). Both estimates suggest that the northwest Atlantic clade is more ancient than the Mediterranean. The Bayesian tree obtained for the northwest Atlantic Ocean suggested that there are two separate American lineages that are composed of a mixture of all localities. The TMRCA of the two American lineages was estimated to be approximately 175,000 and 131,000 years. The substructure detected with the Bayesian methods and not with F ST is likely due to the more sophisticated approaches of Bayesian inference. F ST -based approaches are well understood, widely used, and easily applied, but Bayesian-based analyses allow   . Each vertical bar represents one individual. Acronyms of sampling localities are below the plot and described in Figure 1 and Table 1. Our results support two genetic discontinuities for Bluefish: one in the middle of the Atlantic Ocean and the other in the Mediterranean Sea. However, regional permeability and migration can be inferred to occur across both of them. The Mediterranean barrier could be at the level of the Siculo-Tunisian Strait and may be due to the hydrographic isolation of the Aegean and Ionian and Adriatic seas, or either of them, but not in the Gibraltar Strait or Alboran Sea. The Siculo-Tunisian Strait is a barrier to gene flow for other species (Stefanni and Thorley 2003;Zardoya et al. 2004), and in most cases the major genetic break or limitation to larval dispersal between the Eastern and Western Mediterranean occurs in the straits separating the Adriatic, Aegean, and/or Black seas (Nikula and Vainola 2003;Costagliola et al. 2004;Domingues et al. 2005;Peijnenburg et al. 2006;Perez-Losada et al. 2007;Zulliger et al. 2009). The accurate locations of the barriers for this species cannot be deduced from the present results. more precise estimates (see Pearse and Crandall [2004] for a review).

DISCUSSION
This study provides a genetic analysis of North Atlantic and Mediterranean Bluefish populations. Three clear genetic units were identified and may follow a hierarchical island model. The North American samples were clustered together, without any significant difference between sampling sites, suggesting that Bluefish is panmictic in the northwest Atlantic Ocean. Many studies based on distribution and morphological characteristics had investigated the number of stocks in the East Coast of the United States. In previous studies, the number of stocks identified ranged from six to two (Lund 1961;Lassiter 1962;Lund and Maltezos 1970), but other investigations concluded that either two distinct spawning groups or one stock with two distinct spawning periods occurred in the region (Norcross et al. 1974;Kendall and Walford 1979;Hare and Cowen 1993;Smith et al. 1994). Later, Graves et al. (1993) concluded that there was a single population based on mitochondrial DNA, and a morphometric study (Austin et al. 1999) corroborated the one-stock hypothesis despite the evidence of phenotypic plasticity. Our study confirms the one-stock hypothesis based on both nuclear and mitochondrial DNA, but the Bayesian mitochondrial DNA results suggested that there are two different lineages within the American panmictic population that may have evolved at different times.
Northwest Atlantic Bluefish exhibited higher diversity than European populations, and west Mediterranean populations were more diverse than east Mediterranean populations. Hewitt (1996Hewitt ( , 2000 demonstrated that the last glaciations had a profound effect on genetic diversity. If genetic diversity remains high or increases at the limits of the range, then these populations may be composed of two or more lineages that differentiated in distinct glacial refugia, which is consistent with the two U.S. lineages detected in our study. Therefore, the occurrence of high levels of genetic diversity in previously glaciated areas would suggest colonization from multiple refugia (Griswold and Baker 2002). This could be the case of the northwest Atlantic Bluefish population, which has the highest genetic diversity and was dated as the oldest population, having necessarily passed through several glacial cycles (Riss, Würm). On the opposite side of the Atlantic, the emergence of the Mediterranean clade occurred during the Riss glacial cycle, a harsh period when the Gibraltar Strait was closed due to the decrease in sea level (Loget and van den Driessche 2006). Later, the Mediterranean Sea was closed and part of the Bluefish population was probably isolated for a long period. Further radical hydrographical changes that occurred during the course of the late Pleistocene climatic cycles, such as fluctuations in water level that repeatedly isolated totally or partially the Black Sea, Aegean Sea, and Eastern Mediterranean (Svitoch et al. 2000) may have contributed to differentiate those populations. As reported for other marine species (e.g., Borrero-Perez et al. 2011), this isolation may have been maintained until the present due to the anticyclonic front located in the Peloponnese peninsula (Millot 2005). The Atlantic European clade may have colonized the western Mediterranean after the opening of the Strait of Gibraltar and continued its expansion until the present. The paleogenetic history of Bluefish seems therefore associated with glacial-interglacial cycles, as suggested by the estimated TMRCA of the different clades. TMRCA for the Bluefish analyzed here was dated to the Aftonian II, an interglacial period of favorable conditions with growth and expansion of species, and the TMRCA of both south European and east Mediterranean clades corresponded to glacial periods (fourth cycle of Mindel and Riss glaciation, respectively). Bluefish populations may have migrated to different refuges and stayed isolated for long periods until the climatic conditions became favorable again.  The most probable factor influencing Bluefish population genetic structure seems to be geographical, such as the distance to continents and shallow waters. Small fish are commonly found in shallow coastal waters (e.g., Juanes et al. 1996), and apparent population genetic isolation associated with different continents has been reported before (Goodbred and Graves 1996). Bluefish is a highly migratory species (Juanes et al. 1996). Water temperature and photoperiod are major factors described to influence movement patterns in Bluefish (Olla and Studholme 1972;Wilk 1977;Juanes et al. 1996), but other factors, like salinity, could contribute to shape their population structure. We have compared annual sea surface temperature as well as annual temperature at 200 m depth (maximum depth for Pomatomus saltatrix) with the Bluefish genetic structure detected in this study. In this case, changes in temperature might not explain the barriers detected here or the genetic structure of different Bluefish populations. Similarly, salinity gradients at the surface and at 200 m depth did not explain the genetic structure and barriers detected. Neither factor showed differences along the genetic units detected, suggesting that those factors are not a barrier in our sampling area. In contrast, both annual salinity and water temperature along the North Atlantic Ocean and Mediterranean Sea could explain the permeability across the two boundaries detected and could allow Bluefish to migrate. This connection across barriers could prevent vicariant speciation, because Pomatomus saltatrix is the sole member of its genus and family.
In conclusion, the present study suggests that the paleoecological history (e.g., glaciations and interglacial periods in the late Pleistocene) has been crucial for shaping the present genetic variability and population structure of Pomatomus saltatrix in the North Atlantic Ocean and in the Mediterranean basin.

FUNDING
This work was supported by the Spanish National Grant CGL2009-08279. U.S. sampling was funded through the Bluefish/Striped Bass Research Program (NOAA). Laura Miralles holds a PCTI Grant from the Asturias Regional Government, referenced BP 10-004.

SUPPLEMENTAL MATERIAL
Supplemental data for this article can be accessed on the publisher's website at www://dx.doi.org/10.1080/03632415 . 2014.976701 Figure 5. Maximum population growth rates for the studied Bluefish during the last 50,000 years. Population growth rates (in percentage) estimated with the extended Bayesian skyline of BEAST software for the European (East Atlantic and Mediterranean) and American (West Atlantic) Bluefish clades.