Genetic connectivity between Atlantic bluefin tuna larvae spawned in the Gulf of Mexico and in the Mediterranean Sea

Field collection of ABFT larvae

We conducted larval ABFT collections during the peak of the reproductive season in the two main spawning areas for ABFT, the GOM and the MED (Fig. 1). The GOM was divided into subregions separated at 90°W meridian and designated as eastern (eGOM) and western (wGOM) (Muller-Karger et al., 2015). A total of 76 stations were sampled in the GOM cruise “WS1405” (approved by the National Oceanic and Atmospheric Administration, Washington, D.C., USA with Permit Number TUNA-SRP-14-02) carried out from 28 April to 20 May 2014 aboard the R/V F.G. Walton Smith. The “BLUEFIN14” cruise (approved by the Spanish Institute of Oceanography, Madrid, Spain (Instituto Español de Oceanografía) “BLUEFIN TUNA” project) explored a total of 123 stations and took place from 13 June to 3 July 2014 in the western MED on board the R/V SOCIB, of the Balearic Islands Coastal Observing and Forecasting System. Both cruises used similar standardized methodologies for field collection of fish larvae as described previously (Laiz-Carrión et al., 2015). Larvae were either frozen in liquid nitrogen or preserved in 96% ethanol immediately upon retrieval.

Larvae were identified as ABFT following morphological, meristic and pigmentation characters. Standard length (SL) was measured to the nearest 0.01 mm using Image J 1.44a (National Institute of Health, Stapleton, NY, USA). Developmental stage was determined following Richards (2005). ABFT-positive stations were first determined for each survey (31 stations in the GOM and 63 in the MED), and a subset of stations was selected for genetic analysis from each spawning ground (15 stations in the GOM and 14 in the MED, see Fig. 1).

Genetic analysis of ABFT larvae

A total of 112 ABFT larvae spawned contemporarily in the GOM (62 larvae, 30 from eGOM and 32 from wGOM) and MED (50 larvae) were genotyped (Data S1). DNA extraction from ABFT larval tissue was performed with the NucleoSpin® Tissue XS (MACHEREY-NAGEL, Dueren, Germany), and PCR amplification of microsatellite loci alleles was carried out as described in Uriarte et al. (2019) with modifications (Table S1). Guts were removed prior to analysis to avoid interference from potential piscivory or cannibalism as observed by Uriarte et al. (2019) (except for 18 individuals for which low standard length excluded the possibility of piscivory).

We analyzed eight highly polymorphic microsatellite loci selected according to their moderate-to-high allele count and compatible amplification in two multiplex PCR reactions (Table S1). Allele calling after fragment analysis (STAB VIDA, Caparica, Portugal) was performed with GeneMapper® software v4.0 (Applied Biosystems, Foster City, CA, USA) by two independent readers. We used GeneScan™ 500 LIZ® (Applied Biosystems, Foster City, CA, USA) as size standard. The internal threshold for missing data per individual was 25% (six genotyped loci), which allows for an overall maximum missing data of 13% for locus Tth16-2. A fragment of ~450 bp of the mitochondrial control region (Viñas & Tudela, 2009), was sequenced in a subset of genotyped larvae (22 from GOM and 21 from MED). Partial mtDNA control region sequences were assigned to T. thynnus based on sequence similarity searches performed with the nucleotide Basic Local Alignment Search Tool (BLAST) (https://blast.ncbi.nlm.nih.gov/). Quality trimmed sequences left a fragment of 361 bp aligned detecting 53 haploid binary sites with SNPs. All sequences have been submitted to GenBank (accession numbers MT912036–MT912078).

Statistical analysis of ABFT genetic diversity

Examination of microsatellite genetic data was based on allele frequency estimation with GenAlEx software v6.5 (Peakall & Smouse, 2012) to assess diversity through observed (H_O) and expected (H_E) heterozygosity. Software Cervus v3.0.7 (Kalinowski, Taper & Marshall, 2007) was used to calculate the polymorphic information content (PIC) and to obtain null allele frequencies using the Maximum Likelihood estimator of Summers & Amos (1997). MICRO-CHECKER (Van Oosterhout et al., 2004) was used to further assess null alleles and genotyping errors. GENEPOP 4.7 (Raymond & Rousset, 1995; Rousset, 2008) was used to calculate the inbreeding coefficient F_IS (Weir & Cockerham, 1984) and compliance to equilibrium of Hardy–Weinberg (HW). We followed Weir (1996) correcting probability values due to low counts for certain genotypes with the Markov Chain Monte Carlo (MCMC) approximation (involving 10,000 dememorization steps, 1,000 batches and 10,000 iterations per batch). The presence of full siblings was screened through replicate sibship analyses with multilocus genotypes excluding individuals with missing data. Analyses were run with COLONY (Wang, 2004; Jones & Wang, 2010), assuming an inbreeding polygamous model, performing ten very long runs to assess maximum likelihood sibship. Genotype accumulation curves were explored in R free software (R Team, 2018) in RStudio (RStudio Team, 2019) with the “poppr” package (Kamvar, Tabima & Gṙunwald, 2014).

The software GenAlEx v6.5 (Peakall & Smouse, 2012) calculated fixation indexes to assess partitioning of genetic diversity at microsatellite loci and at the mtDNA control region within subpopulations relative to the total population. The Wright’s F statistic F_ST, G″_ST (Meirmans & Hedrick’s (2011) standardized G_ST further corrected for bias for small k populations), and Jost’s (2008) D estimate of differentiation were calculated, and the associated probability (p) was obtained based on a 999 data permutation test. For ABFT mtDNA haplotypes the F_ST analogous Phi_PT was obtained through Analysis of Molecular Variance (AMOVA) of genetic distances. Spatial genetic analysis of both nuclear and mitochondrial diversity was studied performing a Mantel test with GenAlEx software v6.5 (Peakall & Smouse, 2012). Mantel correlation coefficient Rxy for matrix and spatial autocorrelation within distance classes were performed calculating r by bootstrapping.

Two approaches, Bayesian algorithms and Discriminant Analysis of Principal Components (DAPC), proportionally assigned individuals to population clusters inferred from the microsatellite genetic data. In the Bayesian approach, STRUCTURE v2.3.4 software (Pritchard, Stephens & Donnelly, 2000) was used to analyze ABFT larvae genotyped at six microsatellite loci assuming an ancestry admixture correlated allele frequency model (Falush, Stephens & Pritchard, 2003), and considering prior sampling location information. Running parameters were set to a burn-in of 2.5 × 10⁴ followed by a MCMC simulation of 5 × 10⁴ runs simulating K 1 to 3 populations in 25 iterations. STRUCTURE results were explored with Structure Harvester (Earl & vonHoldt, 2012). Clumpak (Kopelman et al., 2015) was used to detect the consensus solutions for two K clusters that best fit the Bayesian algorithm used by STRUCTURE. In the DAPC approach we used the R package adegenet (Jombart, 2008) to study differences among clusters identified from ABFT larvae genotyped at six microsatellite loci, and also from polymorphisms in mtDNA.

Field collection and genetic characterization of ABFT larvae

During the 2014 spawning season, approximately half of the stations explored were positive for ABFT larvae (Fig. 1), 43% for GOM and 55% for MED. We genetically characterized a total of 112 ABFT larvae spawned contemporarily in the GOM (62 larvae) and MED (50 larvae), excluding overlap in generations. In the GOM, two sample subsets were separated east (30 larvae) and west (32 larvae) of 90°W according to oceanographic features known to influence ABFT catches (Teo & Block, 2010; Muller-Karger et al., 2015), with the eGOM dominated by the Loop Current shedding large anti-cyclonic eddies, which generate mesoscale cyclonic and anti-cyclonic eddies that are the key features of the wGOM, preferred by breeding ABFT, and influencing larval fish distribution (Lindo-Atichati et al., 2012). Mean SL and standard deviation in mm were 4.94 ± 1.05 (MED) and 4.55 ± 0.64 (GOM) (wGOM 4.49 ± 0.49 and eGOM 4.62 ± 0.77).

The mean number of larvae genotyped at eight microsatellite loci are shown in Table S2, for each of the two spawning areas (GOM and MED), or for the three separate geographical regions (wGOM, eGOM and MED again). The number of individuals was sufficient to adequately quantify allele frequencies (Hale, Burg & Steeves, 2012) and further analyze diversity through heterozygosity. The overall number of alleles ranged from six for locus Tth157 to 25 for locus Tth208 (see Table S2), with the highest polymorphic information content (PIC) score (PIC 0.9) obtained for loci Tth208 and Tth1-31. Scoring errors due to stuttering or large allele dropout were excluded with MICRO-CHECKER (Van Oosterhout et al., 2004), whereas null alleles may be present at loci Ttho1 and Tth16-2 according to homozygote excess and null allele frequency (Table S2). Deviations from expected proportions in large populations in HW equilibrium imply non-random mating, selection for certain genotypes, mutations, or small population sizes. Under the assumptions of HW equilibrium, the inbreeding coefficient (F_IS) varies between −1 and +1, measuring the difference between expected heterozygosity (H_E) and observed heterozygosity (H_O) according to allele frequencies. Significant deviations from HW proportions and homozygote excess (positive F_IS) were obtained for loci Tth16-2 and Ttho1 (Table 1), as expected from their null allele frequencies (Table S2). Locus Tth16-2 is the microsatellite that deviated most from HW proportions in all populations and areas (Table 1). Close to zero F_IS and low null allele frequency did not explain deviation from proportions expected in HW equilibrium for loci Tth157 and Ttho7 in the GOM dataset (Table 1).

Characterization of genetic diversity in larval ABFT may be biased if closely related individuals are collected during plankton sampling. The full-likelihood method implemented in COLONY (Wang, 2004; Jones & Wang, 2010) can perform sibship analysis to detect full siblings, and is better than other software for parentage analysis in natural populations (Harrison et al., 2013). In both spawning areas, sibship was inferred from ~50 replicate multilocus genotype datasets and the resulting parameters averaged (Table S3). The hypothesized numbers of families were high and very similar between spawning areas, 45 ± 4 in the GOM and 47 ± 1 in the MED, supporting unbiased characterization of genetic diversity.

Spatial structuring of ABFT genetic diversity

Structuring of genetic diversity between ABFT spawning areas was assessed from multilocus genotypes considering all genotyped microsatellite loci or excluding those with null alleles (Tth16-2 and Ttho1). Partitioning of the total expected genetic diversity between subpopulations or groups under the assumption of equilibrium was quantified through several statistics ranging from 0 to +1. We calculated the general or most frequently used fixation index F_ST, originally established by Wright for two allele systems, and we also obtained Meirmans & Hedrick’s (2011) standardized G″_ST and Jost’s (2008) D estimate of differentiation, which are standardized relative to marker heterozygosity and are thus more appropriate estimators for highly polymorphic microsatellite loci. When we only consider genotypes at six microsatellite loci, excluding the two loci with null alleles, genetic differentiation is captured by the value of G″_ST and D. Close to zero values with significant associated probabilities for deviation of homogeneity or equilibrium were obtained for GOM versus MED comparisons (Table 2). Genetic diversity partitioning between MED and the two subregions wGOM and eGOM, resulted in higher estimator values and levels of significance for MED versus the wGOM, in support of greater connectivity between ABFT breeding in the MED and in the eGOM during 2014. Pairwise comparisons of wGOM and eGOM genotypes supported genetic homogeneity, with non-significant probabilities associated with close to zero statistics F_ST, G″_ST and D, excluding structuring of genetic diversity within the GOM spawning area.

Unique individuals were completely discriminated with four loci according to a genotype accumulation curve (Fig. S1). We obtained diversity estimators from ABFT larvae genotypes at the four most informative loci (PIC ≥ 0.7) that were neutral or in accordance with HW proportions, thus informing of homogeneity of diversity between larval groups. Table 2 shows that four neutral microsatellite loci: Tth208, Tth1-31, Ttho4 and Tth34; detected structuring between larvae spawned in breeding grounds of each ABFT stock, underlining the power of resolution of microsatellite loci as markers for this tuna species. Even multilocus genotypes for only the two loci with the highest PIC were able to significantly detect structuring between GOM and MED larvae (data not shown).

Pairwise comparison of MED and GOM larvae haplotypes at 53 mtDNA control region variable sites (SNPs) through AMOVA resulted in a PhiPT value of 0.029 (with a significant associated data randomization probability of 0.019 with 999 permutations). Phylogenetic analysis of mtDNA sequences to infer evolutionary relationships from a common ancestor failed to separate larvae spawned in the MED and in the GOM (data not shown), as previously reported (Alvarado Bremer et al., 2005; Carlsson et al., 2007).

ABFT genetic diversity partitioning, assessed herein exclusively through larvae resulted in fixation indexes within the ranges reported thus far, analyzing polymorphism in ABFT mtDNA sequences, microsatellite loci or SNPs, in samples of certain (larvae) or almost certain (juveniles) origin (Carlsson et al., 2007; Puncher et al., 2018; Rodríguez-Ezpeleta et al., 2019) (Table 3). We collected ABFT larvae in a precise spatial and temporal frame, and our results are consistent with other genetically characterized juvenile ABFT and larvae collected in different years and areas compiled in Table 3. Besides F_ST, we report other genetic diversity estimators, G″_ST (Meirmans & Hedrick’s (2011) standardized G_ST further corrected for bias for small k populations) and D (Jost’s et al. (2018) estimate of differentiation), as they show a higher variance, and are complementary measures of fixation and allelic differentiation (Bird et al., 2011), arguing against the generalized exclusive use of F_ST.

We assessed whether genotypes were more similar when comparing larvae collected in the same spawning area than when comparing larvae from different spawning areas. Spatial genetic analyses were performed by analyzing the autocorrelation between geographic and genetic data through a Mantel test. Close to zero autocorrelation statistics resulted from analyzing genotypes for six loci (R 0.127, p 0.001) and mtDNA control region binary haplotypes (R 0.051, p 0.022). Genetic structuring within different distance classes (within the same or different spawning areas) was explored through multivariate spatial autocorrelation analysis. Autocorrelation for microsatellite genotypes resulted in a positive r value of 0.024 (with confidence limits of 0.005 and −0.004) for distance classes within the same spawning area and a negative r value of −0.020 (with U of 0.003 and L of −0.004) for distance classes belonging to different spawning areas, in agreement with a weak but significant spatial structuring of genetic diversity.

We investigated connectivity and structuring of genetic diversity through assignment of individual ABFT larvae to two K clusters through a Bayesian analysis of multilocus genotypes in an admixture model considering prior location information (Figs. 2A–2C). In two out of three consensus solutions, representing 32% of all simulations, a gradient from wGOM to MED was observed in the proportion of assignment of each larva to two clusters, with eGOM appearing as an intermediate or mixing area (Figs. 2B and 2C). Figure S2 shows the statistic ∆K (deltaK) proposed by Evanno, Regnaut & Goudet (2005) to detect the real number of clusters based on the rate of change of probabilities between successive K values, and that for our data is necessarily only obtained for K = 2. A different clustering approach is multivariate DAPC analysis that is applied to discriminate clusters through discriminant analysis performed on data previously transformed through principal components analysis (Jombart, Devillard & Balloux, 2010). Assignment of individual ABFT larvae to clusters performing DAPC with one discriminant function (DF) clearly separated larvae spawned in the GOM and in the MED (Fig. 2D). DAPC analysis of variable sites found in the mtDNA control region sequences also separated larvae spawned at Eastern and Western stock reproduction areas (Fig. S3). In Fig. 2D, the similar proportion of assigned individuals to the wGOM and the eGOM supports absence of structuring within the GOM spawning area. Retaining 0.925 variance in two DFs separated ABFT spawned in wGOM, eGOM and MED in three overlapping clusters (Fig. 2E). DAPC is thus in accordance with the wGOM to MED gradient observed through Bayesian clustering, and supports apparent connectivity between the ABFT Western and Eastern stocks. Separate density plots (Figs. 2F and 2G) for DF1 (eigen value 64.72) and DF2 (eigen value 29.34) also illustrate homogeneity between larvae spawned in the GOM.