Phylogeography and Historical Demography of Two Sympatric Atlantic Snappers: Lutjanus analis and L. jocu

Lutjanus analis (mutton snapper) and L. jocu (dog snapper) are mesopredator species with extensive geographic distribution in the Atlantic Ocean. Although historically overfished, their genetic diversity, population structuring, and historical demography along the Brazilian coast are unknown. Here we present genetic data for the hypervariable region 1 (HVR1) of the mtDNA control region of both L. jocu and L. analis, and for cytB of L. analis from distinct geographic regions. Phylogeographic analyses based on HVR1 sequences revealed no geographic structuring of mtDNA lineages for either species. The bimodal mismatch distribution plots of mutton and dog snapper populations implied that they might have experienced secondary contact. Historical demographic parameters estimated that population expansions ranged from 13,631 to 61,171 years before present (ybp) for L. analis and 36,783-55,577 ybp for L. jocu, associated with events that occurred at the end of the last glaciation period. Estimates of the average effective population size for L. jocu were higher than those for L. analis, with the largest population occupying the Brazilian northeastern region coast. High migration rates are maintained among the three northernmost locations, with a substantial decrease to the further southern region. Our study suggests that a tripartite interaction of larval dispersion and interregional adult movement (aggregate spawning), allied to historical contingencies, contributed to contemporary population genetic patterns of these species, and adds relevant information for conservation management of such vulnerable - and valuable - marine resources.


INTRODUCTION
Genetic data are critical for developing strategies to delimit populations and conserve marine species. Identifying geographic borders or diagnostic characteristics for distinct genetic subunits play decisive roles in defining threatened species and proposing conservation units. This information has a direct impact on conservation and management policies (Morin et al., 2010). Nevertheless, spatial distribution patterns of genetic diversity remain scarce for most commercially exploited species (Reiss et al., 2009). These include important marine resources, such as fishes of the family Lutjanidae, which are captured over the entire area of its geographic distribution (Allen, 1985;Burton, 2002;Graham et al., 2008).
Lutjanidae species, known as snappers, occur in the Eastern and Indo-Western Pacific, Eastern, and Western Atlantic. Many species have an extensive geographical distribution, although a few species have restricted ones (Allen, 1985;Lindeman et al., 2001;Moura and Lindeman, 2007). They are demersal fishes of medium to large size, with some species reaching more than one meter, occupying very shallow environments up to 500 m deep (Allen, 1985). Lutjanid fishes constitute one of the most valuable fisheries categories, which contributes to their considerable exploitation. Additionally, they also represent an important food resource for communities dependent on artisanal fishing (Rezende et al., 2003).
Twelve species of the genus Lutjanus have been intensely exploited along the Brazilian coast (Rezende et al., 2003;Klippel et al., 2005;Frédou et al., 2009). Among these fishing targets are Lutjanus analis (Cuvier, 1828) and Lutjanus jocu (Bloch and Schneider, 1801), mainly found in coral reefs, and on the sandy bottoms of bays and estuaries, from Massachusetts (United States) to southeastern Brazil, including the Caribbean and Gulf of Mexico (Allen, 1985).
Accordingly, the Brazilian Lutjanus' genetic patterns are still unknown, hindering appropriate conservation strategies for populations of this genus. This is a particularly necessary and urgent task, since it is estimated that the capture level of the species exceeds nearly 90% of that considered adequate (Frédou et al., 2009). In fact, L. analis is being overfished on the northeast coast of Brazil and is considered a near-threatened species (Frédou et al., 2009;IUCN, 2019).
In order to ascertain the levels of connectivity between populations of L. analis and L. jocu along the Brazilian coast, sequences of the hypervariable region 1 (HVR1) of the mtDNA control region of both species, and also of the cytochrome b (cytB) gene of L. analis were analyzed, comprising wide areas from northeastern to southeastern Brazilian coast. These data proved to be useful for assessing the stock and evolutionary dynamics of these particularly vulnerable species.

Biological Material
Lutjanus analis and L. jocu individuals were collected through artisanal fishing in four areas of the Brazilian coast; three in the NE region (Ceará -CE, Rio Grande do Norte -RN and Bahia -BA) and one in the SE region (Espírito Santo -ES), covering a representative zone of the species' distribution in the South Western Atlantic. Eighty one L. analis individualsthe mutton snapper -and one hundred L. jocu -the dog snapper -were used in phylogeographic analyses ( Table 1).
Fragments of muscle or liver tissues were placed in 2.0 ml microtubes containing 95% ethyl alcohol and stored at a temperature of −20 • C.
Dambe software (Xia and Xie, 2001) was used to determine the possible existence of saturation involving transitions and transversions. Haplotype (h) and nucleotide (π) diversities (Nei, 1987), as well as genetic differentiation between localities, using pairwise F-statistics were estimated using Arlequin 3.5.1.2 (Excoffier and Lischer, 2010). F st values were tested for significance with 10,000 permutations. Analyses of molecular variance (AMOVA) (Excoffier et al., 1992) were also implemented to determine possible structuring between the different geographic regions analyzed. In addition, Spatial Analysis of Molecular Variance (SAMOVA) was used as an approach to define groups of populations that were geographically homogeneous and also maximally differentiated from each other (Dupanloup et al., 2002).
The distribution of haplotype diversity in each population was assessed with a haplotype network for each DNA marker built in TCS (Clement et al., 2000). We used TCSbu (Santos et al., 2015) for graphical representation of the obtained networks.
Mismatch analyses (frequency of pairwise nucleotide site differences between sequences) were carried out to estimate the historical demographic parameters τ, θ 0 and θ 1 , according to the sudden population expansion model (Rogers and Harpending, 1992). Adjustments to an unimodal distribution of pairwise nucleotide differences among sequences, as expected for expanding populations, was tested using Harpending's raggedness index (r). Non-significant r indices suggest a good fit to the growth-decline model, while significant indexes are indicative of a stable population (Harpending, 1994). The sum of square deviations (SSD) between the observed and expected distributions was also considered (Schneider and Excoffier, 1999). Additionally, Fu's F S (Fu, 1997) and Tajima's D neutrality tests (Tajima, 1989), with 10,000 permutations, were used to infer historical demographic variations.
Expansion times expressed in units of mutational time τ were translated to time since expansion as τ = 2ut, where 2u is the mutational rate per site per million years. Similarly, θ estimates (θ 0 = 2N 0 µ and θ 1 = 2N 1 µ were used to estimate the female effective population sizes before (θ 0 ) and after (θ 1 ) expansion, where µ is the mutation rate (Rogers and Harpending, 1992;Rogers, 1995). Confidence intervals of parameters were calculated by applying a parametric bootstrap approach (Schneider and Excoffier, 1999). Values of τ were transformed in real-time by using the estimated expansion time, through the equation u = µmt, where u is the mutational rate by generation, µ (mi) the estimated mutation rate for the sequence analyzed, and mt the number of nucleotide bases involved. The expansion time interval was calculated as τ = 2ut, where t is the estimated time from the occurrence of expansion. To that end, a generation time of 4 years was used for L. analis (Burton, 2002) and 4.6 years for L. jocu, established as means of different estimates ranging from 2 to 6 years (Ault et al., 1998, Froese andBinohlan, 2000).
Estimates of contemporary migration rates among populations for each species were obtained using Migrate-n 4.2.14 (Beerli, 2006;Beerli and Palczewski, 2010). Based on preliminary runs, all priors were set as default, except for theta in L. analis, which was draw from a uniform distribution with mean 0.01, minimum 0.00 and maximum 0.10; and migration rate in L. jocu, that was set with mean 500.0, minimum 0.0 and maximum 5000.0. The analyses were conducted with two independent runs of three long chains over 30,000,000 iterations each, after discarding the first 50,000 steps as burn-in. Posterior samples were recorded every 500 steps. We used a static heating scheme with temperatures 1.0; 1.5; 3.0 and 10,000.0. Convergence for each parameter was evaluated with effective sample size (ESS) values.

Molecular Diversity -cytB and HVR1
Sequencing of the cytB gene from 73 L. analis individuals resulted in the production of a 667 bp sequence. We identified 19 distinct haplotypes differentiated by 20 polymorphic sites, all involving DNA sequence transitions ( Table 2). The sequences showed no evidence of saturation between transitions and transversions (ts/tv). The haplotype and nucleotide diversity varied between 0.88 and 0.61, and 0.001 and 0.005, respectively ( Table 2).
PCR amplification of the HVR1 of L. analis and L. jocu produced 493 bp and 394 bp, respectively, analyzable sequences. This marker in both species also showed no evidence of saturation involving t s /t v . Only five L. analis and nine L. jocu haplotypes were shared within and between the different sites. The amount of nucleotide bases followed the A > T > C > G pattern (Sbisà et al., 1997). Geographic samples showed high genetic variability in terms of haplotype (h) and nucleotide (π) diversity, whose values varied from 0.985 to 1.0, and 0.029 to 0.033, respectively, in L. analis, and from 0.991 to 1.0, and 0.035 to 0.038, respectively, in L. jocu ( Table 2). The L. analis populations of CE, BA, and RN were more polymorphic than the ES population, while nucleotide diversity was similar in all L. jocu populations.
The reduced number of haplotypes shared within and between populations reflects the high h and π values observed, which is expected for populations without genetic depression.

Population Genetic Structure
F st indexes based on HVR1 sequences indicated panmixia between the distribution areas for both L. analis and L. jocu (Table 3). However, the cytB sequences suggested a possible population structure of L. analis in relation to RN and ES (F st = 0.19, p = 0.005). Given the absence of population structure suggested by F st estimates concerning HVR1 sequences, AMOVA of samples from different localities were assessed as a single population. AMOVA analyses in both species indicated that genetic variation is majorly found within populations, with no variance among groups ( Table 4). The absence of valid clusters was also supported by SAMOVA.

Haplotype Networks
Haplotype networks for HVR1 in both species were similar with several singletons present in all populations, thus reflecting the high variability of this sequence and the absence of population differentiation. L. analis cytB network showed a lower number of haplotypes, with a high-frequency central haplotype present in all populations, and with other haplotypes with intermediate frequencies also appearing in several populations and in a few singletons located at the tips of the network (Figure 1).

Migration Rates
Estimated θ values (an estimator of effective population sizes, scaled by mutation) for L. analis, showed higher values in CE and BA populations, while RN and ES showed much lower estimates. In L. jocu, higher values were found in CE, RN and BA and a smaller sizes in ES (Figure 2). Migration rates among L. analis populations presented fewer immigrating individuals per generation in ES, and higher immigration in BA and CE. L. jocu populations presented similar gene flow patterns, with lower immigration rates observed in ES, while a high number of migrants per generation was detected from BA to CE, from RN to BA, and between CE and RN (Figure 2).

Demographic History
Analysis of deviations from neutrality based on Tajima (Table 4). Raggedness (r) indices and the sum of square deviations (SSD) performed on L. analis data were not significant (P > 0.05), also supporting population expansion for   this species. Mismatch distribution graphs (Rogers, 1995) showed a typical bimodal pattern in both species. Associated with other indicators, the bimodal mismatch distribution can be indicative of secondary contact during the period of population expansion of the species (Figure 3).
Coalescence data based on t indicated a period of population expansion between 13,795 and 61,100 ybp, for L. analis and between 36,783 and 55,577 ybp for L. jocu ( Table 5). This parameter was not calculated for ES due to no indication of expansion.

Population Structure of the Mutton and Dog Snappers in Brazilian Coastal Areas
Genetic congeneric comparisons provide an adequate phylogenetic model for estimating contemporary and historical events on population dynamics and genetic population structure. On the Brazilian coast, the species L. analis and L. jocu, present sympatric distribution, and share several biological and ecological characteristics (Allen, 1985), offering an excellent model for comparisons of spatial genetic patterns from historical and phylogenetic perspectives.
Along the Brazilian coast, L. analis and L. jocu can be considered as a unique panmictic unit, suggesting shared evolutionary responses to similar environmental pressures and historical events to which they were subjected. In fact, the F st values, based on the comparative analysis of the mitochondrial DNA HVR1 control region sequences, do not indicated genetic structuring in L. analis and L. jocu. Analyzes based on single locus markers can be prone to peculiarities of the selected marker (Hurst and Jiggins, 2005;Allio et al., 2017), and some patterns need be assessed by using other unlinked markers (Godinho et al., 2008). Nevertheless, analyzes with cytB sequences suggested some genetic structuring of the ES population of L. analis, located in the southernmost sampled site.
The level of genetic structuring in the marine environments seems to be the result of complex historical interactions due to interregional connectivity, physical barriers, restricted migratory behavior, short pelagic larval phase, N ef size, species ecology, and ocean current actions (Planes et al., 2001;Stepien et al., 2001). In contrast, panmixia is the predominant pattern challenging these factors in a large number of marine species (Liu et al., 2016(Liu et al., , 2017, such as other lutjanids of the Brazilian province as L. purpureus (Salles et al., 2006;Gomes et al., 2008) and Ocyurus chrysurus (Vasconcellos et al., 2008).
Ocean currents have important biogeographic implications in the process of genetic differentiation of marine populations. In the Western Atlantic, two main ocean currents derived from the South Equatorial Current, the North Brazil one flowing in the northeastern/northern direction (RN -CE) and the Brazil Current flowing in a northeastern/southeastern direction (RN -BA -ES) (Lumpkin and Garzoli, 2005), dominate the areas of L. analis and L. jocu occurrence. In particular, the Brazil Current, with its northern and southern branches and extensive spatial trajectory, covers vast areas of both species' distributions, favoring gene flow among populations.
Besides the physical factors acting in the Western Atlantic, the absence of genetic structuring is associated with biological and ecological characteristics of the species involved (Lecchini and Galzin, 2003). In this respect, L. analis and L. jocu share biological features that favor long-distance dispersal with others lutjanid species. Among these are included the relatively long pelagic larval phase of approximately one month (Lindeman et al., 2001), high dispersive potential of adult individuals, ontogenetic habitat variation with a direct relationship between body size and depth (Frédou and Ferreira, 2005), formation of reproductive aggregations (Paris et al., 2005, França andOlavo, 2015), long life, and delayed sexual development (Claro, 1981;Allen, 1985;Graham et al., 2008).
Several marine organisms synchronize the release of eggs and larvae by forming spawning aggregations during specific periods. Like others, lutjanids, L. analis and L. jocu are r-strategists that form reproductive aggregations with large densities of individuals in spawning areas (Carter and Perrine, 1994;Lindeman et al., 2001;Paris et al., 2005). This reproductive strategy, although presenting considerable levels of self-recruitment in some areas, contributes larvae to distant populations, promoting increases in the genetic variability and homogeneity of the populations (Paris et al., 2005) presented by both species.
Ecological abilities to explore and colonize habitats can also increase genetic cohesion among L. analis e L. jocu populations along the Brazilian coast. In this respect, continuous environments (environmental corridors) favorable for the maintenance of different ontogenetic phases of the species are found throughout their distribution areas, such as rocky and sandy bottoms and reefs which extend for approximately 3000 km along the NE coast (Maida and Ferreira, 1997), as well as mangroves and estuaries (Cocheret de la Morinière et al., 2003), generating favorable habitats for both species. Additionally, the role of deep reef environments in the formation of effective corridors cannot be ruled out when establishing large population interchanges, even during periods of pronounced marine regression such as during the last glacial maximum (LGM).
The indication of population structure for L. analis between RN-ES based on cytB sequences may likely be associated with historical-geographic conditions (Pinheiro et al., 2015), and low contemporary gene flow. These regions have particular geological and geographic aspects that, under historical ocean level changes, could be responsible for particular demographic and genetic lineage divergences. The Brazilian continental shelf adjacent to the Rio Grande do Norte State has a reduced width and shallow depths as compared with other parts of the Brazilian shelf (Vital et al., 2010). In front of the RN coast extends a submarine volcanic chain, paralleling the northern equatorial coast (east-west direction), and known as the Fernando de Noronha Ridge. This geomorphologic formation consists of an alignment of at least, six volcanic seamounts, which extends in a direction toward emersed areas of the Rocas Atoll and Fernando de Noronha Archipelago. Historical asynchronic demographic patterns involving reef fishes related to these continental and insular areas have been attributed to changes in eustatic sea levels (Souza et al., 2015). During transgression periods, the insular areas were reduced, while demographic expansion in the continental populations had occurred. In contrast, during historical regressions, the insular areas were increased, favoring demographic expansions. In fact, despite the limited information on genetic structuring of reef fish groups in these insular and continental areas, population genetic data for some reef species have shown evidence of differentiated genetic structures among them (Cunha et al., 2014;Souza et al., 2015).
The higher haplotype diversity and the occurrence of a more frequent haplotype in RN and, to a minor degree, to northeastern populations support the possible role of Fernando de Noronha Ridge as an historical functional refuge zone for genetic lineages of L. analis during the sea level changes that occurred in the Pleistocene epoch. Similarly, the continental shelf in the ES region is narrow, also with a sequential chain of underwater seamounts denominating the Vitória-Trindade Seamount Chain (VTC), where the Trindade Island and Martin Vaz Archipelago are the unique emersed areas. During the last oceanic transgression, the VTC also acted as refuge areas, promoting genetic diversification of evolutionary lineages (Pinheiro et al., 2015). Likely, temporary isolation of L. analis lineages in both aforementioned north and south seamount regions created conditions contributing to the genetic divergences detected by cytB analysis.
Other lutjanids species, such the Atlantic species Ocyurus chrysurus, represent a large genetic unit (Vasconcellos et al., 2008), similarly to L. campechanus in the Caribbean, Gulf of Mexico and Florida coast (Garber et al., 2004), or L. kasmira in Indo-Pacific areas, where no signs of genetic structuring along the more than 12,000 km distribution exist, except for the Marquesas Archipelago (Gaither et al., 2010). In fact, panmictic genetic populations are mainly recurrent in lutjanids from other marine regions in the absence of apparent physical barriers.

DEMOGRAPHIC HISTORY
Both species presented values of N ef high enough to maintain a balance between the loss of variability by genetic drift and its replacement by new mutations. However, was evidenced a marked discrepancy between the species, in which L. jocu presented higher N ef values than L. analis.
Differences in contemporary or historically available habitats of L. analis and L. jocu, affected by recurrent cycles of retraction and sea level rise, may be linked to this demographic divergence pattern. Indeed, the widespread effects of such events have been perceived in several marine fishes (Grant and Bowen, 1998;Benzie, 1999;Lourie and Vincent, 2004), including lutjanid species such as L. kasmira, L. fulvus (Gaither et al., 2010) and L. erythropterus (Zhang et al., 2006). The demographic history of L. analis and L. jocu seems to be related to changes during the Pleistocene involving marine environments. Both species present evidence of population expansion along the Brazilian coast, as indicated by results of the Tajima (D) and Fu (Fs) tests. The measure of Fu's Fs statistic, based on the HVR1 control region sequences, revealed negative values significantly different from zero in all localities, thus supporting recent sudden population expansions (Ramos-Onsins and Rozas, 2002). Raggedness indices and SSD applied to searches for additional demographic signs of population expansion in L. analis were also in concordance with a historically expanding pattern.
The L. analis and L. jocu populations presented bimodal mismatch distributions that are usually associated with constant population size. However, this may also indicate the presence of two distinct lineages (e.g., Alvarado-Bremer et al., 2005). In this sense, the networks indicated more than one distinct lineage for bimodal mismatches (L. analis and L. jocu). The data suggested that past fragmentation was followed by population expansions resulting in secondary contacts of genetically diverged populations that had been separated by historical sea level changes during the Pleistocene. The population expansion period was significantly asynchronous between the two species (L. analis -∼ =17.1-19.4 Kybp; L. jocu -∼ =43.1-48.7 Kybp), but consistent with changes which occurred during the last glacial period between 120 and 12 Kybp (Martinson et al., 1987). Events in this period are considered important agents in the demographic patterns of several other lutjanids (Zhang et al., 2006;Gaither et al., 2010).
Due to exploration of different environment depths, the rise and retraction of sea level certainly promoted differential historical habitat conditions for L. analis and L. jocu, respectively. The population expansion of L. analis occurred during the beginning of the LGM period, in which ice sheets reached their maximum volume, between 19 and 26.5 Kybp (Martinson et al., 1987;Clark et al., 2009). This event possibly supplied a geographic expansion of L. analis by providing access to new available ocean areas. In contrast, as L. jocu occupies varied shallow habitats (mangroves, estuaries, and reefs) during its ontogeny, with eventual migration offshore with age , it must have undergone a differentiated demographic impact by changes in these ecologically complex habitats.
High N ef values were found in L. analis and especially in L. jocu populations. This population size discrepancy is an indication of considerably different niches available for these species. In fact, the vast and structurally complex shallow coastal habitats occupied by L. jocu have been shown to favor greater effective population sizes. Both species have their larger population effectives in northeast populations (CE, RN, and BA), with a drastic reduction in the ES region, near the southern limit of their geographic distribution.
The estimation of a large female contingent connected with reproduction, in association with high haplotype diversity, suggested potential resilience to present over-exploitation of both species. In fact, a high N ef value contributes to maintaining the balance between loss of adaptive genetic variation due to genetic drift and its restoration by mutation (Schultz and Lynch, 1997), and this is linked to the evolutionary feasibility of species. In this way, it is a useful indicator of population response to evolutionary and ecological forces (Waples, 2010), and environmental disturbances (Anderson, 2005). Besides providing useful genetic comparative indications for marine populations (Hare et al., 2011), the present N ef estimates are similar to those of the Indo-Pacific lutjanids (Gaither et al., 2010).

CONCLUSION
Genetic patterns indicate that the mutton and dog snappers exhibit panmictic populations along the Brazilian coast with a high level of genetic variability and large effective population sizes. Despite the absence of genetic structuring, as revealed by HVR1 region sequence analysis, for all populations of both species, there are indications of genetic structuring concerning L. analis populations of the RN and ES related to cytB sequences. This is likely associated to genetic divergence in refuge regions during Pleistocene glaciations. Historical contingencies in this period were also responsible for demographic expansions events in these species. In summary, the importance of spawning aggregations in maintaining high genetic variability along large areas of distribution, particular attention should be given to protecting reproductive areas against overfishing and other anthropogenic disturbance. The present data provide useful information for the sustainable exploitation of these overfished species, and create parameters for future monitoring of natural stocks.

DATA AVAILABILITY
The datasets generated for this study can be found in the GenBank access numbers JQ727916-JQ72996.

ETHICS STATEMENT
The animal study was reviewed and approved by Committee of Ethics in the Use of Animals of the Federal University of Rio Grande (# 044/2015).

AUTHOR CONTRIBUTIONS
ED and WM contributed to the conception and design of the study. AS, ED, and MP organized the database and performed the statistical analysis. AS, ED, and WM wrote the first draft of the manuscript. MP, EG-M, MV, MC, LB, and PG wrote sections of the manuscript. All authors contributed to the manuscript revision, and read and approved the submitted version.