Genetic structuring and migration patterns of Atlantic bigeye tuna, Thunnus obesus (Lowe, 1839)

Background Large pelagic fishes are generally thought to have little population genetic structuring based on their cosmopolitan distribution, large population sizes and high dispersal capacities. However, gene flow can be influenced by ecological (e.g. homing behaviour) and physical (e.g. present-day ocean currents, past changes in sea temperature and levels) factors. In this regard, Atlantic bigeye tuna shows an interesting genetic structuring pattern with two highly divergent mitochondrial clades (Clades I and II), which are assumed to have been originated during the last Pleistocene glacial maxima. We assess genetic structure patterns of Atlantic bigeye tuna at the nuclear level, and compare them with mitochondrial evidence. Results We examined allele size variation of nine microsatellite loci in 380 individuals from the Gulf of Guinea, Canary, Azores, Canada, Indian Ocean, and Pacific Ocean. To investigate temporal stability of genetic structure, three Atlantic Ocean sites were re-sampled a second year. Hierarchical AMOVA tests, RST pairwise comparisons, isolation by distance (Mantel) tests, Bayesian clustering analyses, and coalescence-based migration rate inferences supported unrestricted gene flow within the Atlantic Ocean at the nuclear level, and therefore interbreeding between individuals belonging to both mitochondrial clades. Moreover, departures from HWE in several loci were inferred for the samples of Guinea, and attributed to a Wahlund effect supporting the role of this region as a spawning and nursery area. Our microsatellite data supported a single worldwide panmictic unit for bigeye tunas. Despite the strong Agulhas Current, immigration rates seem to be higher from the Atlantic Ocean into the Indo-Pacific Ocean, but the actual number of individuals moving per generation is relatively low compared to the large population sizes inhabiting each ocean basin. Conclusion Lack of congruence between mt and nuclear evidences, which is also found in other species, most likely reflects past events of isolation and secondary contact. Given the inferred relatively low number of immigrants per generation around the Cape of Good Hope, the proportions of the mitochondrial clades in the different oceans may keep stable, and it seems plausible that the presence of individuals belonging to the mt Clade I in the Atlantic Ocean may be due to extensive migrations that predated the last glaciation.


Background
Marine pelagic fishes show broad geographic distribution, large population sizes, and highly migratory movements that are thought to ultimately result in little genetic structuring (e.g. [1][2][3][4]). The above-mentioned biological peculiarities result in complex phylogeographic patterns (e.g. [5,6]), which do not often meet the assumptions (e.g. uniform effective population sizes or symmetric effective migration rates) of classic statistics that measure gene flow (F ST ) [7][8][9], and thus make the study of population genetic variation of marine pelagic fishes particularly challenging. In addition, philopatric behavior, local larval retention mechanism, and historical processes such as e.g. the effect of Pleistocene sea level changes also complicate estimation of genetic differentiation for many marine pelagic fishes [2,10,11].
Given this complexity, accurate determination of marine pelagic fish genetic structuring would ideally require both temporal (different years) and spatial (covering large areas) sampling, analysis of multiple loci (including both mitochondrial, mt, and nuclear data) that improves the estimates considerably [12,13], as well as complementation of classic summary statistics with more recent coalescent-based approaches [14] that convey statistical efficiency and flexibility to the estimation of population genetic parameters by using maximum likelihood (ML) [15][16][17][18] and Bayesian inference (BI) [19][20][21].
Bigeye tuna, Thunnus obesus (Lowe, 1839) is a marine pelagic fish species characterized by large populations and a worldwide distribution (Atlantic Ocean, Indian Ocean and Pacific Ocean) restricted to tropical and subtropical waters (except the Mediterranean Sea) [22,23]. Although capable of long distance movements, conventional and archival tagging indicate regional fidelity for bigeye tuna to geographical points of attraction [24,25]. Catch data from surface gears [22] indicate that the main breeding and nursery area of Atlantic bigeye tuna is located in the Gulf of Guinea, whereas adult feeding grounds are placed at both northern and southern temperate areas (Fig. 1). The exact timing of spawning migrations towards equatorial areas, and whether the adults return to their original northern and southern feeding zones after spawning remain open questions.
Several studies [10,[26][27][28], reported geographic heterogeneity of bigeye tuna mt lineages within and among oceans. Both mt sequence and RFLP data revealed the existence of two highly divergent groups: Clade I that is present both in the Atlantic Ocean and Indo-Pacific Ocean, and Clade II that is almost exclusive to the Atlantic Ocean [10,[26][27][28]. The origin of the two mt clades has been related to temperature fluctuations during the Pleistocene [29] that temporarily impeded bigeye tuna migra-tion around the Cape of Good Hope [26]. Afterwards, unidirectional gene flow of mt Clade I from the Indo-Pacific Ocean into the Atlantic Ocean (promoted by the strong Agulhas Current that flows westward off southern Africa) during an interglacial period would have resulted in secondary contact and the contemporary asymmetrical distribution of the mt clades [10,[26][27][28].
The contemporary presence of individuals of both mt clades in sympatry across the Atlantic Ocean but not in the Indo-Pacific Ocean poses several interesting questions that require further investigation using nuclear data. If individuals of mt Clade I and Clade II show high genetic divergence at the nuclear level, the two lineages could represent distinct cryptic species. Alternatively, individuals from both mt clades could interbreed in the Atlantic Ocean, a process that would render homogeneity at the nuclear level. On the other hand, the presence of individuals of mt Clade I in both the Atlantic Ocean and Indo-Pacific Ocean allows testing whether gene flow around the Cape of Good Hope is still occurring at present within this clade, as well as the directionality of the process.
Besides the above-mentioned historical factors, contemporary ecological factors such as e.g. phylopatric behaviour of adults could play an important role in shaping the current genetic structure of bigeye tuna. In order to discern the relative role of historical and ecological factors, the comparison of molecular markers of mt and nuclear origin is particularly informative. For instance, discordance between mt and nuclear markers in describing population structure have been used to asses sex-biased dispersal or philopatry [30][31][32]. In the case of the bigeye tuna, the comparison of allele-size frequencies at one mt and four nuclear loci in individuals from the Atlantic Ocean and Indo-Pacific Ocean tentatively rejected the hypothesis that the two mt clades could correspond to two cryptic species [33]. However, the population genetic analyses were not conclusive regarding gene flow around the Cape of Good Hope since allele frequencies of two nuclear loci supported the hypothesis of gene flow disruption between the Atlantic Ocean and the Indo-Pacific Ocean whereas those of another two nuclear loci failed to support the potential genetic break [33].
In this study, we gathered the most comprehensive genetic data set on Atlantic bigeye tuna to date. We analyze allele size variation of nine polymorphic microsatellite loci in four Atlantic locations that were surveyed in a previous study [10], and analyzed using mitochondrial control region sequence data. For a comparison of gene flow between ocean basins, we added one sample from the Pacific Ocean and another from the Indian Ocean. Microsatellite allele size data were analyzed with a variety of population genetic methods including hierarchical F-sta-tistics, Mantel tests to determine for isolation by distance, Bayesian estimation of the number of populations, and coalescence-based estimation of the magnitude of gene flow. The main goal of these population genetic analyses was to determine genetic structure of bigeye tuna within the Atlantic Ocean based on nuclear data, and to compare it with mitochondrial evidence. We also tested the null hypothesis of genetic homogeneity (i.e. panmixia) worldwide, and estimated the magnitude and direction of gene flow between the Atlantic Ocean, and the Indo-Pacific region.

Temporal genetic stability within Atlantic bigeye tuna populations
Three locations in the Atlantic Ocean, namely Canada, Canary islands, and Guinea were sampled in two different years (2001 and 2003) in order to test temporal stability of genetic variability. After Bonferroni correction, 12 out of 54 exact tests remained significant, which indicates their departure from Hardy-Weinberg equilibrium (HWE) ( Table 1). Of these, only one locus exact test (TA161 at Guinea2) showed heterozygote excess characterized by a negative F IS value. According to the geographical origin of the samples, most (83%) of F IS significant values were found to belong to Guinea and Guinea 2 sampling sites (Table 1). To test for differences of F IS values between different years and loci, we performed a repeated measure Locations of bigeye tuna samples in the Atlantic Ocean, Indian Ocean and Pacific Ocean (solid circles) Figure 1 Locations of bigeye tuna samples in the Atlantic Ocean, Indian Ocean and Pacific Ocean (solid circles). The textured grey areas show potential feeding grounds at the North and the South, respectively. The arrows indicates the Agulhas Current and its retroflection at the Cape of Good Hope (modified from [29]). Dashed dark lines show the limits of bigeye tuna geographic distribution. The relative proportions of mitochondrial control region Clade I (black) and Clade II (white) for each population of T. obesus at the Atlantic Ocean, the Indian Ocean and the Pacific Ocean [10] are also indicated. ANOVA with year and locus as repeated measures. Both year and locus were modeled as random factors. The selection of the final model was based on the Akaike Information Criterion (AIC) [34]. Two models were considered significantly different when their AIC difference was larger than 2 (AIC ≥ 2.0) [35]. The best fitting model obtained (lowest AIC and significant ΔAIC ≥ 2.5) included the factor locus only, whereas year and the interaction between year and locus were not significant (ΔAIC ≤ 1.7). Furthermore, the analyses of the F ST and R ST pairwise comparisons between first and second year replicate comparisons at Guinea, Canary Islands, and Canada were each not significant ( Table  2).

Genetic variability and population structure among Atlantic bigeye tuna population samples
In order to characterize population genetic variation among Atlantic bigeye tuna, analyses were performed based on four Atlantic Ocean locations (the above mentioned three plus Azores), and included also one Indian Ocean, and one Pacific Ocean locations as outgroups.
Genetic analyses were based only on those individuals that were correctly genotyped for at least seven loci (only nine out of 380 samples did not met this requirement and were excluded). The amount of genetic variability, in terms of average number of alleles (N S ) and observed heterozygosity was similar among sampling sites for the same microsatellite locus (Table 1). However, there were large differences among loci, with values of average observed heterozygosity per locus ranging from H T = 0.675 at locus TTH217 to H T = 0.908 at locus TA121. Overall, the average observed heterozygosity per locus and population was high (H T = 0.742). The total N S varied from 7.16 (for TTH217) to 18.77 (TTH4) ( Table 1). After Bonferroni correction, 7 out of 54 locus exact tests remained significant, which supports their departure from HWE. All showed significant heterozygote deficits in three (locus TTH208), and one (loci TA102, TA113, TA117, and TA161) sampling sites. Only locus TA161 deviated from HWE due to null alleles as detected using a null allele test based on expected homozygote and heterozygote allele size difference frequencies [36]. Since taking this locus out of the analyses did not qualitatively affect population comparisons (data not shown), it was included in further analyses. Tests for linkage disequilibrium did not show any significant value for any of the comparisons.
The allele size permutation test rendered non-significant differences between F ST and R ST estimates (P = 0.3). Multilocus pairwise estimates of F ST showed nine significant comparisons after Bonferroni correction (Table 2). Pacific and Indian Ocean pairwise comparisons with Guinea and Canada sampling sites, but not with Azores and Canary Islands, were each significant. Within Atlantic Ocean comparisons were generally not significant. None of the R ST pairwise comparisons were significant after Bonferroni correction (Table 2). According to a Mantel test, we found no significant correlation between genetic and geographical distances for both Atlantic Ocean samples or across the entire studied geographic range when using F ST (R 2 = 0.23, P = 0.71; R 2 = 0.66, P = 0.99, respectively) or R ST (R 2 = 0.23, P = 0.71; R 2 = 0.66, P = 0.99) estimates ( Fig. 2).   The hierarchical AMOVA revealed overall significant genetic structuring of the analyzed samples ( Table 3).
Most of the total genetic variance was found within populations. The hypothesis that nuclear variation could be structured according to mt clades was clearly rejected (Table 3). Geographic structuring of nuclear variation according to ocean basins (Atlantic Ocean, Indian Ocean, and Pacific Ocean) was also not supported by the AMOVA (Table 3), even when the Indian Ocean and the Pacific Ocean were grouped together (Indo-Pacific region) an analyzed against the Atlantic Ocean (result not shown).
The number of populations (and the assignment of individuals to each population) was estimated using Bayesian inferences. Although the prior parameters for the F model (gamma distribution with mean 0.01 and standard deviation 0.05) were chosen to allow the existence of two possible populations with very similar allele frequencies [37], the highest posterior probability value was found at K = 1 (Fig. 3). If the latent number of populations was not prespecified [38], the highest posterior probability was also found for a partition of one. A Bayesian inference of jointly the probability of assignment of individuals to populations and the number of populations [39,40] estimated that the number of populations with the highest posterior probability was also for K = 1 (P = 0.93).

Effective population size and migration rates
Migration rates and effective population size of bigeye tuna populations were estimated based on geographic Genetic isolation by distance in bigeye tuna inferred from multilocus estimates of F ST and R ST versus geographical distance  regions (ocean basins) of origin using both ML and BI methods [20]. The ML inferences did not return reliable confidence intervals using profile likelihoods even for very long runs (not shown), but none of the runs contradicted the Bayesian analysis. The results from the Bayesian inference were more stable and reliable. Estimates of the mutation scaled population size parameter were higher in the Atlantic Ocean (Θ = 6.07) than in the Indian Ocean (Θ = 2.33) and the Pacific Ocean (Θ = 1.67). These estimates were translated to an average effective population size (Ne) of 15175, 5825, and 4175 bigeye tuna individuals for the Atlantic Ocean, the Indian Ocean, and the Pacific Ocean, respectively (assuming a microsatellite mutation rate of 10 -4 per locus per generation; [41]). Marked differences in ratios of effective population size (mean Ne = 8391.67 ± 5932.19) to census population size (as derived from annual catch data, and averaging 403.000 tones for bigeye tuna; [42]) are commonplace in marine fishes [43,44], and can be explained either by historical events such as e.g. past population bottlenecks or by life history traits such as e.g. strong bias in reproductive success or size-dependent fecundity [44][45][46]. In order to discern among both competing hypotheses, we calculated values of M (a statistic that estimates population bottleneck history using the ratio of the number of alleles to the range in allele size; [47]), which were above the theoretical threshold of 0.68, indicating that none of the Atlantic, Indian and Pacific Ocean populations has experienced recent population reduction. Hence, it is likely that the low Ne observed in bigeye tuna may be due to life-history traits as previously reported for other marine fishes [44].
Migration rates seem to be symmetric between the Indian Ocean and the Pacific Ocean (Fig. 4A). However, immigration from and into the Atlantic Ocean seems to be highly asymmetrical with the Atlantic population providing twice as many immigrants into the Pacific Ocean and the Indian Ocean than those inferred between the Pacific Ocean and the Indian Ocean (Fig. 4A). Furthermore, immigration into the Atlantic Ocean showed the lowest rates. A likelihood ratio test revealed that these results were incompatible with symmetrical scaled immigration rates or with a common rate for all directions (P < 0.00001). Inferred numbers of immigrants per generation between ocean basins ranged between 11.5 and 27.9, Fig.  4B). tive population size (N e ) and time since divergence (t) showed that these loci provided enough statistical power [48] for the detection of genetic variability in Atlantic bigeye tuna. In fact, previously loci TA102, TA113, TA121, and TA161 were applied to analyze population genetic structuring of bigeye tuna in the Indian Ocean [49], and locus TA161 was able to detect weak genetic differentiation in the Pacific Ocean [28]. Significant deviations from HWE were found in few populations at five loci, which were different to those showing departures from HWE in Thunnus albacares [50]. The majority of the detected departures from HWE were due to significant heterozygote deficiency. Most of the significant F IS locus-population comparisons involved Guinea samples of different years versus other Atlantic populations, being deviations attributed to homozygote excess. According to the ANOVA test only the genetic factor (the locus) and neither the time nor the combination of both factors was responsible for the deviation from HWE. On the other hand, the only case of departure from HWE due to null alleles was associated to the Guinea location. All these results together highlight the singularity of the Guinea location, and could be explained in terms of admixture of genetically distinct cohorts (Whalund effect) within this site, owing to the fact that the Gulf of Guinea is known as the main spawning and nursery area for bigeye tuna in the Atlantic Ocean [22,51]. Similar patterns of HWE departures have been associated with spawning or nursery populations in the Atlantic bluefin tuna (Thunnus thynnus) [52], and lemon shark (Negaprion brevirostris) [53].

Genetic variation and population structure of Atlantic bigeye tuna
High values of genetic variability (mean number of alleles per locus and heterozygosity) characterized all studied locations. These values are similar to those previously reported for other scombroid fishes, such as the bluefin tuna [52,54], yellowfin tuna [50] and bigeye tuna [49,55], and consistent with the average for marine fishes [56].
Within the Atlantic Ocean, no significant genetic differentiation between temporal replicates or geographic locations was observed (all F ST and R ST pairwise comparisons rendered low and non-significant values except a F ST pairwise comparison between Canada2 and Canary). Additionally, the Mantel test failed to detect significant correlation between geographical and genetic distances, which discarded isolation by distance among populations within the Atlantic Ocean. Overall, these results support unrestricted gene flow among bigeye tuna populations within the Atlantic Ocean. Similarly, no genetic structuring was reported for bigeye tuna within the Pacific Ocean [28] nor within the Indian Ocean [49]. The lack of genetic differentiation within each ocean basin may not be that surprising given the relatively high dispersal capability of bigeye tunas [24]. Nevertheless, some statistical tests and analyzed loci have shown weak evidence for small genetic differences between most separated locations in the Pacific, suggesting that isolation by distance might be promoting subtle differentiation among populations, which could be at present difficult to document with the current loci and sample sizes [28] The hypothesis of a single panmictic unit for bigeye tuna worldwide is supported by microsatellite data based on the lack of significant differentiation in between-ocean R ST comparisons, and the comparatively low F ST and R ST values obtained. These values are by far lower than the mean F ST considered for marine fishes (0.062, [57][58][59]), but on the same order as those reported in other studies where range-wide panmixia was suggested [60][61][62][63]. Hierarchical AMOVA, Mantel test, as well as clustering and assignment tests rendered consistent results, and also supported the hypothesis of no genetic structure for tuna populations at an inter-oceanic scale.
Therefore, nuclear and mt evidence are clearly discordant with respect to inter-oceanic genetic differentiation. While the former, as mentioned above, fails to support any genetic structuring, the latter supports genetic differentiation between Atlantic Ocean and Indo-Pacific region bigeye tuna populations based on pairwise Φ ST comparisons [10]. Strong phylogeographic association of mtDNA (two clades), compared to microsatellites, has also been previously documented in e.g. blue marlin [64]. In the bigeye tuna, the main source of disagreement largely arises from the existence in the Atlantic Ocean of two highly divergent mt clades [10,26,27], which were not recovered by AMOVA tests based on nuclear data. The origin of the two mt clades (I and II) has been proposed to be associated with past genetic isolation of bigeye tuna Atlantic Ocean populations from those of the Indo-Pacific region during Pleistocene glacial maxima [10,[26][27][28]. Once temperatures rose during inter-glacial periods, gene flow could be re-established through the Cape of Good Hope. The present-day observed pattern would have been achieved when individuals belonging to mt Clade I entered the Atlantic Ocean from the Indo-Pacific region, whereas individuals belonging to mt Clade II apparently remained within the Atlantic Ocean. Microsatellite data indicate that despite the accumulated high sequence mt divergences, individuals from both clades were not reproductively isolated and thus, secondary contact in the Atlantic Ocean led to interbreeding, and present-day genetic homogenization at the nuclear level [33].

Migration patterns of bigeye tuna populations
Bayesian inference of migration rates indicates gene flow rates of the same order of magnitude inferred for other pelagic fishes such as e.g. sardines [65]. There is migration in both directions between the Pacific Ocean and the Indian Ocean, likely facilitated by both the Indonesian Ocean Current that flows from the Pacific Ocean into the Indian Ocean, as well as the seasonal reversal of this current [66]. However, according to our results, the main trans-oceanic migratory activity is occurring around the Cape of Good Hope. The region off South Africa is both topographically and hydrologically complex [29], which complicates understanding of bigeye tuna migration patterns. At both sides of the Cape of Good Hope, the South Atlantic Ocean [67] and the southwest Indian Ocean [68] gyres are connected by the strong Agulhas Current that flows from the Indian Ocean into the Atlantic Ocean [69].
As the Agulhas Current reaches the southern tip of the continental shelf of Africa, it turns almost completely back on itself, and flows eastward as the Agulhas retroflection or return current [69,70]. According to our results (Fig. 4A), immigration from the Atlantic Ocean into the Indo-Pacific region doubles that inferred between the Indian Ocean and the Pacific Ocean whereas scaled immigration rates into the Atlantic Ocean are the lowest, suggesting that immigration is less important for the maintenance of variability in the Atlantic Ocean than in other ocean basins.
Our results indicate that the actual number of individuals moving per generation (Fig. 4B) is several orders of magnitude lower than bigeye population sizes in each ocean (millions of individuals). Thus, the possibility of detecting changes in mt clade proportions in each ocean through catch statistics is extremely small. In fact, the proportions of individuals of mt Clade I (about 25%) and mt Clade II (about 75%) in the studied populations of Atlantic bigeye tuna seemed to remain constant over the years (Fig. 1, [10]) and virtually no individuals of mt Clade II have been reported from the Indian Ocean [49]. Moreover, significant allele-frequency differences at one intronic locus and one microsatellite locus were found between Atlantic and Indo-Pacific samples, supporting simple mixture but little movement of bigeye tuna individuals around the Cape of Good Hope [33]. Thus, it seems that the main invasion of individuals from the Indo-Pacific region carrying mt Clade I haplotypes into the Atlantic Ocean may have occurred during an earlier warmer interglacial period rather than after the last glacial maximum.
Alternatively, lack of nuclear differentiation coupled with strong mt divergences could also indicate instances of male-biased dispersal (females would exhibit phylopatric behavior and thus would be considered sedentary with regards to the ocean basin) [33]. Although, the existence of independent spawning areas in each ocean basin [26] may support the existence of homing behavior in bigeye tuna, and there are reports of regional fidelity in bigeye tuna [24], the existence of homing behavior remains an open question, and tagging data showing male trans-oceanic displacements is required to understand current gene flow among bigeye tuna populations.
Finally, the inferred asymmetrical immigration pattern out and into the Atlantic Ocean could simply reflect the presence in this ocean basin of Clade II, which is not shared with the Indo-Pacific region. The number of immigrants into the Atlantic Ocean could be less extreme because of the presence of the two clades in the Atlantic Ocean, which would produce larger variability, and thus, larger effective population sizes. However, it is important to note here that the inferred migration rates may also reflect a sampling bias since the data set is allele-rich and the present study only includes samples from one Indian Ocean and one Pacific Ocean location. Therefore, these two ocean basins may be underrepresented, and the estimates of the number of immigrants per generation (N m ) have a large variance associated (Fig. 4B).

F ST as a measure of genetic structure and gene flow
Most population genetic analyses performed in this study suggest lack of population structure, and relatively high levels of gene flow among bigeye tuna populations at a global scale in agreement with other population genetic analyses. However, nine out of 36 F ST pairwise comparisons detected significant population structure, mostly involving comparisons between populations of the Atlantic Ocean and those of the Indo-Pacific region ( Table 2). In order to understand this discrepancy, F ST analyses were performed using simulated data of three populations (Atlantic Ocean, Indian Ocean, and Pacific Ocean) (Figure 5). Estimates of F ST values from data with large number of immigrants (like the one estimated from bigeye tuna microsatellite data) showed large 95% confidence intervals that include F ST values up to about 0.2. The presumed differentiation was most likely an artifact of the large number of alleles present in the sample. Although F ST statistics [71] are widely used for describing population genetic structure, they present some limitations, including the implicit assumption of uniform effective population sizes, and symmetric migration rates. Violation of these assumptions is particularly worrisome when using highly polymorphic molecular markers with high mutation rates (e.g. microsatellites) to analyze weakly structured populations [72] with large effective population sizes [56,73,74]. F ST estimations assume that there is no error in the sample frequencies, but with many alleles per locus the precision of sample frequencies is questionable. In such cases, the use of both R ST , which is independent of the mutation rate [75,76], and coalescent methods, which use sample frequencies and not the population frequencies render more accurate and reliable estimates of genetic differentiation than F ST .

Conclusion
Population genetic analyses based on microsatellite data of bigeye tuna populations showed unrestricted gene flow within the Atlantic Ocean, and supported the existence of a spawning and nursery region in the Gulf of Guinea. At a global scale, R ST statistics, Bayesian cluster analyses, and coalescence-based migration rate inferences based on nuclear markers supported lack of trans-oceanic genetic structure, which contrasts with previously detected significant mt divergence of bigeye tuna between Atlantic Ocean and Indo-Pacific region due to the existence of two clades, one restricted to the Atlantic Ocean. Our results support interbreeding between individuals belonging to different clades in the Atlantic Ocean. Given the current distribution of bigeye tuna individuals belonging to the two clades, and the inferred asymmetrical nuclear migration rates between the Atlantic Ocean and the Indo-Pacific region, it is likely that little number of individuals (compare to the actual population sizes in each ocean) may currently move between the Atlantic Ocean and the Indo-Pacific region, maintaining the proportions of the mt clades in both oceans. Therefore, the current presence of individuals belonging to mt Clade I would have originated in an earlier interglacial period.  (Fig. 1). Three Atlantic Ocean sites (Guinea2, N = 50; Canary2, N = 32 and Canada 2, N = 23) were re-sampled a second year to investigate temporal stability of genetic structure in these regions (Fig. 1). Collections included exclusively juveniles except in Canada where adults were also captured (See Table 1

Descriptive statistics
For each locus and sampling site (including second-year samples), observed and expected heterozygosities [77], number of alleles (N A ), and number of alleles standardized to those of the population with the smallest sample size (N S ) [78] were calculated using both GENETIX 4.02 [79] and FSTAT 2.9.3 [80]. Similarly, F IS statistic estimations that detect deviations from HWE, and the linkage disequilibrium test were performed for each locus and sampling site using the program GENEPOP version3.3 [81]. Significance of both analyses was tested with a Markov chain Monte Carlo (MCMC) that was run for 1000 batches of 2000 iterations each, with the first 500 iterations discarded before sampling [82]. P values from multiple comparisons were Bonferroni corrected [83]. MICRO-CHECKER v2.23 [36] was used to explore the  [48] was used to estimate whether the dataset used for genetic population analysis provided enough statistical power for detecting significant genetic differentiation. The nine population sampled were used for testing allele frequency homogeneity at nine locus separately and combining information for the multiple loci using Fisher's exact and traditional chi-squared tests. Simulations were run using various combinations of N e and t (where N e is the effective population size and t is the time since divergence, respectively), leading to F ST of 0.001 and 0.0025, which reflects the magnitude of F ST and N e values estimated from our empirical data. We also estimated the α error (type I) by performing a simulation of no divergence among samples (i.e. setting t = 0 that leads a value of F ST = 0). Results indicated that the probability for detecting population structure were high and statistically significant, corresponding to 76% and 72% for the N e /t combinations of 9,000/18 and 10,000/20 and 100% for both the N e /t combinations of 9,000/45 and 10,000/45, respectively. When F ST was set to zero, the proportion of false significances (α) was 1%, which is lower to the intended value of 5%. Hierarchical genetic structuring of microsatellite data (only first year sampling) was analyzed by assessing the relative contribution among groups, within groups and within population components for partitions of total molecular variance (AMOVA) [89] using ARLEQUIN v3.0 [86] and 20,000 permutations. We specifically tested the null hypothesis of panmixia, as well as of alternative structuring by mt clade (Clade I and Clade II) and by geographical region (Atlantic Ocean, Indian Ocean and Pacific Ocean).

Analyses to establish the number of populations
In order to complement and contrast the results obtained with the classical standard F-statistics [90] we also inferred population structure using a Bayesian approach. The number of populations (K) with the highest posterior probability given the data (only first year sampling) was estimated using STRUCTURE [91], BAPS v3.1 [38] and STRUCTURAMA [92]. In STRUCTURE, we selected the admixture model and the option of correlated allele frequencies between populations (also called the F-model) because it is considered the superior model for detecting structure among closely related populations [37]. MCMC consisted on 1 × 10 5 burn-in iterations followed by 1,5 × 10 6 iterations. We explored K in the range of one and nine and we performed 20 runs for each K value. In BAPS [38], the Bayesian inference does not need the number of populations is pre-specified. We used the option "cluster groups of individuals" to run the program with default conditions. The program STRUCTURAMA [92] infers population genetic structure from genetic data by allowing the number of populations to be a random variable that follows a Dirichlet process prior [39,40]. We run 1 × 10 6 MCMC cycles, discarding the first 1 × 10 5 cycles as burn-in.

Estimation of migration rates and effective population sizes
The program MIGRATE v 2.4 [93] was used to infer the population size parameter Θ (i.e. 4N e μ, where N e is the effective population size and μ is the mutation rate per site) and the migration rate, M (M = m/μ, where m is the immigration rate per generation) among bigeye tuna populations. These analyses used the Brownian mutation model [93] and mutation was considered to be constant for all loci. We used the Bayesian inference [20] and the maximum likelihood [17,18], modes. F ST estimates and a UPGMA tree were used as starting parameters for the estimation of Θ and M. For each locus in the data set, the ML was run for ten short and three long chains with 10,000 and 100,000-recorded genealogies, respectively, after discarding the first 10,000 genealogies (burn-in) for each chain. One of every 20 reconstructed genealogies was sampled for both the short and long chains. The Bayesian run consisted of one long chain with 50 millions-recorded parameter and genealogy changes after discarding the first 10,000 genealogies as burn-in for each locus. For all the analyses we used an adaptive heating scheme with 4 concurrent chains; the analyses were run on a cluster computer using up to 80 compute nodes. We conducted the analyses over two sets of data, the first one with microsat-ellite data structured into two groups equivalent to mt Clades I and II; and the second set of data structured according to geographical regions (Atlantic Ocean, Indian Ocean, and Pacific Ocean).
In order to test whether past populations bottlenecks could be responsible of the observed difference between census and effective population sizes, we used the M_P_val program to calculate M, a statistic that estimates population bottleneck history using the ratio of the number of alleles to the range in allele size [47].

Simulated migration rates and F ST estimations
In order to better understand the limitations of F ST estimation, we explored the range of estimated F ST values for a range of migration rates. We simulated 20 times 100 datasets over the range of the number of immigrants per generation, N e m between 4 and 40 and report the 95% confidence intervals. The migration rate range covers the estimated migration rates of the real data using the coalescent inference method.