Wide-range genetic connectivity of Coney , Cephalopholis fulva ( Epinephelidae ) , through oceanic islands and continental Brazilian coast

The Epinephelidae form a group of species of high biological and economical interests. It ́s phylogeographic patterns are not well known especially the distributed populations in the western region of the Atlantic Ocean. Among the representatives is a small species called Cephalopholis fulva, Coney, which presents a wide geographical distribution, polychromia, hermaphroditism and is quickly becoming a large target for the exploration of commercial fishing. The genetic and historical demography were obtained through the partial sequence analysis of Control Region from six locations on the coastline of Brazil from the northeast coast to the southwest coast, including the oceanic islands of Rocas Atoll and Fernando de Noronha Archipelago. The analyzed samples revealed a high genetic variability and a strong gene flow among the sampled locations. Additionally, the genetic data revealed that population expansions probably occurred due to the changes in the sea levels that occurred during the Pleistocene. The large population connectivity found in Coney constitutes relevant conditions for their biological conservation.


INTRODUCTION
Some species of the Epinephelidae family show peculiar biological characteristics which make them particularly vulnerable to climate changes, fishing and environmental degradation, such as high longevity, late sexual maturation, hermaphroditism, reproductive aggregation and need of nurseries in estuarine regions (Coleman 1999).In fact, around 40% of this family´s species are considered under some level of anthropogenic threat (Morris 2000).
The biological characteristics of the Coney such as their small size, abundance and broad geographic distribution through Atlantic´s Western coast makes this species particularly indicated for phylogeographic analysis.In fact, its geographic distribution extends from the Bermudas and South Carolina (USA), to the Southeast of Brazil, including the oceanic islands in Rocas Atoll, and Fernando de Noronha Archipelago (Heemstra and Randall 1993, Freitas et al. 2003, Sazima et al. 2005) and Trindade and Martim Vaz (Gasparini and Floeter 2001).Polychromatism is a characteristic of this species, not entirely known in all its extension, although it has indications of possible adaptive value, related to the aggressive mimicry of some colors (Sazima et al. 2005).
Some reproductive aspects of marine fishes interfere with the genetic diversity pattern; among them are the sex-determining mechanisms.Like many other Epinephelidae, Coney is a protogynous hermaphrodite, in which the female reaches sexual maturity at 16 cm and begins their sexual reversion into a male when they reach approximately 20 cm, starting to exhibit territoriality and harem formation (Heemstra and Randall 1993, Coelho 2001, Araújo and Martins 2006).According to the evolutionary point of view, hermaphroditism can be an advantage, in relation to gonochorism, among other aspects, when presented in small populations with great fluctuations between sexes (Borgia and Blick 1981).
Despite Coney not being considered an endangered species, its fishing is important in the Caribbean (Trott 2006), and seems to have continually increased due to a displacement of fishing pressure on smaller size species, having in mind the decrease of populations of large fish and the most commercially valued fish (IUCN 2010).In fact, the Coney capture volume has increased since the beginning of the 1990's (Martins et al. 2005).Between the years of 1996 and 2000, the Coney represented 2.4% of the total fishing production in the Northeast of Brazil (Frédou et al. 2006).With an average of 1,116 tons/year, the Coney represented the most abundant species of collected samples in the Central coast of Brazil (Klippel et al. 2005, Martins et al. 2007).As a result of this exploitation, in some areas the effects of fishing were detected through the abundance and space distribution of the individuals in relation to the depth (Coelho 2001).
Signs of overexploitation of C. fulva, as well as its extensive geographical distribution in the Atlantic Ocean and the presence of many color morphs, make studies that identify its phylogeographic aspects and genetic variability levels particularly favorable.This data will provide a comparative base to establish future interpopulational surveys on this species.Here we present genetic data obtained from partial sequences from the D-loop mitochondrial region from samples collected in a broad ocean strip of the Brazilian coast, covering continental and insular coast environments.The elevated genetic diversity and absence of structuring along all of the sampled areas suggests the existence of a single panmictic population along the Brazilian coast.

SPECIMENS COLLECTION
At the conclusion of this research there was no Animal Ethics Committee at the research institute where this study was conducted.Thus, this study was carried out in accordance with Brazilian law regarding the use of laboratory animals Archipelago and Rocas Atoll, and also the coastal regions of the Northeast and Southeast of Brazilian coasts, most precisely in the coast of the states of Ceará (CE), Rio Grande do Norte (RN), Bahia (BA) and Espírito Santo (ES) (Fig. 1) (Table I).The most extreme collection points covered approximately 2,100 km of distance.The samples were taxonomically identified according to Heemstra and Randall (1993).Fragments of hepatic tissue and/or fin tissue were conditioned in micro tubes (1.5 ml) with ethylic alcohol and methyl alcohol (1:1) and stored in temperature of -20°C.ALLYSON S. DE SOUZA et al.

EXTRACTION, AMPLIFICATION AND DNA SEQUENCING
The total DNA extraction was performed according to Sambrook et al. (1989).The hypervariable region 1 (HVR-1) of control region (D-loop) of mtDNA was amplified using primers CR-A and CR-E described by Lee et al. (1995).Each reaction contained 4 μl of dNTP (1.25 mM), 2.5 μl of tampon (10x), 1 μl of MgCl 2 (50 mM), 2.0 μl of each primer (10 pmol/μl), 1-1.5 μl of total DNA, 0.3 μl of Taq DNA Polimerase (5U/μl) (Invitrogen, USA) and water Milli-Q for a final volume of reaction of 25 μl.The used cycling conditions consisted of initial denaturation at 94°C for 2 minutes; 30 denaturation cycles at 94°C for 45 seconds, hybridization at 52°C for 45 seconds, extension at 72°C for 1 minute; and final extension at 72°C for 2 minutes.The PCR products (5 μl) were purified with enzymes exonuclease I (3.3 U/reaction) and shrimp´s alkaline phosphatase (0.66 U/reaction) (GE Healthcare), in thermal cycler, submitting the mix to a cycle of 30 min in 37°C and 15 min in 80°C.The samples were sequenced in EMBRAPA (Brazilian Company of Agricultural Development) Sequencing Platform through the automatic DNA sequencer ABI Prism 3700 (Applied Biosystems).

mtDNA ANALYSIS
The obtained electropherograms were checked using Bioedit (Hall 1999) and an automatic multiple alignment was performed in the Clustal-X application (Thompson et al. 1997).
The saturation between transitions and transversions in the sequences was checked by the software DAMBE (Xia and Xie 2001).The software jModelTest 2.1.4(Guindon andGascuel 2003, Darriba et al. 2012) was used to select the best-fit model of nucleotide substitution.The genetic relationships among samples were evaluated through neighbor-joining method using MEGA5 (Tamura et al. 2011).The significance of the groupings in all generated trees was estimated using bootstrap analysis, based on 1,000 pseudo-replicates.
A haplotype network was set up using the median vector method available in the Network 4.5 program (Bandelt et al. 1999).
The levels of genetic diversity, population structuring, through the estimative of F ST (Wright 1978), Mantel test (Mantel 1967) and AMOVA (Excoffier et al. 1992) were obtained by Arlequin 3.5 (Excoffier and Lischer 2010).The historical demography was estimated by statistic Fs (Fu 1997) and D (Tajima 1989).Complementary to these, the population expansion condition was investigated through analysis of frequencies of the pair-to-pair differences (mismatch distribution), based on three parameters, θ 0 , θ 1 (frequencies of the pair-to-pair differences between the sequences, before and after the population expansion) and τ (time of expansion expressed in mutational time unities) (Rogers and Harpending 1992, Rogers 1995).These estimates were obtained using Arlequin 3.5 (Excoffier and   Lischer 2010).Values of τ were transformed in real time, since the expansion estimate time, from u = μmt, where u is the mutational rate by segment of mtDNA, μ is the mutational rate estimated for the analyzed segment; mt is the amount of the analyzed nucleotide bases.Subsequently, τ = 2ut, where t is the estimated time of expansion occurrence since then.The generation time adopted for the species Coney was one year (Heemstra and Handall 1993, Araújo and Martins 2006, 2009).The mutation rates (μ) adopted for the HVR-1 control region was 8.24 x 10 -8 and 9.30 x 10 -8 (Domingues et al. 2005).In addition, the raggedness index (Harpending 1994) and the sum of the squared deviations (SSD) between observed and expected mismatch distribution were calculated to validate the estimated expansion model (Schneider and Excoffier 1999).

SEQUENCE ANALYSIS
A total of 388 bp of HVR-1 sequence was resolved on 204 specimens (Genbank accession numbers of haplotypes: KC831794 -KC831953) from six geographical localities revealing 94 parsimony informative sites and 153 nucleotide substitutions (23 transversions and 130 transitions).The plot of transition/transversion vs. genetic distance did not indicate the presence of saturation (Fig. 2).The JC Model (Jukes and Cantor 1969) was selected as the best-fit model of nucleotide substitution for the sequences.The content of A/T was clearly higher (69.3%) than G/C (G=11.46%,C=19.28%,T=30.14% e A=39.18%).The haplotype diversity (h) and nucleotide diversity (π) were high and very similar among geographical localities, varying respectively from 0.977 to 0.998 and 0.022 to 0.024 (Table II).These high values of diversity were reflected both in the neighbor-joining analysis, which revealed a large polyatomic tree (there were no significant clusters of samples corresponding to sampling localities; Fig. 3) and in the haplotype analysis, which revealed a star-shaped haplotype network with numerous connections between the median vectors (red circles, Fig. 4), suggesting the existence of a single panmictic population of Coney in Brazilian waters.Furthermore, no genetic differentiation was found among the different color morphs found in Coney.

TABLE III F ST estimatives between geographical populations based on HVR-1 of coney.
Upper diagonal = p values ± confidence interval; below diagonal = Fst values; *p<0.05.

POPULATION STRUCTURING
The pairwise population comparison by F ST (Wright 1978) showed, in most cases, negative or not significant values.Just three situations with significant F ST were found, which were between Bahia and Espírito Santo (F ST =0.020), Fernando de Noronha Archipelago and Espírito Santo (F ST =0.025) and between Rocas Atoll and Bahia (F ST =0.035) (Table III).The Mantel test found no association between the F ST values and geographic distances (p > 0.05).
AMOVA analyses using different population groupings indicate that the higher variation   percentage happens within populations (97.27% -99.54%) and no significant differentiation between the groupings coast-to-coast, coast-to-island or island-to-island (Table IV).

NEUTRALITY AND DEMOGRAPHIC HISTORY
The neutrality tests showed negative values for D and Fs for all sampled geographical localities and groupings (Table V), which indicates the non-neutral nature of the Coney´s HVR-1, and supports the hypothesis of demographic expansion.The mismatch distribution revealed in a general unimodal pattern, especially for the groupings "Insular populations", "Coastal populations" and "All populations", supporting a model of sudden expansion (Fig. 5).In addition to these tests, raggedness index and SSD (Harpending 1994) also supports the idea that Brazilian Coney population is undergoing expansion (Table V).The estimate τ values indicate that demographic expansions in this species happened between 98 and 111 thousand years ago in insular populations (Rocas Atoll and Fernando de Noronha Archipelago) and between 131 and 148 thousand years ago in coastal populations (CE, RN, BA, and ES).Considering a panmixia framework and analyzing all sampled populations as one, the expansion would have happened between 130 -135 thousand years ago (Table V).

POPULATION
The control region mtDNA of the geographic samples of Coney show high values of haplotype diversity (0.977-0.998) and nucleotide diversity (2.2% -2.4%), when compared with another Atlantic grouper, Epinephelus itajara (Silva-Oliveira et al. 2008).Silva-Oliveira et al. (2008) have analyzed populations of E. itajara in the Northern Brazilian coast and found values of haplotype and nucleotide diversity considerably lower, varying from 0.86 -0.53 and 0.5% -0.1%, respectively.E. itajara is a critically endangered species, and according to the authors, its low diversity values were related to high fishing pressure and loss of habitat.
The low F ST values suggest extensive gene flow in almost all sampled regions, besides collected individuals located more to the South (ES) showed little genetic differentiation with those in Fernando de Noronha Archipelago (FNA) and Bahia (BA), more to the North.Likewise, the sample of BA also showed a small differentiation with samples from Rocas Atoll (RA).Curiously, the population of RA shows more similarity with the population of ES, farther south.This condition could suggest recent restrictions between these places; however, a biased sample can not be discarded.Incongruities between genetic patterns and geographical distribution along the Brazilian coast were also observed between populations of the Lutjanidae Ocyurus chrysurus (Vasconcellos et al. 2008).
The estimate of F ST , haplotype network and neighbor-joining analysis can not indicate genetic differentiation along the coast and ocean islands suggesting a genetic connectivity among the sites.In fact, the AMOVA analysis indicated that the highest genetic variation happens within populations.A low differentiation, next to the critical significance limit was shown between populations within groups ((RA, FNA) x (CE, RN, BA, ES)) (ɸ SC = 0.007, p = 0.041).However, when analyzing insular locations (RA, FNA) versus coastal locations (CE, RN, BA, ES), the differences found do not reveal real groupings (p = 0.057 and p = 0.126, respectively).Furthermore, the presence of different color patterns found in Coney seems to be related to environmental or biological factors, because there is no genetic differentiation among the sampled color morphs.However, further analysis using nuclear genes may be used to confirm this condition.
The genetic homogeneity observed in Coney appears very different from other representatives of Epinephilinae, in which a pronounced genetic structuring has been identified in some grouper species such as Epinephelus akaara (F ST ≤ 0.379; p<0.001) (Chen et al. 2008) and Epinephelus labriformis (F ST ≤ 0.661; p = 0.001) (Craig et al. 2006).Although other evidence is needed, the low level of genetic differentiation among the sampled locations could be related to the absence of dependence from estuarine regions in any life stage of the Coney, the adding of continuity of rocky substrates along the Brazilian coast and a directional regime of North and South currents.

COAST
The glaciations of Pleistocene, have often been appointed as being relevant events for the establishment of phylogeographic and population PANMITIC POPULATION OF Cepholopholis fulva ALONG THE BRAZILIAN COAST patterns currently shown by marine fish (Domingues et al. 2005, 2006, 2007, Santos et al. 2006, Craig et al. 2006, Zhang et al. 2006, He et al. 2010, Mobley et al. 2010, Viñas et al. 2010).In this sense, indicative indexes of the samples´ demographic past for each region were used, not only individually but also through geomorphological criteria for samples jointly placed in coastal and insular regions.Both neutrality tests Fs and D, showed negative and highly significant results, like mismatch distribution which appears unimodal, used in insular and coast clusters which appeared compatible with the populations that suffered expansion (Slatkin and Hudson 1991, Rogers and Harpending 1992, Ray et al. 2003, Excoffier 2004).
Dating by oxygen isotopes in two marine terraces of RN´s coast indicated that the Atlantic´s level in this region reached 8 ± 2m above the current sea level around 123.5 ± 5.7 thousand years ago (Fig. 6), period of the Pleistocene´s penultimate sea transgression (Barreto et al. 2002).Similar ocean elevation levels were determined for marine terraces located in BA´s coast (Bernat et al. 1983).This period corresponds to approximately the same period of the last population expansion of Coney in areas of the continental platform.
Distinct events seem to have influenced the asynchronic population expansion in insular and coastal clusters.The demographic expansion in coast cluster (CE, RN, BA, ES) occurred with a smaller number of individuals (θ 0 =0.021) than in insular regions (θ 0 =2.376).In fact, in coastal regions, the data shows a demographic expansion of Coney that is coincidental with the penultimate sea transgression which reached its top around 123.5 ± 5.7 thousand years ago, when the continental shelf was submerged again, making new habitats available to be taken, consequently causing the population expansion, emphasized in coast populations (see Fig. 7B).
This transgression event (during the Glacial Minimun) of Pleistocene is preceded by the maximum sea regression, occurred around 140 thousand years ago, which in some areas reached up to 100 meters below the current level (Hearty 1998) making practically the whole Brazilian continental shelf emerge.In this occasion the loss of broad natural habitat areas in the Brazilian coast possibly resulted in a deep decrease of population effectives in the Brazilian coast (Fig. 7C), like suggested by the low values of θ 0 (Table V).In this period, as the continental shelf was emerged, Coney came to occupy seamounts that were not accessible due to depth, and a narrow area of the continental slope, which the depth has become the limiting factor for the occupation of Coney in the continental slope (Fig. 7C).
More recently, in last glaciation period, in the course of a growing decrease of the Atlantic Ocean's level (Fig. 6), a population expansion involving only contingents of Coney located in insular environments was noticed.Values of θ 0 found for insular populations were much higher than those found in coast populations, indicating population expansion from a higher number of Coney individuals than those in the continent.
Coney inhabits areas in depths from 2 to 35 meters (Gasparini and Floeter 2001), although the species have already been registered in 150 meters of depth (Smith 1997), as represented in Fig. 7. So, its populations seem intensely prone to suffering interference by loss of areas in the continental platform during sea regressions (Fig. 7C).Besides, in general, ocean islands are considered more sensitive environments to environmental effects, they seem to have kept more population effectives during sea Among the conditions that can be suggested for population expansions in Coney, some environmental characteristics of geologic/geographic nature seem to have contributed to a higher θ 0 and asynchrony of the expansion period for insular populations.
The geological conformation and geographic position in which the oceanic islands of Fernando de Noronha Archipelago and Rocas Atoll are located, could have positively contributed to the maintenance of genetic homogeneity.These regions form, equally, conical mounts supported by the oceanic floor around 4,000 meters of depth (as represented in the Fig. 7A), depth much higher than that shown by the continental platform, but also easily impacted by sea regressions.For having conical shape, retractions of sea level could increase the available of rocky areas and consequently the areas of habitats for Coney individuals (Fig. 7C).This condition can explain the fact that the population expansion in Fernando de Noronha Archipelago and Rocas Atoll happened in a period of regression of the Atlantic Ocean´s levels.
According to the geographic view, Fernando de Noronha Archipelago and Rocas Atoll are part of an alignment of, at least, six more submarine mounts of volcanic origin, the Fernando de Noronha Range which extends in continent direction.From this range, currently only these two insular regions are immersed (Almeida 2002, Kikuchi 2002), but under a scenario of decrease in the sea level, the insular environments could be wider and shelter more effectives of Coney from other historically immersed areas as the sea mountains (Fig. 7C).
Although evidences are preliminary, the low indexes of structuring of Coney between insular and continent regions or between continent regions can be explained by biological factors of the species and PANMITIC POPULATION OF Cepholopholis fulva ALONG THE BRAZILIAN COAST geological conditions of the environment in which the species inhabits.This species shows high larval production, with release of 150.000 to 282.000 eggs per female (Heemstra and Randall 1993).Despite the absence of data about the duration of the larval pelagic, in some representatives of the Epinephelidae family, this stage can reach 80 days (Lara et al. 2009), a period which is considered to be long comparing it to other Perciformes.In fact, the broad geographic distribution of Coney (Nelson 2006, Heemstra andRandall 1993), suggests an elevated level of dispersion and colonization of the species.
Associated to these biological characteristics, a system of ocean currents in the Brazilian Atlantic coast (Fig. 1) favors the larval dispersion and the subsequent population homogeneity.The main ocean current present in the Brazilian coast is the South Equatorial Current (SEC), which flows in East/West direction and forks between latitudes 12 and 14, forming the North Brazil Current (NBC) and the Brazil Current (BC).Brazil Current flows parallel to the continental platform, in Northeast/ Southeast direction while North Brazil Current flows in Northeast/North direction (Castro andMiranda 1998, Lumpkin andGarzoli 2005).In the region of Fernando de Noronha´s range of mountains, the Atlantic Equatorial Current (AEC), which originates between the North and Northeast coasts of Brazil, flows under and in opposite direction of the South Equatorial Current (Mendes 2006).The action of single currents on a vast geographic area seems to perform an important role in the larval dispersion and genetic homogeneity of Coney.
Despite Coney being an increasing target of fishing for several years, the genetic profile obtained for the species does not show damaging effects coming from overfishing up until the present.Other than this, evidence of genetic connectivity has not been found among continent/ continent, continent/island and island/island areas.
The broad geographic distribution, the small size of sexual maturation and the decrease in competition by capturing bigger size predators can be positively influencing its abundance and maintenance of genetic variability.So, even though Coney in Brazilian coast does not suggest specific handling actions, the genetic data which was obtained, allow monitoring and accessing the conservation status of the species in the future.

Figure 2 -
Figure 2 -Plot of transition/transversion rates vs. genetic distance of the HVR-1 sequences of C. fulva.

Figure 3 -
Figure 3 -Neighbor-joining tree constructed using Jukes and Cantor distances and inferred from 204 samples of Coney.Each color represents where each sample was collected.CE = Ceará, RN = Rio Grande do Norte, RA = Rocas Atoll, FNA = Fernando de Noronha Archipelago, BA = Bahia, ES = Espírito Santo.Bootstrap values greater than 70 are shown in the branches.

Figure 6 -
Figure 6 -Variations of Atlantic Ocean´s level and periods of population expansion of Coney in Brazilian coast (between 148,000 and 131,000 years; in yellow) and islands (between 111,000 and 98,000 years; in red).Modified from Barreto et al. 2002.

Figure 7 -
Figure 7 -Illustration representing the occupation areas of Coney along the variations of the sea levels during the Pleistocene; A = Submarine relief of the Brazilian coast and oceanic islands and the current sea level; Red areas represent the habitat occupied by Coney during periods of transgression (B) and marine regression (C).

TABLE I Collection points and analyzed samples of coney through the coast of Brazil. n
= Amostral number.PANMITIC POPULATION OF Cepholopholis fulva ALONG THE BRAZILIAN COAST

TABLE IV Hierarchical analysis of molecular variance (AMOVA) based on HVR-1 of coney.
Significant at the 0.05 level. *