Genetic divergence at species boundaries of the dolphinfish (Coryphaena hippurus) in the Tropical Eastern Pacific

Sample collection and DNA isolation

Muscle tissue samples from dolphinfish individuals were opportunistically collected from 2003 to 2006, from artisanal fishing boats operating in the TEP, including the Gulf of California (Table 1; Fig. 1). Samples collected in Mexico included the Pacific coast of Baja California Sur (Bahía Magdalena (BM), Punta Lobos (PL), and Cabo San Lucas (CSL)); within the Gulf of California (Guaymas (GY)); and in the eastern central Pacific (Mazatlán (MAZ), and Chiapas (CH)). In South America, sample locations comprised Ecuador (EC) and Peru (PE). We also included a sample from an oceanic area (OC) located at coordinates 11.886755 N, 122.880087 W (Fig. 1). Total genomic DNA was isolated using the proteinase K-lysis-buffer protocol (Laird et al., 1991) or the Wizard Genomic kit (Promega Cat. No. A1125); DNA was rehydrated in 50–100 μl of TE buffer.

Microsatellite genotyping and quality control

Nineteen microsatellites were genotyped using primers for five loci (Chi002, Chi008, Chi008A, Chi023, and Chi037) designed by Robert Chapmann et al. 2005 (personal communication; GenBank Accession Numbers: AY135025–AY135028, and AY189832), 13 described by Bayona-Vásquez, Díaz-Jaimes & Uribe-Alcocer (2015), and one tetra-nucleotide locus (Chi284) derived from Next Generation Sequencing by Bayona-Vásquez, Díaz-Jaimes & Uribe-Alcocer (2015): with repetitive motive TTCC7: forward primer (5′ - 3′): TGTGGGAGAGGTCATTGCC; reverse primer (5′- 3′): AGAGGCAAGAGTGATGGTGC) (available at 10.5281/zenodo.6787993).

PCR amplification reactions were performed in 10 μl containing 10–100 ng of DNA, 1X PCR buffer (20 mM Tris-HCl, pH 8.4, 50 mM KCl), 0.2 mM of dNTPs, 0.2 μM of each primer, 2.5 mM of MgCl₂, and 0.25 U of Taq DNA polymerase. All forward primers had an M13 tail incorporated (5′-GTAAAACGACGGCCAGT-3′), and a third M13 primer labeled with a fluorophore (6-FAM, VIC, NED or PET) was included in the reaction to bind to the forward primers. Forward, reverse, and florescent-M13 primers were used in 1:10:10 proportions, respectively, and were added to the reaction mix according to the protocol described by Boutin-Ganache et al. (2001). The PCR conditions consisted of 5 min at 95 °C for denaturation, followed by 25 cycles of 30 s at 94 °C, aligning at 56 °C for 45 s, and extension at 72 °C for 45 s; followed by eight additional cycles at 94 °C for 30 s, 50 °C for 45 s, 72 °C for 45 s, and a final extension step at 72 °C for 10 min. Reactions were multiplexed for fragment analysis in an ABI-Prism 3100® sequencer. Genotypes were assigned using GeneMapper® software (Cat. No. 4475073; Thermo Fisher Scientific, Waltham, MA, USA) using GeneScan™ 500 LIZ™ (Cat. No. 4322682; Thermo Fisher Scientific, Waltham, MA, USA) as a size standard.

The existence of null alleles in the microsatellite data was evaluated with MICROCHECKER 2.2.3 (Van Oosterhout et al., 2004) using 1 × 10⁴ randomizations and a confidence interval (CI) of 95%. Departures from Hardy-Weinberg equilibrium (HWE) for each locus were examined by the exact test method. We also tested for linkage disequilibrium (LD) between pairs of loci. For both tests we used GENEPOP 4.7.5 software (Raymond & Rousset, 1995; Rousset, 2008), implementing 1 × 10⁴ dememorizations and iterations with 500 batches. For HWE and LD we applied a Bonferroni correction with the standard procedure (Miller, 1980; Lessios, 1992) using an initial critical value of 0.05. The final data set was used to assess the genetic diversity and structure.

Mitochondrial DNA sequencing

We amplified fragments of mitochondrial gene ND1 and cytochrome B (CYTB) using the primers reported by Díaz-Jaimes et al. (2006) and Hyde et al. (2005), respectively. PCR reactions for sequencing were carried out in 50 μl containing 10–100 ng of DNA, 1X amplification buffer, 10 mM TRIS-HCl (pH 8.4), 50 mM KCl, 1.5 mM MgCl₂, 0.2 mM of each dNTP, 0.1 mM of each primer and 2.5 units of Platinum® Taq DNA polymerase (Cat. No. 10966018; Thermo Fisher Scientific, Waltham, MA, USA). PCR amplifications for ND1 consisted of 95 °C for 5 min followed by 35 cycles of 1 min at 95 °C for denaturation, 1 min at 58 °C for annealing, and a final extension at 65 °C for 3 min, with a final extension at 72 °C for 10 min. Amplicons were purified and sequenced on an ABI 3730xl® automated sequencer by Macrogen Inc. (Seoul, South Korea). The thermal cycling profile for CYTB consisted of an initial denaturation step at 95 °C for 2 min, followed by 35 cycles at 95 °C for 1 min, 52 °C for 1 min, and 65 °C for 1 min, with a final extension of 3 min at 65 °C. Mitochondrial sequences were edited and aligned using ClustalX 1.8 (Thompson et al., 1997). The complete ND1 data sets from Díaz-Jaimes et al. (2006) and Díaz-Jaimes et al. (2010) from the eastern Pacific were reanalyzed, and new sequences from Bahía Magdalena, Punta Lobos, and the Oceanic area were added. Sequences of CYTB reported in this study are available from GenBank (accession numbers: MZ725384–MZ725448). The two loci were analyzed independently because of inconsistency in the amplification of many individuals. Because of this, we could not concatenate the sequences from the two loci.

Data analyses

Gene diversity

High levels of genetic diversity were observed in microsatellite loci across all studied locations. We obtained 218 alleles with a mean of 15.57 ± 6.3 alleles per locus. The total number of alleles (Na) ranged from 75 to 150 (Table 1), and all localities except Chiapas (CH) presented private alleles (Table 1). Chiapas (CH) exhibited the highest mean observed heterozygosity (H_o), which over all localities ranged from 0.64 to 0.80 (Fig. 2). The mean gene diversity (Hs) across localities ranged from 0.66 to 0.81. Mean allelic richness (Ar) ranged from 5.2 to 10.3, averaging 6.4. Although Chiapas (CH) had the highest H_o and Ar values, H_o, Hs, and Ar exhibited wide variability among sampling locations (Fig. 2); therefore, it was unclear whether Chiapas was the most diverse location. Most sites showed positive F_IS values, however, four sites (PL, GY, MAZ, and EC) exhibited a range of plausible values that suggested some degree of heterozygote deficiency, especially for Guaymas (GY) (Fig. 2).

For the mitochondrial analyses, we used the 750 base pair (bp) mtDNA-ND1 sequences from 364 individuals from the Díaz-Jaimes et al. (2006 and 2010) study and added 96 new individuals from the TEP to this dataset, resulting in a total of 117 haplotypes and 106 polymorphic loci. We also obtained a 506 pb fragment from the CYTB locus from 213 individuals, resulting in a total of 65 haplotypes and 72 polymorphic loci. Haplotype (h) diversity was higher than nucleotide diversity (π) for both fragments; h ranged from 0.81 (PL) to 0.96 (CH) for ND1 and from 0.85 to 0.94 for CYTB; π ranged from 0.003 (PL) to 0.005 (CH) for ND1 and from 0.004 to 0.005 for CYTB. Bahía Magdalena (BM) exhibited the lowest number of haplotypes for both fragments (n_h = 12 and 8), whereas Guaymas (n_h = 33) and the Oceanic sample (n_h = 21) had the highest estimations for ND1 and CYTB, respectively (Table 1). The 65 haplotypes for mtDNA-CYTB are available in GenBank with accession numbers: MZ725384–MZ725448.

Genetic structure

For microsatellite data, we were unable to detect a significant correlation between geographic and genetic distances (r = −0.16, P = 0.66). The mean Jost’s D between localities (0.06, P = 0.018) suggested low allelic differentiation, whereas values of G_ST (0.03, P = 0.0) indicated very low and significant differences in allele frequencies across sample sites. Nonetheless, comparisons involving Bahía Magdalena (BM), Peru (PE), and Ecuador (EC) (Table 2) exhibited the highest levels of allelic differentiation. In contrast, none of the mitochondrial pairwise φ_ST results were significant after Bonferroni corrections, with values ranging from 0 to 0.037 for ND1 and 0 to 0.014 for CYTB.

The cross-validation criterion from TESS3 for microsatellite data suggested an optimal K value of 6–7. However, as the variance notably increased for K > 4 (Fig. S1A) and cross-validation scores decrease noticeably with K = 4, we considered K = 4 to be the best value to describe the genetic structure pattern of C. hippurus (Fig. 1B). Nonetheless, both K = 2 and K = 3 were also used (Fig. S1B). The spatial pattern of genetic variation obtained by TESS3 suggested that Bahía Magdalena (BM) was an isolated population, with a second well-defined cluster consisting of Punta Lobos (PL), Cabo San Lucas (CSL), Guaymas (GY), Mazatlán (MAZ), and the Oceanic samples (OC) (Fig. 1B). Further south, TESS3 grouped Chiapas (CH) with Ecuador (EC), and Peru (PE) as a separate distinct population (Fig. 1B). These same four well-defined genetic clusters were corroborated by the DAPC analysis (Fig. S1C, Figs. 1C and 1D). The DAPC results showed almost no overlap and supported Bahía Magdalena (BM) and Peru (PE) as distinct populations despite, their geographic proximity to other sampled locations (Figs. 1C and 1D).

The Ln P(D) values from the STRUCTURE analysis increased when increasing K from K = 1 to K = 2, and reached a plateau with K = 3 (Fig. S1D). The ΔK method suggested K = 2 as the highest level of genetic partitioning in C. hippurus samples (Fig. S1D). Because Ln P(D) values differed slightly between K = 2 and K = 3, we chose K = 2 as the optimal number of clusters given the genetic data. The population structure inferred by the Bayesian algorithm using these values for K, assigned individuals to one of two clusters with high probability values (Fig. 1B). Individuals with high probability of membership to the first cluster (black bars) included samples from Bahía Magdalena (BM), Punta Lobos (PL), Chiapas (CH), and Ecuador (EC). The second cluster (pink bars) contained Cabo San Lucas (CSL), Guaymas (GY), Mazatlán (MAZ), the Oceanic sample (OC), and Peru (PE). Although the three analyses inferred slightly different population structures, all three consistently reflected a genetic divergence at the TEP’s northern and southern transitional areas, which correspond with the boundaries of the species’ distribution range (Fig. 1).

The hierarchical AMOVA grouping by biogeographic provinces using the microsatellite results revealed that ~93% of the genetic variation was explained by differences within populations (Table 3). However, this analysis also indicated a low but significant structure among biogeographic provinces (F_CT = 0.03; P = 0.004). This contrasts with the lack of significant differences of groupings based on dominant oceanographic processes and coastal regions (see Methods) (F_CT = 0.02; P = 0.11). It is noteworthy that the clustering results from TESS3 coincided fairly well with the biogeographic provinces proposed for the TEP, suggesting the influence of differences in environmental conditions among provinces in shaping the genetic structure.

Analyses of the genetic relationships between mitochondrial haplotypes from the two loci revealed two main phylogroups and no evidence of correlations between the haplotype differences within the network and different frequencies among the sampled populations (Figs. 3A and 4A). In both haplotype networks, two dominant haplotypes were found, with similar abundances for each at all locations. Likewise, the Bayesian assignment analyses for the two mitochondrial loci suggested that the most probable number of clusters (K) given the data was K = 2. However, no geographic structuring was detected for either of the mitochondrial fragments (Fig. S2).

Barriers to gene flow and recent migration rates

Using microsatellite data, BARRIER analysis identified genetic boundaries with good statistical support (>65%; Fig. S3), which coincided with the population groups resolved by TESS3 and DAPC. Interestingly, BARRIER also identified a sharp genetic discontinuity surrounding Bahía Magdalena (BM) and Peru (PE). Nonetheless, contrary to the other analyses, BARRIER further subdivided Chiapas (CH) and Ecuador (EC) as two isolated populations. Migration rates from BAYESASS largely agreed with the barriers to gene flow. The migration patterns were asymmetric, with overall low estimates of migration rate (≤0.01), which ranged from 0.2765 to 0.0056 (Table 4). We detected that populations from Guaymas (GY) and the Oceanic sample (OC) exhibited the highest values and thus may connect numerous populations in southern regions of the TEP (e.g., GY to PL, MAZ, and CH; and OC to CSL and PE, respectively).

Effective population size

Estimates of contemporary Ne suggested, in general, low values of Ne for C. hippurus within the TEP. Estimates ranged from 77.9 to 496.4 (Table 5). The genetic groups composed of Chiapas and Ecuador (genetic group C in Table 5) and the one corresponding to Peru (genetic group D) had the lowest values of Ne.

Demographic history

Both mitochondrial and nuclear DNA genetic markers provide evidence of ancient and recent changes in effective population size. For instance, based on the IAM and TPM mutation models, Wilcoxon rank tests suggested a recent contraction in population size for two genetic groups (Table 5). The first one largely corresponded to localities from Mexico (genetic group B in Table 5); and the second included Peru (genetic group D).

The mitochondrial sequences showed a strong signal of population expansion. The distribution of mismatches for the whole dolphinfish population in the eastern Pacific for both fragments was clearly unimodal (Figs. 3B and 4B). For mtDNA-ND1, Harpending’s raggedness index (Hri) and Sum of Squared deviations (SSD) had low and non-significant values (Hri = 0.026, P = 0.057; SSD = 0.004, p = 0.362), supporting the occurrence of a sudden demographic expansion and past fluctuations in population size. Tajima’s D and Fu’s Fs both suggested departures from neutrality for ND1 (Tajima’s D = −2.44, P = 0.0; Fs = −25.96, P = 0.0, respectively) and CYTB (Tajima’s D = −2.46, P = 0.0; Fs = −26.85, P = 0.0, respectively). Demographic Bayesian reconstructions from the mtDNA-ND1 gene suggested a growth of the dolphinfish population that started approximately 126,000 years before present (Fig. 3C), with the population size scaled by a generational time increasing by approximately threefold from 0.73 (CI [10.20–0.27]) to 12.63 (CI [33.85–4.69]). Consistent with this, CYTB suggested a population growth approximately 92,838 years before present, with an increase in population size from 0.50 (CI [5.41–0.18]) to 8.88 (CI [36.62–2.68]) (Fig. 4C).

Genetic differentiation of dolphinfish in the TEP

Without apparent barriers to dispersal for marine species having a pelagic larval stage, we would expect genetically homogeneous populations. However, the low differentiation that has been reported does not necessarily imply high levels of gene flow among populations (Whitlock & Mccauley, 1998), even for species with high dispersal capabilities or those with an extended larval stage (Handal et al., 2020).

An interesting result from this study is the remarkable genetic differentiation between the two most distant locations (Bahía Magdalena (BM) and Peru (PE), Fig. 1) at the northern and southern extremes of the species range boundaries, which separate the tropical and subtropical waters in the TEP. These transitional areas are associated with seasonally complex oceanographic features, such as a convergence of different water masses and strong gradients with seasonal variation in physical (e.g., temperature and salinity) and chemical properties (e.g., dissolved oxygen, and nutrients) (Lluch-Cota et al., 2019). Considering that dolphinfish distribution is directly dependent on sea surface temperature (Zúñiga-Flores, Ortega-García & Klett-Traulsen, 2008), it is possible that transitional areas such as these play a significant role in originating cycles of range expansions-contractions in populations inhabiting them and thus, resulting in genetic differentiation. There is evidence of dolphinfish range expansion during warm ENSO phases, when an abnormal increase of sea surface temperature causes a poleward extension to these fresh and warm water areas usually characterized by relatively cooler temperatures. Studies on the spatio-temporal distribution of dolphinfish within the TEP have suggested that this species expands its habitat in response to increased ocean temperature (Norton, 1999; Torrejón-Magallanes, Grados & Lau-Medrano, 2019), boosting population genetic differentiation by modifying their distribution range during ENSO events. During warm phases, the distribution expands, while during cold events, it contracts. The migratory behavior of dolphinfish also reinforces this hypothesis. There are spatial changes in abundance that reflect a preference for warm temperatures (Zúñiga-Flores, Ortega-García & Klett-Traulsen, 2008; Farrell et al., 2014; Brodie, Hobday & Smith, 2015; Torrejón-Magallanes, Grados & Lau-Medrano, 2019).

The northern boundary of the TEP is a region with a seasonal convergence of three different water masses: Gulf of California, tropical surface, and transitional tropical. Additionally, the equatorial front at the TEP’s southern boundary is a convergence zone that separates the salty and cold equatorial surface water from the fresh and warm tropical surface water. Transitional zones are notably productive (Hill et al., 1998). For both transitional areas, the California and Humboldt currents are driven by winds that blow along the shore toward the equator for a substantial part of the year (Smith, 1995), generating upwelling events and driving nutrient-rich waters from deeper layers to superficial waters, promoting dense phytoplankton blooms that are the base of ocean food webs (Chan, 2019). This high productivity generates feeding grounds for many pelagic species that aggregate for various functions, including spawning (Dingle, 2014). The dolphinfish is an opportunistic feeder that consumes easily caught prey and shows seasonal changes in diet (Tripp-Valdez, Galván-Magaña & Ortega-García, 2010). The TEP boundaries are attractive for dolphinfish because of the abundance of prey species (Olson, 2001; Acha et al., 2004) and floating debris (e.g., driftwood, detritus, and other flotsam), with which dolphinfish are strongly associated, showing both homing and high levels of site fidelity (Girard et al., 2007; Whitney et al., 2016; Perle et al., 2020). The specific behavior in which dolphinfish tend to stay and return after foraging excursion, and their ability to orientate towards a fish aggregating device (Girard et al., 2007), support restricted gene flow within the boundaries of the study area and isolation of populations.

These seasonal genetic shifts are also supported by the barriers to gene flow identified in the region and the migration rates found in this study. The suggested differences between Chiapas with northernmost locations (Fig. 1) may be related to oceanographic conditions in the Gulf of Tehuantepec and oceanic circulation. The Gulf of Tehuantepec is influenced by intense seasonal wind events during winter, producing strong upwelling, low sea-surface temperatures, mesoscale eddies, and shallow thermocline (González-Silvera et al., 2004; Fiedler & Lavín, 2017). In addition, wind stress across Central America gives rise the poleward Costa Rica current which flows to the north reaching the Gulf of Tehuantepec where it is interrupted and restricts larval dispersal northward. Therefore, the oceanography of the Gulf of Tehuantepec may influence the genetic interchange of the Chiapas population in relation to those from northern Mexico. Although mark-recapture studies have shown that dolphinfish tend to remain residents near the areas where they were originally tagged (Kingsford & Defries, 1999; Merten, Appeldoorn & Hammond, 2014; Perle et al., 2020), some individuals have also performed longer distance migrations from northern to southern Mexico in response to seasonal temperature changes (Perle et al., 2020). This agrees with the overall low migration rates found and the asymmetric connectivity between distant regions in this study (e.g., from Guaymas to Chiapas), which could be facilitated by possible annual migration route during winter (Perle et al., 2020). Further, the geographic distribution of genetic variation found in this study indicated that the transoceanic connectivity between individuals in the Oceanic (OC) and Peru (PE) sites (depicted by K = 2 from STRUCTURE and the migration rate analysis) reflects that the deep-water zone does not constitute insurmountable barriers separating populations (Robertson & Cramer, 2009). Hence, larval dispersal promoted by oceanic currents might still eventually occur, meaning that dolphinfish populations are not entirely isolated.

Finally, the highly dynamic oceanographic processes observed in the transitional areas of the TEP may also contribute to the genetic variation observed. The mixing of water masses that differ in temperature and salinity increases the mesoscale dynamics, promoting small and long-lived energetic eddies and/or permanent upwellings (Montecino & Lange, 2009; Pantoja et al., 2012). Oceanographic conditions at transitional areas are characterized by high seasonal variations in temperature, salinity and primary productivity. North of the equator, off the cost of Mexico, the transitional area between tropical and subtropical waters known as transitional tropical water (TTW), is part of Eastern Tropical-Subtropical Convergence, where the Gulf of California water and the Tropical Surface Waters also converge (Fiedler & Talley, 2006). Strong gradients in temperature and oxygen have been reported in the TTW, originating from the thermocline and oxycline sinking in the frontal area (30 m) towards the California Current waters (90 m) (Cruz-Hernández et al., 2019). South of the equator, the equatorial cold tongue converges with the cooler water of the Peru Current, generating a front whose seasonal northern and southern limits depend on the intensity of the Peru Current (Fiedler & Lavín, 2017). The equatorial cold tongue clearly marks the limits of the tropical surface water at Central America and Mexican latitudes. This may prevent the dispersal of larvae and limit adult movements, contributing to seasonal genetic subdivisions observed in the transitional areas.

Conversely, the merging of these water masses, the strong environmental gradients, and their interaction with feeding areas and other fine-scale mechanisms of the transitional areas could promote local adaptations for individuals inhabiting them. This hypothesis agrees with previous findings of genetic differentiation in areas that exhibit environmental gradients or patchiness for other species such as Clupea harengus in the Baltic Sea; Amphiprion bicinctus in the Red Sea; and Merluccius productus in the Eastern Pacific (Jørgensen et al., 2005; Nanninga et al., 2014; García-De León et al., 2018). Furthermore, our AMOVA analysis reflects significant differentiation when considering the environmental differences that distinguish the biogeographic provinces. However, additional studies are required to demonstrate a link between the heterogeneous environmental conditions in these areas and the genetic differentiation found to verify that it is indeed environmental factors that drive genetic structure.

Genetic diversity and effective population size

Historically, dolphinfish is one of the main fishery resources of many countries within the TEP, especially for Ecuador and Peru (Aires-da-Silva et al., 2016). Estimates of heterozygosity and the number of alleles in dolphinfish were in agreement with other studies of large pelagic fishes, such as Thunnus thynnus thynnus, T. albacares, T. orientalis, and Xiphias gladius (Carlsson et al., 2004; Díaz-Jaimes & Uribe-Alcocer, 2006; Nomura et al., 2014; Aguila et al., 2015; Riguik et al., 2020), as well as in a comparative study that highlighted the considerably high heterozygosity levels of marine fishes (H_O: 0.67–0.70; Martinez, Willoughby & Christie, 2018).

Our estimates of contemporary Ne suggest that dolphinfish populations are below the minimum threshold of 1,000 for retaining their evolutionary potential in the long term (Frankham, Bradshaw & Brook, 2014), raising concerns that the survival of dolphinfish populations in the TEP could be in jeopardy. Deviations from the Hardy-Weinberg equilibrium were detected in some populations due to a large number of loci with heterozygote deficit, leading to positive values of inbreeding (F_IS). This pattern is unlikely to be attributed to a high frequency of null alleles because loci were extensively tested for genotyping artifacts that could produce null alleles. Loci Chi389, Chi853, and Chi967 showed null allele frequencies of about 25% in only one population (one of each); however, loci had a minor effect on G_ST and Jost’s D estimates. Also, no evidence for significant allele drop-out was found using MICROCHECKER. Likewise, a Wahlund effect also seems unlikely because no cryptic genetic structure was detected using clustering analysis. Therefore, we suggest that the positive values of F_IS and the low effective population size constitute a warning that the dolphinfish populations are vulnerable to fishing pressure.

Low Ne estimations have been documented for other similar species such as X. gladius in the TEP. Estimates of Ne for X. gladius ranged from 100 (or lower) to 1,000 individuals and are comparable to those reported for C. hippurus in this study. The notably decreased Ne for C. hippurus parallels the remarkable decline in X. gladius fishery (Yüncü et al., 2020). Harvesting may increase mortality and cause deviations in drift-mutation equilibrium, which in species with short generational time may increase the loss of genetic variation over relatively few generations, which in turn is enhanced by the presence of genetic structure (Allendorf et al., 2008). Consequently, a divergent population subjected to prolonged mahi-mahi fishing practices, such as in Peru, may easily cause very low values of Ne. Although C. hippurus in Mexico is restricted to sport-recreational fishing within 50 nautical miles from the coastline (Diario Oficial de la Federación, 2013), it is common to find dolphinfish specimens within commercial landings recorded as by-catch with no official reporting, which prevents estimates of the total by-catch of this species in Mexico (Alejo-Plata, 2012). Dolphinfish show early sexual maturity, high fecundity, and an asynchronous spawning occurring in waters relatively close to the coast throughout the year in the tropics (Cheung, Lam & Pauly, 2008; Moltó et al., 2020). In the southern Gulf of California, dolphinfish reproduction occurs mainly during the warm months of summer-autumn (Zuñiga-Flores et al., 2011) whereas in southern Mexico reproduction occurs in May–July, and November–January (Alejo-Plata, Díaz-Jaimes & Salgado-Ugarte, 2011). In Costa Rica, spawning occurs in January and February (Campos et al., 1993). In Ecuador, dolphinfish spawn throughout the year, although the months of maximum recorded reproduction are October-December in the north, and from January to March in central and southern waters (Zuñiga-Flores & Lavayen, 2014). In the Economic Exclusive Zone of Peru spawning records are during austral summer (Solano et al., 2015). These studies show that C. hippurus spawn mainly in tropical waters, with females using coastal habitats as feeding areas more frequently than males (Alejo-Plata, 2012). As a consequence, females are more vulnerable to be capture near to the coast during the spawning season. The resulting alteration of the sex ratio would potentially reduce the effective population size and thus increase the probability of mating with relatives. Thus, it is crucial to consider that dolphinfish may be prone to genetic drift and inbreeding, which would compromise their long-term evolutionary potential.

Mitochondrial patterns

Temporal population reductions and subsequent expansions influence genetic structures, especially for haploid and uniparental genomes. Population expansions to new areas are usually initiated by few individuals at the extremes of the distribution, followed by a rapid and large increase in the population size. We found signals of rapid population expansion of dolphinfish within the TEP, evidenced by the marked difference between haplotype (very high) and nucleotide (very low) diversities estimated for the mtDNA and the star-like phylogeny. Previous studies suggested that habitat reductions or bottlenecks are associated with low nucleotide diversity resulting from a limited number of shared sequences and numerous low-frequency haplotypes separated by few mutations (Díaz-Jaimes et al., 2006). This arrangement was found in the haplotype networks of both mitochondrial sequences and the observed pattern of nucleotide and haplotype diversity. In turn, the Bayesian Skyline plot, the unimodal mismatch distribution, and the rejection of the neutrality test clearly supported a scenario of population expansion after a period of low effective population size in our samples of dolphinfish. One possible explanation is that fragments of populations of dolphinfish when faced with unfavorable environmental conditions may survive in refugia. Upon restoration of favorable environmental conditions, these populations could rapidly rebound, colonize new sites and continue to grow over subsequent generations. The timing of expansion as estimated from the Bayesian Skyline plot was about 126,000 years before present for ND1 and 92,838 for CYTB; these periods coincide with the penultimate interglacial period during the late Pleistocene. This period was characterized by sea temperature two to three degrees warmer than the current values and a sea level five to six meters higher than present (reviewed by Bergoeing, 2017). This is not the first study that associates mitochondrial genetic variation in C. hippurus with population expansion after cooling episodes during the Pleistocene, and cycles of population contraction and expansion following changes in sea surface temperature and habitat availability in the Eastern Pacific (Díaz-Jaimes et al., 2006).

It has been previously hypothesized that a lack of genetic divergence among C. hippurus populations could be due to extensive region-wide gene flow (Díaz-Jaimes et al., 2006; Tripp-Valdez et al., 2010) or because genetic drift could not impact genetic variation due to large population size (Díaz-Jaimes et al., 2010). Even if oceanographic features such as eddies, gyres, and fronts with substantial environmental variability can limit the dispersal of larvae or more mobile life stages, a large effective population size counteracts the effects of genetic drift on genetic variation, blurring evidence of population subdivision (Díaz-Jaimes et al., 2010). Thus, we propose that the low genetic differentiation in microsatellite markers of C. hippurus indicates a recent genetic divergence with enhanced genetic drift in populations at the latitudinal extremes of the species range, rather than the previous postulation of high gene flow among all populations.