High Levels of Genetic Connectivity among Populations of Yellowtail Snapper, Ocyurus chrysurus (Lutjanidae – Perciformes), in the Western South Atlantic Revealed through Multilocus Analysis

In the present study, five loci (mitochondrial and nuclear) were sequenced to determine the genetic diversity, population structure, and demographic history of populations of the yellowtail snapper, Ocyurus chrysurus, found along the coast of the western South Atlantic. O. chrysurus is a lutjanid species that is commonly associated with coral reefs and exhibits an ample geographic distribution, and it can therefore be considered a good model for the investigation of phylogeographic patterns and genetic connectivity in marine environments. The results reflected a marked congruence between the mitochondrial and nuclear markers as well as intense gene flow among the analyzed populations, which represent a single genetic stock along the entire coast of Brazil between the states of Pará and Espírito Santo. Our data also showed high levels of genetic diversity in the species (mainly mtDNA), as well a major historic population expansion, which most likely coincided with the sea level oscillations at the end of the Pleistocene. In addition, this species is intensively exploited by commercial fisheries, and data on the genetic structure of its populations will be essential for the development of effective conservation and management plans.


Introduction
Preservation of the biological diversity of any ecosystem is essential for evaluation of the distribution and connectivity of its populations [1] and the factors that determine these patterns. Considering the marine environment, opportunities for isolation to occur between populations are rare [2][3][4]. Many marine fish species tend to present a high degree of genetic connectivity, despite being distributed over thousands of kilometers of ocean, although this is often attributed to the from commercial fishing ports (obtained from the near-shore artisanal fishery) and the fishing with rafts, performed by local fisherman. The specimens were identified based on the specific literature [10,18].

Ethics Statement
All specimens were obtained from dead individuals, procured through direct purchase from commercial landings in the localities mentioned above. Ocyurus chrysurus is not endangered or protected along the Brazilian Coast. Therefore, there was no need to apply for a license for collection or approval by the Animal Ethics Committee. The specimens were transported with the authorization of the Brazilian Environment Ministry.

Laboratory procedures
In the laboratory, tissues samples (from the muscle, fin or tongue) were extracted from each specimen and frozen until analysis. All of the specimens were included in the Lutjanidae tissue bank of the Applied Genetics Laboratory at the Bragança, Campus of the Federal University of Pará (UFPA). Their genomic material was obtained using a phenol-chloroform and enzymatic extraction protocol and was precipitated with sodium acetate, isopropanol and ethanol [19].
The amplicons were purified with PEG 8000 (polyethylene glycol) following the protocol [25] and sequenced via the dideoxy-terminal method [26] using the reagents of the Big Dye kit (ABI Prism Dye Terminator Cycle Sequencing Reading Reaction-PE Applied Biosystems, Carlsbad, CA, USA). The precipitate was sequenced through capillary electrophoresis in an ABI 3500 automatic sequencer (Applied Biosystems).
Initially, only one of the strands of each genomic region was sequenced (see Table 1). When ambiguities were observed in the chromatograms, the sample was sequenced in both directions and/or twice in the same direction (especially in the intragenic regions) to avoid errors in the identification of heterozygous individuals. For the intron of Insulin-like Growth Factor (IGF), it was necessary to design an internal primer (IGF B-5'-CATTGATATTCCTGNTCGTTCA-3') to obtain sequences in both directions. All haplotypes were deposited in Genbank under the accession numbers KM596919 to KM507050, except for GH-5, because that fragments smaller than 200 bp are not accepted in Genbank (these sequences are listen in S1 Supporting Information).

Database
The DNA sequences were edited and aligned in BIOEDIT v. 7.1.3.0 [27]. In the case of the nuclear loci, the individual heterozygotes were detected when double peaks were observed at the same position in both directions in the chromatograms. Heterozygotic events caused by indels were diagnosed through visual analysis of the chromatograms, with the alleles being reconstructed in INDELLIGENT v. 1.2. (http://imperialis.inhs.illinois.edu/dmitriev/indel.asp) [28].
The gametic phase of each nuclear marker was defined using the PHASE algorithm [29], available in DNAsp v 5. 10.01 [30]. The runs consisted of 1,000 burn-in iterations and 1,000 principal iterations, with a thinning interval of 1. The algorithm was applied five times, with the fifth chain being ten times longer than the others. The haplotypes that returned a probability of less than 0.8 were excluded from the analyses.
For the nuclear data, the minimum number of recombination events was estimated via the Rm method [31], available in DNAsp v 5. 10.01 [30]. As the results of this analysis are strongly affected by homoplasy [32], the significance of the number of recombination events was evaluated using the Ф W test [32], available in SPLITS TREE v. 4.12.6 [33]. Linkage disequilibrium was analyzed the EM algorithm in ARLEQUIN v. 3.5.1.3 [34], which was run five times, with 20,000 permutations.

Characterization of genetic diversity and levels of gene flow
Determination of the number of polymorphic sites and the identification of possible stop codons (in the case of codifying regions) were performed in MEGA 6 [35]. The identification, quantification, and distribution of the haplotypes were determined in DNAsp v 5. 10.01 [30]. The genetic variability of the populations was evaluated based on the haplotype (h) and nucleotide (π) diversity indices [36] obtained from ARLEQUIN v. 3.5.1.3 [34]. The haplotype network was generated using HAPLOVIEWER [37] based on a maximum parsimony tree produced in DNAPARS, available in the package PHYLIP v. 3. 6 [38], in accordance with Salzburger et al. [37].
The genetic homogeneity of the O. chrysurus populations was initially evaluated through an analysis of molecular variance (AMOVA) [39] for each marker individually and subsequently through a multilocus approach, both with 10,000 permutations. This analysis permitted partitioning of the results into within-and between-population variation. In addition, F st values [40] were used to evaluate the gene flow between pairs of populations. These analyses were run in ARLEQUIN v. 3.5.1.3 [34], with a subsequent adjustment of the p values using the False Discovery Rate test [41]. To verify the existence of isolation by distance, Mantel tests were performed using a matrix of genetic (Fst/(1-Fst)) and geographic (km converted to Ln) distances [42]. Negative Fst values were expressed as zero. These analyses were conducted in IBDWS (http://ibdws.sdsu.edu/*ibdws) [43], with 10,000 permutations.
For comparison between Brazilian and Caribbean populations, the control region sequences used by Vasconcellos et al. [16] (accession numbers EF624354-EF624359; EF624361-EF624373) were included in the network, AMOVA, and pairwise Fst analyses as well the Mantel test.
Bayesian methods, using STRUCTURE v. 2.3.4 [44], were applied to assign individuals to populations. This procedure places individuals into K clusters, where K is chosen in advance but can be experimentally varied throughout independent runs. K values between 1 and 8 were tested, using a model with admixture and no locprior (i. e., only genetic data is used for the assignment of individuals to a given K). For this analysis, only nuclear data were employed. Each run consisted of 1,000,000 steps (burn-in = 10%, and each value of K was implemented 10 times). The number of K was inferred by comparing the mean values of Ln Prob obtained in Structure Harvester (http://taylor0.biology.ucla.edu/structureHarvester/) [45].
Cluster analyses were also conducted in STRUCTURAMA [46]. For this analysis, the mitochondrial and nuclear data were grouped. The runs consisted of 2,000,000 generations, (burnin = 20%). For the values of K, we employed the following distribution (K = expk (2)). The runs were summarized using the "showtogetherness" command.
To check the fit of the historical population dynamics to a model of exponential growth, we used a mismatch distribution [47] together with the SSD and raggedness indices. Mismatch analyses were conducted in DNAsp v 5 10 01 [30], rates of SSD and raggedness, were implemented in Arlequin v. 3.5.1.3 [34] based on 10, 000 permutations.
Historic fluctuations in the demography of O. chrysurus were visualized using a Bayesian Skyline Plot (BSP) [48] and an extended Bayesian Skyline Plot (EBSP) [49]. These procedures were run in BEAST v.1.7.4 [50], based on evolutionary models suggested by JMODELTEST 2.1.1 [51,52] (HKY for Cytb, IGF 2, GH 5, ANT 1 and HKY+ I + G for CR). The analyses were based on the strict molecular clock used for the teleost control region, with a substitution rate of 10% per million years [53,54]. Two runs were performed using different random seeds, including 200 million generations for each BSP run and 400 million for each EBSP run, with samples taken at intervals of 10,000 generations, 10% of which were discarded as burn-in. The convergence and mixing of the chains were inspected visually in TRACER v.1.5 [55]. The convergence and mixing were considered to be appropriate when all of the ESS values for each of the parameters analyzed were above 200.
Tajima's D [56] and Fu's Fs [57] were also calculated, given that in addition to the detection of possible deviations from neutrality, these values may be used to evaluate demographic patterns, such as population expansion. These analyses were run in ARLEQUIN v. 3.5.1.3 [34], with their statistical significance being assessed using 10,000 permutations.

Mitochondrial DNA
A total of 602 base pairs (bp) was sequenced from the Control Region in 152 O. chrysurus specimens, and 645 bp of the Cytochrome b gene was sequenced in 170 specimens ( Table 2). Considering the respective evolutionary rates of the two markers, similar patterns concerning their distribution and haplotype frequencies were observed in the different populations examined. In the CR, which includes 93 polymorphic sites, a total of 91 haplotypes were identified, the most common of which was shared by 27 specimens and was present at all localities except Bahia. All other CR haplotypes were either unique or occurred at low frequencies and were distinguished by a small number of mutations. Only 12 polymorphic sites were found in Cytochrome b; however, a total of 12 haplotypes were identified, two of which were very common, being shared by 74 (44%) and 67 (39%) of the specimens and being found in all of the populations analyzed.
The indices of haplotypic diversity were high in the CR (h = 0.963±0.010) and lower for Cytb (h = 0.653±0.002), and the same pattern was observed in the case of nucleotide diversity, with π = 1.7% for the CR, but only π = 0.15% for Cytb (Table 2). AMOVA indicated that most of the variation in both markers occurs within populations, rather than between them, with low and non-significant Ф ST values being obtained (Table 3), and this finding was further corroborated by the non-significant F st values obtained in the pairwise comparisons between populations (Table 4). However, the comparison between the populations from the Brazilian coast and Belize revealed that approximately 30% of the variance is explained by differences between these two regions ( Table 3). The Fst values between the Brazil and Belize populations were greater than 0.20 and were highly significant for all comparisons (Table 4), indicating particularly high differentiation between these stocks.
The genetic homogeneity of the populations from the Brazilian coast was also emphasized by the distribution of the haplotypes, given the lack of any clear geographic pattern in the network (Fig. 2). This genetic homogeneity over a broad geographic scale was further supported by the Mantel test, which rejected the scenario of isolation by distance (IBD), although when the comparisons included Belize there was some evidence in favor of a scenario with IBD, (Mantel test; r = 0.7354, p = 0.05) (S1 Fig.). The haplotypes identified in Belize were not shared by any locality analyzed along the Brazilian coast, supporting the significance of the Mantel test.
With regard to the neutrality of the data, the obtained Fs values were significant in some populations, or when the specimens were grouped in a single population for both mitochondrial markers ( Table 2). The D values were not significant for any population ( Table 2).
The Bayesian Skyline Plot (Fig. 3) indicated the historic occurrence of an increase in the effective size of the O. chrysurus populations, dated to the end of the Pleistocene. These results are consistent with a process of historical expansion of yellowtail snapper populations, as indicated by significant negative Fs values [57] and by adjusting the mismatch distributions to model population growth (Fig. 4).  (Rm = 8; p = 0.03). When the indel region was removed, recombination was not detected (p = 0. 25), and all of the analyses for this marker were performed excluding this region. The highest diversity was recorded in the IGF 2 sequences (h = 0.785 ± 0.017; π = 1.4%), with intermediate values (h = 0.510±0.031; π = 0.57%) being recorded for GH 5 and considerably lower values being obtained for ANT 1 (h = 0.112±0.024; π = 0.03%) ( Table 2).

Nuclear DNA
The results of the AMOVA revealed that most of the genetic variance was present within, rather than between the populations, which is consistent with high levels of genetic similarity between populations (Table 3). This finding was reinforced by the frequency and distribution  of the identified haplotypes (Fig. 2) as well as by non-significant values of the pairwise Fst (Table 4) and Mantel tests (S1 Fig.). Likewise, the Bayesian clustering analysis showed that a scenario with only one group is the most likely (S2 and S3 Figs.). With regard to the neutrality indices (Table 2), the obtained values were negative and not significant. The mismatch distributions were unimodal for all markers, except for IGF 2, which presented a clearly bimodal distribution (Fig. 4). This pattern can be observed in populations that have experienced long intervals with stable sizes or that have undergone a subtle decrease in their effective size, followed by an event of expansion of range, or in populations who have experienced a weak bottleneck [58]. The EBSP analysis indicated the occurrence of an increase in the effective size of the O. chrysurus populations (Fig. 3), supporting the results of the BSP for the mitochondrial Control Region.

Genetic structure in the yellowtail snapper
The pattern of genetic homogeneity observed in O. chrysurus in the present study is similar to that found in other lutjanids, such as Lutjanus kasmira [59], Pristipomoides filamentosus [60], L. campechanus [6], and L. purpureus [5]. All of these species exhibit a high dispersal capacity during at least one phase of their life cycle. However, a pattern of significant genetic structure has been observed for other species of this Family, such as Lutjanus erythropterus [53], Lutjanus fulvus [59], and Lutjanus synagris [61].
Differences in larval behavior as well as the characteristics of the environment where these species live may explain these discordant patterns of genetic connectivity and appear to directly influence the evolutionary history of these fish. However, it is important to note that O. chrysurus is not genetically homogeneous throughout the whole of its geographic distribution. Vasconcellos et al. [16] found marked differences between the Caribbean and Brazilian populations of O. chrysurus, while Saillant et al. [15] detected genetic sub-structuring in the populations of these species located in neighboring areas of the Caribbean and Gulf of Mexico. In general, these data are consistent with the existence of effective barriers in the ocean between northern Brazil and the Caribbean, (e.g., circulation patterns and ocean currents), as identified in previous phylogeographic studies [62][63][64], as well as between the Caribbean and the Gulf of Mexico [15,65]. Nevertheless, regarding genetic differentiation involving Brazilian populations, our results demonstrate that the model of isolation by distance cannot be discarded.
Juvenile and adult yellowtail snappers tend to remain within their area of settlement over the course of their lives [13,14]. This characteristic, in addition to specific habitat features, such as the pattern of currents, may account for the genetic differentiation observed in previous studies, with limited gene flow occurring between populations in some areas [15].
However, the larvae of O. chrysurus display a strong swimming ability [66], and they are most likely able to travel long distances, which could lead to intensive mixture of individuals. Moreover, these larvae present a planktonic phase of approximately one month and an offshore distribution, and during this larval stage, they preferably inhabit surface waters, which are more subject to the influence of ocean currents [12]. The dynamics of ocean surface circulation in the south of the western Atlantic region is primarily a result of bifurcation of the South Equatorial current [67]. Some studies, however, demonstrate that the bifurcation of the south equatorial current is not an effective barrier to dispersal for many marine taxa [68,69], furthermore, the main currents parallel to the Brazilian coast (i.e. Brazil Current / North Brazil Current) show slight seasonal changes in direction (see http://www.aoml.noaa.gov/phod/graphics/ dacdata/seasonal_brazil.gif) that could allow some connectivity between the locations analyzed. Thus, all of these features are consistent with high genetic similarity along the Brazilian coast (coast Pará to Espírito Santo), as demonstrated in this study, and they may be the main factors regulating the genetic connectivity between the Brazilian populations of O. chrysurus, as observed in another lutjanid, L. purpureus [5].

Genetic Diversity
The results of the present analysis revealed high levels of genetic variation in the investigated yellowtail snapper populations, especially in the mitochondrial markers, which presented a large number of haplotypes, similar to that recorded for L. purpureus [5] and L. campechanus [6], with similar genetic variability demonstrated between localities, over a distance of approximately 3,000 km along the Brazilian coast. This outcome would be expected for populations linked by intense gene flow, a scenario commonly observed in a number of marine species [70]. The high indices of genetic diversity recorded in the present study, especially for the mitochondrial markers, and particularly the Control Region, are a common feature of marine teleosts, including lutjanids [5,6,53,71,72], and have been recorded previously in O. chrysurus [16].
This high genetic variability cannot be interpreted as an absence of an impact of fishing on these populations because the commercial exploitation of this species is relatively recent (approximately 50 years and less than 20 generations). Examples from the literature reveal several situations where overfishing apparently shows no direct relationship with genetic diversity indices. For example, Hoarau et al. [73] did not report decreased levels of genetic variability in microsatellite markers in a temporal analysis of Pleuronectes platissa covering almost 100 years, even though this species has been heavily exploited since the XIX century. Additionally, Pinsky & Palumbi [74] recently performed a meta-analysis involving hundreds of fish species and observed high levels of genetic diversity for snappers of the genus Lutjanus, even for species that are considered to be overfished.
Distinct levels of genetic diversity were observed at the three nuclear loci analyzed in the present study. The most variable nuclear locus was the intron of the IGF 2 gene, which presented a degree of polymorphism comparable to that observed in the mitochondrial Control Region, indicating that it is a potentially useful marker for population-level and phylogeographic studies in lutjanids, as observed for Centropomus [75].

Demographic History
The yellowtail snapper inhabits coastal waters [10] and is therefore relatively susceptible to sea level oscillations [76]. During the last glacial maximum, for example, approximately 90% of the present-day continental shelf of the Caribbean region was above sea level, due to a decrease in the sea level of 120 meters [77]. This situation almost certainly led to the extinction of species or lineages, as well as contraction processes and population expansion [78,79,80].
The temperature fluctuations that occurred in this period may have had a great influence on the historical demography of a number of marine species from this region. For example, Rocha et al. [81] lists temperature fluctuations that occurred in the Atlantic (see Sachs et al. [82]) as one of the events responsible for population growth in species of the Chromis genus. O. chrysurus generally prefers water at temperatures between 24-30°C, with temperatures above 34°C being lethal for this species [83]. Thus, temperature changes beyond these thresholds could also lead to strong fluctuations in the effective size of O. chrysurus populations.
Historic events of population expansion are commonly identified in demographic studies on marine fish species. In the western Atlantic, events of this type have been recorded in Cynoscion guatucupa [84] and Chromis multilineata [81] as well as in lutjanids, such as the snappers L. purpureus [5] and L. campechanus [6].
The results obtained for the O. chrysurus stocks analyzed in the present study clearly indicate population expansion along the Brazilian coast. This process is similar to that observed in Caribbean O. chrysurus, although Vasconcellos et al. [16] concluded that the Brazilian populations were demographically stable because all tests neutrality ("Fs" and "D") were unable to reject a neutral model of evolution. In this case, it is possible that the discrepancies in relation to the previous study are related to intrinsically different sensitivities of the applied neutrality tests. Finally, Fs test, was significant only for the total population for the mitochondrial markers examined in the present study. Fu's test shows greater statistical power and sensitivity for detecting events of geographical expansion [57] and is therefore more suitable for supporting the expansion scenario indicated by the mitochondrial markers.
The expansion of the Brazilian populations of O. chrysurus is also clear from the BSP and EBSP analyses, which identified an abrupt increase in population size approximately 25,000 years ago, which would coincide with the last glacial maximum in this region (see Barreto et al.) [85]. In fact, the recent Pleistocene glaciation events are well documented in the marine environment [82,86], and the marked changes in water temperatures and sea levels almost certainly had great effects on the historical dynamics of the population size of O. chrysurus as well as other marine organisms

Final considerations
The present study is the most comprehensive investigation of the Brazilian populations of Ocyurus chrysurus undertaken to date. The results clearly revealed the existence of a single panmictic population along the Brazilian coast (sensu Selkoe et al.) [87], and as reported by Vasconcellos et al. [16], there is a high level of substructure between the Brazilian and Caribbean populations. However, it is worth mentioning that the absence of individuals sampled between Brazil and Belize makes it difficult to perform further analyses of the historical connectivity between populations in this geographic range.
In relation to the high levels of genetic diversity observed, mainly among mitochondrial markers, these findings cannot be interpreted as a lack of any impact of fisheries on these populations, given that the commercial exploitation of this species is relatively recent, and the markers used in the present study would not be appropriate for the detection of alterations on this time scale (see Gaither et al. [88]). Furthermore, the deterioration of genetic diversity due to overfishing is more clearly demonstrated in populations that exhibit a low effective size (e.g., a few hundred individuals) and low migration rates [74], which may not be the scenario for Ocyurus chrysurus along the Brazilian coast.
The data presented herein suggest that the populations of yellowtail snapper along the Brazilian coast may represent a single genetic stock, and if so, they should be managed as a unit. This conclusion has strong implications for the management of this species. For example, for various snappers, including O. chrysurus, the formation of shoals for spawning is reported to occur at a number of locations [89]. Thus, the protection of these areas should be given a higher priority because they apparently have a greater impact on the population structure over a vast area.
Further studies involving other classes of molecular markers (e.g., microsatellite markers, SNPs, adaptive loci) as well as studies tracking adults and larval dispersal should be conducted to obtain a better understanding of the structure of O. chrysurus populations.