Natalia Diaz-Arce

Pierre-Alexandre Gagnaire

David E Richardson

John Walter

Sophie Arnaud-Haond

Jean-Marc Fromentin

Deirdre Brophy

Molly Lutcavage

Piero Addis

Francisco Alemany

Robert Allman

Simeon Deguara

Igaratza Fraile

Nicolas Goñi

Alex R Hanke

F Saadet Karakulak

Ashley Pacicco

Joseph M Quattro

Jay R Rooker

Haritz Arrizabalaga

Naiara Rodriguez-Ezpeleta

Natalia Díaz-Arce

John F Walter Iii

Saadet Karakulak

Naiara Rodríguez-Ezpeleta

Unidirectional trans-Atlantic gene flow and a mixed spawning area shape the genetic connectivity of Atlantic bluefin tuna

Keywords: [4-6]: Large migratory fish, genetic connectivity, Atlantic bluefin tuna, introgression, single nucleotide polymorphisms

The commercially important Atlantic bluefin tuna (Thunnus thynnus), a large migratory fish, has experienced notable recovery aided by accurate resource assessment and effective fisheries management efforts. Traditionally, this species has been perceived as consisting of eastern and western populations, spawning respectively in the Mediterranean Sea and the Gulf of Mexico, with mixing occurring throughout the Atlantic. However, recent studies have emerged challenging this assumption by revealing weak genetic differentiation and identifying a previously unknown spawning ground in the Slope Sea used by Atlantic bluefin tuna of uncertain origin. To further understand the current and past population structure and connectivity of Atlantic bluefin tuna, we have assembled a unique dataset including thousands of genome-wide Single Nucleotide Polymorphisms (SNPs) from five hundred larvae, young of the year and spawning adult samples covering the three spawning grounds and including individuals of other Thunnus species. Our analyses support two weakly differentiated but demographically connected ancestral populations that interbreed in the Slope Sea. Moreover, we also identified signatures of introgression from albacore into the Atlantic bluefin tuna genome, exhibiting varied frequencies across spawning areas, indicating strong gene flow from the Mediterranean Sea towards the Slope Sea. We hypothesize that the observed genetic differentiation may be attributed to increased gene flow caused by a recent intensification of westward migration by the eastern population, which could have implications for the genetic diversity and conservation of western populations. Future conservation efforts should consider these findings to address

Introduction

Conservation of fisheries resources relies on the assessment and management of self-sustaining units called stocks, whose delimitation often oversimplifies species population dynamics [START_REF] Begg | Stock identification and its role in stock assessment and fisheries management: an overview[END_REF][START_REF] Stephenson | Stock complexity in fisheries management: a perspective of emerging issues related to population sub-units[END_REF][START_REF] Reiss | Genetic population structure of marine fish: mismatch between biological and fisheries management units[END_REF]). Yet, failing to account for stock complexity can induce overfishing and ultimately result in fisheries collapse [START_REF] Hutchinson | The dangers of ignoring stock complexity in fishery management: the case of the North Sea cod[END_REF], highlighting the importance of integrating knowledge on spatial structure and connectivity into management plan development processes [START_REF] Kerr | Lessons learned from practical approaches to reconcile mismatches between biological population structure and stock units of marine fish[END_REF]. In this context, the potential of genomic approaches is increasingly being harnessed to tackle a diverse range of fisheries management related questions, such as assessment of population structure, connectivity and adaptation to local environments [START_REF] Ovenden | Ocean's eleven: a critical evaluation of the role of population, evolutionary and molecular genetics in the management of wild fisheries[END_REF][START_REF] Bernatchez | Harnessing the Power of Genomics to Secure the Future of Seafood[END_REF], even when genetic differentiation is low, as observed in highly migratory fish like striped marlin [START_REF] Mamoozadeh | Genome-wide SNPs resolve spatiotemporal patterns of connectivity within striped marlin (Kajikia audax), a broadly distributed and highly migratory pelagic species[END_REF], blue shark [START_REF] Nikolic | Stepping up to genome scan allows stock differentiation in the worldwide distributed blue shark Prionace glauca[END_REF]) and yellowfin tuna (Barth et al. 2017b). In particular, the preservation of fish genetic diversity and conservation of locally adapted populations has gained importance in the face of rapid environmental changes and increasing fishing pressure [START_REF] Bonanomi | Archived DNA reveals fisheries and climate induced collapse of a major fishery[END_REF]. Species resilience may depend on adaptive evolution capacities [START_REF] Hoffmann | Climate change and evolutionary adaptation[END_REF]Sgrò 2011, Bernatchez 2016), making essential the inclusion of adaptive variation in genomic studies focusing on managed species [START_REF] Fraser | Adaptive evolutionary conservation: towards a unified concept for defining conservation units[END_REF][START_REF] Valenzuela-Quiñonez | How fisheries management can benefit from genomics?[END_REF][START_REF] Xuereb | Incorporating putatively neutral and adaptive genomic data into marine conservation planning[END_REF].

The Atlantic bluefin tuna (ABFT, Thunnus thynnus ) is a large and emblematic highly migratory species that inhabits waters of the North Atlantic Ocean and adjacent seas [START_REF] Fromentin | Lessons from the past: investigating historical data from bluefin tuna fisheries[END_REF]Powers 2005, Collette et al. 2011). ABFT has been heavily exploited for millennia and the emergence of the sushi-sashimi market in the 1980s turned it into one of the most valuable tuna species in the international fish trade (Fromentin et al. 2014a). This high value, coupled with poor governance, led to three decades of high fishing pressure and ultimately to overexploitation. By 2011, ABFT was considered endangered by the IUCN [START_REF] Collette | High Value and Long Life -Double Jeopardy for Tunas and Billfishes[END_REF]. Following the implementation of a strict management plan in the late 1990's, signs of population rebuilding have been documented [START_REF] Iccat | SRCS Advice to the Comission[END_REF], but uncertainties around ABFT biology suggest that an overly simplistic management paradigm could compromise the long-term conservation of the species (Fromentin et al. 2014a[START_REF] Brophy | Combining genetic markers with stable isotopes in otoliths reveals complexity in the stock structure of Atlantic bluefin tuna (Thunnus thynnus)[END_REF]. ABFT has been managed as two separate units since 1981: the Western and Eastern stocks, which are separated by the 45°W meridian and are assumed to originate from the two spawning areas located in the Gulf of Mexico and the Mediterranean Sea respectively [START_REF] Iccat | Report of the Standing Committee on Resarch and Statistics[END_REF]. Several studies on the population structure and stock dynamics support two reproductively isolated spawning components (Gulf of Mexico and the Mediterranean Sea): electronic tagging studies [START_REF] Block | Electronic tagging and population structure of Atlantic bluefin tuna[END_REF] have found no individual visiting both spawning areas, and otolith chemical signatures [START_REF] Rooker | Crossing the line: migratory and homing behaviors of Atlantic bluefin tuna[END_REF] and genetic data [START_REF] Rodríguez-Ezpeleta | Determining natal origin for improved management of Atlantic bluefin tuna[END_REF]) support natal homing. Nevertheless, numerous studies also detected evidence of regular trans-Atlantic movements across the 45°W meridian boundary line and of mixed foraging grounds along the North Atlantic [START_REF] Block | Electronic tagging and population structure of Atlantic bluefin tuna[END_REF][START_REF] Rooker | Crossing the line: migratory and homing behaviors of Atlantic bluefin tuna[END_REF][START_REF] Arregui | Movements and geographic distribution of juvenile bluefin tuna in the Northeast Atlantic, described through internal and satellite archival tags[END_REF][START_REF] Rodríguez-Ezpeleta | Determining natal origin for improved management of Atlantic bluefin tuna[END_REF]. In response to these findings, the International Commission for the Conservation of Atlantic Tunas (ICCAT) recently adopted a management procedure for ABFT that accounts for mixing between the two stocks (ICCAT 2023). Given recent advancements in stock of origin assignment and increased samples from the mixing areas, it is important to determine if the modeled dynamics are consistent with the new data. Specifically, when applying individual origin assignment based on subsets of informative genetic markers of ABFT captured in the North Atlantic Ocean [START_REF] Rodríguez-Ezpeleta | Determining natal origin for improved management of Atlantic bluefin tuna[END_REF][START_REF] Puncher | Individual assignment of Atlantic bluefin tuna in the northwestern Atlantic Ocean using single nucleotide polymorphisms reveals an increasing proportion of migrants from the eastern Atlantic Ocean[END_REF], it was observed that 10-25% of individuals could not be clearly assigned to either spawning ground. Moreover, a combined analysis of genetic and otolith microchemistry data resulted in contrasting or unresolved origin assignments [START_REF] Brophy | Combining genetic markers with stable isotopes in otoliths reveals complexity in the stock structure of Atlantic bluefin tuna (Thunnus thynnus)[END_REF].

Amidst uncertainty surrounding ABFT stock dynamics, the recent discovery of ABFT larvae in the Slope Sea (Richardson et al. 2016a) adds another layer of complexity to our knowledge of the reproductive ecology of the species. Subsequent oceanographic studies [START_REF] Rypina | Investigating the suitability of the Slope Sea for Atlantic bluefin tuna spawning using a high-resolution ocean circulation model[END_REF]) and larval collections [START_REF] Hernández | Support for the Slope Sea as a major spawning ground for Atlantic bluefin tuna: evidence from larval abundance, growth rates, and particle-tracking simulations[END_REF] provide additional evidence of spawning activity in this area. Tagging information further revealed that mature size fish occurred in the Slope Sea in spring and summer coinciding with the spawning season estimated for the found larvae [START_REF] Galuardi | Complex migration routes of Atlantic bluefin tuna (Thunnus thynnus) question current population structure paradigm[END_REF], supporting the hypothesis of spawning strategy in the western Atlantic. The implications of Slope Sea spawning generated debate and controversy (Richardson et al. 2016b[START_REF] Safina | Data do not support new claims about bluefin tuna spawning or abundance[END_REF], Walter et al. 2016), with one of the key unknowns being the connectivity between the Slope Sea and the other two spawning grounds. In addition, some studies found evidence of migratory changes in ABFT, including the return to [START_REF] Horton | Tracking Atlantic bluefin tuna from foraging grounds off the west coast of Ireland[END_REF][START_REF] Nøttestad | The comeback of Atlantic bluefin tuna (Thunnus thynnus) to Norwegian waters[END_REF][START_REF] Aarestrup | First tagging data on large Atlantic bluefin tuna returning to Nordic waters suggest repeated behaviour and skipped spawning[END_REF]) and even expansion [START_REF] Jansen | Atlantic bluefin tuna (Thunnus thynnus) in Greenland -mixed-stock origin, diet, hydrographic conditions, and repeated catches in this new fringe area[END_REF] of its geographic range. These changes coincided with a strong recovery of the Mediterranean Sea spawning biomass during the last two decades and the increased presence of eastern origin fish in the western Atlantic [START_REF] Aalto | Stock-of-origin catch estimation of Atlantic bluefin tuna (Thunnus thynnus) based on observed spatial distributions[END_REF]).

To disentangle the population structure and connectivity of ABFT, we genotyped and analyzed thousands of genome-wide single nucleotide polymorphism (SNPs) from a total of five hundred ABFT larvae, young of the year (YoY) and adults from the two well-known spawning grounds (Gulf of Mexico and Mediterranean Sea) as well as the recently discovered Slope Sea spawning ground. We studied individual genomic diversity, tested for admixture between spawning grounds and inferred the demographic history of ABFT for the first time. We screened for adaptive genomic variation, incorporating samples from other Thunnusspecies to evaluate the impact of gene flow between species as an additional contribution to adaptive genomic diversity. Finally, we integrated information obtained from neutral, potentially adaptive and introgressed genetic markers to reconstruct the connectivity patterns of the ABFT across its entire distribution.

Methods

A summarized schematic view of samples and methods is shown in Figure S1.

Sampling, DNA extraction and additional data collection

Larvae, YoY and adult samples of ABFT from the Mediterranean Sea (n= 260), the Gulf of Mexico (n=210) and the Slope Sea (n=49) were obtained from scientific surveys and commercial fisheries from these three spawning grounds (Table S1; Figure 1A). From each adult and YoY, a ˜1cm 3 piece of muscle or fin tissue sample was excised and immediately stored in RNA-later or 96% molecular grade ethanol at -20ºC until DNA extraction. Larvae were collected with a 60 cm diameter bongo net or a 2 x 1-meter frame net and immediately preserved in 96% molecular grade ethanol. Genomic DNA was extracted from about 20 mg of tissue or from whole or partial larvae (eyeballs or tails) using the Wizard® Genomic DNA Purification kit (Promega, WI, USA), following manufacturer's instructions for "Isolating Genomic DNA from Tissue Culture Cells and Animal Tissue". Extracted DNA was suspended in Milli-Q water and concentration was determined with the Quant-iT dsDNA HS assay kit using a Qubit® 2.0 Fluorometer (Life Technologies). DNA integrity was assessed by electrophoresis, migrating about 100 ng of GelRed-stained DNA on a 1.0% agarose gel. For selected specimens, spawning capability was assessed by histologic inspection following the criteria described in [START_REF] Brown-Peterson | A Standardized Terminology for Describing Reproductive Development in Fishes[END_REF]. Additionally, for females, ovulation within the past two days was determined from identification of postovulatory follicle complexes, which are assumed to degrade within 24-48 h (McPherson 1991[START_REF] Schaefer | Spawning time, frequency, and batch fecundity of yellowfin tuna, Thunnus albacares, near Clipperton Atoll in the eastern Pacific ocien[END_REF][START_REF] Aranda | GnRHa-induced spawning in cage-reared Atlantic bluefin tuna: An evaluation using stereological quantification of ovarian post-ovulatory follicles[END_REF]. For selected adult specimens, sagittal otoliths were prepared for analysis of stable isotope signatures of the otolith core (yearling period) according to the protocol described in [START_REF] Rooker | Evidence of trans-Atlantic movement and natal homing of bluefin tuna from stable isotopes in otoliths[END_REF] and analyzed for δ 13 C and δ 18 O on an automated carbonate preparation device (KIEL-III Thermo Fisher Scientific, Inc., Waltham, Mass.) coupled to a gas-ratio mass spectrometer (Finnigan MAT 252 Thermo Fisher Scientific, Inc.) at the University of Arizona. Stable isotopes of carbon and oxygen (δ 13 C and δ 18 O) are reported relative to the PeeDee Belemnite (PDB) scale after comparison to an in-house laboratory standard calibrated to PDB.

Cytochrome oxidase subunit I gene fragment amplification and sequencing

A fragment of the mitochondrial cytochrome oxidase subunit I (COI) gene was amplified for a representative subset of 86 individuals using the FishF1 (5'-244 TCAACCAACCACAAAGACATTGGCAC-3') and FishR1 (5'-TAGACTTCTGGGTGGCCAAAGAATCA-3') primers (Ward et al. 2005) in a total volume of 20 μl with 0.2 μl of Dream Taq Polymerase (Thermo Fisher Scientific), 2 μl of Dream Taq Buffer 10X (Thermo Fisher Scientific), 0.4 μl of each primer and 50 ng of total DNA using the following profile: an initial denaturation step at 95°C for 3 min, 35 cycles of 30 sec at 98°C, 30 sec at 54°C and 60 sec at 72°C, and a final extension of 72°C for 10 min. Products were visualized on 1.7% agarose gels, purified with GE Healthcare Illustra ExoProStar (ref. US77705) and Sanger sequenced. The newly generated 86 sequences were edited using SeqTrace 0.9.0, submitted to Genbank (Accession numbers MT037084-MT037149, MT037151-MT037170) and aligned with BioEdit (v7.2.5) together with other publicly available representative COI sequences of albacore (accession number KT074102) and ABFT (accession number DQ107585), including the alalungalike (accession number GQ414567) haplotypes (Table S2). Diagnostic positions between ABFT and albacore haplotypes were used to detect mitochondrial introgression from albacore to ABFT samples.

RAD-seq library preparation, sequencing and read filtering

Restriction-site-associated DNA libraries of 519 ABFT individuals were prepared following [START_REF] Etter | SNP discovery and genotyping for evolutionary genetics using RAD sequencing[END_REF] . Input DNA (ranging from 50 to 500 ng) was digested with the SbfI restriction enzyme and ligated to modified Illumina P1 adapters containing 5bp unique barcodes. Pooled DNA of 32 individuals was sheared using the Covaris(r) M220 Focused-ultrasonicator Instrument (Life Technologies) and size selected to 300-500 pb on agarose gel. After Illumina P2 adaptor ligation, each library was amplified using 14 PCR cycles. Each pool was paired-end sequenced (100 pb) on an Illumina HiSeq2000. Demultiplexing, quality filtering (removing reads whose average Phred score is lower than 20 and truncating them to 90 nucleotides to remove low-quality bases at the end) and PCR duplicate removal were preformed using theprocess radtags and clone filter modules ofStacks version 2.3e [START_REF] Rochette | Stacks 2: Analytical methods for paired-end sequencing improve RADseq-based population genomics[END_REF].

RAD-tag assembly and SNP calling

Five RAD-seq derived catalogs (Figure S1; Table S3) were generated. Three of them included ABFT individuals and were either de novo assembled ('nuclear de novo') or mapped to the Pacific bluefin tuna nuclear (PBFT, Thunnus orientalis ) [START_REF] Suda | Improvement of the Pacific bluefin tuna (Thunnus orientalis) reference genome and development of male-specific DNA markers[END_REF] or ABFT mitochondrial (accession number NC -014052) genomes ('nuclear mapped' and 'mito'). The two others were mapped to the PBFT genome and included ABFT individuals and 4 Southern bluefin tuna (Thunnus maccoyii ), 4 albacore (Thunnus alalunga ) and 5 PBFT individuals (Diaz-Arce et al. 2016) ('nuclear mapped + others') or only ABFT larvae and 4 albacore individuals ('nuclear mapped + ALB'). Both reference-mapped and de novo assembled catalogs were generated for testing possible bias introduced by the use of the reference genome from a closely related species, which is less fragmented than the one available for ABFT (Accession GCA 003231725). Three nuclear mapped catalogs were generated including different groups of individuals to maximize the number of informative markers included for each type of analysis. In order to avoid inclusion of kins in the resulting datasets, which could bias some population structure results, a genetic relatedness matrix using the GCTA toolbox (Yang et al. 2011) was generated using the genotypes obtained from the 'nuclear mapped' catalog (generated as described below), and only one individual (the one with the highest number of assembled RAD tags) of each resulting pair with relatedness higher than 0.1 was included in subsequent analyses. Referencebased assemblies were performed by mapping the quality-filtered reads to the corresponding reference genome using the BWA-MEM algorithm [START_REF] Li | Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM[END_REF], converting the resulting SAM files to sorted and indexed BAM files using SAMTOOLS [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF]) and filtering the mapped reads to include only primary alignments and correctly mate mapped reads. De novo assembly was performed using the ustacks, cstacks, sstacks and tsv2bam modules of Stacks version 2.3e with a minimum coverage depth of 3 reads per allele (this is, each of the two possible versions of one biallelic SNP variant), a maximum of 2 nucleotide mismatches between two alleles at a same locus and a maximum of 6 mismatches between loci (Rodriguez-Ezpeleta et al. 2019). For all mapped and de novo catalogs, SNPs were called using information from paired-end reads with the gstacks module of Stacks version 2.3e. For the 'mito' catalog only samples with no missing data for the three diagnostic positions used for detecting introgression were kept, and the heterozygous genotypes, considered to be related to sequencing or assembly errors, were removed. For the rest of the RAD catalogs, only samples with more than 25,000 RAD loci and SNPs contained in RAD-loci present in at least 75% of the ABFT ('nuclear mapped' and 'de novo') or in 75% of the individuals from each of the species included ('nuclear mapped + others' and 'nuclear mapped + ALB') were kept and exported into PLINK [START_REF] Purcell | PLINK: a tool set for whole-genome association and populationbased linkage analyses[END_REF]) using thepopulations module of Stacks version 2.3e. Using only SNPs derived from read 1, increasing threshold values for minimum genotyping rate for individuals and SNPs were applied to obtain a final genotype table with a minimum genotyping rate of 0.95 and 0.85 per SNP and individual, respectively (except for the 'nuclear mapped + ALB' catalog for which thresholds were 0.95 and 0.90, respectively). SNPs were filtered using different minor allele frequency (MAF) thresholds considering sample sizes of the different datasets to exclude from the analysis rare non-informative variants that are susceptible to being derived from sequencing or assembly errors. For the 'nuclear mapped' and 'nuclear de novo' catalogs, SNPs with a MAF < 0.05 were removed, for the 'nuclear mapped + others' catalog, SNPs with MAF < 0.05 in ABFT and MAF < 0.25 in each of the other species were removed, and for the 'nuclear mapped + ALB' catalog, SNPs with a minimum allele count of two in ABFT and those variable within albacore were removed. For all nuclear catalogs, SNPs failing Hardy Weinberg equilibrium test at a p-value threshold of 0.05 in Mediterranean Sea larvae or Gulf of Mexico larvae groups were removed. Resulting genotype tables including all SNPs or only the first SNP per tag were converted to genepop, structure, PLINK, BayeScan, immanc, VCF and treemix formats using populations or PGDSpider version 2.0.8.3 [START_REF] Lischer | PGDSpider: an automated data conversion tool for connecting population genetics and genomics programs[END_REF] . From the 'mito' catalog, genotypes for the three diagnostic positions used for detecting introgression were extracted using PLINK [START_REF] Purcell | PLINK: a tool set for whole-genome association and populationbased linkage analyses[END_REF]).

Genetic diversity and population structure estimates

The following analyses were performed on the 'nuclear mapped' and 'nuclear de novo' datasets including only the first SNP per tag. Genome-wide average and per-SNP pairwise F ST values were calculated using GENEPOP [START_REF] Raymond | GENEPOP (version 1.2): population genetics software for exact tests and ecumenicism[END_REF] both including all individuals or only larvae. Significance (p < 0.05) of F ST values was estimated by performing 10,000 permutations. Principal Component Analysis (PCA) was then performed using the adegenet R package [START_REF] Jombart | adegenet 1.3-1: new tools for the analysis of genome-wide SNP data[END_REF] to illustrate the main axes of genetic variation among individuals. The number and nature of distinct genetic clusters was investigated using the model based clustering method implemented in ADMIXTURE [START_REF] Alexander | Fast model-based estimation of ancestry in unrelated individuals[END_REF]) assuming from 2 to 5 ancestral populations (K) and setting 5000 bootstrap runs. A first ADMIXTURE run was launched for each value of K to check the number of steps necessary to reach the default 0.001 likelihood value during the first run. This information was used to set the "-c" parameter (steps to be fulfilled in each bootstrapped run) that would assure convergence for each analysis (from 20 to 100 steps) for the bootstrapped runs. The value of K (ranging from 2 to 10) with lowest associated error value was identified using ADMIXTURE's cross-validation procedure. The convertf function from ADMIXTOOLS software [START_REF] Patterson | Ancient admixture in human history[END_REF] was used to convert from PLINK to eigenstrat format and the qp3Pop function was used to calculate F3 statistic and Z-score associated values [START_REF] Patterson | Ancient admixture in human history[END_REF], testing for all possible admixture scenarios grouping separately samples from different locations and age classes (Table S5) on the 'nuclear mapped' catalog dataset.

Demographic History

We used the unfolded three-dimensional joint Site Frequency Spectrum (3D-JSFS) to infer the ABFT demographic history. The 3D-JSFS was constructed for Mediterranean Sea, Slope Sea and Gulf of Mexico populations using the allele counts of biallelic variants included in the VCF file obtained from the 'nuclear mapped + ALB' catalog, which included 4 albacore samples for variant orientation. Derived allele counts were averaged over all possible resampling of 20 genotypes within each of the three ABFT locations and singletons were excluded. We performed historical demographic model comparison by fitting separately 10 candidate models (Table S6) to the observed JSFS using a diffusion approximation approach implemented in δaδi v1.7.0 [START_REF] Gutenkunst | Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data[END_REF]) and an optimization routine based on consecutive rounds of optimizations [START_REF] Portik | Evaluating mechanisms of diversification in a Guineo-Congolian tropical forest frog using demographic model selection[END_REF]. We adapted existing divergence models to include the three different possible dichotomous branching of the three populations involving two splits, a simultaneous split of the three populations from an ancestral populations and a scenario of split between the Mediterranean Sea and Gulf of Mexico populations followed by an admixed origin of the Slope Sea. We fitted each of these divergence scenarios with or without allowing constant migration rates between populations from split to present. Ancestral effective population size (N A ), migration rates and time estimates scaled to theta (4N A μ ) and the percentage of variable sites correctly oriented with respect to the ancestral state were estimated for all models. Model selection was performed using the Akaike information criterion and goodness-of-fit was assessed by generating 100 Poisson-simulated SFS from the model SFS, fitting the model to each simulated SFS, and using the log-likelihood and log-transformed chi-squared test statistic to generate a distribution of simulated data values against which the empirical values can be compared [START_REF] Portik | Evaluating mechanisms of diversification in a Guineo-Congolian tropical forest frog using demographic model selection[END_REF].

Loci under selection

Loci potentially influenced by selection were screened from the 'nuclear mapped' catalog considering all SNPs using two approaches. The reversible jump Markov chain Monte Carlo approach implemented in BAYESCAN 2.1 [START_REF] Foll | A Genome-Scan Method to Identify Selected Loci Appropriate for Both Dominant and Codominant Markers: A Bayesian Perspective[END_REF] was applied by grouping samples per location, setting default parameters of 50000 burn-in steps, 5000 iterations, 10 thinning interval size and 20 pilot runs of size 5000. Candidate loci under selection with a posterior probability higher than 0.76 (considered as strong according to the Jeffery's interpretation in the software manual) and a false discovery rate (FDR) lower than 0.05 were selected. We then used the multivariate analysis method implemented in the pcadapt R package, which does not require a prior grouping of the samples, following [START_REF] Luu | pcadapt: an R package to perform genome scans for selection based on principal component analysis[END_REF] recommendations and selected outlier SNPs following the Benjamini-Hochberg procedure. Pairwise linkage disequilibrium between all filtered SNPs obtained from those scaffolds which contained candidate SNPs under selection was measured using the R package LDheatmap. PCAs were performed using the adegenet R package [START_REF] Jombart | adegenet 1.3-1: new tools for the analysis of genome-wide SNP data[END_REF] based on outlier SNPs, and variants obtained from one genomic region found to be under high linkage disequilibrium (scaffolds BKCK01000075 (partially) and BKCK01000111) from the 'nuclear mapped' and the 'nuclear mapped + others' catalogs. Individual heterozygosity values based on SNPs within this region from the 'nuclear mapped' were calculated using PLINK [START_REF] Purcell | PLINK: a tool set for whole-genome association and populationbased linkage analyses[END_REF].

Tests for nuclear introgression

Nuclear introgression from albacore to ABFT was tested by applying the statistical model implemented in TreeMix [START_REF] Pickrell | Inference of Population Splits and Mixtures from Genome-Wide Allele Frequency Data[END_REF] and ABBA/BABA analyses [START_REF] Kulathinal | The genomics of speciation in Drosophila: diversity, divergence, and introgression estimated using low-coverage genome sequencing[END_REF][START_REF] Green | A draft sequence of the Neandertal genome[END_REF][START_REF] Durand | Testing for ancient admixture between closely related populations[END_REF] on 'nuclear mapped + other' dataset with only one SNP per tag. The latter test was also performed excluding or including only those SNPs located within the genomic region found under high linkage disequilibrium (scaffolds BKCK01000075 (partially) and BKCK01000111). TreeMix was used to estimate historical relationships among populations and species by estimating the maximum likelihood tree for a set of populations allowing historical gene flow events. TreeMix was run allowing from 0 to 10 migration events, obtaining an increasing number of possible gene flow events and associated likelihood values. We followed the author's recommendations [START_REF] Pickrell | Inference of Population Splits and Mixtures from Genome-Wide Allele Frequency Data[END_REF] to select the most probable number of migration events as that showing best associated likelihood values after stopping adding additional migration events so that the results remained interpretable. The ABBA/BABA test, which measures the excess of derived alleles shared between a candidate donor species and one of two tested groups (in this case, one ABFT group) compared with the other group taken as a reference (a different ABFT group), was performed on the allele frequencies of the derived allele in albacore and ABFT locations, based on the ancestral state defined by the Southern bluefin tuna taken as an outgroup. Patterson's D statistic was calculated for all possible combinations of target and reference groups of ABFT, always considering albacore as the candidate donor species. Additionally, inter-species absolute divergence (dxy) between Mediterranean ABFT larvae and albacore individuals was estimated at each polymorphic position from the 'nuclear mapped + other' catalog. PCAs were performed using the adegenet R package [START_REF] Jombart | adegenet 1.3-1: new tools for the analysis of genome-wide SNP data[END_REF] based on all filtered SNPs and only those SNPs from the region under high linkage disequilibrium from the 'nuclear mapped + others' catalog.

Results

Genetic differentiation between Atlantic bluefin tuna spawning components

We studied the population genetic structure and connectivity of Atlantic bluefin tuna using a genome-wide single nucleotide polymorphisms (SNP) dataset. Our study includes reference samples (i.e., larvae and YoY ABFT captured at or close to where they were hatched and adults caught on the spawning grounds during the spawning season) from the Gulf of Mexico and Mediterranean Sea, and from a more recently discovered spawning ground in the Slope Sea used by ABFT of unknown origin (Figure 1A, Table S1). Significant genetic differentiation between samples from these three spawning grounds was revealed by both an unsupervised clustering analysis of genetic ancestry (ADMIXTURE) (Figures 1B andS2A) and a principal component analysis (PCA) (Figure 1C and S2B), with significant pairwise genetic differentiation (F ST ) between reference larvae from different spawning grounds ranging from 0.0007 to 0.003 (Table S4). This contemporary genetic structure was associated with a mixture of two genetic ancestries (Figure S3), hereafter called GOM-like, predominant in the Gulf of Mexico (average GOM-like proportion across Gulf of Mexico individuals was 0.81 SD ± 0.22), and MED-like, predominant in the Mediterranean Sea (average MED-like proportion across Mediterranean individuals was 0.82 SD ± 0.11). Compared to the Gulf of Mexico and the Mediterranean Sea, the Slope Sea showed a large variance in individual ancestries ranging from GOM-like to MED-like, with a high proportion of admixed ancestries (average GOM-like proportion across Slope Sea individuals was 0.68 SD ± 0.22) (Figure 1D). This supports the Slope Sea as a mixed spawning ground where genetic admixture occurs. Additionally, whereas all Mediterranean individuals had a homogeneous MED-like genetic origin, ancestry profiles of Gulf of Mexico individuals were more variable, including 15 adults (out of 156) (Table S1) with a clear MED-like genetic profile (average GOM-like proportion across Gulf of Mexico individuals excluding 15 MED-like was 0.86 SD ± 0.14) (Figure 1D,C). Otolith microchemistry composition available for 6 of these 15 MED-like adults is compatible with Mediterranean Sea origin (Figure S4). In agreement with these results, admixture tests (F3 statistics) showed that the Slope Sea component is the result of admixture between the two other components (Figure 1E, S5 and Table S5). No admixture was found in the Mediterranean Sea, neither in larvae from the Gulf of Mexico. In contrast, admixture was detected in adult samples from the Gulf of Mexico (Figure 1E). Demographic history inferences (mami) (Figures 1F,S6 and Table S6) supported that the Slope Sea and the Gulf of Mexico spawning components share a recent common ancestry, and that there is strong contemporary migration from the Mediterranean Sea and the Gulf of Mexico towards the Slope Sea. Migration rates in all other directions are much weaker, the strongest being the migration from the Slope Sea back to the Gulf of Mexico, which is three times lower than in the opposite direction.

Overall, these results support a historical scenario of two ancestrally differentiated populations that interbred in the Slope Sea and to a lesser extent in the Gulf of Mexico. This later finding is supported by the presence of spawning capable MED-like individuals in the Gulf of Mexico, including one female that had ovulated less than 48 hours prior to capture (Table S1). A high level of gene flow between GOM-like and MEDlike ancestries, as evidenced by the detection of frequent admixture, could imply either recent demographic changes (e.g. in migration rates or effective population sizes), or the existence of mechanisms that maintain genetic differentiation through time (e. g. cryptic barriers to gene flow or introgression from another source lineage).

Gene flow from Mediterranean Sea towards the Slope Sea revealed by interspecific introgression

To understand why genetic differentiation is maintained despite presumable ongoing gene flow, we studied the effect of interspecific introgression. Mitochondrial introgression from albacore (Thunnus alalunga ) into ABFT has been previously reported, with all introgressed individuals detected so far found in the Mediterranean Sea and the Slope Sea [START_REF] Alvarado Bremer | Comparative phylogeography of Atlantic bluefin tuna and swordfish: the combined effects of vicariance, secondary contact, introgression, and population expansion on the regional phylogenies of two highly migratory pelagic fishes[END_REF][START_REF] Viñas | Facts and uncertainties about the genetic population structure of Atlantic bluefin tuna (Thunnus thynnus) in the Mediterranean. Implications for fishery management[END_REF], but not in the Gulf of Mexico. According to three diagnostic positions for mitochondrial ancestry (Table S2), we found albacore origin introgressed mitochondria in individuals of all age classes in both the Mediterranean Sea (4%) and the Slope Sea (6%), but also to a lower extent in Gulf of Mexico adults (1%) (Figure 2A and Table S1). These results were confirmed at the nuclear level by a tree-based analysis of population splits and admixture using allele frequency data (TreeMix), which supported an introgression event from albacore into the Mediterranean Sea ABFT (Figure 2A andS7). In accordance with this deviation from a strict bifurcating evolutionary history, we also found an excess of derived allele sharing between albacore and both the Slope Sea and the Mediterranean Sea with respect to the Gulf of Mexico (ABBA/BABA test, Figure 2B). However, the contribution of albacore introgression to genetic differentiation between Mediterranean Sea and Gulf of Mexico ABFT populations appears to be limited at best (Figure S8) and therefore it could contribute but not explain the maintenance of genetic differentiation between MED-like and GOM-like individuals.

Our results confirm the introgression of mitochondrial DNA and show for the first-time traces of introgression at the nuclear level from albacore to the ABFT. The gradient of albacore ancestry further suggests that this introgression occurred (or has been more intense) in the Mediterranean Sea where the signal is strongest, and then diffused towards the Slope Sea and, to a lesser extent the Gulf of Mexico through multigenerational gene flow. These different introgression signal intensities from east to west support gene flow between the Mediterranean Sea and Slope Sea spawning components, which is in accordance with the admixed nature of Slope Sea individuals. On the other hand, the differences in the introgression signal intensities also imply reduced gene flow from the Mediterranean Sea to the Gulf of Mexico, which contrasts with the frequency of MED-like spawners observed in the Gulf of Mexico. This, together with the fact that introgression is not the main contributor to genetic differentiation between spawning components, strongly suggests that other mechanisms, such as local adaptation, maintain genetic differentiation in the presence of gene flow, or that migration towards the Gulf of Mexico has increased recently.

Observed genetic differentiation between Atlantic bluefin tuna spawning components cannot be attributed to local adaptation acting on few loci of large effect

To better understand the evolutionary processes behind genetic differentiation in ABFT, we separately studied genetic diversity at neutral (i.e. those that are mostly influenced by genetic drift and migration) and outlier SNPs markers (i.e. those that are potentially under selection or in tight linkage with selected loci). Removing the 123 identified outlier markers did not change the overall population structure pattern nor differentiation values (Figure S9), suggesting that observed genetic differentiation cannot be explained by local adaptation only. On the other hand, analyses based on the 123 markers identified as potentially under selection provided higher genetic differentiation values among spawning grounds (Figure S9) but revealed three groups of samples that do not correspond to the overall population structure (Figure 3A andS10) and that are neither related to laboratory nor phenotypic sex effects (Figure S11). The 20% of the SNP markers that contribute the most to this grouping are located within the same region of the genome (mapping on two scaffolds spanning 2.63 Mb region of the PBFT reference genome) (Table S7) and show strong pairwise linkage disequilibrium (LD) across the whole region (meaning that variant versions of SNP pairs are nonrandomly associated and the same combination is often found among individuals haplotypes) (Figure 3B). The SNPs located within this high-LD region support a three-grouping pattern (Figure S12), with the intermediate group of individuals in the PCA presenting increased heterozygosity values (Figure 3C). This suggests the existence of two main haplotypes (unique allelic combinations across multiple SNPs) in this region combined into three possible genotypes (e.g. AA, AB, BB), which shows characteristics typical of a chromosomal inversion. These two haplotypes, presumably the inverted and collinear versions, are present at different frequencies among spawning grounds, the rarest found to be homozygous only in the Mediterranean Sea, where it is more frequent, and the alternative being more abundant in the Gulf of Mexico and Slope Sea (Figure 3D). A PCA based on genetic markers from this genomic region including otherThunnus species (Figure 4A) showed that homozygous individuals for the most abundant haplotype group associated with the PBFT, whereas those homozygous for the rarest variant were close to albacore. By contrast, the PCA based on the genome-wide SNP dataset showed the expected grouping pattern reflecting species membership, where PBFT and ABFT cluster together and were separated from albacore and Southern bluefin tunas (Figure S13). Test for deviation from a strict bifurcating evolutionary history (ABBA/BABA) showed a much more pronounced signal of introgression from albacore into the Mediterranean and Slope Sea spawning ground samples in the high-LD region (Figure 4B) than when considering the overall genome (Figure 2B). Yet, this pattern remained when removing all SNPs from the high-LD region (Figure S14), indicating that the signal of introgression is present genome-wide.

These results suggest that the genetic differentiation observed between ABFT from different spawning grounds is maintained despite gene flow between the Mediterranean Sea and the Slope Sea and cannot be explained by local adaptation acting on a few loci of large effect. Additionally, a large genomic region of albacore ancestry, introgressed into the ABFT genome in the Mediterranean Sea, has retained high LD while expanding towards the western Atlantic, following the previously detected genome-wide signal of albacore ancestry. Altogether, our results point towards a situation where the two ancestral genetic components of ABFT (Western Atlantic and Mediterranean) have initially diverged in isolation, independently experiencing genetic drift combined with introgression of genetic material from albacore in the Mediterranean Sea. More recently, homogenization between western Atlantic and Mediterranean components could have been initiated by the intensification of gene flow, without completely eroding existing genetic differentiation.

Discussion

Our results, based on a comprehensive ABFT genome-wide SNPs dataset (including larvae from the Slope Sea and spawning adults for the Mediterranean Sea and the Gulf of Mexico), confirmed that current ABFT populations originated from two ancestral populations as previously hypothesized [START_REF] Rodríguez-Ezpeleta | Determining natal origin for improved management of Atlantic bluefin tuna[END_REF]). Yet, they also revealed interbreeding in the Slope Sea and an eastern-western unidirectional gene flow that challenges the assumption of two isolated spawning areas. Moreover, the identification of previously unreported inter-specific introgressed regions in the ABFT nuclear genome and potentially adaptive markers within a newly discovered putative chromosomal inversion provided evidence to suggest that there have been recent changes in ABFT connectivity which holds significant implications for the conservation of the species.

Strong admixture in the Slope Sea as the result of a potential source-sink dynamic

The observed heterogeneously admixed genetic profiles of Slope Sea larvae and YoY ABFT support recurrent interbreeding between migrants from the Gulf of Mexico and the Mediterranean Sea in the Slope Sea, which contributes to the admixed genetic background of this spawning area. This observation is compatible with tagging data, which shows adult individuals that enter the Gulf of Mexico or the Mediterranean Sea also visit the Slope Sea spawning area [START_REF] Block | Electronic tagging and population structure of Atlantic bluefin tuna[END_REF][START_REF] Aalto | Evidence of bluefin tuna (Thunnus thynnus) spawning in the Slope Sea region of the Northwest Atlantic from electronic tags[END_REF]. Otolith microchemistry data provides evidence of individuals with Mediterranean Sea and Gulf of Mexico origin compatible profiles in this area [START_REF] Siskey | Forty years of fishing: changes in age structure and stock mixing in northwestern Atlantic bluefin tuna (Thunnus thynnus) associated with size-selective and long-term exploitation[END_REF].

Our analyses support that the Slope Sea component originated from the Gulf of Mexico population, and that mixing with the Mediterranean population started later. Thus, even if evidence of spawning activity in the Slope Sea dates back to the 1950s [START_REF] Baglin | A preliminary study of the gonadal development and fecundity of the western Atlantic bluefin tuna[END_REF][START_REF] Mather | Historical document : life history and fisheries of Atlantic bluefin tuna[END_REF] and could have started much earlier, it is most likely that the now observed genetic differentiation of the Slope Sea is due to an increase in the immigration rates from the Mediterranean component towards the Slope Sea. In fact, heterogeneous genetic profiles of individual ABFT from the Slope Sea indicate a diverse genetic composition of spawners, a situation at odds with the scenario of an exclusively self-sustained population at equilibrium. Moreover, previous studies using otoliths have reported highly variable proportions of Mediterranean origin individuals in the western Atlantic across the last five decades [START_REF] Secor | Depressed resilience of bluefin tuna in the western atlantic and age truncation[END_REF][START_REF] Siskey | Forty years of fishing: changes in age structure and stock mixing in northwestern Atlantic bluefin tuna (Thunnus thynnus) associated with size-selective and long-term exploitation[END_REF], Rooker et al. 2019[START_REF] Kerr | Mixed stock origin of Atlantic bluefin tuna in the U.S. rod and reel fishery (Gulf of Maine) and implications for fisheries management[END_REF] and [START_REF] Puncher | Individual assignment of Atlantic bluefin tuna in the northwestern Atlantic Ocean using single nucleotide polymorphisms reveals an increasing proportion of migrants from the eastern Atlantic Ocean[END_REF] detected that the proportion of individuals genetically assigned to Mediterranean origin increased over the past two decades at some northwestern Atlantic areas, particularly among individuals younger than 15 years, which is compatible with dynamically changing migratory trends. Demographic connectivity is of major importance for fisheries management, as it directly affects productivity and a stocks recruitment. Despite the limited knowledge about the spawning dynamics in the Slope Sea, our data suggest asymmetrical genetic connectivity towards the Slope Sea, possibly acting as a sink spawning area which is receiving rather than exporting individuals, although its admixed nature could hamper the detection of gene flow from the Slope Sea towards the Mediterranean Sea and the Gulf of Mexico. This highlights the importance of understanding the demographic interdependence of the Slope Sea with the other components, especially in view of recent studies proposing the Slope Sea as a major spawning ground [START_REF] Hernández | Support for the Slope Sea as a major spawning ground for Atlantic bluefin tuna: evidence from larval abundance, growth rates, and particle-tracking simulations[END_REF]. One important knowledge gap is thus the understanding of Slope Sea born individuals' life cycle. The currently observed genetic profiles are compatible with Slope Sea born individuals showing i) Slope Sea spawning-site fidelity, ii) limited spawning, iii) spawning in the Gulf of Mexico and iv) MED-like individuals born in the Slope Sea spawning in the Mediterranean Sea. Unfortunately, weak genetic differentiation, typical in marine fishes with large population sizes, together with the presence of intermediate and heterogeneous (and presumably temporally variable in proportions) genetic profiles hamper the clear identification of Slope Sea born individuals based solely on genetic markers. Thus, we suggest that exploration of the dynamics of these individuals may require the use of integrated methods, such as the combination of genetic markers with otolith microchemistry [START_REF] Brophy | Combining genetic markers with stable isotopes in otoliths reveals complexity in the stock structure of Atlantic bluefin tuna (Thunnus thynnus)[END_REF]. The capability of identifying Slope Sea born individuals and monitoring their presence across the ABFT distribution range, together with an increase in larval sampling efforts in this spawning area, would allow us to obtain and analyze temporal samples to understand their life cycle and estimate admixture rates in the Slope Sea across generations.

Evidence of presence of Mediterranean origin fish in the Gulf of Mexico could suggest recent changes in Atlantic bluefin tuna connectivity patterns

Overall, our results support a historical split between the western Atlantic and Mediterranean spawning grounds followed by a subsequent split between the Gulf of Mexico and Slope Sea and trans-Atlantic unidirectional gene flow from the Mediterranean into western Atlantic spawning grounds. While admixture in the Slope Sea is reflected in the larval and juvenile individual genetic profiles, larvae captured in the Gulf of Mexico were pure GOM-like. Previous work suggested weak input of Mediterranean alleles in the larvae collected from the western Gulf of Mexico in the year 2014 [START_REF] Johnstone | Genetic connectivity between Atlantic bluefin tuna larvae spawned in the Gulf of Mexico and in the Mediterranean Sea[END_REF]. However, we have not detected evidence of such genetic connectivity in larval samples from the western and eastern sides of the Gulf of Mexico collected before this date (from years 2007 to 2010) despite using thousands of SNP genetic markers. Nevertheless, the number of larvae collected in the Gulf of Mexico with individual genetic profile available for this study remains limited (n=27) and the presence of MED-like spawning adults suggests potential genetic connectivity between the Mediterranean Sea and the Gulf of Mexico. Given that Slope Sea individuals' ancestry proportions cover the range of MED-like individual genetic profiles, it would also be possible that these MED-like individuals have their origin in the Slope Sea. The detection of MED-like individuals in the Gulf of Mexico originating from the Slope Sea is made likely due to its proximity.

Due to the heterogeneous profile of the Slope Sea, the observed proportion of MED-like individuals in the Gulf of Mexico originated in the Slope Sea could only be explained by a high number of Slope Sea MEDlike individuals entering the Gulf of Mexico, unless MED-like individuals originated in the Slope Sea under a scenario of even higher inflow from the eastern Atlantic. Otolith microchemistry analyses revealed that genetically MED-like individuals captured in the Gulf of Mexico showed an otolith isotopic composition of oxygen (δ 18 O) intermediate between the Gulf of Mexico and the Mediterranean spawning areas, suggesting that these were probably not born in the Gulf of Mexico. Interestingly, these intermediate values are consistent with the proposed signature range of a potential third contingent, compatible with a Slope Sea or Mediterranean origin of individuals showing early and/or more intense migratory behavior [START_REF] Brophy | Combining genetic markers with stable isotopes in otoliths reveals complexity in the stock structure of Atlantic bluefin tuna (Thunnus thynnus)[END_REF]. These observations allow for different possible origins of the MED-like individuals captured in the Gulf of Mexico. Additional observations of the genetic composition of adult ABFT spawning in the Slope Sea coupled with further knowledge about the migratory behavior of Mediterranean ABFT would be needed to assign the origin of these MED-like individuals more accurately.

Regardless of the origin of the MED-like individuals in the Gulf of Mexico, a few dozen migrants exchanged per generation is theoretically sufficient to erase genetic differentiation between populations [START_REF] Waples | Separating the wheat from the chaff: patterns of genetic differentiation in high gene flow species[END_REF][START_REF] Lowe | What can genetics tell us about population connectivity?[END_REF][START_REF] Gagnaire | Using neutral, selected, and hitchhiker loci to assess connectivity of marine populations in the genomic era[END_REF]. The low F ST values reported in this study are common among marine fishes with large population sizes, high dispersal rates and wide-ranging distributions [START_REF] Hauser | Paradigm shifts in marine fisheries genetics: ugly hypotheses slain by beautiful facts[END_REF][START_REF] Da Fonseca | Population genomics reveals the underlying structure of the small pelagic European sardine and suggests low connectivity within Macaronesia[END_REF][START_REF] Fuentes-Pardo | The genomic basis and environmental correlates of local adaptation in the Atlantic horse mackerel[END_REF]. While genetic differentiation of ABFT between the Mediterranean Sea and the Gulf of Mexico could persist despite admixed individuals in the Slope Sea, the number of migrants detected in the Gulf of Mexico is theoretically expected to lead to genetic homogeneity between eastern and western born ABFT and is thus not easily compatible with significant F ST values. Histological inspection confirmed the presence of at least one MED-like female which had spawned less than 48 hours before capture in the Gulf of Mexico, suggesting that MED-like individuals spawn in the Gulf of Mexico. One possible explanation for the observed levels of genetic differentiation would be negative selection against Mediterranean genes preventing successful gene flow. Hence, we have explored different sources of genetic variation, which could help to explain the maintenance of genetic differentiation between highly demographically connected populations. More specifically, we have explored the effect of locally adaptive alleles, which could maintain genetic differentiation in the presence of gene flow [START_REF] Tigano | Genomics of local adaptation with gene flow[END_REF] and inter-specific introgression, which can trigger different evolutionary processes, such as the input of adaptive alleles [START_REF] Huerta-Sánchez | Altitude adaptation in Tibetans caused by introgression of Denisovan-like DNA[END_REF] or the contribution to reproductive isolation [START_REF] Duranton | The contribution of ancient admixture to reproductive isolation between European sea bass lineages[END_REF]). However, we did not find evidence to confirm the maintenance of genetic differentiation through either local adaptation or reproductive isolation. These would lead to much higher levels of F ST at loci involved in incompatibilities and/or selection than the ones observed in this work. Besides, genome-wide and homogeneously distributed genetic differentiation between Mediterranean Sea and Gulf of Mexico reference individuals at neutral alleles reflects that differentiation is not primarily driven by introgression or adaptation, but by the effect of historical long-term isolation between Mediterranean and Atlantic populations, as indicated by their inferred demographic history. Moreover, interbreeding in the Slope Sea implies genetic compatibility between GOM-like and MED-like individuals, which makes the barrier to gene flow hypothesis unlikely to explain the maintenance of genetic differentiation. Another possibility is that selection undetected in this study impedes incorporation of MED-like alleles in the Gulf of Mexico. On one hand, widespread polygenic selection remains difficult to reject based on outlier detection tests as the ones used here, as these are underpowered to detect weakly selected loci, although genetic differences are unlikely maintained by weak polygenic selection in the presence of gene flow. On the other hand, the use of reduced representation sequencing could lead to missed localized selective sweeps. Analyses based on whole genome sequencing data would allow one to detect adaptation signals missed in our study. Alternatively, we propose that the currently observed ancestry patterns could be explained by recent secondary contact following genetic divergence of both ancestral populations after long-term isolation with reduced or no migration, and that the high observed migration rates, in the absence of barriers to gene flow and if sustained over time, could ultimately lead to genetic homogenization and the consequently the loss of genetic differentiation.

Possible drivers and implications of changes in inter-spawning area connectivity

The observed levels of genetic differentiation, hardly compatible with the high level of gene flow suggested by our results in constant equilibrium, suggest that connectivity patterns between ABFT spawning grounds could be subjected to temporal changes and that an increase in migration from the Mediterranean Sea towards the known western Atlantic spawning areas could have a genetic homogenizing effect. Such a homogenizing effect is expected to be correlated with migration rates and the effective population size (Ne) of the recipient population [START_REF] Lowe | What can genetics tell us about population connectivity?[END_REF]Allendorf 2010, Gagnaire et al. 2015), which in turn relates to the number of adult individuals among other factors [START_REF] Waples | What Is Ne, Anyway?[END_REF]) and would consequently be affected by fluctuations in the population's abundance. The abundance of ABFT stocks have undergone strong changes during the last ˜60 years. After both the western and the eastern Atlantic stocks reached a critical status, including the collapse of several fisheries around the 1960-1970s (Fromentin 2009[START_REF] Porch | The journey from overfishing to sustainability for Atlantic bluefin tuna, Thunnus thynnus[END_REF], the western Atlantic stock has not recovered as rapidly as the eastern Atlantic stock (of Mediterranean origin), whose estimated abundance has been one order of magnitude larger for several decades (ICCAT 2017), despite decades of conservation efforts. This slow recovery could result from a regime shift over the last decades, due to the combination of oceanographic changes in the equatorial Atlantic and overfishing (Fromentin et al. 2014b) possibly affecting both migration rates and effective population sizes. Fluctuations in the abundance and distribution of eastern Atlantic bluefin tuna during the last century were largely explained by the Atlantic Multidecadal Oscillation (AMO) [START_REF] Faillettaz | Atlantic Multidecadal Oscillations drive the basin-scale distribution of Atlantic bluefin tuna[END_REF]) and long-term trends in temperature [START_REF] Ravier | Are the long-term fluctuations in Atlantic bluefin tuna (Thunnus thynnus) population related to environmental changes?[END_REF]. Moreover, coinciding with the last negative AMO period starting in the 1960s, ABFT had disappeared from several North-East Atlantic areas where it is reappearing during the current positive AMO phase starting in the mid-1990s [START_REF] Horton | Tracking Atlantic bluefin tuna from foraging grounds off the west coast of Ireland[END_REF][START_REF] Nøttestad | The comeback of Atlantic bluefin tuna (Thunnus thynnus) to Norwegian waters[END_REF][START_REF] Aarestrup | First tagging data on large Atlantic bluefin tuna returning to Nordic waters suggest repeated behaviour and skipped spawning[END_REF]. Furthermore, increasing catches of ABFT in Greenland waters show that the northern limit of ABFT distribution was expanded northwards during the last decade by mostly individuals of Mediterranean genetic origin [START_REF] Jansen | Atlantic bluefin tuna (Thunnus thynnus) in Greenland -mixed-stock origin, diet, hydrographic conditions, and repeated catches in this new fringe area[END_REF]. Based on electronic tagging data, the proportion of individuals of eastern origin present in the western Atlantic has also increased during the last two decades [START_REF] Aalto | Stock-of-origin catch estimation of Atlantic bluefin tuna (Thunnus thynnus) based on observed spatial distributions[END_REF]. Interestingly, AMO positive or warm phases as well as current global warming involve an increase in habitat suitability in most of these northern areas (Fromentin et al. 2014b[START_REF] Faillettaz | Atlantic Multidecadal Oscillations drive the basin-scale distribution of Atlantic bluefin tuna[END_REF]). In summary, population size, distribution, and migratory behavior of ABFT has been undergoing changes during the last decades, probably due to changes in environmental conditions, fishing pressure, conservation efforts or combined effects of these. These changes could explain migration intensification from the recently expanded eastern stock towards the more slowly recovering western Atlantic stock. Our analyses support that the recent short-term genetic effects of immigration on the variance in ancestries are much stronger in the Slope Sea than in the Gulf of Mexico. This could be due to behavioral preferences or favorable conditions which would make it easier for Mediterranean individuals to reach and reproduce in the Slope Sea, or due to a smaller effective population size compared to the Gulf of Mexico. The homogenizing effect of this unidirectional migration would ultimately lead to the genetic swamping of the GOM-like genetic component. It is thus possible that genetic differentiation is detected in the existing samples because gene flow is relatively recent relative to mean generation times estimated to average 9.6 years for the Western stock [START_REF] Collette | High Value and Long Life -Double Jeopardy for Tunas and Billfishes[END_REF], and that this genetic divergence will be attenuated across generations in the future. While genetic connectivity does not necessarily equate with demographic dependence, genetic erosion would not necessarily imply demographic decline either, but could have more unpredictable demographic consequences.

Inter-specific gene introgression from albacore to Atlantic bluefin tuna

We detected for the first-time signatures of introgression from albacore tuna in the nuclear genome of the ABFT, which contradicts a previous report [START_REF] Ciezarek | Phylotranscriptomic Insights into the Diversification of Endothermic Thunnus Tunas[END_REF]. The most probable inter-lineage gene flow event estimated by the TreeMix analysis happened between the albacore and the ABFT Mediterranean population. The presence of mitochondrial introgression has also been reported in PBFT [START_REF] Chow | Phylogenetic relationships between tuna species of the genus Thunnus (Scombridae: Teleostei): Inconsistent implications from morphology, nuclear and mitochondrial genomes[END_REF]; however, although we included very few PBFT samples, we did not find any sign of nuclear introgression in this species. Considering that ABFT and PBFT evolved from a recent common ancestor [START_REF] Díaz-Arce | RAD-seq derived genome-wide nuclear markers resolve the phylogeny of tunas[END_REF]) and that they show very little genetic divergence, our results suggests that introgression between albacore tuna and ABFT happened after the split between the ABFT and PBFT lineages. Thus, to explain the presence of albacore-like mitochondrial genomes in PBFT we hypothesize either parallel introgression events between albacore tuna and both ABFT and PBFT, or that mitochondrial introgressed genomes present in PBFT have been introgressed through genetic exchanges with the ABFT. Likewise, among ABFT individuals, the signal of introgression from albacore is stronger in the Mediterranean and the Slope Sea and nearly absent in the Gulf of Mexico. Thus, we conclude that introgression from albacore most likely occurred in the Mediterranean ABFT population, where spawning areas for ABFT and albacore overlap [START_REF] Alemany | Characterization of the spawning habitat of Atlantic bluefin tuna and related species in the Balearic Sea (western Mediterranean)[END_REF], and was subsequently transmitted to the Slope Sea. Besides, the nearly complete absence of both nuclear and mitochondrial introgression in Gulf of Mexico individuals suggests that introgression happened after the split between MED-like and GOM-like ancestral lineages which is consistent with the inferred scenario of historically restricted gene flow between east Atlantic and the Gulf of Mexico spawning components. Overall, the east-west gradient of introgressed albacore alleles confirms the inferred connectivity patterns presented in this study.

We identified a particular genomic region with characteristics typical of a chromosomal inversion, such as high linkage disequilibrium and increased heterozygosity values in samples occupying intermediate positions in a local PCA [START_REF] Barth | Disentangling structural genomic and behavioural barriers in a sea of connectivity[END_REF][START_REF] Puncher | Chromosomal inversions in the Atlantic cod genome: Implications for management of Canada's Northern cod stock[END_REF][START_REF] Jiménez-Mena | Weak genetic structure despite strong genomic signal in lesser sandeel in the North Sea[END_REF]. Our analysis supports that, as reported for other species [START_REF] Jay | Supergene Evolution Triggered by the Introgression of a Chromosomal Inversion[END_REF], the origin of this inversion was introduced into the ABFT genome as the result of a past introgression event from albacore tuna. This region aggregates a high number of outlier genetic markers. However, high linkage-disequilibrium could bias towards the detection of false positives within this region due to the synergic signal of dozens of variants. We could not associate the presence of this introgression nor chromosomal inversion with ecological or environmental factors; yet, introgression represents an important source of genetic adaptive variation playing an important role favoring speciation through processes such as introgression of favored alleles [START_REF] Arnold | Adaptation by introgression[END_REF][START_REF] Hedrick | Adaptive introgression in animals: examples and comparison to new mutation and standing variation as sources of adaptive variation[END_REF][START_REF] Clarkson | Adaptive introgression between Anopheles sibling species eliminates a major genomic island but not reproductive isolation[END_REF] or reproductive isolation [START_REF] Abbott | Hybridization and speciation[END_REF][START_REF] Duranton | The contribution of ancient admixture to reproductive isolation between European sea bass lineages[END_REF]. Moreover, the literature abounds with examples showing that chromosomal inversions are associated with local adaptation in the presence of gene flow [START_REF] Berg | Three chromosomal rearrangements promote genomic divergence between migratory and stationary ecotypes of Atlantic cod[END_REF], Barth et al. 2017a[START_REF] Wellenreuther | Eco-Evolutionary Genomics of Chromosomal Inversions[END_REF][START_REF] Huang | Multiple chromosomal inversions contribute to adaptive divergence of a dune sunflower ecotype[END_REF], Le Moan et al. 2020[START_REF] Mérot | Locally Adaptive Inversions Modulate Genetic Variation at Different Geographic Scales in a Seaweed Fly[END_REF][START_REF] Thorstensen | A chromosomal inversion may facilitate adaptation despite periodic gene flow in a freshwater fish[END_REF]. Given that allele frequency differences between the Mediterranean and the Gulf of Mexico components are stronger in the candidate chromosomal inversion than the mean genome-wide differentiation, ascertaining its role in the ABFT adaptation could be of great relevance to understand the species resilience to changes in environmental conditions [START_REF] Muhling | Predicting the effects of climate change on bluefin tuna (Thunnus thynnus) spawning habitat in the Gulf of Mexico[END_REF][START_REF] Erauskin-Extramiana | Large-scale distribution of tuna species in a warming ocean[END_REF].

Implications for Atlantic Bluefin tuna conservation and management

Conservation of ABFT is challenged by past and future fishing pressure (Fromentin et al. 2014a[START_REF] Secor | Depressed resilience of bluefin tuna in the western atlantic and age truncation[END_REF], which has sharply increased over the last five years following the rebuilding of the Mediterranean ABFT population (ICCAT 2023), and by changes in environmental conditions (often interacting with fishing pressure), which have been shown to alter population size and productivity, migratory behavior and spatial distribution [START_REF] Ravier | Are the long-term fluctuations in Atlantic bluefin tuna (Thunnus thynnus) population related to environmental changes?[END_REF]. From a conservation perspective, hybridization between genetically differentiated lineages, in this case between GOM-like and MED-like individuals, could increase each population's genetic diversity, leading to the incorporation of potentially adaptive genomic variation and reducing vulnerability to environmental changes [START_REF] Brauer | Natural hybridization reduces vulnerability to climate change[END_REF]. However, strong unidirectional gene flow could provoke genetic swamping of the western Atlantic spawning areas jeopardizing ABFT genetic diversity [START_REF] Roberts | Marine genetic swamping: hybrids replace an obligately estuarine fish[END_REF]. In this sense, large effective population sizes, which could increase following a rebuilding of abundance at the different spawning areas [START_REF] Hoey | Genetic decline and recovery of a demographically rebuilt fishery species[END_REF]) would counteract the homogenizing effect of genetic drift. In the absence of accurate estimations of ABFT effective populations sizes [START_REF] Puncher | Spatial dynamics and mixing of bluefin tuna in the Atlantic Ocean and Mediterranean Sea revealed using next-generation sequencing[END_REF], further genetic monitoring of temporal samples could help to understand potential ongoing trends in genetic diversity conservation [START_REF] Hoban | Comparative evaluation of potential indicators and temporal sampling protocols for monitoring genetic erosion[END_REF][START_REF] Oosting | Unlocking the potential of ancient fish DNA in the genomic era[END_REF].

From a fisheries management perspective, the confirmation of ongoing admixture in the Slope Sea challenges the paradigm of two isolated ABFT stocks. However, large knowledge gaps related to the dynamics of Slope Sea individuals, the magnitude of the Slope Sea spawning in terms of recruitment, and its demographic connectivity with other components hinders explicit modelling of it as a distinct stock. Nonetheless, the recently adopted management procedure (ICCAT 2023) does explicitly consider spawning in the Slope Sea. Our study highlights the need for further monitoring combining multidisciplinary data such as larval sampling, tagging, otolith microchemistry signature and genetic origin to understand the Slope Sea population dynamics and the relevance of this spawning area in demographic and evolutionary terms. Yang, J., S. H. Lee, M. E. Goddard, and P. M. Visscher. 2011. GCTA: A Tool for Genome-wide Complex Trait Analysis. The American Journal of Human Genetics 88 :76-82. S4 and Figure S5). (f) Visual representation of the best-fit demographic model, where arrow and branch widths are proportional to directional migration rates (m) and effective population sizes (n) respectively, and where T represents the duration of population splits. Estimated parameter values are given in units of 2nA, where nA is the effective size of the ancestral population, related to the population-scaled mutation rate parameter of the ancestral populations by θ=4nAμ. Numbers indicate the percentage of individuals (from those included in the tree) showing the introgressed mitochondrial haplotype for each location and age class (abbreviations as in Figure 1). On the upper right, zoom on the phylogenetic relationships among Atlantic bluefin tuna groups. (b) D statistical values estimated from the ABBA/BABA test used to detect introgression from albacore to different targets (rows) using different references (colours). The higher the value, the more introgressed is the target group respect to the reference. values from an ABBA/BABA test based on variants located within the genomic region of high-linkage disequilibrium, using Southern bluefin tuna as an outgroup, albacore as a donor species and all different groups of Atlantic bluefin tuna as alternative targets ordered along the y-axis.

Figure 1 .

 1 Figure 1. Population structure and connectivity of Atlantic bluefin tuna. (a) Map showing capture location and life stage of Atlantic bluefin tuna samples included in this study. Capture location of adults from the Gulf of Mexico are enclosed within the purple rounded polygon to fulfil confidentiality requirements. (b) Estimated individual ancestry proportions assuming two ancestral populations. (c) Principal Component Analysis (PCA) of genetic variability among Atlantic bluefin tuna samples, following colour codes identical to (b). (d) Density distribution of individual MED-like ancestry proportions per spawning ground. (e) F3statistics for each combination of sources and target populations, where the Slope Sea and Mediterranean Sea contain larvae and young of the year and larvae, young of the year and adults, respectively (see detailed results in TableS4and FigureS5). (f) Visual representation of the best-fit demographic model, where arrow and branch widths are proportional to directional migration rates (m) and effective population sizes (n) respectively, and where T represents the duration of population splits. Estimated parameter values are given in units of 2nA, where nA is the effective size of the ancestral population, related to the population-scaled mutation rate parameter of the ancestral populations by θ=4nAμ.

Figure 2 .

 2 Figure 2. Interspecific introgression between albacore and Atlantic bluefin tuna. (a) Phylogenetic tree estimated by TreeMix based on nuclear data allowing one migration event (the arrow indicates migration direction and rate).Numbers indicate the percentage of individuals (from those included in the tree) showing the introgressed mitochondrial haplotype for each location and age class (abbreviations as in Figure1). On the upper right, zoom on the phylogenetic relationships among Atlantic bluefin tuna groups. (b) D statistical values estimated from the ABBA/BABA test used to detect introgression from albacore to different targets (rows) using different references (colours). The higher the value, the more introgressed is the target group respect to the reference.

Figure 3 .

 3 Figure 3. Outlier markers in Atlantic bluefin tuna cluster within one 2.63Mb genomic regions showing high long-distance linkage disequilibrium. (a) PCA performed using the 123 outlier SNPs showing the threecluster grouping (shades of blue) where shapes and colours of samples are those indicated in Figure 1. (b) SNP pairwise linkage disequilibrium plot among the 110 SNPs found within a high linkage region covering scaffolds BKCK01000075 (partially) and BKCK01000111 of the reference genome where most of the SNPs contributing to PC1 from (a) are located. (c) Boxplot showing heterozygosity values (y axis) at the three sample groups show in (a), represented by the same blue colour , and based on the 110 SNPs within the genomic region shown in (b). (d) Proportion of samples from each location and age class assigned to each of the three groups shown in (c).

Figure 4 .

 4 Figure 4. Evolutionary origin of Atlantic bluefin tuna within the region of high-linkage disequilibrium. () Principal component analysis (PCA) including other Thunnus species (albacore in blue, Southern bluefin tuna in green and Pacific bluefin tuna in yellow) performed using 156 genetic variants located within the genomic region under high linkage disequilibrium hosting a candidate structural variant (b) Estimated Dvalues from an ABBA/BABA test based on variants located within the genomic region of high-linkage disequilibrium, using Southern bluefin tuna as an outgroup, albacore as a donor species and all different groups of Atlantic bluefin tuna as alternative targets ordered along the y-axis.

Posted on 22 Jun 2023 -The copyright holder is the author/funder. All rights reserved. No reuse without permission. -https://doi.org/10.22541/au.168742629.93690601/v1 -This a preprint and has not been peer reviewed. Data may be preliminary.

Acknowledgements

We wish to thank all the participants of the GBYP biological sampling program, in particular Naiara Serrano, Xiker Salaberria and Inma Martín for tissue sample handling and storage, Iraide Artetxe-Arrate and Ai Kimoto for database management and Iñaki Mendibil and Elisabete Bilbao for excellent laboratory work. We also would like to thank Marty Kardos for his valuable comments, contributing to the quality improvement of this manuscript. This study has been carried out under the provision of the ICCAT Atlantic Wide Research Program for Bluefin Tuna (GBYP), funded by the European Community (Grant SI2/542789), Canada, Croatia, Japan, Norway, Turkey, United States (NMFS NA11NMF4720107), Chinese Taipei, and the ICCAT Secretariat. This work was also supported by a grant (project GENGES) and a predoctoral

fellowship to ND-A, from the Department of Agriculture and Fisheries of the Basque Government. The contents of the paper do not necessarily reflect the point of view of ICCAT or of the other funders. This manuscript is contribution number X [to be added upon acceptance] from the Marine Research Division of AZTI.

Data Accessibility Statement

Demultiplexed sequences are deposited in the SRA (Bioproject PRJNA804694).

Scripts used to perform the analyses described in this manuscript can be found athttps://github.com/rodriguez-ezpeleta/abftconnect.

Benefit-Sharing Statement

This research is a result of a collaborative agreement between partners which are all included as co-authors. Each Partner respected the Nagoya Protocol on Access and Benefit-sharing entered into force on the 12 October 2014 produced by the United Nations Convention on Biological Diversity to carry out its activities under this Agreement and performed the necessary formalities to the competent authorities.

Author contributions

NDA, HA and NRE designed research. NDA, PAG, SAH, AP and NRE contributed analytical tools. DER, JFW, PA, FA, RA, SD, ARH, FSK, JMQ and JRR contributed samples. NDA analyzed data. NDA, PAG, DER, JFW, SAH, JMF, DB, ML, IF, NG, JR, HA and NRE interpreted data. NDA wrote the paper, with insightful contributions from all authors. All authors revised the manuscript and agreed with its publication.

Declaration of Interests

The authors declare no competing interests.