Population differentiation in Haemulon plumieri juveniles across the northern coast of the Yucatan Peninsula

The white grunt Haemulon plumieri is an important local artisanal fishery resource along the northern coast of the Yucatan Peninsula. In the present study, the population structure of recently post-settled H. plumieri was assessed in 4 locations (Celestun, Dzilam, Yalahau and Chacmochuch) from the northern coast of the Yucatan Peninsula by using a phenotypic approach based on otolith shape analysis and by a genotypic approach using 4 microsatellite DNA markers. Both approaches indicated that the 4 location samples were different from each other: significant differences in otolith mean shapes were observed (PERMANOVA, F3,196 = 10.9879, p = 0.0001), and significant genetic differentiation was also observed (global FST = 0.0560, p < 0.001). The isolation-by-distance test revealed a significant correlation between the FST values and the geographic distance in the subpopulations (Z = 98.1407, R2 = 0.8680, p = 0.001). The greatest genetic differences were observed between the Celestun and Chacmochuch samples (FST = 0.069, p < 0.001), and the lowest genetic differences were observed between the Yalahau and Chacmochuch samples (FST = 0.029, p < 0.001). Our data indicate that dispersal may be limited in this species, and this could make white grunts vulnerable to local overfishing.


INTRODUCTION
The white grunt Haemulon plumieri (Lacepède, 1801) is a small coastal fish inhabiting the Western Atlantic, from Chesapeake Bay south to the Gulf of Mexico and the coast of Brazil (Lindeman & Toxey 2002).This species supports an important local artisanal fishery in the northern coast of the Yucatan Peninsula, with some other species of higher commercial importance such as snappers (Lutjanidae) and groupers (Serranidae).However, in this area, the white grunt population remains as an unregulated resource with no maximum capture fees and unspecified catch size limits.For this reason, it is important to generate baseline knowledge of H. plumieri's population structure in southern Mexico to provide information for future conservation and management efforts.
In this sense, the most successful way to evaluate the population structure of a given species is to use multiple approaches to maximize the probability of defin-ing its population structure (Begg & Waldman 1999, Abaunza et al. 2008).Thus, a recommended protocol is the complementary application of genotypic and phenotypic-based approaches (Ihssen et al. 1981).The phenotypic approach based on analysis of the shape of fish otoliths has proven to be very useful (Friedland & Reddin 1994, Stransky et al. 2008), as has the use of DNA microsatellite markers (Cadrin et al. 2005, Purcell et al. 2006, 2009).Microsatellite loci are highly polymorphic markers (Jarne & Lagoda 1996) and are able to detect low but significant levels of differentiation among populations, even in the presence of high levels of gene flow (Hauser & Ward 1998), with the capability to detect differences on short spatial scales (up to a few hundred kilometers) (Knutsen et al. 2003).
In the present study, we used otolith shape analysis and DNA microsatellite markers to assess the population structure of post-settled white grunt Haemulon plumieri juveniles (within the same cohort or breeding season) from the northern coast of the Yucatan Peninsula.

Sampling
Fishes were collected during 2 wk with a trawling epibenthic net at depths less than 1.5 m in January 2010 at 4 locations along the northern coast of the Yucatan Peninsula: Celestun (CEL), Dzilam (DZI), Yalahau (YAL), and Chacmochuch (CHA) (Fig. 1).These localities have abundant prairies of the seagrass Thalassia testidinum.In the field, fishes were identified as Haemulon spp.and were fixed in 97% ethanol.They were then confirmed in the laboratory as H. plumieri according to the identification keys for the family Haemulidae (Lindeman 1986, Lindeman & Toxey 2002).Fifty H. plumieri (TL ranged from 4.0 to 11.0 cm) from each sampling site were chosen according to Ruzzante (1998), who suggested that at least 50 individuals are necessary for accurate estimation of genetic distances.Fish were measured for total length and standard length.Then, ~50 mg of gill tissue was removed from each individual, placed in vials with 90% ethanol, and stored at −20°C for the molecular analyses.

Otolith shape analysis
Sagitta otoliths were removed from each fish under the stereomicroscope (2×) (Olympus SZ51), using the 'open the hatch method' described by Secor et al. (1995).They were cleaned and stored dry in vials.Subsequently, each otolith was placed over a black surface with the sulcus side placed downward and the rostrum placed forward.All analyses were performed with the otolith pointing in the same direction to minimize distortion errors.A white transmitted light, high-contrast photo was taken with the aid of a stereomicroscope (2×) (Olympus SZ51).Each photo was properly labeled to avoid confusion and analyzed with the OPTIMAS image-analysis system v.6.1 (Bioscan).This analysis uses the elliptical Fourier shape analysis (Kuhl & Giardina 1982), which fits an arbitrary closed curve to an ordered set of data points in a 2-dimensional plane and generates a set of shape-representative variables subsequently used for statistical comparisons.For more details, see Kuhl & Giardina (1982), González-Salas & Lenfant (2007), and Villegas-Hernández et al. (2008).Briefly, the data were termed elliptic Fourier descriptors using the Shape 1.2 software (Iwata & Ukai 2002).First, we calculated the elliptic Fourier descriptors and the amplitudes of each harmonic, then calculated the average Fourier power spectrum (FP): this is the number of harmonics needed to describe effectively the shape of the otolith (Renaud et al. 1996).Finally, the total number of harmonics necessary for the statistical analysis was calculated as the total number of harmonics that were above 99% of the cumulative average power spectrum of Fourier (FPC), defined by FPC = Σ n Fig. 1.Location of sampling sites on the northern coast of the Yucatan Peninsula, Mexico

Genetic analysis
DNA extraction was carried out with a commercial purification kit (Invitrogen).Four microsatellite loci (HfAAC54, HfAAC31, HfAAC46, and HfAAC10) from the white grunt genome were amplified by PCR using a modification of Williams et al. (2004).The working volume was of 20 µl and contained 1.2 µl of DNA (80 to 100 ng µl −1 ), 11.6 µl of ultrapure sterile distilled water, 2.5 µl of 10× loading buffer TermoPol (New England BioLabs), 1.5 µl MgCl 2 (50 mM), 2.0 µl of dNTPs (10 mM), 0.5 µl of each primer (forward and reverse) (10 mM), and 0.2 µl of DNA Taq polymerase (5000 U µl −1 ) (New England BioLabs).PCR amplifications were performed on a TC-412 thermocycler (Techne) programmed for 5 min at 94°C for initial DNA denaturation, followed by 30 cycles at 94°C for 30 s, 30 s for annealing of each primer (50°C for HfAAC54; 55°C for HfAAC31, HfAAC46, and HfAAC10), and 30 s at 72°C for extension followed by a final extension step at 72°C for 5 min.All PCR products were resolved on 6% denaturing polyacrylamide gels with a vertical electrophoresis Sequi GenII R-GT Sequencing Cell (BIO-RAD) and stained with silver dye (Bassam et al. 1991).The resulting alleles from each specimen were identified and compared for relative size with a 50 bp sequence ladder (Promega).

Data analysis
Otolith shape data Data were analyzed with a PERMANOVA analysis (Anderson 2005), based on the Bray-Curtis dissimilarity measure (4999 random permutations and 95% confidence) that allows testing and partitioning the multivariate variation (defined by the distance measure used) according to individual factors in a fully balanced multi-way ANOVA (Anderson 2005).When the PERMANOVA was significant at the 95% confidence level, the levels of each factor were investigated through a posteriori pairwise comparisons using 4999 random permutations to obtain p-values.Furthermore, to support the PERMANOVA analysis, a non-parametric generalized discriminant analysis using a data array of the harmonic coefficients (4 coefficients × n harmonics) was also carried out to test the classification success of a number of discriminant functions to predict the identity of the white grunt population (Anderson & Robinson 2003).All otolith shape variables were examined for normality and homogeneity of variances prior to the generalized discriminant function analysis (Zar 1999).Then, the performance of the discriminant functions was evaluated using Cohen's Kappa statistic (κ), which provides an objective means of calculating the chancecorrected percentage of agreement between actual and predicted groups.The κ-values range from 0 to 1; 0 indicates that the discriminant analysis yields no improvement over chance, and 1 indicates a perfect agreement (Titus et al. 1984).

Microsatellite data
Measures of genetic variability including the number of alleles per locus, allele frequencies, observed heterozygosity, and gene diversity (heterozygosity expected from Hardy-Weinberg equilibrium [HWE] assumptions) were calculated using Genetix v.4.04 (Belkhir et al. 2004) and Arlequin v.3.1 (Excoffier & Schneider 2005). Wright's (1978) inbreeding coefficient (F IS ) was calculated for each locus and sample to evaluate if the white grunt samples were in HWE.The significance of F IS differing from zero was assessed with 10 000 permutations (Weir & Cockerham 1984).Departures from HWE can be caused by biological processes, such as inbreeding or population substructure (i.e. the Wahlund effect), or by technical issues, such as null alleles.Therefore, the software Micro-Checker v.2.3.3 (van Oosterhout et al. 2004) was used to infer the most probable technical cause of HWE departures, including null alleles and/or misscoring due to stuttering (Brookfield 1996).
Inter-population differentiation was evaluated with Slatkin's (1995) linearized F ST (a measure of genetic distance among populations).The statistical significance of these estimates was tested by a re-sampling jacknife permutation method (10 000 permutations) over loci and samples (Weir 1990) using Genetix v.4.04.The values of probability were modified with a sequential Bonferroni correction (Rice 1989).To assess the relative partitioning of genetic variation among sample sites, an analysis of molecular variance (AMOVA) (Cockerham 1969(Cockerham , 1973) ) was carried out both on multilocus genotypes and locus by locus with Arlequin v.3.1 with 10 000 permutations.
The effect of the geographic distance on the genetic differentiation was assessed by the isolationby-distance test (Mantel pro cedure) using the software IBD v.1.53(Bo honak 2002) to obtain the correlation be tween the geographic distance matrix (km, previously measured by shorted distance by sea) and the F ST values from the subpopulations.Also, a pairwise analysis of genetic distances among subpopula-tions was evaluated with Nei's index (Ds) (Nei 1978) with Genetix v.4.04 (Belkhir et al. 2004).
Finally, to identify any potential cryptic genetic variation within and among samples, a Bayesian clustering analysis was performed in STRUCTURE v.2.2 (Pritchard et al. 2000, Falush et al. 2003), assigning individual fish to groups without using prior information about their origin.The software was run for 10 interactions (repeat runs) per K cluster, and each run consisted of 10 000 bur-in steps and 10 000 Markov-chain Monte Carlo (MCMC), assuming independent allele frequencies among populations and the admixture ancestry model.Our analysis considered numbers of conceptual populations from 1 to 8; however, we performed an analysis to determine the optimal number of populations using the methodology of Evanno et al. ( 2005) using the STRUCTURE HARVESTED website program (Earl & von Holdt 2012) to estimate the probability that individual genotypes fall into k clusters and to assign individuals to specific clusters.Then, we plotted the genotypic proportions of each white grunt specimen within each predefined sample.

Otolith shape analysis
According to the FPC analysis, the first 12 harmonics were enough to describe the shape of the sagitta otolith at 99.9% of the total variation in the otolith shape.The first harmonic was eliminated, and only 11 harmonics were used in the Fou rier coefficients and amplitudes.The PERMANOVA analysis indicated differences in otolith mean shapes in the 4 sites (F 3,196 = 10.9879,p = 0.0001), and the subsequent pairwise comparisons showed that all samples were significantly different from each other (Table 1).The generalized discriminant analysis generated 3 statistically significant (p < 0.05) functions, indicating differences among the samples (Fig. 2).The 2 first discriminant functions accounted for > 87% of the total variation among the samples.In the classification matrix, a total of 94.0% of the individuals were classified correctly within each sample (94% in CEL, 92% in DZI, 96% in YAL, and 94% in CHA).Cohen's κ co efficient (κ = 0.92) confirmed the high success rate of classification described above.The classification efficiency was 92% higher than if it occurred by chance.

Genetic variation within samples
Overall, a total of 28 alleles were observed in the analyzed samples, displaying a similar degree of allelic variation at all loci, ranging from 6 alleles in locus HfAAC54 to 8 alleles in locus HfAAC46.Expected heterozygosity (over all 4 loci) ranged from 0.81 in CEL to 0.85 in YAL ( ferroni's correction (p > 0.0125).Null alleles were not detected in MICROCHECKER (data not shown).Thus, all subsequent analyses were performed using a multilocus approach.

Genetic differentiation between samples
Independently, the loci displayed very similar global linearized F ST values: HfAAC46 (F ST = 0.057, p < 0.001), HfAAC54 (F ST = 0.052, p < 0.001), HfAAC31 (F ST = 0.059, p < 0.001), and HfAAC10 (F ST = 0.051, p < 0.001).The multilocus genetic analysis revealed significant differentiation between sites (global F ST = 0.0560, p < 0.001), with 5.6% of the total genetic variation among subpopulations.The pairwise comparisons using heterozygote F ST and linearized F ST revealed significant differences among all sites (Tables 3 & 4).The greatest differentiation was observed in CEL vs. CHA (F ST = 0.069) and the lowest differentiation in YAL vs. CHA (F ST = 0.029) (Table 4).These results were reinforced with the multilocus AMOVA, which identified significant genetic structure among samples (global F ST = 0.0560, p < 0.001), with 5.5% of the total genetic variation explained by the variation among sampling sites (Table 4).The isolation-by-distance test revealed a significant correlation between the F ST values and the geographic distance in the subpopulations (Z = 98.1407,R 2 = 0.8680, p = 0.001) (Table 3).Nei's (1978) genetic distances (Ds) showed that CEL and CHA were very dissimilar sites (Ds = 0.432), while YAL and CHA were the most similar sites (Ds = 0.197).Finally, the Bayesian clustering analysis computed in STRUCTURE indicated that the 4 populations were sufficient as the optimal number of populations was K = 4, inferred from the modal value of ΔK = 8.03 at K = 4 (Fig. 3), corresponding to the highest level of structure or real cluster detected, in which 89.8% of the individuals were assigned to the site where they were sampled (Fig. 4).

DISCUSSION
To our knowledge, the present work is the first study using both phenotypic and genotypic approaches to explore the structure of juveniles of white grunt Haemulon plumieri at a regional level in southern Mexico.The application of the otolith shape analysis and DNA microsatellite markers indicated that the 4 sampling sites (CEL, DZI, YAL, and CHA) were significantly dif ferent from each other.The PERMANOVA and the generalized discriminant analysis, computed with the  4. Pairwise comparisons using linearized F ST differentiation (above diagonal) and multilocus Nei's genetic distances (Ds) and in parentheses, the real shortest distance by sea (km) (below diagonal) among sampling locations.Linearized F ST values in bold were significant at p = 0.001 otolith shape data, showed phenotypic differ entia tion among the samples.According to Friedland & Reddin (1994), a classification success of > 75% is considered acceptable to use otolith shape analysis as a discrimination tool.In the present study, 92 to 96% of the specimens were correctly classified within each sample, suggesting that the otolith shape-variation analysis can be a feasible tool to evaluate structuration of this species at a regional level.This type of geographic variation in otolith shape has also been observed at different spatial scales among different fish populations, such as the Atlantic mackerel Scomber scombrus (Caston guay et al. 1991), cod Gadus morhua (Campana & Casselman 1993), Atlantic salmon Salmo salar (Friedland & Reddin 1994), haddock Melanogrammus aeglefinus (Begg & Brown 2000), and horse mackerel Trachurus trachurus (Murta et al. 1996, Stransky et al. 2008).Variations in the shape of the otolith in fishes are related in part to genetic differences (Castonguay et al. 1991, Friedland & Reddin 1994) as well as to the influence of the environment.For instance, it has been shown that in similar habitats, fishes produce otoliths that have similar patterns in growth and shape (Campana & Casselman 1993, Villegas-Hernández et al. 2008).Thus, evaluating otoliths of the same species of fishes in different landscapes and at different distances can provide reliable results to evaluate any population structure.
On the other hand, with data from the microsatellite, we detected significant genetic structure with only 4 sampling sites (F ST = 0.0560, p < 0.001).The linearized F ST values ranged from 0.029 to 0.069.These values were higher than those found on French grunts Haemulon flavolineatum (F ST = 0.003), bicolour damselfish Stegastes partitus (F ST = 0.003), and bluehead wrasse Thalassoma bifasciatum (F ST = 0.0002) collected in the whole Caribbean Sea (Purcell et al. 2006(Purcell et al. , 2009)).Similar to French grunts, white grunts have an isolationby-distance pattern of genetic differentiation, suggesting that their dispersal distance is limited.Grunts, including H. plumieri, have relatively short planktonic larval durations (14 to 15 d on average) (Lindeman et al. 2001).This may explain the shorter dispersal distances for these species (Purcell et al. 2006).According to Palumbi (2003), evidence of genetic structuring on small spatial scales could indicate low dispersal among populations.
Despite the genetic differentiation through isolationby-distance over hundreds of kilometers between the YAL-CHA and CEL-DZI samples, in future studies, it could be necessary to evaluate if the circulation of ocean currents influences larval transport and recruitment with the genetic consequences of large variance in individual reproductive success to generate patterns of patchy genetic structure at fine scales (Cowen et al. 2000).In contrast to terrestrial environments, where topographic barriers occur, discon - tinuities in marine populations are caused mainly by meso-scale and submeso-scale oceanic circulation that minimizes larval dispersal over large distances (Cowen et al. 2000).This causes genetic differences among subpopulations (Carvalho & Hauser 1994).
Both Doherty et al. (1995) and Riginos & Victor (2001) found that there is a significant positive correlation between genetic differentiation and planktonic larval duration and concluded that marine species with long pelagic larval stages (> 40 d) usually have low levels of genetic differentiation among populations because their larvae have a greater chance of migration with a major chance of gene flow, contributing to a healthy gene pool.Thus, the relatively short planktonic larval duration of Haemulon plumieri and the population structure observed in the present study by using genotypic and phenotypic approaches may indicate that low gene flow and differentiation between the studied sub-populations could likely be caused by the species low larval dispersal and by the currents.Marine populations may exhibit chaotic genetic heterogeneity at small spatial scales due to variation among recruiting cohorts (i.e.several pulses of recruits) (Hedgecock 1994a,b).A cohort of progeny produced by a single episode of spawning could have less allelic diversity than the adult population, causing larger allele-frequency differences among cohorts of juveniles than the corresponding allele-frequency differences among adult populations (Hedgecock & Pudovkin 2011).As we only sampled a single cohort of recently settled juveniles in the present study, we cannot explore this possibility; however, with H. flavolineatum which have a similar life history, no evidence of significant temporal variability or differences between juveniles and adults were found (Purcell et al. 2006).

CONCLUSIONS
We found similar patterns by genotypic and phenotypic spatial heterogeneity in juveniles of Haemulon plumieri along the northern coast of the Yucatan Peninsula.This species supports an important local fishery in southern Mexico, and although currently, there are no data of its overfishing status on the North coast of the Yucatan Peninsula, its long life cycle characterized by slow growth rates and late maturity (Murie & Parkyn 2002) may make this species vulnerable to local overfishing.Their structuration at a relatively small scale might be related to their low larval dispersal as well as the ocean currents.In future studies, it is necessary to analyze juvenile and adult samples collected over a period of several years and over a larger spatial scale to corroborate our findings.Also, it is important to evaluate the fishing effort and mortality of H. plumieri and to improve suitable management strategies to preserve the whole gene pool variability of the population in this study area and on a large scale.Results from the present study provide baseline information on the genetic status of H. plumieri populations that are important for southern Mexico.

Fig. 3 .
Fig. 3. Graphical method allowing detection of the true number of groups K*.Mean log probability values [lnP(K )] for each run from K = 1 to 8. The mean likelihood L(K) (±1 SD) over 10 runs for each K value (black dots) and the difference of the rate of change of the likelihood function with respect to K is also shown as L'(K ).The modal value of ΔK (gray bars) calculated as ΔK = m|L′′(K)|/s[L(K)], where m is the mean and s is the SD, represents the true K* or the uppermost level of structure (4 clusters)

Table 2 .
Number of alleles (Na), observed heterozygosity (Ho), expected heterozygosity (He), and inbreeding coefficient (F IS ) across all samples and all loci (see Table1for abbreviations)