Genomic assessment reveals signal of adaptive selection in populations of the Spotted rose snapper Lutjanus guttatus from the Tropical Eastern Pacific

Sample collection

A total of 192 individuals were collected from 17 locations spanning the species shoreline range (Fig. 1, Table 1). Samples were collected and preserved following procedures described in (Torres-Hernández et al., 2022). Tissue samples were stored at the Laboratorio de Genética at the Instituto de Ciencias del Mar y Limnología, Universidad Nacional Autónoma de México, and in the Icthyological Collection of the Laboratorio de Biología Acuática, of Universidad Michoacana de San Nicolas de Hidalgo, México (SEMARNAT registration number CPUM-PEC-227-07-09). Permits for specimen collection were provided by the following institutions in each country: Mexico (PPF/DGOPA-035/15; SGPA-DGVS-02920/15 and F00.DRPBCPN.DIR.RBAR-100/2015-CONANP), El Salvador (MARN-AIMA-004-2013), Panama (SC/A-17-19), Costa Rica (007-2013-SINAC and R-056-2105-OT-CONAGEBIO), and Ecuador (013/2012 PNG; N° 21-2017-EXP-CM-2016-DNB/MA and MAAE-DBI-CM-2021-0152).

Library preparation

The RADcap protocol (Hoffberg et al., 2016) was used to obtain genomic data. For the first step, fifteen individuals from nine sampling sites across the L. guttatus distribution range were selected to design baits for capture sequences using the 3RAD protocol to obtain genomic data (Bayona-Vásquez et al., 2019). For the 3RAD library, DNA was digested using three restriction enzymes (BamHI, ClaI and MspI) which were ligated to adapters with internal iTru indexed primer (Glenn et al., 2019). Libraries were size-selected for 500 bp ± 10% using Pippin Prep and pooled in equimolar ratios and sequenced on a HiSeq PE150 at the Oklahoma Medical Research Foundation. The reads obtained were processed in Stacks v1.42 (Catchen et al., 2011; Catchen et al., 2013), and assembled de novo to generate a FASTA file catalog of the polymorphic loci that were present in at least five of the nine populations to design the SR-Snapper-Bait. We obtained 1973 RAD-loci from which were designed baits at Arbor Biosciences. Using custom filter criteria, we selected two baits—R1 and R2—for 1000 loci which served as the pseudo-reference genome.

For the second step, we generated 3RAD libraries for 192 individuals, which were normalized to 20 ng/uL, using 5 uL per sample for a final concentration of 100 ng/uL. The samples were digested with the enzymes described above and ligated to 96 internal index adapters. Two independent pools with 96 samples per capture were constructed. For each pool a single-cycle PCR was carried out using the Kapa HiFi HotStart kit. A degenerate primer (8N) was used to index individuals for downstream filtering of PCR duplicates. After cleaning the PCR products with a 1.2:1 volume ratio of SpeedBeads:DNA, three limited-cycle PCR replicates were performed to increase the concentration (for detailed description see Hoffberg et al., 2016). This final product was run on a 1.5% agarose gel to validate the library construction. The capture process was done following the manufacturer’s protocol, using temperatures between 62.5 °C and 65 °C, selecting baits with GC content between 30% and 60%. We kept one bait for Read1 and Read2, and those that passed the filters. This set of baits was synthesized as a custom RNA myBaits^® kit (for details of the RADcap protocol, see File S1).

Data processing

Reads of 150 bp were obtained for each pool. In Stacks v2.0 (Rochette, Rivera-Colón & Catchen, 2019) clone_filter and process_radtags modules were used to remove PCR duplicates, demultiplex samples, and trim the ligated adapters. This yielded final reads with the same length (139 bp). From this initial set of data, we retained only individuals that had at least 1000 reads and at least 4x of read depth (Hoffberg et al., 2016). Reads were mapped with the mem algorithm in BWA v 0.7.7 (Li & Durbin, 2009), using the designed baits as the reference (pseudo-genome). Finally, using the BAM files we ran the ref_ map.pl pipeline in stacks v2.0 for locus-building and SNP-calling. Finally, in Stacks v2.0, the populations module was used for retained loci presenting at least 50% of the individuals within a locality (r = 0.5), and with a maximum observed heterozygosity higher than 70% (h = 70) to exclude potential paralogues (Pedraza-Marrón et al., 2019). To filter the data we used VCFtools v0.1.15 (Danecek et al., 2011) to remove individuals with more than 20% missing data (an exception was made for localities with ≤ 5 samples, where we kept all individuals; Table S1). After the filter steps obtained a total dataset of 123 individuals which were used in the downstream analysis. From this final dataset, we removed data with allele frequencies less than 1%. We avoided local linkage disequilibrium by considering only one SNP every 1000 bases, using PLINK v. 1.9 (Purcell et al., 2007). To corroborate the absence of linkage disequilibrium, we calculated the Index of Association (ia) in R package Poppr (Kamvar, Tabima & Grünwald, 2014).

Outlier detection

We identified outlier loci using two different approaches. The first, BayeScan v.2.1 (Foll & Gaggiotti, 2008), is a Bayesian approach that identifies loci under selection using differences in allele frequencies. This method uses a linear regression to decompose F_ST coefficients into population-specific and locus-specific components. The BayeScan analysis was conducted using default settings, 20 pilot runs of 5,000 iterations each, burn-in of 50,000 iterations and sampling every 5,000 generations with a thinning interval of 10; individuals were grouped into localities that were considered “populations”. The second method, pcadapt was used to define loci under selection, this approach is based on Principal Component Analysis (PCA). Markers were considered to be under selection if they had a false recovery rate (FDR) correction of 0.01 in PCAdapt v.4.0.2 (Privé et al., 2020). We applied the default parameters, including α value of 0.1 and the first two PCs. All SNPs that were identified as outliers by both methods were separated, and for downstream analyses we considered three datasets: neutral loci (NL), outlier loci (OL), and all loci (combined [CL]). Sequences from the OL were blasted in the NCBI dataset (https://blast.ncbi.nlm.nih.gov/Blast.cgi) with the blastn option. Additionally, we scanned the OL dataset with the genome from L. erythropterus to evaluate if loci under selection were related with a protein (BioProjects number: PRJNA662638).

Estimates of genetic diversity

The functions basic.stats and allelic.richness from the R package hierfstat v.0.04.22 (Goudet & Jombart, 2015) were used to obtain population genetic statistics for the three datasets (CL, NL and OL). For each location we estimated the expected (H_e) and observed (H_O) heterozygosity, inbreeding coefficient (F_IS) and allelic richness (Ar). Confidence intervals for F_IS values were calculated with 1000 bootstraps in the boot.ppfis function on hierfstat. Paired F_ST were calculated following the Weir and Cockerham method (1984) in Arlequin version 3.5.2.2 (Excoffier & Lischer, 2010). The significance level (α) was adjusted by the Benjamini & Yekutieli (2001) correction (B-Y), as proposed by Narum (2006), by dividing the critical value of α by the sum of the numbers of tests.

Assessment of population structure

We conducted a spatial analysis of molecular variance (SAMOVA) in SAMOVA 2.0 (Dupanloup, Schneider & Excoffier, 2002) considering K values from 1 to 17 (we used as the maximum K the total number of localities as recommended by Pritchard, Stephens & Donnelly, 2000), and using 10, 000 iterations without any previous grouping. The groups that best explained the genetic structure were tested by AMOVA (Analysis of molecular variance), allowing 10% missing data in ARLEQUIN version 3.5.2.2 (Excoffier & Lischer, 2010). The Discriminant Analysis of Principal Components (DAPC) was made with adegenet 2.0.1 (Jombart & Ahmed, 2011). For all datasets, this analysis was carried out in two ways—without grouping localities and grouping localities based on the SAMOVA results. The find.cluster function was used to identify the K value. The best-supported number of clusters was identified by comparison of the Bayesian Information Criterion (BIC) for the different values of K. After running the DAPC, 70 PCs, which accounted for approximately 80% of the total variation in the data set, were retained using the dapc function. Finally, the DAPC scatterplots were obtained using the ggplot2 package (Wickham, 2016) in R. To further evaluate the population structure, we used a maximum likelihood approach in ADMIXTURE 1.3.0 (Alexander & Lange, 2011). This software calculates ancestry of individuals that are unrelated using a large SNP genotype dataset. We ran ADMIXTURE under the EM algorithm as optimization method, with 2000 for bootstrap, the cross-validation procedure was used to define the most suitable K-value, and the best value of the coancestry cluster was selected based on the lowest cross-validated error (CVE).

IBD, IBE, IBR

Isolation by distance (IBD) is widely used to explain differences among populations in species with a broad geographic range (Wright, 1978). To test the IBD hypothesis we performed a Mantel test using a correlation between F_ST/1-F_ST values and the geographic distance among localities (transformed to natural logarithm (LN)). Geographic distances were calculated with marmap v1.0.4 package (Pante & Simon-Bouhet, 2013) in R. This approach considers the bathymetry of the marine environment, which can accurately determine distance while considering vertical movements of marine organisms. The values of 5 to 120 m (Robertson & Allen, 2015), were used to constrain the least cost path (LCP) algorithm. This corresponds with the depth range of habitat. Besides, mantel test was carried out using the Euclidian geographic distance, because the current could have an effect in the dispersal of pelagic larval in Lutjanid as reported in Pedraza-Marrón et al. (2019).

To test isolation by environment (IBE) and isolation by resistance (IBR), we obtained 13 oceanographic variables from Bio-ORACLE (Tyberghein et al., 2012; Assis et al., 2017). The IBE was first assessed testing the collinearity between 12 variables considering values obtained from the sampling location. A Pearson’s correlation was carried out with the R package GGally (Schloerke et al., 2018), and variables strongly correlated (>0.8) were discarded. Five oceanographic variables (temperature, salinity, primary production, phytoplankton and chlorophyll) were retained to generate the environmental matrix between localities. Environmental distances were obtained with the Canberra distance method, using the “havegdist” function from the R package vegan (Oksanen et al., 2019). For the IBR the current velocity variable was used, layer was crop to the TEP polygon size and converted to raster format using the R Package Raster (Hijmans, 2020). Resistance pairwise distance matrix was generated with Circuitscape v. 4.0.5 (McRae, 2006). The circuitscape algorithm predicts patterns of gene flow from a heterogeneous landscape, considering possible pathways between a particular region and calculating the average cumulative resistance between localities sampled (McRae & Beier, 2007). Finally, a Mantel test was performed in vegan with the Spearman’s correlation using 9999 permutations. Genetic distances were used as dependent variable and the environment and resistance were considered as independent variables.

Spatial migration

We evaluated the spatial population structure as a function of potential barriers to gene flow using EEMS (estimated effective migration surface) software (Petkova, Novembre & Stephens, 2015). This approach allowed us to differentiate between regions with strong migration rates and regions where possible barriers to gene flow exist. This method employs a matrix of genetic distances and geo-referenced samples, and a distribution polygon to determine the migration surface. For the genetic distance, a dissimilarity matrix was calculated using a bed file with the bed2diffs module from EEMS (Petkova, Novembre & Stephens, 2015). The distribution polygon was drawn following the potential habitat of L. guttatus and obtained with the polyline method in the Google Maps API v3 tool (http://www.birdtheme.org/useful/v3tool.html). For the EEMS model we used three independent runs considering 200, 300, and 450 demes, a burn-in of 1,000,000 steps and MCMC length of 5,000,000 iterations. The convergence of runs and the spatial surface (maps of effective migration rate (m)) and effective diversity (q) was visualized using the Reemsplots function in R.

SNP calling and outlier detection

We obtained 24 million reads with the RADcap protocol. After filtering data and removing duplicates, 10 million paired reads were retained. The first evaluation of the quantity of reads per sample left 154 samples, which recovered 10814 SNPs with 39.64% missing data. This dataset was filtered for the missing data, generating a final dataset consisting of 123 individuals from 17 localities (Fig. 1, Table S1). From those 123 samples we retained 2003 bi-allelic SNPs in 942 loci (loci were constructed by Stacks using the baits as pseudo-genome, (as they do not have biology response, we used only to name the regions containing SNPs), which correspond to the combined loci dataset (CL; neutral and outlier). The two approaches used to detect loci putatively under selection obtained as a consensus panel of 145 candidate outlier SNPs, which were separated into the OL dataset, while another 1858 SNPs remained as neutral loci dataset (NL). Contrasting results obtained with the OL and NL datasets will be shown in the main manuscript, while results for CL are reported in File S2.

Genetic diversity

The ia indicate a lack of disequilibrium linkage in NL and OL (${\overline{\text{ṛ}}}_{\mathrm{d}}=0.0242$, p = 0.433 and ${\overline{\text{ṛ}}}_{\mathrm{d}}=0.0163$, p = 0.84 respectively) (Fig. S1) Values of F_IS for NL were almost all positive; only TSA (F_IS = −0.483) had a negative value (Table 1). For the OL dataset, all F_IS values were positive (ranging between 0.130 for TSA and 0.572 for CLB). Confidence intervals of F_IS were above zero (Fig. S2), suggesting some degree of heterozygote deficiency (Table 1). For LPA and SRO samples, it was not possible to calculate the F_IS (Table 1). The observed heterozygosity (H_O) values were similar between the NL dataset (0.047–0.081) and the OL (0.053–0.078) (Table 1). For expected heterozygosity (H_e), NL had lower values (0.047–0.127) than OL (0.033–0.175) (Table 1). Allele richness (Ar) for NL ranged from 1.05 to 1.12, and the OL ranged from 1.06 to 1.11 (Table 1).

Population genetic structure

For the NL dataset, the results evaluating different groups were non-significant and failed to recover any geographic group (results for two groups are showed in Table S2). However, when using the outlier loci (OL) and combined loci (CL) (File S2), results based on two groups showed higher and significant F_CT (0.4860; p = 0.001 and 0.0653; p = 0.001). The first group included all of the localities for the northern distribution of the species (SRO, LPA, LOR, TSA, SON, SIN, NAY, COL, MCH, PTO, GRO, OAX, SLV and CRI, hereafter known as Northern Group), and the second cluster recovered three localities from the Southern distribution (PAN, CLB and ECU, hereafter known as Southern Group) (Table S2).

Pairwise-sample F_ST estimates using the NL dataset were not significant for almost all the comparisons (Table 2). After adjusting the significance level for multiple testing (p = 0.015), only the comparison between SLV vs PAN was significant (F_ST = 0.031, p = 0.004). In contrast, the OL dataset presented higher estimations of F_ST, resulting in highly significant differences involving the three locations from the Southern group (Table 2): PAN (F_ST ranging from 0.306 to 0.498 with corresponding p-values below to 0.015), CLB (F_ST = 0.199–0.688, p < 0.015) and ECU (F_ST = 0.306–0.644, p < 0.015). None of the comparisons between these southern locations and OAX were significant, although they showed high pairwise F_ST; this is likely due to the fact that the OAX location consisted of a single individual.

The hierarchical AMOVAs tested the hypothesis of Northern and Southern Groups. AMOVA using the NL dataset failed to detect significant genetic differences for locations among groups (F_CT = 0.005; p = 0.553). In contrast, the AMOVA using the OL dataset showed significant genetic differences for variance among Northern and Southern groups (F_CT = 0.4860; p =  < 0.001) (Table 3).

DAPC analysis recovered the Northern Group formed with the localities SRO, LPA, LOR, TSA, SON, SIN, NAY, COL, MCH, PTO, GRO, OAX, SLV and CRI, while the Southern Group was formed with localities PAN, CLB, and ECU (Fig. 2). The Admixture analysis using NL recovered K = 2 (Fig. S3A); while using OL, K = 3 was obtained (Fig. S3B) as the number of clusters that best explained the genetic subdivision. Although the better explanation for NL was 2 clusters, the plot showed no evidence of population structure (Fig. 3A). For OL the better result was the formation of three clusters, nevertheless, the results show a clear separation into of two mains groups, Northern and Southern groups remain as the most differentiated (Fig. 3B).

IBD, IBE and IBR

The Mantel test to assess the isolation by distance model failed to detect a significant correlation between genetic and geographical distances using either dataset (NL: r = 0.0587; p = 0.278 and OL: r = 0.0970; p = 0.174; Fig. S4). Second Mantel test was carried out evaluating the IBD in each of the two main groups, results do not recover any pattern of relation between the geographic distance and the genetic variation (Fig. S5).

Mantel test for the IBE and IBR for the NL, failed to detect correlation between the environment variables tested (IBE) and currents velocity (IBR) (r = 0.1241; p = 0.1158 and r =  − 0.05345; p = 0.5898 respectively) (Figs. 4A and 4C). On the other hand, mantel test for OL recovered a correlation between the genetic distance versus the oceanographic environment variables and currents velocity (r = 0.3403; p = 0.0013 and r = 0.5207; p = 0.002) (Figs. 4B and 4D).

Spatial migration

Evaluation of the spatial genetic structure using EEMS with the NL dataset recovered patches of migration throughout the distribution area of L. guttatus in the Gulf of California (LPA and LOR) and between NAY and COL (Fig. 5A), and a migration corridor was found off the coast of CRI, which flows northward and could be associated with the Costa Rica Coastal Current (CRCC) (Fig. 5A). There was no clear barrier to migration between the Northern and Southern Groups (Fig. 5A), although migration between CRI and PAN was restricted. A barrier for gene flow was inferred off the Mexican coast (Fig. 5A).

Results for the EEMS using the OL dataset contrast with the NL, showing a possible zone that generate local adaptation between the Northern and the Southern Groups which concur with the heterogeneous environment found in the Panama Bight Gyre (Fig. 5B, Fig. S6). Like the NL results, there was a barrier between the Gulf of California localities and the Mexican Pacific shoreline (Fig. 5B). Corridors for migration were found between LOR and LPA in the Gulf of California, NAY and COL on the Mexican Pacific coast, and among PAN, CLB and ECU (Southern group). As in the NL dataset, a corridor associated with the CRCC that flows northward was inferred with the EEMS (Fig. 5B).

For the genetic diversity inferred with EEMS was similar between the NL and OL datasets. High diversity rates were observed in SRO, LOR, and LPA (Fig. S7), while there was lower diversity in CRI and GRO-MCH in the Northern group (Fig. S7). In the Southern group, PAN had the highest diversity rates, and CLB and ECU showed the lower diversity (Fig. S7).

Genetic variation in Lutjanus guttatus

Lutjanus guttatus is a representative species of the Lutjanidae family, known as snappers, which includes a large number of commercially important fishes (FAO, 2016). Genetic characterization of populations is essential to preserve the genetic diversity of this fishing resource. Most evaluations of population structure and phylogeography of snapper species have used mitochondrial data and have suggested panmictic populations (Garber, Tringali & Stuck, 2004; Zhang, Cai & Huang, 2006; Gomes, Sampaio & Schneider, 2012; Reguera-Rouzaud et al., 2020). Hernández-Álvarez et al. (2020) recovered a panmictic pattern for L. Peru and L. guttatus, showing the limitations of using a single genetic marker to resolve genetic structure for species with large effective population sizes and constant but probably very small gene flow, like the Snappers. The low resolution of the mitochondrial data for this group was recently evaluated in the delimitation of L. campechanus and L. purpureus, which were first considered a single species based on mitochondrial data (L. campechanus) (Gomes, Sampaio & Schneider, 2012). Nevertheless, SNP data shown to clearly segregate these two taxa into two well-differentiated species (Pedraza-Marrón et al., 2019).

Contrasting patterns of genetic differentiation that we found between neutral and outlier loci in L. guttatus are consistent with the hypothesis that the genomic data, particularly polymorphic loci as SNPs, are better for resolving genetic questions in L. guttatus where a model of selection under migration was recovered. Loci showing adaptive selection clearly separated the populations into Northern and Southern clusters (F_CT = 0.4861; p < 0.001; Table 3), while no differences were found when using neutral loci. The lack of differences for neutral loci, in both mitochondrial and nuclear DNA, suggests the existence of sufficient gene flow to maintain the connectivity among populations at a large spatial scale, probably through larval dispersal as observed in pelagic species (Ward et al., 1994; Díaz-Jaimes et al., 2010). Meanwhile, the observed differences in outlier loci may be related to adaptations to environmental factors that promote genetic differences in particular loci (Graham et al., 2018). Since selection is acting in the presence of gene flow, the observed pattern should be consistent with the migration-selection equilibrium model (Graham et al., 2018; Bernardi, 2022). The migration-selection balance plays an important role in the genetic divergence of populations for marine species (Bernardi, 2022). If selection is strong enough to overcome the effects of migration, it results in differentiation at loci subjected to adaptive processes but not at neutral loci, which results in a distinct pattern determined by equilibrium between migration and genetic drift (Ribeiro, Lloyd & Bowie, 2011; Graham et al., 2018; Bernardi, 2022).

Outlier loci have been shown to be informative in detecting patterns of shallow divergence in marine species (Milano et al., 2013; Torrado et al., 2020). The evaluation of markers affected by natural selection are useful for determining the effect of oceanography on populations of marine species (Milano et al., 2013; Pecoraro et al., 2018; Torrado et al., 2020). Unfortunately, for non-model species the lack of a reference genome hinders the evaluation of the effect that selection may have on genetic structure (Milano et al., 2013). Our blast evaluation for OL failed to identify any particular protein that could be affected by the selection that we detected, although a reference genome was used. This could be due to the lack of a reference genome for the species or the method employed in the present work, as we selected a polymorphic region of the genome generating 1973 RAD loci of 280 pb long, which could have resulted in a most reduced representation of the L. guttatus genome.

The efficiency of SNPs for detecting genetic variation in marine organisms has been widely confirmed, but it is important to highlight the advantages of SNP data to analyze separately outlier and neutral loci, which allowed us to identify genetic patterns at different scales. As an example, genomic data in yellowfin tuna Thunnus albacares populations, found genetic structure among the Atlantic, Indian and Pacific oceans, but by using a set of 33 outlier loci, additional differences were detected at the intra-oceanic level for the Pacific and Indian populations (Pecoraro et al., 2018). In addition to the usefulness of SNPs for detection of genetic structure at small scale in marine organisms, outlier loci allow understand how habitat heterogeneity promotes the adaptation of organisms to environmental factors. On the other hand, demographic process determined with neutral loci allow to identify barriers for gene flow between populations (Pecoraro et al., 2018). Nevertheless, in most of the studies using neutral loci it is usually revealed weak or undetectable genetic structure (Mamoozadeh, Graves & McDowell, 2019; Maroso et al., 2021), which contrasts with the well-defined genetic differences obtained with the loci under selection. This was reported in Sparus aurata (Gilthead seabream), where neutral loci detected three genetic clusters in the Mediterranean, “but this variation was not clear and low supported”, while by using outlier SNPs there were revealed two additional populations in the eastern Mediterranean related with adaptations to environmental factors (Maroso et al., 2021). However, differences resulting from adaptive selection have limitations on elucidating the effect of demographic processes or other evolutionary forces as genetic drift and therefore on their usefulness to reconstruct the evolutionary history of populations. The contrasting results obtained with neutral and outlier loci for Lutjanus guttatus in this study reinforce the idea that demographic process as migration may be influencing our ability to detect genetic structure for the species, while there is a region of the genome influenced by adaptive selection that could be related with environmental factors. The challenge in the identification of factors contributing in the adaptation of organisms to environment lies in understanding how multigenic features are related to environment variables (Moore et al., 2014).

Oceanographic conditions influencing genetic variation

Historically, genetic structure for rocky reef fishes distributed in the TEP has been found to be strongly associated with the presence of sandy gaps in the TEP that act as barriers (Sandoval-Huerta et al., 2019; Torres-Hernández et al., 2022). However, snapper species are characterized by high gene flow mediated by pelagic larval dispersal (Hernández-Álvarez et al., 2020), and they inhabit different habitats during different life stages (juveniles individuals are found in mangroves or estuaries and adults inhabit rock substrate), making these sandy gaps less important as barriers (Reguera-Rouzaud et al., 2021). This is corroborated by our results, since the limit between the Northern and Southern genetic groups was associated with oceanographic conditions and the currents velocity (Fig. 4) more than with the presence of the sandy gaps (Fig. 5). The complex oceanographic conditions found in the Gulf of Panama have promoted endemism in the region (Robertson & Allen, 2015). The seasonal winds that generate a strong surface upwelling variation between the eastern and western Gulf of Panama have allowed variation in the growth rates of coral reef (Randall et al., 2019). Furthermore, the wind system has an effect in temperature, salinity, oxygen among other factors, promoting heterogeneous environment conditions, which have promoted genetic structure in fishes as the red-head goby (Elacatinus punticulatus) (Sandoval-Huerta et al., 2019), and in species of Lutjanus such as the red (L. guttatus) and yellow (L. argentiventris) snappers (Reguera-Rouzaud et al., 2021). Correlation found with our OL and the oceanographic condition found in the Gulf of Panama (Fig. 4B) are congruent with the hypothesis that the environment in the region is playing an important role in adaptive selection.

The barrier for the two main groups was found at the Gulf of Panama, an area where the current system and temperature are influenced by the Inter Tropical Current Zone (ITCZ). The ITCZ is located north of the equator (between 4–8°N), representing the limit between the northern and southern eastern Pacific areas, and the effect of the temperature as one of the variables that explain the genetic variation found with the OL in L. guttatus (Fig. 4) corroborated the important role for selection that the environment plays in the Gulf of Panama. Also, the ITCZ represents the southern boundary of the eastern Pacific warm pool (EPWP) located along the coast of southwestern Mexico and Guatemala, where warm, low-salinity surface water and a shallow thermocline predominate (Fiedler & Lavín, 2017), which is in accordance with the five environment variables that found to be correlated with the populations structure (Fig. 5). South of the ITCZ, the Equatorial Cold Tongue predominates, which is largely influenced by cold waters from the Peru Current and equatorial upwelling (Wyrtky, 1981). These contrasting oceanographic conditions may influence the observed differences due to adaptive selection, and gene flow through the ITCZ may be not sufficient to homogenize allele frequencies for outliers subjected to selection. The Pacific shelf of Panama is characterized by complex oceanographic patterns caused by the existence of two semi-open ocean basins with contrasting hydrography in terms of temperature, salinity and primary productivity (D’Croz & O’Dea, 2007). In this shelf, the Northern part of the Eastern Pacific Countercurrent ends, splitting into a northern and a southern branch, increasing differences in environmental factors. The southern branch flow reaches South American coasts, while the northern branch is diminished by the prominence of the Azuero Peninsula, which may limit larval dispersal toward northern areas (Fig. S7). Limited gene flow may result in the disruption of the selection-migration balance, generating the pattern of differentiation observed in outlier loci (Table 2 and Table 3).

The effect of the currents were also found with EEMS results, which recovered a barrier for the gene flow located in front of the Mexican coast (Fig. 5), Although any of the other results indicate a genetic differentiation from this locations and considering that must be taken with caution because locations from the Gulf of California have small sample size, correlation found between the genetic distance for OL and the current velocity (Figs. 4B and 4D) reinforced the idea that the current could contributed in limit movement of individuals between locations in the Northern group. In previous work, genetic structure was found in the region of the Gulf of California and Sinaloa, this genetic differentiation are related with the current system (Sandoval-Huerta et al., 2019) and habitat discontinuity (Reguera-Rouzaud et al., 2021). This suggests that future exploration with more samples and using another genomic markers are necessary to determine whether oceanographic conditions in the Northern locations are acting as a barrier for L. guttatus.

Geographic distance and migration

Geographic distance did not show an effect on genetic differences for L. guttatus populations using either NL or OL, so IBD is not contributing to any pattern of genetic structure in L. guttatus (Figs. S4 and S5). The pattern of migration recovered with results from EEMS indicates that the currents facilitate the dispersal of larvae, and that the migration is lower near to the shoreline but wider in offshore waters (Fig. 5). The high residence of adult individuals along the coast supports the dispersal through larval drift (Hernández-Álvarez et al., 2020) and broadly coincides with the migration model recovered, where corridors for migration are consistent with the trajectory of the Costa Rica Coastal Current (CRCC) (Fig. 5).

For Lutjanus guttatus, larval dispersal seems to be sufficient to homogenize populations throughout the species range for neutral loci, but not sufficient to counteract the genetic differences resulting from adaptive selection in outlier loci. Our results coincide with those observed in the blue-banded goby, where the current system from the Galapagos Archipelago allowed gene flow, recovering a panmictic pattern for the neutral loci, while results using outlier loci indicated local adaptation associated with oceanographic conditions like temperature (Bernardi, 2022).

Some studies using microsatellites in snappers have assessed the effect of oceanographic gradients and/or environmental factors on snapper species distributed in the TEP (Reguera-Rouzaud et al., 2021). Movements of adults for lutjanids are limited to the feeding areas (Hernández-Álvarez et al., 2020), while juvenile could move through the mangroves or estuaries. This general pattern of movement influenced by the current system in the TEP is not observed in the Gulf of Panama, where the currents limit dispersal for reef fishes (Sandoval-Huerta et al., 2019; Pedraza-Marrón, 2014). Currents also have showed an effect in genetic variation of lutjanids (Reguera-Rouzaud et al., 2021). Genetic studies in L. peru based on microsatellites rejected the panmixia model after analyzing the role of larval dispersal on the species’ genetic structure (Munguia-Vega et al., 2018). Likewise, recent studies with microsatellites found that oceanographic gradients and the currents have contributed to patterns of genetic structure in L. peru and L. argentiventris (Reguera-Rouzaud et al., 2021) in the TEP, this idea is according with the correlation found with the genetic variation of the OL and the current velocity of the TEP (Fig. 4D; Fig. S6). Therefore, the heterogeneous environment of the TEP could be acting as a selective pressure in snapper species.

Conservation implications

Lutjanus guttatus is an important fishery resource throughout the TEP (Hernández-Álvarez et al., 2020). In countries such as Costa Rica, Colombia and Mexico, several studies have been done to develop technologies to establish aquaculture of the species (Rojas et al., 2009; Instituto Nacional de Pesca, 2013; Abdo-De la Parra et al., 2015; Chacón-Guzmán et al., 2020). Although the genetic information used for the establishment of management units is limited, a population study carried out in the Olympia oyster concluded that the use of the SNP analysis evaluating neutral and outlier loci separately have implications in the management of fisheries species (Silliman, 2019). The overall genetic structure found in L. guttatus across the TEP seems to be related with a signal of selection, and, as reported by Silliman (2019), adaptive loci could be used for genetic monitoring protocols.

Finally, the genetic structure separating the spotted rose snapper into two main clusters (Northern and Southern groups), is important information for the conservation of this important fishery resource. The use of SNPs highlights the importance of selecting adequate markers for evaluating genetic structure in snapper, and particularly in L. guttatus, to avoid erroneous conclusions. Although the species is cataloged as least concern by the IUCN, there is an evident lack of information. The general idea of panmixia for L. guttatus remains controversial in considering that the species maintains constant gene flow. Nevertheless, the selection pattern found herein indicates that in fact the oceanographic conditions in the TEP are influencing the genetic structuring of L. guttatus. Although to date an accurate assessment for fishery status has not been made, we consider necessary to reevaluate the effect that the fishery could have on the species employing additional genomic data. With our results two clusters that are affected by selection should highlight the need of genetic data for being considered for management purposes. Further studies expanding the genome coverage would improve the genetic information for Lutjanus guttatus in order to confirm the absence of genetic structure in neutral regions and the relevance for the species of inhabiting in a heterogeneous environment, and contribute to the species conservation.