Lack of genetic differentiation in yellowfin tuna has conservation implications in the Eastern Pacific Ocean.

Yellowfin tuna, Thunnus albacares, is an important global fishery and of particular importance in the Eastern Pacific Ocean (EPO). According to the 2019 Inter-American Tropical Tuna Commission (IATTC) assessment, yellowfin tuna within the EPO is a single stock, and is being managed as one stock. However, previous studies indicate site fidelity, or limited home ranges, of yellowfin tuna which suggests the potential for multiple yellowfin tuna stocks within the EPO, which was supported by a population genetic study using microsatellites. If numerous stocks are present, management at the wrong spatial scales could cause the loss of minor yellowfin tuna populations in the EPO. In this study we used double digestion RADseq to assess the genetic structure of yellowfin tuna in the EPO. A total of 164 yellowfin tuna from Cabo San Lucas, México, and the Galápagos Islands and Santa Elena, Ecuador, were analysed using 18,011 single nucleotide polymorphisms. Limited genetic differentiation (FST = 0.00058-0.00328) observed among the sampling locations (México, Ecuador, Peru, and within Ecuador) is consistent with presence of a single yellowfin tuna population within the EPO. Our findings are consistent with the IATTC assessment and provide further evidence of the need for transboundary cooperation for the successful management of this important fishery throughout the EPO.

Introduction

Overfishing is one of the greatest threats to the yellowfin tuna, Thunnus albacares [1]. The International Union for the Conservation of Nature (IUCN) has categorized this species as least concern, with a declining population trend [2]. Five intergovernmental associations, the Tuna Regional Fisheries Management Organizations (tRFMOs), are dedicated to monitor and manage tuna and tuna-like populations [3]. Some of these organizations have worked on rebuilding stocks, most noticeably those of Atlantic bluefin tunas [4]. Four tRFMOs manage populations of the yellowfin tuna, these are: the International Commission for the Conservation of Atlantic Tunas (ICCAT); the Indian Ocean Tuna Commission (IOTC); the Western and Central Pacific Fisheries Commission (WCPFC); and the Inter-American Tropical Tuna Commission (IATTC) [5]. Efforts to rebuild yellowfin tuna fisheries have been implemented throughout the species range [3], which has included reducing quotas of fleets that surpassed their limits [6].

Various studies suggest that, at a global level, yellowfin tuna is split into several populations. Genetic assessments by tRFMOs have reported the presence of four different stocks, which are located in the Atlantic Ocean, Indian Ocean, Eastern Pacific Ocean (EPO), and West Central Pacific Ocean [5, 7–9]. However, the finer resolution of population structure within oceans shows contradictory results, which could be partially attributable to the lack of a suitable model for the analysis of genetic marker data of marine species. The latest assessment on the status of yellowfin tuna in the EPO assume the existence of a single stock [10]. However, Díaz-Jaimes & Uribe-Alcocer [11] suggest genetic differentiation between yellowfin tuna sampled in the northern (Gulf of California, Mexican coast and southwest Revillagigedo Islands) and southern EPO (Perú). Their study used microsatellite markers (7 loci) to detect genetic differences among subregions; a method that is not always reliable for determining the genetic structure in species like tunas due to their large effective population sizes [12]. Ely et al. [13] estimated the long term effective population size of females with mismatched distributions of D-loop sequences, suggesting a large long-term number of effective female breeders and low levels of genetic differentiation between Atlantic and Pacific yellowfin tuna individuals. Conversely, Barth et al. [8] estimated the demographic history of yellowfin tuna, detecting significant differentiation of Atlantic and Indo-Pacific yellowfin tuna populations and showing that the true values of effective size (Ne) are far below those estimated by Ely [13], which may suggest greater genetic differentiation within populations. However, the sample size used was smaller than 12 individuals per sampling area. Estimations such as the one for Atlantic bluefin tuna with large sample size have not been carried out in yellowfin tuna [14].

In addition to genomic analyses, a number of other studies suggest that yellowfin tuna have a home range of 1,800 km, and exhibit high regional fidelity [15-21]. Studies which support the high fidelity theory include, analyses of yellowfin tuna parasites [22], stable isotopes [23], and pop-up satellite tags [24, 25]. The findings of these studies support what Díaz-Jaimes & Uribe-Alcocer found [11], however, as we previously mentioned, the use of microsatellite in their analyses could have some limitations due to their sample size and the utility of microsatellites. The limitations of microsatellites can be overcome using Next Generation Sequencing (NGS); a meta-analysis of yellowfin tuna population genetic studies emphasised the importance of NGS over other techniques [26]. NGS techniques offer a broader representation of the genome and allows the identification of single nucleotide polymorphisms (SNPs) loci that are possibly undergoing selection. SNPs can be found in DNA fragments less than 100 bp in size versus microsatellites that have larger amplicons [27]. Though more markers are required for SNPs analyses, high numbers of SNPs provide greater power for assessing genetic structure [28]. Additionally, these loci provide a better opportunity to distinguish populations with low differentiation, by discriminating neutral evolutionary processes like genetic drift from those caused by selection that could result in local adaptations [29].

Hauser and Carvalho [30] emphasize the need for regular updates of population studies of commercial fish species because populations vary over time, and there is a risk of losing minor populations when multiple populations are managed as a single stock [31]. Considering that several studies suggest residence and limited movement of yellowfin tuna, as well as structured populations (albeit with a small sample size), multiple yellowfin tuna stocks may be present in the EPO. Currently, yellowfin tuna in the EPO is managed as a single stock, which could cause the loss of minor yellowfin tuna populations in the region, if they exist. NGS approaches in marine organisms provide a finer resolution of population structure using several thousand SNPs [32, 33], and is suggested to be the most appropriate tool to yellowfin tuna populations [26]. Integrating genomic data into fisheries management is a critical step to improve management practices, to ensure the ecological sustainability of fish stocks [34]. By using NGS approaches, this study tested the hypothesis that a single genetic stock of yellowfin tuna exists within the EPO, to determine the most suitable spatial scale for the effective management of this important fishery.

Materials and methods

Sample collection

Yellowfin tuna samples were collected from three artisanal fishing ports each associated with distinct fishing grounds: [1] CSL-MEX: Cabo San Lucas (20.985°N, -109.539°W), La Paz, México, covering the northern tropical EPO; [2] GAL-ECU: Santa Cruz (0.132°S, -91.383°W), Galápagos Islands, Ecuador, covering the central tropical EPO; and [3] STR-ECU: Santa Rosa fishing port (-4.022°S, -82.953°W), Santa Elena, Ecuador, covering the continental coast of Ecuador and northern Perú (Fig 1). Sampling occurred over a three-year period, from February 2015 to October 2017, under permits MAE-DNB-CM-0046 and MAE-DNB-CM-2016-0041. All yellowfin tuna sampled in this study were caught by artisanal fishing fleets for commercial purposes and no fish was euthanized for this research. Samples from a total of 230 individual yellowfin tuna were initially collected, 64 from CSL-MEX, 103 from GAL-ECU, and 63 from STR-ECU. The sample size of this study is comparable with other studies (20–40 samples per location) testing similar hypotheses [7, 14, 35] and enough to provide an accurate genetic assessment of the yellowfin tuna population [36].

Fig 1

Approximate sampling areas of yellowfin tuna in the Eastern Pacific Ocean.

This figure was produced in QGIS.

Fork length (FL) of each individual fish was measured to the nearest centimetre (cm) using a fish board; fish ranged in size from 34–200 cm FL. A sample of red muscle, approximately 1cm3 in size, was removed from individual yellowfin tunas, at fishing landing sites. New disposable scalpels were used to collect each sample, to avoid cross contamination. Muscle samples were stored in individually labelled Eppendorf tubes in 97% ethanol (ensuring the entire sample was immersed) for transportation and further processing. Sexing individuals requires the extraction of their gonads, which in turn reduces the commercial value of the fish. Consequently, we could not determine our samples’ sexes as fish were sampled on the proviso that they retained their market value.

ddRAD genotyping and quality filtering

DNA extraction, double digestion RAD library preparation and sequencing, and SNP genotyping were carried out by Diversity Arrays Technology (DArT, Canberra, Australia). The patented protocol DArTseq™ was used as it is a cost-effective methodology, and has been optimized for the analysis of yellowfin tuna high-throughput genotyping data [7, 37, 38]. Following automated DNA extraction using Tecan Freedom EVO 150 robots (Männedorf, ZH, Switzerland), samples were digested using PstI and SphI methylation-sensitive restriction enzymes. These enzymes were chosen as they avoid the repetitive fractions of the genome, which is uninformative for SNP discovery [26, 38], and also because they have been used in previous NGS studies of yellowfin tuna [7]. PstI and SphI compatible adapters, which included an Illumina flow cell attachment sequence and a PCR primer sequence, were ligated to the restriction fragments. Additionally, the PstI adapter incorporated unique barcode sequences to identify samples within pooled libraries. Only fragments capped with both PstI and SphI adaptors were amplified. PCR conditions consisted of an initial denaturation step at 94°C for 1 min, followed by 30 cycles at 94°C for 20 sec, 58°C for 30 sec, and 72°C for 45 sec, with a final extension step at 72°C for 7 min. Bridge PCR was then applied to further amplify and sequence the libraries on the Illumina HiSeq 2500 (San Diego, CA, USA) using MJ Tetrad PCR thermocyclers (Hercules, CA, USA) [7].

The resulting sequences were processed using DarTseq analytical pipelines DarT-Soft14 version [7]. The primary quality control pipeline consisted of two filtering steps of the FASTQ files, applying decreasingly stringent selection criteria to the barcode region to increase demultiplexing accuracy. The first step stringently filters out barcode sequences with a Phred score threshold of 30, and the second step filters out sequencing reads with Phred score threshold of 10 [39]. Identical sequences were then collapsed into fastqcoll files and “groomed” using DarT’s proprietary algorithm. This algorithm replaces a low-quality base from a singleton tag with the correct base through a template composed of collapsed tags with multiple members. The corrected fastqcoll files were then used in the DArT-Soft14 secondary pipeline for SNP calling. The final dataset produced by these pipelines consisted of individual genotype reports for each processed sample presented as a matrix of SNP loci, their global call rate, polymorphic information content (PIC), and their co-dominant status.

The SNP dataset was further filtered using the package radiator [40] in R v3.5.3 [41] to remove outliers and noise that could interfere with polymorphism estimates. Filter parameters for SNP quality control included DArT reproducibility (proportion of reproducible genotypes calculated through the comparison of technical replicates), monomorphic markers (to remove markers in one unique genotype), common markers (to remove markers not found across all individuals), minor allele count, call rate (proportion of samples successfully genotyped), SNP number (maximum number of SNPs per marker), SNP position on the fragment sequenced, linkage disequilibrium, Hardy-Weinberg Equilibrium, duplicated genomes (potential duplicated samples estimated via pairwise genome similarity) and mixed genomes (potential mixed samples or poor SNP discovery caused by DNA quality, sequencing procedure, among others). Quality control for individual samples with low DNA quality or DNA contamination included filters of missingness (proportion of non-genotyped values), heterozygosity (proportion of heterozygous loci) and total coverage (average sequencing depth across markers) (S1 Table). DArT DS14 pipeline delivered 80,000 filtered markers. After filtering using radiator [40], a total of 164 individuals: 35 from CSL-MEX, 88 from GAL-ECU, and 41 from STR-ECU; and 18,011 loci were retained after quality control (S1 Table) and used for the remaining analyses.

Data analysis

Basic genetic diversity measures of the quality controlled data set were conducted in R v3.5.3 [41]. Population genetic statistics such as Allelic richness (Ra), observed heterozygosity (Ho) and expected heterozygosity (He) were calculated using the adegenet package [42]. A Mann Whitney test was applied to determine the statistical significance of observed heterozygosity. Analyses for Molecular Variance (AMOVA) were run in the R package Poppr [43]. Global and pairwise F values were calculated using DartR [44], and the inbreeding coefficient was calculated using the R package InbreedR [45]. The mean probability of identity by descent across loci was calculated using the R package rrblup [46]. Population genetic structure among the three study areas was assessed by genetic differentiation (F), the absolute allele frequency difference (AFD) [47], Discriminant Analysis of Principal Components (DAPC), Principal Component Analysis (PCA), and non-spatial clustering (Structure). To determine if any observed population structuring was a result of isolation by distance, a Mantel test was performed. All analyses were conducted using the R package DartR [44]. PCA and non-spatial clustering analysis were conducted in STRUCTURE [48] and CLUMPP [49] software packages using the StrataG interface [50]. STRUCTURE was run with and without admixture with prior information on geographic origin included. Run parameters including the range of K (1–10) values assessed, the number of replicates, as well as burning and number of iterations were set on default values. The outputs produced by STRUCTURE were introduced in the Structure Selector [51] and Clumpak [52] to evaluate the optimal value of K, which was K = 1. According to Li & Liu, the Evanno method could underestimate the results [51]. Structure Selector [51] uses four statistical methods (MEDMEDK, MEDMEAK, MAXMEDK and MAXMEAK) to calculate the K number more accurately. A DAPC and PCA were generated in the adegenet package [42]. The DAPC was run without a priori grouping of samples through alpha optimization (via the dapc() and a.score.optim() commands).

Results

Genetic diversity

The global expected heterozygosity (Hs) was 0.14 and the overall observed heterozygosity (Ho) was 0.13 (S1 Fig). No statistical differences among the heterozygosity values (p = 0.493), discarding excess or deficit of heterozygotes were observed. The allelic richness was similar across all three sampling areas, 1.77 for CSL-MEX, 1.78 for GAL-ECU, and 1.75 for STR-ECU (Table 1). The AMOVA results show that the greatest variability occurs within individuals, with a mean value of 96.3%. Variation among individuals within sampling areas was 3.6%, while only 0.1% of the variation was explained by differences among sampling areas (S2 Fig).

Table 1

Summary of genetic diversity statistics for each sampling area of yellowfin tuna.

Locations	Ar	Ho	Hs
CSL-MEX	1.77	0.11 (0.02)	0.11
GAL-ECU	1.78	0.11 (0.01)	0.15
STR-ECU	1.75	0.10 (0.01)	0.10
Total		0.11 (0.02)	0.11

Ar: allelic richness, Ho: observed heterozygosity, Hs: expected heterozygosity. Standard deviation in parentheses.

Genetic structure

Limited genetic differentiation (F = 0.00328) was observed across the 164 individuals. Pairwise F analyses identified similarly low genetic differentiation between the three sampling areas, with the greatest differentiation observed between CSL-MEX and STR-ECU (FST = 0.0005 followed by CSL-MEX and GAL-ECU (F = 0.0004) and the lowest between GAL-ECU and STR-ECU (F = 0.0003) (Table 2). No isolation by distance was observed (Mantel statistic r: 0.9452, p = 0.16). Low allele frequency difference (AFD) values observed for pairwise comparisons of sampling locations ranged from 0.014 to 0.016 and pairwise genetics differentiation analyses indicated low levels of genetic differentiation among the three different sampling areas (Table 2). The inbreeding coefficient was 0.0609. The visualization of broad-scale population structure using a DAPC with all loci did not resolve any structure among individuals from CSL-MEX, GAL-ECU, and STR-ECU (Fig 2). Similarly, no divergence between sampling sites and no variation within sites was observed in the PCA (Fig 2). This lack of separation of individuals from the three sampling areas based on their genetic makeup is particularly evident when the individual density distribution of the first retained principal components from the discriminant function are examined (Fig 2).

Table 2

Pairwise genetic differentiation analyses for the yellowfin tuna among sampling areas.

		CSL-MEX	GAL-ECU	STR-ECU
F st	CSL-MEX (n = 35)	-	0.00006	0.00001
	GAL-ECU (n = 88)	0.0004	-	0.00001
	STR-ECU (n = 41)	0.0005	0.0003	-

Fig 2

(A) Genetic clusters of yellowfin tuna using analysis of Principal Components (PC) identified by Adegenet (B) Discriminant Analysis of Principal Components of yellowfin tuna within the Eastern Pacific Ocean, using 18,011 markers grouping samples by sampling site. CSL-MEX: Cabo San Lucas, México (n = 35), GAL-ECU: Galápagos Islands (n = 88), Ecuador, STR-ECU: Santa Elena, Ecuador (n = 41).

Sample sizes per sampling area are given in parentheses after the sampling area name; values above the diagonal are p-values; values below are F-values.

Stock clustering analyses or structure are largely consistent with those of DAPC and PCA, showing a lack of structure among individuals in each group, and any differences between groups are minimal (S4 Fig). All Structure Selector analyses suggest a single cluster (S3 Fig). Similarly, the absence of groups were observed in the Clumpak cluster analyses, suggesting no structure in the population (S4 Fig).

Discussion

A lack of genetic structure was observed for yellowfin tuna among the three sampling areas and no genetic difference by isolation was observed across spatial distances of 1,000–4,000 km, suggesting that yellowfin tuna within the EPO should currently be managed as a single stock. Pairwise analyses identified lowest genetic differences among individuals sampled at the two closest sampling sites, the Galápagos Islands and Santa Elena, Ecuador, separated by approximately 1,050 km. However, the greatest genetic differentiation was not observed between the two most distal locations (4,372 km) as might be expected if migratory distances played a role in the species differentiation, but rather among individuals from southern México and the Galápagos Islands (3,223 km). Structure analyses support the limited genetic differentiation observed among individuals from the different locations; low F scores and AFD values suggest high levels of genetic homogeneity among yellowfin tuna from the three sampling sites. This is further supported by the Mantel test analysis and the low inbreeding coefficient. The absence of deficit or excess of heterozygotes allowed us to discard the Wahlund effect in the population [53]. The lack of observed genetic differentiation in yellowfin tuna populations within the EPO region, conflicts with a previous study using microsatellite markers [11]. Our findings suggest there is insufficient genetic population structure in this region of the EPO to support multiple yellowfin tuna stocks, and therefore yellowfin tuna within the EPO should be managed as a single transboundary fishery, in accordance with the IATTC. However, future studies that increase the sample size, such as more sampling locations, or utilize more advance genomic techniques (e.g., whole genome sequencing and close kin mark recapture) may more clearly define the limited population structure observed in this study and may challenge these findings.

Our results showed a similar allelic richness among the three locations. This result and the similar patterns found in the PCA of each region, support the idea of limited differentiation among the sampling areas. Nevertheless, allelic richness is only one measure of genetic diversity at a given time [54], so continuous monitoring is crucial for achieving good management plans. In 1961, Orange [55] described several areas of the Pacific Ocean near Central America, including the Revillagigedo archipelago in México, and Cocos Island in Costa Rica, as possible spawning grounds for yellowfin tuna. Recent observations reported within the past ten years also support these findings and have demonstrated high occurrences of yellowfin tuna larvae within this regional area of the EPO [56, 57]. These studies provide evidence for the limited genetic differentiation that was observed in yellowfin tuna across our sampling areas, despite being more than 2,000 km apart. According to a recent study, movements/migration patterns of the yellowfin tuna occurred between close areas and showed higher dispersion rates in the fish released in offshore areas compared to the individuals released close to the coast or around islands [58]. If the region around Cocos Island is a spawning area, native individuals could disperse both to northern and southern areas of the EPO, reaching all three sampling areas of this study. While there is potential for temporally displaced spawning aggregations to produce mixed stocks of yellowfin tuna in the EPO, additional spatially and temporally structured sampling programs will be required to further address this hypothesis.

In addition to overfishing and anthropogenic impacts (e.g. pollution, habitat disruption), tuna populations in the EPO face other important threats like oxygen-minimum zones, which occur at elevated frequencies in this region [59], and have expanded in size over the past 20 years [60]. Likewise, the rise of water temperatures, as a consequence of anthropogenically driven climate change, could also reduce the availability of oxygen for marine organisms like tuna [61]. Due to their high oxygen demand, low dissolved oxygen might be a limiting factor for tunas, restricting their movements to shallow waters [62]. Moreover, the historical warming of the ocean has had a negative effect on marine fisheries production [63] causing changes in species regional distributions and, therefore, affecting the long-term availability of the fisheries [64]. These alterations could influence the migration patterns of yellowfin tuna, disrupting the gene flow, which ultimately could lead to the fragmentation of the population [65-67]. It will be important to track any changes in genetic structure of yellowfin tuna in the EPO, and other commercial fisheries, as adaptive management related to these changes may be required. The limited genetic structure observed in this fishery should be considered as a metric or indicator within fisheries management policies in this area, especially considering that this population is declining [2].

Several countries in the region, especially Ecuador and México, profit from the intense fishing of this stock and, therefore, regional management and conservation strategies are required for this shared resource. Hillborn et al. [68] provide evidence that intense management measures generally result in stock recovering or, at least, staying close to management target levels. Management measures include: adopting total allowable catch limits, establishing seasonal closures, and, setting catch limits for the longline fishery and effort limits for the purse seine fishery [69]. Identify management actions that have worked in other fisheries, particularly tuna fisheries, and the harmonization of fishing regulations for yellowfin tuna throughout the EPO region could promote the sustainability of the fishery. In 2017, IATTC [70] established several management measures that included limiting the number of active fish aggregation devices (FADs) per vessel, a 72-day closure for purse seine vessels greater than 182 tons carrying capacity, and a seasonal closure of the purse seine fishery in “El Corralito”, towards west of the Galápagos Islands, where three species of tunas (yellowfin tuna, bigeye, and skipjack) converge. These measures were planned to be applied between 2022 and 2024 and their overall effect is yet to be evaluated [71]. Our results can be used as steppingstone for the design and implementation of national and regional management policies rigorously complied by all members of IATTC, which should be linked to conservation efforts like the recent expansion of the Galápagos Marine Reserve and the establishment of the Galápagos-Cocos Swimway, a marine corridor that protect species as they migrate between protected areas. Boerder et al. [72] implied that the Galápagos Marine Reserve has a net positive effect on pelagic fisheries associated to the archipelago, emphasizing the importance of large-scale marine protected areas as both fisheries management and biodiversity conservation tools [73]. Additionally, Orange [55] suggested that the area around the Galápagos Islands is a spawning area for yellowfin tuna, which should be investigated further. It is critical to continue establishing management measures focused on the Galápagos Islands, given their importance for biodiversity and the fisheries, using the most accurate data to avoid the fragmentation of yellowfin tuna populations.

To inform adaptive management, it is essential to perform NGS population studies of the yellowfin tuna fishery, in combination with ecological and reproductive biology studies. Moreover, it is important to design management policies that do not only focus on controlling fishing activity, but also manage other threats, such as climate change, illegal and unrecorded catches, and pollution, as part of an ecosystem-based management approach. Despite the existence of international tuna regulation agencies, it is critical that governments and local entities that promote the research collaborate in the development of common regulations that guarantee the survival of this species. Understanding key aspects of the status of the species (e.g. habitat range, abundance, and number of stocks) is critical for stakeholders to develop and establish appropriate management measures that promote the biological sustainability of the stocks, which ultimately guarantees the survival of the species and food security in the region [68].

SNP filtering summary using radiator (Gosselin et al. 2020) detailing the thresholds for each filter parameter.

(TIF)

Click here for additional data file.

Genome-wide individual’s mean observed heterozygosity.

(TIF)

Click here for additional data file.

Analysis of molecular variance (AMOVA) of yellowfin tuna populations amongst basins nested within sampling areas.

(TIF)

Click here for additional data file.

Structure selector analyses.

(PDF)

Click here for additional data file.

Clumpak cluster analyses.

(PDF)

Click here for additional data file.

Metadata explanation of the SNPs single row file.

(XLSX)

Click here for additional data file.

SNPs single row file.

(RAR)

Click here for additional data file.