Phylogeography of weakfish Cynoscion guatucupa (Perciformes: Sciaenidae) from the southwestern Atlantic

The genetic structure of populations/species was established during the Quaternary glaciations. Over the last 250 ka (Pleistocene), the South American marine biogeographic history recorded three main glaciations: the most extensive one between 140 and 180 ka, a minor one between 60 and 70 ka, and the last glaciation approximately between 15 and 35 ka. With the aim of assessing the pattern of molecular diversity and historical demography of weakfish (Cynoscion guatucupa), a 365 bp sequence of the mitochondrial control region was amplified at four coastal sites located in the southwestern Atlantic. Haplotype diversity was high, whereas nucleotide diversity was low and similar at each sample site. AMOVA failed to detect population structure. This lack of differentiation was subsequently observed in the distribution of samples sites in the haplotype network. Fu’s Fs was negative and highly significant while the mismatch analysis yielded a unimodal distribution, indicating a global population expansion. The Bayesian skyline plot revealed a coalescence time of weakfish population at approximately 210 ka, and a very rapid expansion from 180-190 ka, probably caused by a habitat expansion, as these two events coincide in time.


INTRODUCTION
The genetic structure of populations/species was established during the Quaternary glaciations.This genetic imprint of cyclic climate changes was mainly reported in terrestrial and freshwater species (Hewitt, 2000).Fewer studies have characterized the recent biogeographic histories of marine species (Grant andBowen, 1998, 2006;Beheregaray, 2008;Marko et al., 2010).In South America, the marine biogeographic history recorded 3 main glaciations over the last 250 ka (kilo annum = 10 3 years) (Pleistocene): the most extensive one (Oxygen Isotopic Stages -OIS 6) between 140 and 180 ka, a minor glaciation (OIS 4) between 60 and 70 ka and the last one (OIS 2) approximately between 15 and 35 ka (Rabassa et al., 2005;see Ruzzante et al., 2008).Ice reached 53ºS and glacialinterglacial cycles generated significant variations in the environmental conditions affecting the Pampean and Patagonian regions (Rabassa et al., 2005).The glaciations in this region led not only to global and sea temperature reduction, including a sea-level lowering of up to 100-140 m during the full glacial episodes (Fig. 1), but also to changes in the ocean current patterns and to the displacement or eradication of coastal habitats (Rabassa et al., 2005; see Grant and Bowen, 2006, for similar effects at other latitudes).Presently, the greatest diversity of marine species in the southern Atlantic spreads over four marine fronts, depending on the temperature and salinity features thereof, which could act as retention areas for larval marine fish (Acha et al., 2004).Coastal marine areas in South America are primarily affected by the Brazilian warm current and the Malvinas cold current, which converge between 35ºS and 45ºS (Piola and Matano, 2001;Acha et al., 2004).Though the subtropical convergence is at present located at approximately 38ºS (Fig. 1) (Piola and Matano, 2001), during the Pleistocene glaciations it could have been located near Espíritu Santo (20ºS) (Santos et al., 2006).In the present interglacial period, this cold water area moved southwards, resulting in a transition zone and establishing a new coastal area for colonization of some species.
These climate changes (temperature, marine currents and loss of coastal habitats) may have affected the abundance and geographic distribution of marine species of temperate coastal habitats.In addition, Pleistocene climatic fluctuations in South America may be responsible for the origin and distribution of some coastal estuarine fish which found refuge in lagoons.This is the case of Odonthestes perugiae on the Brazilian coast (Beheregaray et al., 2002) and of Brevoortia aurea on the Uruguayan and Argentine coast (García et al., 2008).Moreover, the divergence in two phylogenetic clades of Macrodon ancylodon could be attributed to historical changes in marine currents and water temperatures during the Pleistocene (Santos et al., 2006).The majority of the studies on southwestern Atlantic marine fish (O.argentinensis, M. ancylodon, B. aurea, Micropogonias furnieri) show a star-like phylogeographic pattern and population expansions (Beheregaray and Sunnucks, 2001;Beheregaray et al., 2002;Santos et al., 2006;García et al., 2008;Pereira et al., 2009) linked to the Pleistocene climate changes.Nevertheless, there remains a long way to go before a complete picture of this region biogeography can be drawn.
Marine fish populations may have expanded their distributional range as a consequence of climatic changes, so signatures of this expansion should coincide in time with some drastic climate event.The classical view in this respect is that the species underwent a severe bottleneck at the last glacial maximum (LGM), dated 23-25 ka, and then started growing.Conversely, another approach sustains that the species were not affected by LGM, thereby suggesting a long term population history (Marko et al., 2010).Along these lines, it has been suggested that benthic and pelagic species reacted differently to the Pleistocene ice-sheet expansion, which probably led to a considerable reduction in the suitable habitat for benthic species development (Rabassa et al., 2005;Janko et al., 2007).Therefore, to advance our understanding of coastal marine species distribution, it seems useful to study the population genetics, including aspects of comparative historical demography between species (Marko et al., 2010).
The present study focuses on the striped weakfish (Cynoscion guatucupa Cuvier, 1830) (Perciformes: Sciaenidae), which is a widespread pelagic-demersal fish predominantly found on the coast of South America, ranging from Rio de Janeiro, Brazil (22°S), to Chubut province, Argentina (43°S) (Cousseau and Perrotta, 2004).C. guatucupa belongs to a group of about 20 species corresponding to a multispecific demersal fishery (Carozza et al., 2001).Within this group, C. guatucupa is considered the second species in commercial importance after the whitemouth croaker (Micropogonias furnieri) (Ruarte and Aubone, 2004).The striped weakfish body size is close to 320 mm at maturity; its spawning period ranges from October to early April, with a main peak in October-November (Ruarte and Aubone 2004;Ruarte et al., 2004).It is a long-lived species (20-23 years).The purpose of this study is to provide a population genetic analysis using a partial sequence of the mitochondrial control region to assess the molecular diversity pattern, the historical demography and the estimation of the weakfish C. guatucupa expansion time on the Atlantic coast of South America.
DNA was extracted from small muscle pieces using the Chelex method (Estoup et al., 1996).Amplification of a partial fragment of the mitochondrial control region was conducted using primers A and E designed for teleost fish (Lee et al., 1995).PCR conditions using 50 µl reaction volumes were: 5 µl of 10× buffer, 3.6 µl of Cl 2 Mg (25 mM), 5 µl of each primer (200 mM), 6 µl of dNTPs (100 mM), 1 µl of Taq polymerase (5 U/µl) and 5 µl of genomic DNA, the reaction volume being completed with water.The PCR was performed at 94°C for 3 min, followed by 5 cycles of 94°C for 30 s, 45°C for 1 min and 72°C for 1 min, followed by 34 cycles of 94°C for 30 s, 58°C for 1 min and 72°C for 1 min, with a final extension at 72°C for 5 min (Beheregaray and Sunnucks, 2001).The PCR product was purified and sequenced in MACROGEN (Korea).
Sequences were manually aligned using PROSEQ (Filatov, 2002).For each sampling site, haplotype (h) and nucleotide (π) diversities were estimated using DNAsp v5 (Librado and Rozas, 2009).The genetic structure among Cynoscion guatucupa haplotypes was assessed through the analysis of molecular variance (AMOVA) using 10000 permutations in ARLE-QUIN 3.11 (Excoffier et al., 2005).The spatial differentiation of samples from Argentina was analyzed by estimating the pairwise F ST parameter (10000 permutations) using ARLEQUIN.The genealogical relationships between sequences were inferred by the haplotype network constructed with the method of statistical parsimony by Templeton et al. (1992) using TCS 1.21 (Clement et al., 2000).The substitu-tion model of the control region was established using jModelTest 0.1.1 (Posada, 2008).
Tajima's D (Tajima, 1989) and Fu's Fs (Fu, 1997) tests were performed to discriminate mutation/drift equilibrium and to evaluate the hypothesis of population expansion through the significant excess of low-frequency haplotypes.For neutral markers, significant negative values of D and F can be expected in cases of population expansion.Moreover, the demographic history of the control region was investigated using mismatch distributions, which are the distribution of pairwise differences among haplotypes (Rogers and Harpending, 1992).This method can discriminate whether a population has undergone a sudden population expansion or has remained stable over time.Mismatch distribution was compared with a distribution expected under a model of sudden population expansion and any deviation from the model was evaluated by calculating the raggedness index (r) (Harpending, 1994) and R 2 (Ramos-Onsins and Rozas, 2002).The populations that have undergone large expansion are expected to exhibit unimodal mismatch distributions with low r or R 2 values, while stable populations produce a variety of multimodal distributions with a high r or R 2 index.Though the R 2 index has greater statistical power than r for small sample sizes (n=10), Fu's Fs is more accurate than R 2 for larger sample sizes (n > 20) (Ramos-Onsins and Rozas, 2002).The significance of the indexes (D, Fs, r, R 2 ) was assessed with 10000 coalescent simulations using DNAsp.
In addition, we used two methods to estimate the time of Cynoscion guatucupa population expansion.First, the expansion time was estimated directly from the mismatch distribution with the statistic τ (tau) and translated into absolute time in years (t), using the equation t = τ/2µk, where µ is the per-year mutation rate and k is the number of nucleotides of sequence.Confidence intervals of t were obtained using a parametric bootstrap approach in ARLEQUIN.Second, we dated population size changes over time (Bayesian skyline plot: BSP) using the BEAST program (Drummond and Rambaut, 2007).BSP analysis was performed using a relaxed molecular clock, with three runs of 10 millions steps (MCMC) each, in which trees and parameters were sampled every 1000 steps.To estimate the expansion time by the mismatch and BSP analyses, a molecular clock of 5%/ma (mega annum =10 6 years) was used for the control region, also considering the 4-6% range (Bowen et al., 2006;Ruzzante et al., 2008).

RESULTS
The 365 bp segment of the control region amplified in 68 individuals allowed us to identify 62 haplotypes, 61 substitutions (47 transitions and 14 transversions) and 52 polymorphic sites.Haplotype sequences were deposited in GenBank (accession numbers JF423124-JF423185).For all the species (Table 1) the mean number of differences between pairs of sequences was 6.688 (sd = 3.194), haplotype diversity was 0.997 and nucleotide diversity 0.018.The two indexes were very similar at each sample site from Argentina (Table 1).
AMOVA results (Table 2) indicated that most genetic variation among the control region haplotypes was distributed within populations (f ST =-0.011,P=0.763).Accordingly, no significant differences were obtained from the analysis between sampling sites from Argentina (Ubatuba was not included due to the small sample size), i.e.Samborombón-Mar Chiquita (F ST =0.003, P=0.330), Samborombón-El Rincón (F ST =-0.017,P=0.904) and Mar Chiquita-El Rincón (F ST =0.006, P=0.273), respec-tively.The haplotype network showed a complex phylogeographic pattern (Fig. 2).The network displayed a geographically unstructured pattern with few common central haplotypes shared by at least two coastal areas from Argentina and many unique haplotypes; the three haplotypes from Brazil are mixed among the haplotypes from Argentina (Fig. 2).
Tajima's D was negative but only significant for the species.However, Fu's Fs, which are more sensitive to demographic changes, were negative and highly significant for the species and each sampling site (Table 1).The mismatch analysis yielded a unimodal distribution (Fig. 3) and a low and significant raggedness index (r=0.007,P=0.005), though a low and marginally significant R 2 index (R 2 =0.059,P=0.064) (Table 1).At the species level, not only Fs and D, but also r indicated population expansion.R 2 in contrast was not significant, but the power of Fs is greater than the R 2 index when the sample size is greater than 20 and the number of polymorphic sites greater than 50 (Fig. 3A y B in Ramos-Onsins and Rozas, 2002).According to that the above, both conditions are met by the species level sample (Table 1).
Since AMOVA failed to detect genetic differentiation, we pooled all the samples for the history demographic analysis.The τ parameter was 5.914 (confidence intervals, a=0.05: 5.049-6.629),estimating the expansion time with a molecular clock of 5% at 160 ka (140-180 ka).The BSP analysis, using a substitution pattern adjusted to a model of GTR + I + G, revealed that the haplotypes coalesced nearly 210 ka ago (5%/ma rate), when the species suffered a significant decrease in population size, and that the expansion began about 180-190 ka ago (Fig. 4).The BSP pattern, considering a mutation rate of 6-4%, indicated a coalescence time of 175 to 250 ka, respectively.In spite of these different coalescence times, the skyline trend is similar: a strong population expansion took place between 180 and 80 ka, remaining stable until the present time.

DISCUSSION
This study provides the first insight into the mitochondrial DNA variability of a coastal fish, Cynoscion   guatucupa, denoting an extraordinarily high polymorphism which shows 62 haplotypes (57 occurring only once) out of 68 individuals.C. guatucupa displayed high haplotype and low nucleotide diversity, which could be ascribed to rapid population growth and accumulation of mutations after a period of low effective population size.This conclusion seems to be supported by the unimodal mismatch distribution (low r and R 2 ) as well as by the significant excess of low-frequency haplotypes (negative D and Fs), both interpreted as sudden expansion indicators.The population expansion of C. guatucupa in the southwestern Atlantic, which was estimated directly from the mismatch distribution, started 160 ka (140-180) ago.In this respect, BSP analysis reported a rapid population expansion, starting between 180 and 190 ka.From a palaeo-climatic viewpoint, these periods coincide with the beginning of the largest cooling phase in the last 250 ka (Fig. 4).During this period, cold marine water moved northwards, changing the Malvinas-Brazil convergence (38º to 20º) (Fig. 1) and setting a new coastal area for C. guatucupa colonization (expansion) from southernmost distributions.Furthermore, pre-LGM coalescence times have recently been observed in marine fishes in the southern Atlantic at the Antarctic Convergence (Zane et al., 2006;Matschiner et al., 2009), in the northeastern Atlantic (Larmuseau et al., 2009) and in the northeastern Pacific (Wilson, 2006;Marko et al., 2010), suggesting that some marine species were not affected by a severe bottleneck at LGM.
However, the pattern is not necessarily reflected in all coastal fish species, as they could differ in their adaptive response to climate changes (mainly temperature, marine currents and loss of coastal habitats), and therefore could display a different pattern of genetic di-versity.Moreover, if demographic expansion is linked to a habitat distribution expansion, we could expect a different pattern for organism inhabiting pelagic and benthic coastal habitats.Eradication in coastal habitats affects most probably benthic rather than pelagic marine species (Janko et al., 2007).
A similar, though slightly posterior (120 ka ago), pattern of demographic expansion is found in the pelagic Antarctic fish Pleuragramma antarcticum (Zane et al., 2006), though not in the benthic Trematomus bernacchi and T. pennelli, in which population expansions correlate with the last glacial retreat (Janko et al., 2009).In addition, the demersal-benthic species M. furnieri (in two samples from the Argentine coast at approximately 36°S, 56°W and 39°S, 61°W) shows no signs of population expansions (Pereira et al., 2009).The pattern is probably explained by the recent population bottleneck or founder event undergone by a single or few mtDNA lineages (Pereira et al., 2009).M. furnieri suffered a bottleneck in the last glacial period, which could be interpreted in the light of coastal habitat loss with post LGM recolonization.Thus, the feeding habits of M. furnieri rely more on coastal benthic invertebrates (Giberto et al., 2007) than those of C. guatucupa, which is more fish species-dependent (Lopez Cazorla, 1996).
This study clearly provides evidence of an older demographic expansion and no evidence of population phylogeographic structure of C. guatucupa from the southwestern Atlantic coast.This conclusion seems to be supported by the AMOVA and by the haplotype network.Moreover, hypervariable nuclear markers (microsatellite loci), which can increase the discrimination power, also point to the lack of restrictions to present-day gene flow in this species in geographical Argentine locations (Sabadin et al., 2010).This outcome is consistent with the apparent lack of barriers in open marine environments.
While the mutation rates considered allowed for a coarse approach to the real expansion time, further research should be conducted in order to include the calibration date of molecular clocks (cytochrome b is the best candidate), more samples sites (especially from the Brazilian coast) and organisms with different life histories, so as to assess the pattern and ecological and evolutionary processes determining the genetic diversity of populations/species in this part of the world.

Fig. 2 .
Fig. 2. -Mitochondrial DNA network of control region sequences for Cynoscion guatucupa.The size of the ovals is proportional to the haplotype frequency in each population: Ubatuba (black), Samborombón (dark gray), Mar Chiquita (gray) and El Rincón (light gray), respectively.Each single line indicates one mutation between haplotypes, small circles (white) dividing single lines represent missing haplotypes.

Fig. 4 .
Fig. 4. -Climate change and the demographic histories of Cynoscion guatucupa over the last 250 ka.The most significant glaciations of southern South America are shaded (see text).The time series on the x axis are derived from mtDNA control region sequences, and the y axis represents the product of effective population size and generation length.The thick solid lines are median estimates under the assumption of per site mutation rate of 5%/ma, and the dotted lines indicate 95% highest posterior density (HPD) regions.The thin top solid line shows the median obtained under the assumption of per site mutation rate of 4%/ma, while the thin bottom line depicts the median obtained under the assumption of per site mutation rate of 6%/ma.

Table 1 .
-Genetic diversity (N = sample size; Nh = number of haplotypes; Np = number of polymorphic sites; k = average number of nucleotide differences between pairs of sequences; h = haplotype diversity; π = nucleotide diversity), neutrality and demographic test (Tajima's D; Fu's Fs; r and R 2 ) for the 365 bp of the control region from Cynoscion guatucupa in El Rincón (RIN), Mar Chiquita (MCH), Samborombón (SAM) and Ubatuba (UBA).