Genetic connectivity in the spotted rose snapper Lutjanus guttatus (Lutjaniformes: Lutjanidae) between Mexico and Panama throughout the Tropical Eastern Pacific

Abstract The spotted rose snapper, Lutjanus guttatus, is an important fishery species with high potential for aquaculture. Genetic characterization of its natural populations is necessary to avoid stock collapse and loss of genetic diversity. Previous studies carried out in the Tropical Eastern Pacific (TEP), however, have shown contrasting results in the genetic structure of fish populations, particularly in species of Lutjanidae. Therefore, to understand the genetic structure of spotted rose snapper in the TEP, twelve microsatellite loci were used to assess the genetic diversity and explore the hypothesis of population genetic structure in samples of the species collected throughout the TEP. Fin clips from 186 sampled individuals (27 to 49 per site) were analyzed from five sites in the three regional biogeographic provinces, delimited by shoreline reef habitat breaks: La Paz (Cortez province), Colima and Oaxaca (Mexican province), Chiriqui and Port of Panama (Panamic province). Results of global Analysis of Molecular Variance (AMOVA), population pairwise FST, hierarchical AMOVA, and a discriminant analysis of principal components (DAPC) reflected a panmictic population involving the entire set of sampled sites. The role of larval dispersal, post-recruitment migration, and marine current dynamics as drivers of genetic connectivity in this species is discussed.


INTRODUCTION
The spotted rose snapper, Lutjanus guttatus (Steindachner, 1869), is distributed throughout the Tropical Eastern Pacific (TEP) from the Gulf of California, Mexico, to Peru, including oceanic islands.It is an important food and recreational fishery species with a high market price of US$6-8 per kg (Sarabia-Méndez et al., 2010;Ibarra-Castro et al., 2012) and a high potential for aquaculture in Mexico (García-Ortega et al., 2005).Despite its great economic importance, there are few studies that use molecular techniques to identify stocks in L. guttatus, critical information for the to managing exploitation of mixed or locally discrete stocks and avoiding loss of genetic diversity and possible stock collapse (Pauly et al., 1996).
The TEP is a region that extends along 2,500 km from the equator north to the southernmost tip of the Baja California Peninsula (Kessler, 2006).The TEP coast is a highly dynamic environment, with sea temperatures ranging from warm to temperate, upwelling systems and various large gyres, alternating currents, and large rocky-habitat discontinuities that may greatly influence the genetic connectivity of populations (Robertson, Cramer, 2009;Sandoval-Huerta et al., 2019).These physical characteristics can affect distributions of species with narrow environmental tolerances and influence the dispersal of pelagic larvae, resulting in variable gene flow (from reproductive isolation to high connectivity) between adjacent populations (García-De León et al., 2018;Sandoval-Huerta et al., 2019).Noé Díaz-Viloria, Adriana Max-Aguilar, Mailin I. Rivera-Lucero, Elaine Espino-Barr, Nicole Reguera-Rouzaud, Andrea Casaucao-Aguilar and Ricardo Perez-Enriquez There are several hypotheses about biogeographic partitioning in the TEP, where environmental and ecological differences have promoted speciation in the absence of isolation of diverging populations (Briggs, Bowen, 2012).Robertson, Cramer (2009) determined that three biogeographic provinces exist: the Cortez (Gulf of California and the southernmost Pacific Baja California), the Panamic (southward) and the Island province.However, Walker (1960) had previously defined two provinces in Mexico based on the distribution of locally endemic reef fishes: a Cortez Province from the Pacific coast of Baja California below 25° N, including the Gulf of California, and a Mexican province for the remainder.Within the TEP, there are also two major breaks in the distribution of shoreline reef habitats, consistent with shoreline extensions of sand and mud, named the Sinaloan Gap (370 km in the SE Gulf of California), and the Central American Gap (around 1000 km, from the Gulf of Tehuantepec, Mexico, to El Salvador).These gaps separated three mainland provinces (Cortez, Mexican, and Panamic) defined by Hastings (2000).Finally, the TEP was split into two provinces, the Galapagos and the remainder of the region, by Spalding et al. (2007).
Many studies carried out in the TEP have obtained measures of gene flow to explore levels of population genetic differentiation and evaluate the influence of habitat gaps and oceanographic processes, with contrasting results.Gene flow rates among coral and North Pacific hake (Merluccius productus (Ayres, 1855)) populations along coast are high, although populations at the northernmost and the southernmost peripheries appear to be more genetically isolated (Lessios, Baums, 2017;García-De León et al., 2018).Strong subdivisions between populations of the goby Elacatinus puncticulatus (Ginsburg, 1938) were better explained by local oceanographic processes than the largest habitat discontinuities (Sandoval-Huerta et al., 2019).Meanwhile, basin-wide connectivity and shallow population structure in the olive ridley sea turtle (Lepidochelys olivacea (Eschscholtz, 1829)) seems due, in part, to their low nesting site fidelity and broad foraging ranges (Silver-Georges et al., 2020).
Recently, three studies of population genetics in different snapper species (Lutjanidae) in the TEP were completed.The first study carried out using sequencing of the mtDNA control region on Lutjanus peru (Nichols & Murphy, 1922) (~800 bp) and L. guttatus (576 bp) found high overall levels of genetic diversity and a lack of genetic differentiation for both species (Hernández-Álvarez et al., 2020).These results indicate that equatorial and subtropical residents display high levels of connectivity and highlight that no significant effect of environmental differences between Cortez and Panamic provinces exist.In contrast, a study using 13 microsatellite loci on L. peru and 11 microsatellite loci on L. argentiventris (Peters, 1869) evaluated genetic diversity across 10 and five locations, respectively (Reguera-Rouzaud et al., 2021).Significant genetic structure was identified in both species, but the pattern of genetic structure differed between species.These authors suggested two possible drivers, including isolation by distance (IBD) at a spatial scale of more than 2,500 km and the presence of potential barriers to gene flow at smaller scales (< 250 km).The most recent study of L. guttatus, covering nearly all its distributional area, used 2003 single nucleotide polymorphisms (SNPs); including neutral loci (1858 SNPs) and outlier loci (145 SNPs) to assess genetic variation and population genetic structure (Mar-Silva et al., 2023).For neutral loci (NL), no differences were found, but with outlier loci (OL) two clusters were found dividing at the Gulf of Panama, suggesting the role of selection in generating genetic differences in L. guttatus.

Genetic connectivity in Lutjanus guttatus
Because microsatellites are assumed to be neutral markers, codominant with Mendelian inheritance, have higher mutation rates than mtDNA (Liu, Cordes, 2004;Hernández-Álvarez et al., 2020), and higher levels of genetic diversity in number of alleles per locus (N A ) and heterozygosities (H O and H E ) than SNPs (Mar-Silva et al., 2023), they are better suited to studying mutation-drift equilibrium and gene flow among populations.The hypothesis of the existence of population genetic structure in L. guttatus in the TEP was therefore retested using a set of 12 mostly tetranucleotide microsatellite loci to explore genetic diversity and neutral population genetic structure among five sites distributed throughout the TEP.

MATERIAL AND METHODS
Sampling.Fin clips were collected from 30-50 individuals from each of five locations across the three putative mainland provinces in the TEP (three from Mexico and two from Panama) and preserved in ethanol (80%) (Fig. 1; Tab. 1).Note that this study did not include specimens from the Galapagos Islands, an island province where the species is also reported (Robertson, Allen, 2015), because originally was focused on the three mainland provinces scheme.Captures of L. guttatus at every site were supported by commercial fishermen.Aljanabi, Martinez (1997) that contained a modified homogenizing buffer (5M NaCl, 1M Tris-HCl, 0.5M EDTA, 10% SDS, pH 8.0), and DNA extracts were standardized to 30 ng/µl.A set of 12 microsatellite loci were selected based on their numbers of alleles (moderate to high polymorphism) and general conformation to Hardy Weinberg Equilibrium (HWE) expectations (Perez-Enriquez et al., 2020;Tab. 2).The PCR was carried out following Schuelke (2000) in an 11 μl volume containing 1 μl of DNA (30 ng/μl), Taq Buffer (1×), MgCl 2 (1.5 mM), dNTPs (0.25 mM), forward primer (0.1 μM), reverse primer (0.4 μM), M13+dye (0.4 μM; 6-FAM, VIC, NED, or PET (Applied Biosystems; Tab. 2), and Taq polymerase (0.04 U/μl).The PCR thermal cycling conditions were as follows: 94ºC for 5 min; 30 cycles at 94ºC for 30 sec, annealing temperature (Tab.2) for 45 sec, and 72ºC for 45 sec; eight cycles at 94ºC for 30 sec and 53ºC for 45 sec; and a final extension at 72ºC for 10 min.The PCR products were mixed in poolplex (Tab.2), and 2.0 μl of every poolplex were used in fragment analyses on an ABI3130 automated DNA sequencer (Applied Biosystems) at the University of Arizona Genetics Core.Fragments sizes were obtained relative to the GeneScan 500 LIZ Size Standard (Applied Biosystems).
Microsatellite genotyping and data analysis.Alleles were sized using the GeneMarker version 2.4.0 (Softgenetics, 2012).Individuals with more than 15% missing data were removed (two from Colima, one from Oaxaca, and one from Port of Panama), and 186 individual L. guttatus were retained for further analysis (Tab.S1).Binning within each allelic class was carried out with Flexibin version 2.0 (Amos et al., 2007), with all 186 retained individuals successfully genotyped at all 12 loci (2,232 genotypes).Allelic frequencies and null allele frequencies were obtained with Arlequin version 3.5 (Excoffier, Lischer, 2010) and FreeNA (Chapuis, Estoup, 2007), respectively.To test if the alleles were drawn from the same distribution in all populations, Fisher exact tests were carried out using Genepop version 4.7 (Raymond, Rousset, 1995;Rousset, 2008).Genetic diversity indices, including the number of alleles (N A ), effective number of alleles (N EA ), number of private alleles (N PA ), and observed and expected heterozygosities (H O and H E ), were assessed with GenAlEx ver.6.5 (Peakall, Smouse, 2012).Kruskal-Wallis tests for N A , N EA , H O , and H E looking for differences among locations were performed with STATISTICA version 8 (StatSoft.Inc., Tulsa, OK, USA).The HWE and Linkage disequilibrium (LD) were evaluated through exact tests with Arlequin, and sequential Bonferroni adjustment at α = 0.05 was carried out to control the effects of multiple testing (Rice, 1989).Three poolplexes were combined using amplicons from four loci having similar annealing temperatures and tagged with different flouorescent dyes.
The statistical power to find genetic differentiation among populations was tested with POWSIM version 4.1 (Ryman, Palm, 2006) with the following parameters: sample size from 29 to 49 individuals (based on mean n, Tab.S2) and one to 12 loci, with ~20 alleles each.In addition, a predefined F ST = 0.0025 was employed, representing the lower non-significant F ST value found in this study.Although, in L. peru (F ST = 0.0159, P = 0) and L. argentiventris (F ST = 0.019, P = 0) the significant F ST values were higher (Reguera-Rouzaud et al., 2021).Finally, interaction/permutation factors for the Fisher's exact test (10,000, 1,000, 10,000), five populations, an effective population size of 2,000, 10 generations under drift, and 1,000 runs were also included in simulations.
Genetic population structure was evaluated through a global Analysis of Molecular Variance (AMOVA), pairwise F ST , and a hierarchical AMOVA using Arlequin.Three groups were created for the hierarchical AMOVA: 1) La Paz, 2) Colima and Oaxaca, and 3) Chiriqui, and Port of Panama.The significance of the covariance components associated with the different possible levels of genetic structure (within populations, within groups of populations, among groups) was tested using non-parametric permutation procedures (10,100 permutations).These three groups were in agreement with the biogeographic provinces proposed by Hasting (2000).A phylogenetic tree was constructed using Nei genetic distances, the UPGMA method, and 10,000 bootstrap replicates with PHYLIP (Felsenstein, 2005).Statistical support was obtained by bootstrap percentages in every branch.
Discriminant Analysis of Principal Components (DAPC) was performed in R (R Development Core Team, 2011) package adegenet (Jombart, 2008) to explore clustering in the data (Jombart, Collins, 2021).The data were transformed through a PCA and a discriminant analysis was applied to the retained principal components to minimize intra-group variability while maximizing inter-group variability (Jombart, Collins, 2021).First, genotype data were imported using import2genind, using a genetix file (.gtx).Second, the best number of K clusters was determined de novo using find.clusters,we keep all the information, specifying to retain 200 PCs.After obtaining the graph of Bayesian information criterion (BIC), and no clear elbow was observed, the best number of K clusters was obtained by the difference of K i+1 -K i .Third, DAPC was run with dapc(obj, grp$pop), 80 PCs and four discriminant functions retained, accounting for 82% of variance.Each DAPC was cross-validated and rerun with suggested PCs according to the best proportion of successful outcome prediction.Two STRUCTURElike plots with membership of all individuals and with mixed individuals having no more than 0.5 probability of membership were created using compoplot in adegenet.

RESULTS
The allelic frequencies showed two patterns.First, most medium and low frequency alleles (< 10%) were present at most locations (data not shown).Second, there were frequency differences among locations for some of the most common alleles in nearly all loci.However, after Fisher exact tests, such differences in distribution of allelic frequencies were significant only for Lgut38 and Lgut44 (P < 0.05), and none were significant after sequential Bonferroni adjustment (P > 0.006).
Higher values in N A and N EA were found at the center of the species distribution range (Colima and Oaxaca, Mexico) than in locations at the extremes of the sampled range.The opposite was observed in H O , where higher values were found for La Paz and the Port of Panama than for the rest of locations (Fig. 2).Such results were an effect of the sample size per location (Tabs.3, S2).However, after Kruskal-Wallis tests, no differences among locations were observed in N A (P = 0.28), N EA (P = 0.74), H O (P = 0.56), or H E (P = 0.99).Heterozygote deficits were observed but not significant (P > 0.006) in most loci except Lgut37, Lgut38, and Lgut44.These three loci had null allele frequencies greater than 0.05 and deviated from HWE in most locations.After sequential Bonferroni adjustment, however, zero loci showed significant HWE deviations (P < 0.004) (Tab.S2) or signs of linkage disequilibrium (LD) (P > 0.0005), through the five locations.Statistical power testing showed that there was greater than 95% capacity to detect significant genetic structure with ten or more loci (Fig. 3).Yet, global AMOVA (F ST = 0.00012, P = 1), population pairwise F ST (Tab.4), and hierarchical AMOVA (F CT = 0.00126, P = 0.1215) showed no evidence of population genetic structure.Exclusion of the three loci with high frequencies of null alleles (Lgut37, Lgut38, and Lgut44) did not change these results, and subsequent analyses therefore included all 12 loci.The topology of the UPGMA tree did group the sampled populations by their respective biogeographic provinces, but this result is tempered by low bootstrap support for Chiriquí-Port of Panama (Fig. 4).and P values (below diagonal).Note that all comparisons were not significant (P > 0.05).Genetic connectivity in Lutjanus guttatus DAPC showed five mixed genetic groups (Fig. 5), and most individuals in each cluster were found to have a very high (90-100%) membership probability to their own group, despite a variable but lower proportion of samples exhibiting mixed memberships.No transitional zones with higher proportions of putatively admixed individuals were observed between pairs of neighboring sites (Fig. 6).Of 23 individuals with membership probabilities lower than 0.5, most exhibited membership to more than two clusters (Fig. 7).

DISCUSSION
Lutjanus guttatus did not show any evidence for genetic population structure by traditional analysis, such as global AMOVA, population pairwise F ST , hierarchical AMOVA, and IBD.However, DAPC showed distinctive local genetic subunits of a large population but with no clear cuts among them.The presence of admixed individuals, with membership to more than one external clusters, could be outlining contemporary connectivity that both a) connects the closest sampled populations and b) spans the entire geographic range (e.g., the major cluster found in La Paz; Cluster LAP) is also found present in admixed individuals in Panama, where Cluster PP predominates), with c) no transitional zones of admixture.These results as well as the low membership probabilities of 23 individuals reflect panmixia involving the entire set of sampled populations including all three mainland provinces, which implies that long-distance connectivity is as prevalent as short-distance exchanges in L. guttatus.Notably, despite the lack of genetic population structure, the UPGMA tree showed three apparent genetic groups in L. guttatus, including La Paz, Colima with Oaxaca, and Chiriquí with Port of Panama, although this last grouping had low bootstrap support.Such apparent genetic groups were in good agreement with the Cortez, Mexican and Panamic provinces suggested by Hastings (2000).
Juveniles of L. guttatus can use the soft bottoms adjacent to rocky reefs (Mariscal-Romero, Van der Heiden, 2006;Saucedo-Lozano et al., 2006) or mangroves to feed (Allen, 1995;Gutierrez-Barreras, 1999;Vega et al., 2015;2016b;Medina-Contreras et al., 2021).Such ecological adaptations to different environmental conditions in different habitats enable nomadic individuals of L. guttatus to migrate around habitat discontinuities that restrict movement in L argentiventris (absence of mangroves) and L. peru (absence of rocky reefs), possibly resulting in the connectivity between the Gulf of California-Colima and Oaxaca-Chiriquí and Port of Panama regions seen in this study.
Larval dispersal and possible migration of nomadic individuals were mentioned by Mar-Silva et al. (2023) as mechanisms that may have differentially contributed to the high contemporary genetic connectivity seen among locations.Mar-Silva et al. (2023) found no differences using neutral loci (dataset of 1858 SNPs) and genetic differences using outlier loci (dataset of 145 SNPs) suggesting the role of selection in this case.Results of the present study were in agreement with the conclusion of "no differences" or panmictic population, which is explained because microsatellite loci are assumed to be neutral markers.F ST , when based on neutral genetic markers, estimates the degree to which populations have diverged from one another as a result of gene flow and genetic drift, without the selection effect in the equation (Freeland, 2006).
If genetic connectivity is the result of larval dispersal or migration of juveniles, preadults, or adults of L. guttatus every generation, it may result in greater resilience to local extirpation because fishing areas could be recolonized relatively quickly (Craig et al., 2006).On the other hand, local extinction could modify the pattern of connectivity, increasing the relative geographic distance among populations (Saavedra-Sotelo et al., 2010).Nevertheless, additional research that includes more sites within the distribution range of the species, such as the Galapagos islands, which represents a province not assayed in this study, uses potentially adaptive genetic markers (e.g., SNPs), and adequate sample size per site (Flesch et al., 2018) is required to both further improve our understanding of the population dynamics of L. guttatus in the TEP as well as the influence of environmental variables on its genetic makeup.

FIGURE 1 |
FIGURE 1 | Sampling sites for adult Lutjanus guttatus along the Tropical Eastern Pacific.La Paz (PAZ), Colima (COL), and Oaxaca (OAX) are in Mexican waters and Chiriquí (CHI) and Panama Port (PAN) are in Panama.

FIGURE 2 |
FIGURE 2 | Genetic diversity in Lutjanus guttatus from five sampled locations.Numbers of alleles (N A ), effective alleles (N EA ), and private alleles (N PA ) and observed (H O ) and unbiased expected heterozygosities (uH E ).

FIGURE 3 |
FIGURE 3 | Estimates of the statistical power (percent) for finding significant population genetic structure when using different numbers of microsatellite loci.

FIGURE 4 |
FIGURE 4 | UPGMA dendrogram of Lutjanus guttatus from five sampled locations along the Pacific coast of Mexico (La Paz, Colima, Oaxaca) and Panama (Chiriquí, Port of Panama) based on Nei's genetic distance (1972).Numbers on the nodes indicate the percent of times the illustrated topology was found with 10,000 bootstrap replicates.

FIGURE 8 |
FIGURE 8 | Mantel test of correlation between genetic (y-axis) and geographic distances (x-axis) among sampled locations.The black line represents the central tendency among the dots in the scatterplot.

TABLE 1 |
Lutjanus guttatus tissue collections.Samples sizes (n), mean standard lengths (centimeters), standard deviations (inside parentheses), collection dates, and biogeographic province of origin in the Tropical Eastern Pacific.

TABLE 2 |
The twelve microsatellite loci selected for population genetic analysis on Lutjanus guttatus.

TABLE 3 |
Genetic diversity in Lutjanus guttatus from the five sample locations.Sample sizes (n), number of alleles (N A ), observed (H O ), and expected (H E ) heterozygosities.Note that values with asterisk showed significant deviations from Hardy-Weinberg Equilibrium after Bonferroni adjustment (P < 0.004, Tab.S2).

TABLE 4 |
Pairwise F ST values using 12 loci for sampled Lutjanus guttatus populations (above diagonal)