Promiscuous Speciation with Gene Flow in Silverside Fish Genus Odontesthes (Atheriniformes, Atherinopsidae) from South Western Atlantic Ocean Basins

Promiscuous Speciation with Gene Flow in Silverside Fish Genus Odontesthes (Atheriniformes, Atherinopsidae) from South Western Atlantic Ocean Basins Graciela Garcı ´a 1 *,Ne ´stor Rı ´os 1 , Vero ´nica Gutie ´rrez 1 , Jorge Guerra Varela 2 , Carmen Bouza Ferna ´ndez 2 , Bele ´nGo ´mez Pardo 2 , Paulino Martı ´nez Portela 2 1Seccio ´n Gene ´tica Evolutiva, Facultad de Ciencias, UdelaR, Montevideo, Uruguay, 2Departamento de Gene ´tica, Facultad de Veterinaria, Campus de Lugo, Universidad de Santiago de Compostela, Lugo, Spain Abstract The present paper integrates phylogenetic and population genetics analyses based on mitochondrial and nuclear molecular markers in silversides, genus Odontesthes, from a non-sampled area in the SW Atlantic Ocean to address species discrimination and to define Managements Units for sustainable conservation. All phylogenetic analyses based on the COI mitochondrial gene were consistent to support the monophyly of the genus Odontesthes and to include O. argentinensis,O. perugiae-humensis and some O. bonariensis haplotypes in a basal polytomy conforming a major derivative clade. Microsatellites data revealed somewhat higher genetic variability values in the O. argentinensis-perugia populations than in O. bonariensis and O. perugia-humensis taxa. Contrasting population genetics structuring emerged from mitochondrial and microsatellites analyses in these taxa. Whereas mitochondrial data supported two major groups (O. argentinensis-perugiahumensis vs. O. bonariensis-perugiae-humensis populations), microsatellite data detected three major genetic entities represented by O. bonariensis,O. perugiae-humensis and an admixture of populations belonging to O. argentinensis-perugiae respectively. Therefore, the star COI polytomy in the tree topology involving these taxa could be interpreted by several hypothetic scenarios such as the existence of shared ancestral polymorphisms, incomplete lineage sorting in a radiating speciation process and/or reticulation events. Present findings support that promiscuous and recent contact between incipient species sharing asymmetric gene flow exchanges, blurs taxa boundaries yielding complicated taxonomy and Management Units delimitation in silverside genus Odontesthes from SW Atlantic Ocean basins. Citation: Garcı ´aG,Rı ´os N, Gutie ´rrez V, Varela JG, Bouza Ferna ´ndez C, et al. (2014) Promiscuous Speciation with Gene Flow in Silverside Fish Genus Odontesthes (Atheriniformes, Atherinopsidae) from South Western Atlantic Ocean Basins. PLoS ONE 9(8): e104659. doi:10.1371/journal.pone.0104659 Editor: Valerio Ketmaier, Institute of Biochemistry and Biology, Germany Received May 29, 2014; Accepted July 11, 2014; Published August 15, 2014 Copyright: ß2014 Garcia et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: The authors confirm that all data underlying the findings are fully available without restriction. All relevant data are within the paper and its Supporting Information files. Funding: This research received financial support from the project Fondo Marı ´aVin ˜as_2009_1_2793 (FMV_2009_1_2793_Project) granted by the Agencia Nacional de Investigacio ´n e Inovacio ´n (ANII) of Uruguay. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing Interests: The authors have declared that no competing interests exist. * Email: [email protected] Introduction The New World presents multiple examples of atherinid species flocks or adaptive radiations arising from habitat transitions [1,2,3]. Silverside fish from South America constitute a exciting model to understand the scenario of fish speciation driven by divergent natural selection [1,2]. The silverside genus Odontesthes includes 20 nominal species [4] distributed in marine, estuarine and freshwater environments of tropical and temperate regions in South America [5]. Most Odontesthes species co-occur in the same habitats and they are characterized by a great morphological homogeneity [6]. The low morphological divergence between species and the high meristic plasticity within species together with the tendency of local populations to form micro-geographic habitat associations had led to complicated taxonomy among silverside taxa [7]. Among freshwater representative species, two of them O. bonariensis and O. hatcheri are endemic of rivers and lakes located east of the Andes in subtropical and temperate areas [8]. The distribution of these species was originally allopatric: O. hatcheri occuring in the South (Patagonia), whereas O. bonariensis occupying Central and Northern Argentina, South Brazil and Paraguay. The occurrence of the spontaneous hybridization between both species in a communal laboratory tank has been reported [9]. On the other hand, marine silversides generally have similar life history strategies, occurring in large numbers in semi-isolated populations in estuaries and coastal lagoons [10,11,12,1]. Ten species of Odontesthes are endemic of a chain of small shallow lakes spread along the South Western Atlantic Ocean coastal plain [13,14]. Among them, in Patos Lagoon estuary and its adjacent marine coastal area occurs O. argentinensis and O. incisa, whereas PLOS ONE | www.plosone.org 1 August 2014 | Volume 9 | Issue 8 | e104659 in the freshwater habitats of Patos-Mirim lagoon system can be found O. bonariensis, O. humensis,O. retropinnis and O. aff. perugiae. Most Odontesthes species represent economically important resources for artisanal and recreational fisheries in South America and particularly O. bonariensis shows a great potential for aquaculture development [15]. The identification of incipient ecological species represents an opportunity to investigate the current evolutionary process where adaptive divergence and reproductive isolation are associated [2]. Beheregaray and Sunnucks [2] found that niche divergence due to estuarine colonization by marine silverside fish led to isolation by adaptation and speciation in the presence of high gene flow, one of the most convincing reports of parapatric speciation in aquatic organisms from the Southern Hemisphere. Beheregaray et al. [16] explored the role of adaptive diversification and recent sea-level changes as evolutionary drivers in the O. perugiae species complex which comprises several allopatric and sympatric morphotypes found in the lakes and rivers of southern Brazil, Uruguay and northern Argentina [17]. Most morphotypes have uncertain taxonomic status and are endemic to the vast system of lakes of the Coastal Plain of Rio Grande do Sul State (CPRS), southern Brazil [17]. Beheregaray et al. [16] performed a phylogeographic reconstruction of radiations in the South American coastal freshwater O. perugiae species complex, and also reported some of the most rapid speciation rates for a vertebrate group. Gene flow among hybridizing species with incomplete reproductive barriers blurs species boundaries, while selection under heterogeneous local ecological conditions or along strong gradients may counteract this tendency [18]. Thus, phylogeographic approach provides a valuable framework to identify signatures of divergent natural selection associated with ecological divergence and the possible occurrence of reticulation events among incomplete reproductively isolated taxa. In this study we implement a phylogeographic analysis based on mtDNA coding sequences (cytochrome oxidase subunit I, COI) and ten microsatellite loci to access in the species boundaries and to test possible reticulation and introgression events among Odontesthes taxa from the SW Atlantic Ocean, the Rı ´odela Plata estuary and in the Uruguay River basins. At the same time this information will contribute for a long-term success of Management Units for sustainable conservation of these taxa in fisheries and aquaculture. Materials and Methods Sample collection and DNA extraction All sampling protocols for this scientific study were approved by CNEA (Comisio´n Nacional de Experimentacio´n Animal) from Uruguay. A total of 163 individuals of Odontesthes from 20 sampling sites through three major regions, the Rı ´o de la Plata (RP) estuary (N = 45), Lower Uruguay and Negro river (UNR) basins (N = 23) Figure 1. Distribution map of 20 sampling sites through three major areas, lower Uruguay and Negro river (UNR) basins, the Rı ´ode la Plata (RP) estuary, and associated coastal lagoons and sites from SWAtlantic Ocean (AC) in Odontesthes as follows: UNRLas Can ˜as beach (CA B ), Yaguarete stream (YA S ), Pavo ´n stream (PV S ), Baygorria dam (BA D ), Rinco ´n del Bonete dam (RB D ), Ansina town (AN T ); RPBuceo Port (B P ), Hatchery and Carrasco lake (CA L ), Pando stream (PA S ), Pinar beach (PN B ), Solis Chico stream (SC S ); Solis Grande stream (SG S ); Piriapolis beach (PR B ), Sauce Lagoon (S L ), Chascomus Lagoon (CH L ), Argentina; AC-, Garzo ´n Lagoon (G L ), Rocha Lagoon (R L ), Castillos Lagoon (C L ), Valizas stream (VA S ), National Institute for Fisheries Research and Development (INIDEP), Argentina. doi:10.1371/journal.pone.0104659.g001 Promiscuous Speciation Silverside Fish Odontesthes SW Atlantic Ocean PLOS ONE | www.plosone.org 2 August 2014 | Volume 9 | Issue 8 | e104659 and Atlantic coast (AC) sites of SWA Ocean (N = 95) were included in the present study. All these samples were primarily ascribed to O. bonariensis (Ob), O. perugiae species complex (Op), O. humensis (Oh), O. argentinensis (Oa) and only two specimens from O. incisa (Oi) according to Dyer [19] morphological diagnosis. Tissue samples were obtained from artisanal gillnets fisheries operating in these areas during 2006–2012. The sampled areas are shown in Figure 1 and Appendix S1. Sample codes are as follows: collecting site name and the corresponding environments in lowercase (i.e.: S = stream, P = port, B = beach, L = lagoon or lake, D = dam, T = town). Tissues of the voucher specimens were deposited in the collection of the Evolutionary Genetics Section in the Faculty of Sciences, University of the Republic, Montevideo, Uruguay. Genomic DNA of sacrificed specimens was isolated from muscle tissue (fixed in ethanol 95%) using sodium chloride protein precipitation, followed by ethanol precipitation modified from Medrano et al. [20]. PCR amplifications and sequencing of the mitochondrial COI gene A fragment of 650 bp from the COI gene was amplified using FishF2 and FishR1 primers [21]. Reaction volume was 10 mL containing 10X supplied buffer, 0.25 mM MgCl 2 , 0.2 mM of each dNTP (10 mM), 0.25 mM of each primer (10 mM), 0.1 units of Taq DNA polymerase (Invitrogen) and approximately 100 ng/ul of template DNA. Cycling conditions consisted of one initial denaturation at 94uC for 5 min followed by 35 cycles of 94uC for 30 s, 52uC for 30 s, 72uC for 1 min and a final extension of 72uC for 10 min. Amplified COI products were sequenced directly on both strands in a Perkin-Elmer ABI Prism 377 Automated Sequencer (MACROGEN, Seoul, Korea). Sequence alignments were performed using Clustal X 1.8 [22]. Statistical analyses of sequences from COI data set Corrected estimates of pairwise sequence divergence were obtained using Kimura’s [23] two-parameter algorithm (K2P) implemented in MEGA 5.0 [24]. Within a population, DNA polymorphism was measured by calculating the proportion of segregating sites (S), the haplotype diversity (h) [25], and the nucleotide diversity (p) [25] with ARLEQUIN v3.11 [26] and DnaSP version 4.50 [27] programs. Tajima’s [28] and Fu’s [29] tests implemented in DnaSP 4.50 [27] were performed to check the mutation/drift equilibrium and any departure from neutrality. Significance of Fu’s Fs [29] and Tajima’s D[28] values was evaluated using the coalescent algorithm comparing the observed value with a null distribution generated by 10,000 replicates, and giving an empirical population sample size and the observed number of segregating sites. Phylogenetic analysis and divergence time estimates of the mitochondrial gene The phylogeographic relationships among mitochondrial COI haplotypes in Odontesthes populations from the sampled area were assessed by using two different methodologies. A non-model based method (MP, maximum-parsimony) was implemented in PAUP* 4.0b10 [30] following an equally weighted MP analysis using heuristic search (MULPARS option, stepwise addition, treebisection-reconnection [TBR] branch swapping, 100 replicates). A strict consensus between rival trees was computed to reconcile equally parsimonious topologies. The degree of confidence assigned to nodes in the trees was assessed by bootstrapping with 500 replicates. On the other hand, two model based approaches were also used, i.e., maximum-likelihood (ML) and Bayesian inference (BI), implemented in PAUP* 4.0b10 [30] and BEAST v.1.5.4 [31], respectively. In ML and BI analyses, the best-fitted nucleotide substitution model for each data set was determined in Modeltest v.3.7 [32] based on the Akaike information criterion [33], which simultaneously compares multiple nested or non-nested models. In the COI data set among the 56 models of nucleotide substitution, the best fit was the HKY+Cmodel [34] with gamma distribution (C). The gamma distribution shape parameter value was 0.18. The likelihood scores estimated for these models were used as the prior settings for the ML analysis in the data set (2lnL = 21602.50). Heuristic search (again with 100 replicates of stepwise addition and TBR branch swapping) in ML analyses was implemented in PAUP* 4.0b10 [30]. The robustness of the nodes was determined after 1,000 bootstrapping replicates as implemented in PhyML 3.0 (http://atgc.lirmm.fr/phyml), according to the algorithm developed by Guindon et al. [35]. In this case, the NNI (a fast nearest neighbour edge interchange search) swapping algorithm option was implemented. Nonparametric bootstrap values above 75% were considered to be robust support for clades [36]. All trees were rooted by means of an outgroup criterion using sequences of O. regia,O. incisa,O. smitti,O. hatchery and O. platensis and a more distantly taxon Atherina hepsetus retrieved from the GenBank. For the data set, divergence time of nodes and the age of the most recent common ancestor (tMRCA) were estimated with the BEAST v.1.5.4 software [31]. This program performs Bayesian statistical inferences of parameters by using MCMC (Monte Carlo Markov chain) as a framework. Input files were generated with Beauti v.1.5.4 [31] assuming uncorrelated lognormal trees and a Yule speciation process as prior information. The nucleotide substitution model and its parameter values were selected according to the Modeltest v.3.7 [32] results. An uncorrelated lognormal relaxed molecular clock, which allows rate variation among lineages, was implemented using an estimated rate for mitochondrial genome of 0.023 [2]. We carried out two independent runs of 10 million generations. Trees and parameters were sampled every 1,000 iterations, with a burn in of 10%. Results of each run were visualized in the Tracer v.1.5 program [37] to ensure that stationarity has been achieved and that convergence has been reached. Each analysis was repeated many times to optimize the operators of parameters until no suggestion message appeared in the log file. The timing of clade divergence and the tMRCA were estimated in million years ago (Mya) with a mean and a 95% HPD (lower and upper). Posterior probabilities and the maximum credibility tree were calculated using the TreeAnnotator v.1.5.4 software [31]. AMOVA, isolation by distance and historical demography To determine the genetic structure of Odontesthes samples the variance components among hierarchical partitions in the dataset were assessed by Analysis of Molecular Variance (AMOVA) [38]. The Euclidean metric of Excoffier et al. [38] was used to construct the pairwise distances matrix. The genetic variation was partitioned into three components, i.e., among groups (F CT ), among populations within groups (F SC ), and among individuals within populations (F ST ), after disregarding either their original populations or their groups. For both molecular markers, populations were ascribed to three major sampling areas such as Atlantic coast (AC), Rio de la Plata (RP) and Uruguay and Negro Promiscuous Speciation Silverside Fish Odontesthes SW Atlantic Ocean PLOS ONE | www.plosone.org 3 August 2014 | Volume 9 | Issue 8 | e104659 river basins (UNR), and different grouping hypotheses for the populations were tested. The significance of the observed Fstatistics was tested using the null distribution generated from 3,000 non-parametric random permutations of the data matrix variables and P-values were adjusted with sequential Bonferroni corrections for multiple comparisons [39]. Relationships and geographical distribution of the haplotypes were analysed in the haplotype network constructed with NETWORK v. 4.6.0.0 (http://www.fluxus-engineering.com/ sharenet.htm), which implements the median-joining method, in the absence of recombination [40]. The network was optimized using maximum parsimony criterion. Population subdivision and the level of genetic isolation among sampling sites were measured assuming an infinite sites model [41]. Pairwise estimates F-statistics were calculated in ARLEQUIN v3.11 [26]. To determine to what extent the geographic distance could explain the genetic differentiation among locations, a test for isolation by distance was performed using the Mantel test [42]. In this case, this test determines if there is a significant correlation between the geographic distance matrix (represented by the minimum coastline or river contour distance in kilometers) and the pairwise Fst matrix between collecting sites. The significance of the Z value (Mantel coefficient) was calculated using random permutation procedures implemented in the Mantel Nonparametric Test Calculator 2.0 [43]. Statistical significance was accessed through 1,000 permutations. To assess to the historical demography of Odontesthes we compared the observed frequency distribution of pairwise nucleotide differences among haplotypes (i.e., mismatch distribution) in relation to the expected under a sudden population expansion model [44] implemented in ARLEQUIN v3.11 [26] and DnaSP version 4.50 [27] programs. The significance of the assumed model was tested using the sum of squares deviations (SSD) between the observed and expected data by means a parametric bootstrapping approach (1,000 permutations) and considering the Harpending’s raggedness index [45]. The mismatch distribution will be multimodal in stable populations and unimodal in expanding ones. The time of a possible population expansion (t) can be calculated as t=2ut [44], where tis the mode of the mismatch distribution and uis the mutation rate of the sequence (such that u=mm T , where mis the mutation rate/site/generation and m T is the number of nucleotide base pairs). If the sudden expansion model was not rejected, then twas converted to time since expansion (t) in years before present as follows: [YBP (t=t/2u)]. For Odontesthes silversides, the mtDNA substitution rate was estimated in 0.023 [2]. Because time (t)is measured in generations and the age at sexual maturity for Odontesthes was calculated as minimum population doubling time 1.4–4.4 years (http://www.fishbase.org), to convert to time since expansion in years, we have multiply by the generation time of a mean 2.9 years. Analysis of microsatellite markers A total of 120 individuals from 13 populations were analyzed using these nuclear markers (Appendix S1). Ten polymorphic microsatellite loci developed for Odontesthes were amplified: Odon02, Odon09, Odon27, Odon38, Odon39 [2]; and Obo01, Obo26; Obo46; Obo54 and Obo77 [46]. The forward primer of each pair was fluorescently labeled as follows: Odon02, Odon25, Odon39, Obo01, Obo54 with 59-FAM; Odon27, Odon38, Obo26 and Obo77 with 59-HEX; and finally Odon09 with 59NED. PCR amplifications were carried out in a reaction volume of 10 ml (final concentrations in parenthesis) each containing DNA extract (400 ng/ul); dNTPs (0.1 mM each); primers (10 mM each); MgCl2 Invitrogen (0.8–2.5 mM); Taq DNA Polymerase Invitrogen (0.04 U/ml); and Invitrogen buffer (1X). Amplification conditions were those proposed by [2] and [46] respectively. The PCR reactions were carried out in a Verity 96-Well Thermal Cycler (Applied Biosystems) and the PCR products separated on an ABI 377 automated sequencer. The amplified fragments were genotyped using an ABI 3730 DNA Sequencer (Applied Biosystems) and visualization of the results was performed using the program GeneMapper 3.7 software (Applied Biosystems). Alleles were scored using a GeneScan 500 LIZ Size Standard and Genotyper software (Applied Biosystems, Inc.). Statistical and population structure analyses based on nuclear markers Among all populations, only 13 were analyzed with microsatellites. To implement the analysis of the Odontesthes data set, based on biogeographic criteria and to avoid statistical bias due to the low number of samples in some collecting sites, the populations were first ascribed to different taxa as follows: Oa including population from G L collecting site; Ob belonging to populations from CA L ,S L and C L ;Oph populations from CA B ,BA D ,RB D ; finally Oap embracing populations from B P ,PN B ,PA S ,SC S ,SG S , PR B and R L . The number of alleles, the allelic richness, the expected heterozygosity corrected for sampling bias, the observed heterozygosity, the polymorphic information content and the estimated null allele frequency were calculated for each locus in the whole population per taxon using CERVUS version 3.0.3 [47]. GENEPOP 4.0.10 [48] was used to perform the exact test for Hardy-Weinberg (HW) equilibrium by microsatellite loci (test multi-population) and by population (test multi-locus) using the Markov chain method with 1,000 iterations. Linkage disequilibrium between loci and deviations from Hardy-Weinberg equilibrium for each locus were tested by a Markov chain method following the algorithm of Guo and Thompson [49] and using the Bonferroni [50] correction for multiple comparisons (a= 0.05). All the analyses outlined above were implemented in GENEPOP 4.0.10 [48]. Wright’s F-statistics (Fis, Fst, and Fit [51]) over populations and loci were calculated by FSTAT version 2.9.3.2 [52]. To detect the presence of scoring errors or the possible presence of null alleles, we analyzed the genotypic matrices obtained with the Micro-Checker software [53]. Neighbor Joining tree based on D A distance [54] was constructed using Populations, 1.2.30 software package [55]. An analysis of population subdivision and clustering of individual genotypes was implemented with STRUCTURE v. 2.2 [56] by a MCMC method. We considered 1 to 13 different populations (K = 1 to K = 13). Ten independent runs employing an admixture model were implemented with a burn-in period length of 50,000 iterations, followed by 100,000 MCMC replicates. The average of these independent runs was calculated and the true value of ‘‘K’’ was accessed following the approach detailed in the manual of STRUCTUREv. 2.2 (http://pritch.bsd. uchicago.edu/structure.html). Different groups of hypotheses and populations as sources of variation were assessed in the AMOVA considering all ten loci using ARLEQUIN 3.1 software package [26]. Furthermore, F ST values for pairwise comparisons of the 13 Odontesthes populations and their significance level for genetic differentiation (P = 0.05) and Rst were tested additionally with FSTAT [52]. Promiscuous Speciation Silverside Fish Odontesthes SW Atlantic Ocean PLOS ONE | www.plosone.org 4 August 2014 | Volume 9 | Issue 8 | e104659 Population divergence and migration rates from both molecular markers To discriminate between the relative effects of divergence and gene flow on the speciation process, we analyzed our data set under the Isolation with Migration model [57]. The ‘‘isolation with migration’’ model in IMa does not assume gene flow and genetic drift are in equilibrium, making it the most appropriate for recently diverged populations that share haplotypes and alleles due to both gene flow and ancestral polymorphism. The model assumes that an ancestral population splits into two descendant populations that may continue to exchange genes after separation. Following [1] we consider Ob as a freshwater sister taxon of Oa and Oph, sharing a common freshwater ancestor with these taxa. The method estimates posterior probability distributions for both ancestral and actual population sizes, directional migration rates between the two populations, and the time elapsed since population splitting. An MCMC approach is used to draw a sample from the posterior distribution of genealogies and to estimate three types of population parameters: population size (h=4Nu), splitting time (t=Tu, where Tis the time in generations since the common ancestry, and it is of the same order of 4N) and migration rates (2NM =4Nu6m/2). The priors were finally set as follows: the upper bound of population sizes q= 10, splitting times t= 4 and migration rates m= 2, respectively. We run the MCMC simulations with 100,000 burn-in steps and 10,000.000 sampled genealogies. The posterior distributions of migration rates and population sizes are derived analytically from the sampled genealogies. Results Genetic variation in the mitochondrial COI gene in Odontesthes species from SWA Ocean basins This study includes a data set of 655 bp of mitochondrial COI gene from 156 individuals belonging to populations of O. argentinensis, O. perugia, O. humensis and O. bonariensis (GenBank accession numbers: KJ854753–KJ854894, see Appendix S1). Moreover, other sequences from Odontesthes species and one more distantly related genera (Atherina hepsetus) were retrieved from the GenBank and included for both the pairwise distance comparisons and the phylogenetic analyses. Among 36 COI haplotypes initially assigned to O. argentinensis, 30% of them were shared with O. perugiae and 25% with O. humensis respectivelly. Among 7 COI haplotypes initially grouping O. bonariensis sequences, 58% of them were shared with O. perugiae and 8% with O. humensis. Therefore we partitioned the statistical analysis in two different data sets: O. argentinensis-perugiae-humensis (Oaph)andO. bonariensis-perugiae-humensis (Obph). The Oaph populations showed higher haplotype diversity (h) and nucleotide diversity (p)thanObph (Table 1). Thirty six haplotypes were found in Oaph populations whereas only seven in Obph taxa. Except for the three most common haplotypes (H_1 and H_2 in Obph and H_6 in Oaph), most haplotypes represented rare variants that explained the observed haplotype diversity in each taxa (Appendix S1). A significant excess of low-frequency haplotypes and thereby negative and significant values of both Tajima’s and Fu’s neutrality tests were observed in Oaph indicating a departure from neutrality, whereas Obph presented only negative and significant values in Tajima’s D test (Table 1). These values would be consistent with populations that experienced demographic expansion scenarios or alternatively selective sweeps. In Oaph populations the average of the corrected pairwise K2P sequence divergence between COI haplotypes was higher than in Obph (Table 1). The average pairwise distances between haploTable 1. Estimates of DNA polymorphism in COI gene of Odontesthes populations from SW Atlantic Coast, Rı ´o de la Plata estuary and Uruguay-Negro River basins. Base pairs Variable Sites S Number of Haplotypes Haplotype Diversity p Kimura 2P Distance (Tv + Ts) D Fs Oaph 684 40 36 36 0.843 (0.034) 0.003 (0.013) 0.007 (0.001) 22.397 (P,0.01) 254.046 (P,0.00) Obph 684 20 10 7 0.696 (0.058) 0.002 (0.006) 0.004 (0.002) 22.797 (P,0.001) 20.715 (P.0.10) S = Average of polymorphic segregating sites; Haplotype diversity (h = gene) (Nei, 1987); p= Nucleotide diversity (Nei, 1987). Corrected Kimura 2P distances (1980). D = Neutrality test (Tajima, 989). Fs = Neutrality test (Fu, 1997). Standard deviation in brackets (SD). doi:10.1371/journal.pone.0104659.t001 Promiscuous Speciation Silverside Fish Odontesthes SW Atlantic Ocean PLOS ONE | www.plosone.org 5 August 2014 | Volume 9 | Issue 8 | e104659 types of O. argentinensis-perugiae-bonariensis-humensis taxa and other Odontesthes species included in present study (O. regia,O. platensis,O. smitti and O. hatchery) was 0.04760.011 (mean 6 SE), whereas the divergence between the former and O. incisa was 0.07560.018. The average divergence between the ingroup and the outgroup A. hepsetus was 0.61160.153. Phylogenetic analyses Present phylogenetic analyses included 43 haplotypes from Oaph and Obph populations. All performed phylogenetic analyses (ML and BEAST) conducted using the HKY+Gmodel of sequence evolution, clearly identified a major monophyletic and recently derivate clade with a high posterior probability of occurrence, including minor monophyletic clades which collapsed in a basal polytomy joining most of the 41 haplotypes of Oaph and a well supported clade of Obph (Fig. 2). Other minor clade integrated by two O. bonariensis haplotypes collapsed basal to the major clade. Other species from the genus Odontesthes (O. regia,O. smitti,O. hatchery and O. platensis) were reciprocally monophyletic in relation to the minor clade of O. bonariensis and the major derivate clade, whereas O. incisa was the most basal taxon of the genus Odontesthes. Simultaneously we used a reference calibration time for nodes by assuming a substitution rate of conventional rate for mitochondrial genome of 0.023 mutations/site/per million Figure 2. Tree topology generated using the HKY + Cmodel of molecular evolution based on 43 COI gene haplotypes (H) of Odontesthes from lower Uruguay and Negro river basins, the Rı ´o de la Plata estuary, SWA Ocean basins. Bayesian phylogeographic inference framework implemented in BEAST 1.5.4 and the estimated divergence dates. Numbers above branches refer to the Bayesian posterior probability of occurrence for clades while bootstrap support values from ML bootstrap are shown below branches. The bottom bar summarizes the time-scale divergence dates in Mya. doi:10.1371/journal.pone.0104659.g002 Promiscuous Speciation Silverside Fish Odontesthes SW Atlantic Ocean PLOS ONE | www.plosone.org 6 August 2014 | Volume 9 | Issue 8 | e104659 Table 2. Pairwise F ST values based on COI data set of Odontesthes populations from SW Atlantic Coast, Rı ´o de la Plata estuary and Uruguay-Negro River basins. YA S CA B PV S BA D RB D AN T R L SC S SG S PA S PN S G L B P C L VA S S L CA L PR B CH L INIDEP YA S _CA B _PV S (Op)0 BA D _RB D _AN T (Oph) 0.036 0 R L (Oap)0.229 0.085 0 SC S (Oap)0.261 0.092 0.023 0 SG S (Oap)0.315 0.103 20.031 0.013 0 PA S (Oap) 0.191 0.076 20.04 20.043 20.074 0 PN S (Oap) 0.147 0.040 0.024 20.042 0.049 20.013 0 G L (Oa)0.249 0.088 20.026 0.022 20.038 20.035 0.036 0 B P (Oap)0.262 0.086 20.020 20.056 20.065 20.057 0.009 20.026 0 C L _VA S (Obp) 0.206 0.207 0.405 0.582 0.616 0.551 0.483 0.503 0.539 0 S L (Ob) 0.133 0.086 0.259 0.305 0.340 0.301 0.215 0.288 0.302 0.196 0 CA L (Ob)0.480 0.313 0.312 0.400 0.450 0.387 0.283 0.355 0.377 0.709 0.220 0 PR B (Oap)0.333 0.195 0.211 0.253 0.313 0.250 0.026 0.243 0.255 0.667 0.332 0.499 0 CH L (Ob) 0.067 0.021 0.238 0.273 0.318 0.270 0.173 0.260 0.272 0.147 0.087 0.450 0.333 0 INIDEP (Oa)0.195 0.112 0.122 0.092 0.189 0.103 20.051 0.141 0.141 0.503 0.246 0.321 0.100 0.214 0 F ST significant values are in bold (P= 0.05). (See Fig. 1 and Appendix S1). doi:10.1371/journal.pone.0104659.t002 Promiscuous Speciation Silverside Fish Odontesthes SW Atlantic Ocean PLOS ONE | www.plosone.org 7 August 2014 | Volume 9 | Issue 8 | e104659 years [16] to capture a plausible time interval (lower and upper estimate) for clade divergence. All Odontesthes clades showed a high probability to have diverged between 0.1 and 2.5 Mya (Quaternary). The divergence between the genus Atherina and Odontesthes occurred in the Miocene. Population genetic structure, isolation by distance and historical demography Table 2 shows the pairwise F ST values of the COI data set among the 15 collecting sites analyzed with this marker. Low population genetic structure was detected among localities ascribed to Oap from RP estuary and AC areas respectively. Nevertheless, these localities appeared divergent to those ascribed to Obph from some RP estuary sites, and UNR basins. Remarkably the CA L collecting site seems to be the most divergent from all the remaining ones. The most plausible population structuring based on COI data set among tested hypotheses in the AMOVA was addressed following two different grouping criteria: (a) assigning all populations to two-group of samples; (b) forming three groups of populations (Table 3). The two-group hypothesis (a) pointed out that most genetic variation was distributed among groups (W CT ), suggesting a remarkably higher level of genetic structure when samples from marine and estuarine morphs Table 3. Analysis of molecular variance (AMOVA) based on COI gene of Odontesthes populations from SW Atlantic Coast, Rı ´odela Plata estuary and Uruguay-Negro River basins. Hypothesis Source of variation df Sum of squares Variance components Percentage of variation Wstatistics a Among groups 1 8.864 0.10574 Va 21.21 W CT =0.21209 Among population within groups 13 8.690 0.03416 Vb 6.85 W SC = 0.08695 Within populations 139 49.855 0.35867 Vc 71.94 W ST = 0.28060 b Among groups 2 9.588 0.09463 Va 19.44 W CT =0.19441 Among population within groups 11 7.419 0.03316Vb 6.81 W SC = 0.08456 Within populations 138 49.539 0.35897 Vc 73.75 W ST = 0.26253 Two grouping hypotheses among all tested: a) conforming two groups of populations as follows: 1-populations from S L ,C L ,CA L and CH L vs. 2-populations from VA S ,R L , SC S ,SG S ,PA S ,PN B ,G L ,B P ,PR B , INIDEP, CA B ,YA S ,BA D ,RB D and PV S ;b) separating three groups of samples as follows:1-populations from CA B ,PV S ,YA S ,BA D and RB D ;2populations from VA S ,R L ,SC S ,SG S ,PA S ,PN B ,G L ,B P ,PR B , INIDEP and 3populations from S L ,C L ,CA L and CH L . (See Fig. 1 and Appendix S1). doi:10.1371/journal.pone.0104659.t003 Figure 3. Haplotype network (constructed with NETWORK v. 4.6.0.0 software) of Oaph and Obph taxa. Black dots represent missing haplotypes and circle size is proportional to haplotype frequency. Different colours in each circle indicate the collecting sites as described in the Figure 1. doi:10.1371/journal.pone.0104659.g003 Promiscuous Speciation Silverside Fish Odontesthes SW Atlantic Ocean PLOS ONE | www.plosone.org 8 August 2014 | Volume 9 | Issue 8 | e104659 ascribed to Oaph populations was considered as a separate group from the other including freshwater samples from Obph taxa. The haplotype network based on COI gene (Fig. 3) showed a strikingly star-shaped topology including the two most frequent haplotypes with a high proportion of singletons, typical of populations that have suffered a recent demographic expansion. One of the most frequent and central haplotypes (H_6) including samples of Oaph is present in 11 sampling sites and is shortly interconnected by one to three step-mutations to most haplotypes of the network, belonging to Oaph populations. A single step mutation separated H_6 from the other most frequent haplotype (H_2), which included samples belonging to Obph taxa from 6 sampling sites. Remarkably, the network topology showed some loops involving the central H_6 and H_2 haplotypes and their respective derivate ones. Therefore these alternative links may be representing equally good connections due to homoplasy or perhaps the existence of peripheral reticulation events among them. Taking into account all collecting sites, negative values in Mantel test were observed (r = 20.126, p = 0.050), showing a negative correlation between genetic and geographic distances and excluding the isolation by distance model of population differentiation in Odontesthes. Figure 4 shows an unimodal mismatch distribution pattern in the COI data set which adjusted to the distribution predicted by the growth–decline population model [44] in the Oaph populations. The sum of squares deviations was SSD = 0.148 (P.0.06) and Harpending’s Raggedness index was 0.486. The estimate parameter under the model was t= 3.218. The time of expansiondecline in Oaph based on a substitution for this marker was estimated to have started around 227,000 YBP. In the Obph populations data set the sum of squares deviation value was SSD = 0.180 (P.0.05) and Harpending’s Raggedness index was 0.492, therefore population growth-decline hypothesis was accepted for this taxon. The estimate parameter under the model was t= 7.843 and the time of expansion-decline in Obph populations were estimated to have started around 497,000 YBP. Genetic variability in multi-locus nuclear data The measures of microsatellite genetic variation including total number of alleles, allelic richness, heterozygosity and HWE deviation for each locus and the corresponding average across all loci per taxon are showed in Table 4. Significant departures from Hardy-Weinberg equilibrium were found at some loci of the following populations: Oap (Odon27 and Obo54); Oph (Obo26, Odon27, Odon39 and Odon02); Ob (Obo39, Obo54 and Odon02). This could be due to a Wahlund effect, with a reduction Figure 4. Mismatch distribution in Odontesthes species under the growth–decline population model using mtDNA COI data set. (a). Oaph and (b). Obph populations. doi:10.1371/journal.pone.0104659.g004 Promiscuous Speciation Silverside Fish Odontesthes SW Atlantic Ocean PLOS ONE | www.plosone.org 9 August 2014 | Volume 9 | Issue 8 | e104659