Population Genetic Structure of Cnesterodon decemmaculatus (Poeciliidae): A Freshwater Look at the Pampa Biome in Southern South America

The Pampas is a Neotropical biome formed primarily by low altitude grasslands and encompasses the southernmost portion of Brazil, Uruguay, and part of Argentina. Despite the high level of endemism, and its significant environmental heterogeneity, Pampean species are underrepresented in phylogeographic studies, especially aquatic organisms. The Pampean hydrological system resulted from a long history of tectonism, climate, and sea level changes since the Neogene. In this study, we examined the population genetic structure of Cnesterodon decemmaculatus, a freshwater fish species that occurs throughout most of the Pampa biome. We characterized mitochondrial and autosomal genetic lineages in populations sampled from Southern Brazil and Uruguay to investigate (1) the correspondence between current drainage systems and evolutionary lineages, (2) the demographic history for each genetic lineage, and (3) the temporal depth of these lineages. Overall, we found that the major evolutionary lineages in this species are strongly related to the main Pampean drainage systems, even though stream capture events may have affected the distribution of genetic lineages among drainages. There was evidence for recent population growth in the lineages occupying drainages closest to the shore, which may indicate the effect of quaternary sea-level changes. In general, divergence time estimates among evolutionary lineages were shallow, ranging from 20,000 to 800,000 years before present, indicating a geologically recent history for this group, as previously reported in other Pampean species. A Bayesian phylogeographical reconstruction suggested that an ancestral lineage probably colonized the Uruguay River Basin, and then expanded throughout the Pampas. This evolutionary scenario may represent useful starting models for other freshwater species having a similar distribution.

Introduction

The South American Pampa, or “Pampas,” is a Neotropical biome dominated by natural grasslands spreading over plains of Uruguay, Northern Argentina, Southern Brazil, and part of Paraguay (Pallarés et al., 2005). However, from the biological standpoint, the Pampas are far from homogeneous. Indeed, many proposals for small biological “provinces” have been made, especially regarding on the different plant communities (Cabrera and Willink, 1980; Overbeck et al., 2007), which reflects the mosaic of soil types resulted from its complex geological history. At a smaller scale, such heterogeneity leads to high species endemism and significant genetic structure in Pampean species studied so far (e.g., Freitas et al., 2012; Turchetto et al., 2014; Felappi et al., 2015). While some phylogeographic studies in the Pampas have found little geographic structure and shallow gene trees [<0.1 million years ago (mya)] (Speranza et al., 2007; Turchetto et al., 2014), other studies found strong geographic structure and a deep mitochondrial gene tree (>2.5 mya) (Felappi et al., 2015). However, these studies did not include freshwater organisms, and, therefore, it is difficult to predict the level and depth of the genetic structure exhibited by these species. Species occurring in the Pampa biome have been underrepresented in studies of phylogeography and conservation genetics (Lawler et al., 2006; Beheregaray, 2008; Turchetto-Zolet et al., 2013), which further complicates the understanding of general drivers of biological diversification in this biome.

The Pampas hydrological system resulted from a long history of tectonism, climate, and sea level changes since the Neogene (Casciotta et al., 1999). During the Quaternary, soil erosion and other marine regression/transgression cycles have continued to shape hydrological systems (Tonni and Cione, 1997; Quattrocchio et al., 2008) by altering the relationship among tributaries, isolated drainages, lagoons, and estuaries (Martin and Dominguez, 1994; Thomaz et al., 2015). These geomorphological processes promoted successive stream captures events, allowing obligate freshwater species to disperse using temporary connections (e.g., Loureiro et al., 2011; Candela et al., 2012; Grassi et al., 2017). Stream capture (also known as river capture, headwater capture, or drainage rearrangement) is a geomorphological process that consists in the contact between neighboring drainages. It occurs when the tributary of a river basin starts flowing toward a neighbor basin, and results in icthyofauna dispersal between them (Bishop, 1995). Temporary connections between adjacent basins caused by fluctuations in the sea level during glacial–interglacial cycles may have also influenced the current distribution of several freshwater fish genera (e.g., Ketmaier et al., 2004; Swartz et al., 2007, 2009; Thomaz et al., 2015), including the Pampean taxa Australoheros, Cnesterodon, Jenynsia, and Corydoras (e.g., Bruno et al., 2011, 2016; Ponce et al., 2011).

It is well accepted that the distribution of freshwater fish lineages mainly reflects the paleogeography of a specific region (Bermingham and Avise, 1986; Bernatchez and Wilson, 1998; Avise, 2000; Lévêque et al., 2008). Cnesterodon decemmaculatus (Jenyns 1842) is endemic to the Pampa biome and is one of the most widespread freshwater fish species in this region. It is found in the freshwater Ecoregions 332 (lower Uruguay) and 334 (laguna dos Patos basin) described by Abell et al. (2008), but exclusively in grassland environments in the southern part of these ecoregions, associated to the Pampa biome, being absent in the northern portion of these ecoregions located in the Atlantic Forest Biome. It is further found in the freshwater ecoregions 345 (Subtropical Potamic axis) and 347 (Bonaerensean Atlantic), also in the Pampas. Thus, this species constitutes an excellent model to examine historical processes that played a major role in shaping the genetic structure of freshwater fishes in the Pampas. In addition, this species can be found in a wide range of habitats, including rivers, ponds, and shallow wetlands, even though it has a poor swimming capacity in high-speed currents (Trenti et al., 1999). A recent study using C. decemmaculatus populations from the southern Pampean region showed low genetic divergence and evidences of recent colonization of this area (Bruno et al., 2016). However, this study focused on the coastal drainages from Argentina (ecoregion 347 – Bonaerensean Atlantic), which represent a small portion of the natural distribution of C. decemmaculatus in the Pampas.

In this study, we used mitochondrial and nuclear genetic markers to evaluate the genetic structure of C. decemmaculatus over its distribution. More specifically, we were interested in answer the following questions: (1) How strong is the correspondence between current drainage systems and evolutionary lineages? (2) Is there evidence of ancient population growth or reduction for each genetic lineage in this species? (3) What is the temporal depth for the evolutionary divergence among lineages and populations? (4) What is the colonization history for this taxon in the Pampas?

Materials and Methods

We sampled 99 individuals from 37 localities (Figure 1), covering a significant part of the species’ range. Our sampling ranged from one to eight individuals per locality (Supplementary Table S1). All individuals included in our analyses were fixed in 96% ethanol and deposited in the fish collection at the Federal University of Rio Grande do Sul (Universidade Federal do Rio Grande do Sul – UFRGS). All collections were performed under the approval of the government authorities of Brazil and Uruguay (Ministério do Meio Ambiente, Brazil – SISBIO, license number 12038-2, and Dirección de Recursos Naturales del Ministerio de Ganadería, Agricultura y Pesca, Uruguay). The collection and euthanasia of the specimens were approved by the ethics committee of the UFRGS, license number 24434.

FIGURE 1

FIGURE 1. Geographical location of sampled sites showing the limits of the major drainage systems (assumed as different populations). The dashed gray line indicates the distribution of Cnesterodon decemmaculatus according to Lucinda, 2005. Vouchers of the specimens used in the study, and the geographical coordinates for sampling sites are available in Supplementary Table S1.

Total genomic DNA was isolated from muscle tissue with cetyltrimethyl ammonium bromide (CTAB) as described by Doyle and Doyle (1987). We tested two mitochondrial (mtDNA) genes: cytochrome oxidase subunit I (COI) (Herbert et al., 2003) and NADH dehydrogenase 2 (ND2) (Sorenson, 2003) and two nuclear (nDNA) genes: SH3 and PX domain containing 3-like protein (SH3PX3) and myosin heavy chain 6 (Myh6) (Li et al., 2007). However, SH3PX3 and COI had insufficient variation in a preliminary sample of 10 individuals per basin and, therefore, were excluded from further analysis. PCR amplification protocols for ND2 and Myh6 followed Sorenson (2003) and Li et al. (2007), respectively. PCR products were checked on a 1% agarose gel, purified with Exonuclease I and Shrimp Alkaline Phosphatase (GE Healthcare^®) and sequenced using the Sanger method at Macrogen Inc. (Seoul, South Korea). Both DNA strands were sequenced, checked, and aligned with GeneiousPro v.4.8 (Drummond et al., 2007). Haplotypes of the Myh6 gene were estimated in PHASE 2.1 (Stephens et al., 2001) using 10,000 steps sampling every 10 steps and discarding the first 1000 steps as burnin. We ran PHASE several times using different starting seeds to ensure the reliability of the final estimate. Mitochondrial haplotypes were defined using DnaSP 5.10 (Rozas et al., 2003). All sequences obtained in the present study were deposited in GenBank (KU214332–KU214430 and KU214252–KU214331, for ND2 and Myh6 genes, respectively).

Dispersal of obligatory freshwater fish requires connections of aquatic habitats between adjacent basins and, most of taxa have their distribution range narrowly coincident, if not equal, with the hydrographic basins boundaries in which they live (Ward et al., 1994; Albert and Reis, 2011). Based on this, and also because we found a general genetic structure associated to major Pampean drainages, we defined four geographical groups for population analysis: Uruguay River Basin (URU), Negro River Basin (NEG), Mirim Lagoon Basin (MIR), and Southern coastal Basins of Uruguay (SOU) (Figure 1 and Supplementary Table S1). Because independent drainage systems represent plausible boundaries for isolated biological populations for freshwater organisms, individuals collected in different localities within each basin were merged into a single population, resulting in a sample of 62 individuals for URU, 11 for NEG, 15 for MIR, and 11 for SOU. Following Abell et al. (2008), URU and NEG represent the freshwater Ecoregion 332 (lower Uruguay), MIR represents Ecoregion 334 (laguna dos Patos), while SOU have populations located in both Ecoregions 334 and 345 (lower Paraná).

To determine the correspondence between current drainage systems and evolutionary lineages, we inferred the evolutionary relationship among haplotypes using a median-joining network (Bandelt et al., 1999) estimated in NETWORK 4.1.0.9¹. The resulting genetic structure was quantified using an analysis of molecular variance (AMOVA) (Excoffier et al., 1992) and pairwise Φ_STs calculated in Arlequin 3.5 (Excoffier and Lischer, 2010). In addition, we tested for population structure without assuming any hypothesis based on hydrography using the Bayesian clustering algorithm implemented in BAPS (Corander et al., 2008). We used the mixture population model and set the maximum number of populations (K) to 37, the number of sampled local populations. Because stream capture events may affect gene genealogies and molecular diversity patterns within drainages, we decided to run several analyses excluding localities for which we may have evidence for stream capture. We inferred a possible stream capture event whenever we found shared haplotypes between neighbor drainages, especially when populations closely related to watersheds were involved (see the section “Discussion” for further details for each case).

To test for evidence of ancient population growth or reduction for each population (URU, NEG, MIR, and SOU) we first estimated descriptive statistics of genetic diversity, such as nucleotide (π) and haplotype diversity (H), followed by Tajima’s D (Tajima, 1983) and Fu’s F_S (Fu, 1997) neutrality tests, which were calculated in Arlequin 3.5 (Excoffier and Lischer, 2010). We performed neutrality tests under two schemes: considering all individuals for a given basin, and excluding localities for which there was evidence of headwater capture. We also estimated the effective population size (N_E) and the population growth parameter (G) for each population without the localities with evidence for headwater capture (namely URU 7, URU 8, NEG 35, NEG 36, MIR 31), using both markers in the program LAMARC v. 2.1.6 (Kuhner, 2006). We performed a maximum-likelihood search based on three replicates of 10 initial chains and two final chains. The initial chains consisted of 250 samples drawn every 20 steps, and a burn-in of 1000 samples for each chain. The two final chains used the same burn-in and sampling interval, but consisted of 10,000 samples. We assumed a generation time of 1 year based on estimates of sexual maturity for the poeciliid Poecilia reticulata (Reznick and Bryga, 1987).

Finally, we inferred the divergence times among populations and the colonization history for C. decemmaculatus based on two approaches. First, we estimated a time-calibrated gene-tree genealogy for mtDNA lineages using Bayesian inference in BEAST 1.7.5 (Drummond and Rambaut, 2007). We used a coalescent constant-size tree prior with a random starting tree and the TN93+G model of sequences evolution, as determined by the corrected Akaike Information Criterion (AICc) in jModelTest2 (Darriba et al., 2012). We used each drainage as a discrete trait to allow the estimation of the most likely location of all ancestors in the mtDNA tree. While this may be indicative of dispersal events, it should not be viewed as a formal test of stream capture hypotheses. The MCMC was run twice for 10 million generations each, and 10% of the samples was discarded as burn-in. A normal strict molecular clock was assumed, calibrated with a normal distribution for molecular substitution rate of 8.6 ± 0.1 × 10^-9 substitutions per site per year (s/s/y), which has been suggested in the literature for Cyprinodontiformes based on a dataset which included the ND2 gene (Hrbek and Meyer, 2003).

Second, we used a coalescent-based species-tree analysis (Heled and Drummond, 2010) to infer the colonization history of each river basin using the StarBeast model in the program BEAST 1.7.5 (Drummond and Rambaut, 2007). For this analysis, we used populations as terminals, and used both markers with unlinked gene-trees (mitochondrial – ND2 and nuclear – Myh6). We excluded localities showing evidence for headwater capture, as the species-tree method does not allow gene flow among terminals. We used a rate-reference prior to calibrate the molecular clock for Myh6 (Ferreira and Suchard, 2008). Under this approach, a time-calibrated genealogy (i. e., ND2) is taken as a reference to estimate the relative molecular clock associated with alternative gene trees (i.e., Myh6). This is a useful approach when different partitions evolve at different rates, but there is no prior information on the molecular clock for all independent partitions in the matrix (Ferreira and Suchard, 2008). All remaining parameters and priors followed the Bayesian phylogenetic analyses mentioned above. In all Bayesian analysis, sampling sufficiency was evaluated by monitoring effective sample size (ESS) and ensuring that all values were higher than 200.

Results

Ninety-nine sequences were obtained for the mtDNA ND2 gene (980 bp) yielding 31 different haplotypes defined by 42 variable sites (28 parsimoniously informative). For the autosomal Myh6 gene we had 188 sequences (681 bp) resulting in 25 different haplotypes defined by 21 variable sites (16 parsimoniously informative). In general, the mtDNA haplotype network showed little haplotype sharing among populations, and a close relationship among haplotypes from the same population (Figure 2A). Populations of MIR and SOU exhibited related haplotypes, except for two SOU localities that had haplotypes H21 and H30, which are more closely related to haplotypes found in URU. Concerning Myh6 haplotypes (Figure 2B), URU, MIR, and SOU shared the most frequent haplotypes (H5 and H1), and H9 was the only haplotype present in populations of all drainage systems. The close relationship among haplotypes for both markers is suggestive of a shallow genealogical history for this species (see below).

FIGURE 2

FIGURE 2. Median-joining networks for (A) mtDNA ND2 haplotypes, and (B) nDNA Myh6 haplotypes. Circles size are proportional to the observed frequency of each haplotype. Cross marks represent mutational differences between haplotypes. (C) Geographic location of each mtDNA haplotype. Different colors represent the six genetic groups identified in BAPS analysis.

Considering all samples, the Bayesian clustering analysis found K = 6 genetic populations (Figure 2C). There was a marked difference in the number of different genetic population occurring in each major drainage. Five genetic populations occurred in URU, three in NEG, and two in both SOU and MIR. Not considering samples possibly affected by stream capture events, we found four genetic populations in URU, two in SOU, and one in both NEG and MIR, suggesting significant genetic structure within URU mainly due to three genetic groups with a restricted geographical distribution in this drainage (Figure 2C).

The mtDNA genealogy for C. decemmaculatus (Figure 3) suggests the most recent common ancestor (MRCA) of all mitochondrial lineages dates from 0.8 mya [95% credible interval (CI) 0.5–1.2 mya], and were most probably located in URU (PP = 0.58). A well-supported clade (posterior probability PP = 1.00), sister to all other lineages, occurred in isolated populations inhabiting elevated plains from three basins: URU, MIR, and NEG (haplotypes H11, H12, H14, and H19). These haplotypes represented one of the Bayesian genetic populations. However, the location of the MRCA for this clade was uncertain (PP < 0.5 for all locations). Its sister clade, which had moderate support (PP = 0.82) had URU as the most likely location (PP = 0.75), and showed H25 (another Bayesian genetic population) sister to a well-supported clade (PP = 0.99) that also has URU as its most likely location (PP = 0.80). In turn, this clade harbored four Bayesian genetic populations that corresponded to three well-supported clades (PP > 0.90) plus a “mixed” genetic population that included H1, H7, H8, H30, and H21. Regarding the most likely ancestral location, all remaining genetic groups had URU as the most likely location (PP > 0.85), with one exception. The clade, which was almost exclusive to MIR and SOU (though H10 also occurred in NEG), had SOU or MIR as its most likely ancestral location (PP = 0.46, PP = 0.30, respectively).

FIGURE 3

FIGURE 3. Maximum clade credibility tree among ND2 lineages found in Cnesterodon decemmaculatus individuals. Posterior probabilities (PP) are shown above the branches for selected nodes with PP > 0.8. Pie charts beside selected ancestral node represent the PP of ancestral location among each of the four populations according to the color scheme shown in the inlet.

Considering the full dataset, both markers showed high haplotype diversities (H = 0.94 for ND2 and H = 0.79 for Myh6). Both haplotype and nucleotide diversity were higher for URU and SOU compared with NEG and MIR (Table 1). There was significant genetic population structure (ND2 Φ_ST = 0.26, and Myh6 Φ_ST = 0.30, P < 0.001). Excluding localities that may represent stream capture events increased Φ_ST estimates (ND2 Φ_ST = 0.59, and Myh6 Φ_ST = 0.35, P < 0.001). For ND2, this was still lower than the Φ_ST estimate among the genetic populations delimited with BAPS (Φ_ST = 0.75, P < 0.001), which is expected, given that Bayesian populations were defined based on ND2 data. However, this trend was not followed by Myh6 (Φ_ST = 0.19, P < 0.001). Pairwise Φ_ST values indicated that geographic groups were significantly different from each other except MIR vs. SOU (Table 2). Taking into account possible events of stream capture resulted in increased Φ_ST values between NEG and the remaining populations for both genetic markers, and between the other pairs only for ND2 (Table 2). Because localities from URU might represent different genetic populations, we also quantified the genetic structure between them. For both ND2 and Myh6, we found significant genetic structure in URU (0.81 and 0.39, respectively, P < 0.001), which was robust to the exclusion of some localities due to stream capture events (0.64 and 0.31, respectively, P < 0.001).

TABLE 1

TABLE 1. Diversity indexes obtained for the mitochondrial and nuclear markers, respectively (ND2 values/Myh6 values), for the four populations considered in the study.

TABLE 2

TABLE 2. Pairwise Φ_ST values for mtDNA ND2 (lower diagonal) and nDNA Myh6 (upper diagonal) markers.

Neutrality tests provided mixed evidence for population growth in some populations. However, they were dependent on which statistic or genetic marker were used (Table 1). On the other hand, maximum likelihood-based demographic estimates revealed population growth for SOU and MIR (Table 3). Effective population size was also larger for SOU and MIR (∼3 × 10⁶ effective individuals), while NEG had the lowest value (∼180,000 effective individuals; Table 3 and Figure 4). The combined analysis of mtDNA and nDNA markers suggested that the split between NEG and the clade URU+MIR+SOU, which was also the root of the tree, dates back from ∼0.6 mya (95% CI 0.27–0.95 mya). This was followed by the divergence between URU and MIR+SOU around 0.12 mya (95% CI 0.04–0.19 mya), and the subsequent divergence between MIR and SOU <30,000 years ago (95% CI 23,000–57,000) (Figure 4).

TABLE 3

TABLE 3. Estimates for the effective population size (N_E) and growth parameter (G) in all populations.

FIGURE 4

FIGURE 4. Maximum clade credibility tree considering ND2 and Myh6 under a coalescent approach. The thickness of the branches is proportional to effective population size (N_E) estimates for each geographical group (see Table 3 for more details). The color scheme for each branch followed the estimates of the most likely location for ancestral nodes, shown in Figure 3.

Discussion

In this study, we characterized the major genetic structure patterns of populations of C. decemmaculatus in most of its distribution in the Pampa biome. Differently from the star-like pattern found for Argentinean populations of C. decemmaculatus in the southern Pampas (Bruno et al., 2016), which indicates a very recent evolutionary history, our results reveal a significant genetic structure that, in general, paralleled broad drainage systems, especially considering mtDNA lineages (Figures 2A, 3). However, there is no simple relationship between clades and basins, such that a single basin may contain more than one mtDNA clade, and different basins may share the same mtDNA haplotype (Figures 2A, 3). Thus, understanding population structure imposed by the drainage system is necessary, but not sufficient, to account for the genetic patterns observed in our study. This is reinforced by the Bayesian clustering of populations, which suggested K = 6 as the bet number of genetic populations. Three alternative (but not exclusive) scenarios that may explain the discrepancy between drainage and genetic structure are: (1) migration among neighbor populations (active dispersal), (2) headwater capture (geodispersal), or (3) shared ancestral polymorphism. The first scenario seems unlikely because C. decemmaculatus has a poor swimming capacity in high-speed currents (Trenti et al., 1999), being unable to disperse actively across long distances. On the other hand, Cnesterodon species are commonly found in shallow wetland habitats (Maltchik et al., 2014), suggesting that geodispersal may represent a more likely alternative for these species. In this case, wetland habitats close to the watersheds between different drainages may facilitate headwater capture events, as minor environmental changes may reshape the hydrological network at these sites, impeding a clear-cut relationship between genetic lineages and drainages.

However, distinguishing between the geodispersal and ancestral polymorphism hypotheses may be more nuanced. This is further complicated by the fact that it is very difficult to set up an explicit test of stream capture hypotheses. One example of a likely case of stream capture involving closely related haplotypes occurring in different drainages is represented by mtDNA haplotypes H11, H14, and H19 (Figure 3). Even though URU has received higher support as the most probable location of the MRCA, mtDNA network analysis indicates that this clade evolved in NEG with subsequent stream capture events toward URU and MIR. First, there are seven exclusive substitutions in this clade (Figure 2A), suggesting it may have evolved in relative isolation. Second, the sampling points in URU and MIR containing these lineages are adjacent to NEG. Third, H14 and H19 are probably descendent from H11, which occurs in NEG (Figure 2A). Under this interpretation, this clade would reflect three independent stream capture events: NEG would have been colonized from URU between 0.2 and 0.8 mya (see below), while MIR and URU would have received migrants from NEG more recently (<0.2 mya). Even in the alternative scenario in which URU is the true location for the MRCA of this clade, two headwater capture events would have occurred, since this lineage would have reached NEG from URU by ∼0.2 mya, dispersing to MIR later on. Two other examples of haplotype sharing that may reflect stream capture toward NEG include H8 (from URU) and H10 (from SOU). In both cases, these haplotypes are more closely related to other lineages from specific drainages, and their sequence identity suggests a recent arrival in NEG. However, while H10 occurs in a site close to SOU, H8 occurs relatively distant from the watershed with URU. An alternative hypothesis for the occurrence of H8 in NEG would involve upstream active dispersal from lower URU populations despite the low dispersal ability of Cnesterodon (Trenti et al., 1999).

The tectonic activity associated with stream capture events in eastern South America may be as old as 1.6 mya (Saadi et al., 2002). The presence of H11 and H19 in MIR and H8 in NEG is in agreement with a pattern in which more eastern drainages capture upland shield drainages. However, the presence of H14 in URU and H10 in NEG indicates that these events could also occur in the opposite direction, suggesting a highly dynamic scenario in the Pampas, possibly because this biome is largely flat. Indeed, Loureiro et al. (2011) suggested drainage rearrangements in both directions across the MIR/NEG watershed based on the distribution pattern of the killifish Austrolebias in Uruguay. Similarly to Cnesterodon, Austrolebias inhabit shallow wetland habitats, and have little capabilities for active dispersal (Maltchik et al., 2014), reinforcing the relevance of passive geodispersal in the Pampas. In contrast, the most parsimonious explanation for the presence of haplotypes H21 and H30 in SOU is incomplete lineage sorting, given the low phylogenetic signal associated with location and clade structure for these haplotypes, and because the sampling points in which these haplotypes were found are distant from the URU watershed.

Our results suggest that the MRCA of C. decemmaculatus reached the lower Uruguay before 0.6 myr (Figure 4). The mtDNA genealogy suggests that from URU, C. decemmaculatus would have colonized NEG and SOU basins independently (Figure 3). The high genetic diversity in URU is consistent with its role as the ancestral location. Even though we have found a case of headwater capture from NEG to MIR, the most likely colonization route leading to MIR is through SOU, given that most haplotypes found in SOU and MIR belong to the same mtDNA haplogroup. The Bayesian delimitation of genetic populations also suggested that most localities in MIR and SOU belong to the same genetic population, which is also in agreement with pairwise indices of genetic structure (Table 2). These finding also highlight that, for this species, the genetic diversity in the lower Uruguay ecoregion would be much higher than in Laguna dos Patos (sensu Abell et al., 2008). Furthermore, given that URU is neighbor to SOU, but not MIR, a colonization route from URU to MIR and then SOU is less likely. Indeed, Bruno et al. (2016) suggested that the Argentinean coastal populations of C. decemmaculatus descend from upstream URU populations associated with the Rio de la Plata mouth, which highlight the mouth of the Uruguay river as a putative source for populations currently inhabiting isolated drainages flowing to Río de la Plata. Alternatively, SOU/MIR could have been colonized from NEG. This hypothesis would be strengthened if we assume that H10 evolved in situ in NEG before dispersing eastward. Additional sampling efforts, especially in NEG, but also in regions not sampled by this study, such as the Argentinian coast, will be required to refine these scenarios, as well as explicitly testing for hypotheses of stream capture or incomplete lineage sorting.

While the Pampean region has been always dominated by grasslands during the whole Pleistocene (Behling et al., 2005), changes in precipitation regimes and in the sea-level may have affected populations of C. decemmaculatus. For example, it could be expected that periods of increased humidity would have favored the formation of new wetland areas and promoted population expansion (Zemlak et al., 2008; Jones and Johnson, 2009). On the other hand, marine transgressions would have caused local extinctions (Villwock and Tomazelli, 1995). In this regard, the high genetic diversity in SOU was surprising, since this area suffered at least three marine transgression events since the Pleistocene (Villwock and Tomazelli, 1995) that could have led to declines in population size and local extinction. The high N_E and the significant negative values of Fu’s neutrality test for mtDNA indicate that either local extinction did not affect these populations, or that the high genetic diversity resulted from recent and strong population growth (Table 3). We also found evidence for population growth in MIR. However, differences in genetic diversity and N_E estimates between MIR and SOU may reflect different population histories after marine transgressions. Alternatively, some founder effect may have reduced genetic diversity in MIR following its colonization from SOU. Coastal populations from SOU may have also benefited from an extended coastal plain during marine retractions, facilitating population growth during these periods, as have been suggested for other species occurring along the coastal rivers of southern and southeastern Brazil (Thomaz et al., 2015).

Author Contributions

AR-F, LM, and NF designed the study, interpreted the results, and contributed to the final version of the manuscript. AR-F did the laboratory work. AR-F and NF analyzed the data.

Funding

This research was supported by Conselho Nacional de Desenvolvimento Científico (CNPq), Programa de Pós-graduação em Biologia Animal (PPGBAN-UFRGS), and Programa de Pós Graduação em Genética e Biologia Molecular (PPGBM-UFRGS).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to thank Marcelo Loureiro Barrela (Universidad de la República, Uruguay), Juliana Wingert, Alice Hirschmann, and Juliano Ferrer, who provided specimens for study or assistance during field work. They also thank Jéferson Fregonezi, Andréa Thomaz, and three reviewers for contributions in a preliminary version of this manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2017.00214/full#supplementary-material

Footnotes

1. ^ www.fluxus-engineering.com

References

Abell, R., Thieme, M. L., Revenga, C., Bryer, M., Kottelat, M., Bogutskaya, N., et al. (2008). Freshwater ecoregions of the world: a new map of biogeographic units for freshwater biodiversity conservation. BioScience 58, 403–414. doi: 10.1641/B580507

CrossRef Full Text | Google Scholar

Albert, J. S., and Reis, R. (2011). Historical Biogeography of Neotropical Freshwater Fishes. Berkeley, CA: University of California Press. doi: 10.1525/california/9780520268685.001.0001

CrossRef Full Text | Google Scholar

Avise, J. C. (2000). Phylogeography: The History and Formation of Species. Cambridge: Harvard University Press.

Google Scholar

Bandelt, H. J., Forster, P., and Rohl, A. (1999). Median-joining networks for inferring intraspecific phylogenies. Mol. Biol. Evol. 16, 37–48. doi: 10.1093/oxfordjournals.molbev.a026036

PubMed Abstract | CrossRef Full Text | Google Scholar

Beheregaray, L. B. (2008). Twenty years of phylogeography: the state of the field and the challenges for the Southern Hemisphere. Mol. Ecol. 17, 3754–3774. doi: 10.1111/j.1365-294X.2008.03857.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Behling, H., Pillar, V., and Bauermann, S. G. (2005). Late quaternary grassland (Campos), gallery forest, fire and climate dynamics, studied by pollen, charcoal and multivariate analysis of the São Francisco de Assis core in western Rio Grande do Sul (southern Brazil). Rev. Palaeobot. Palynol. 133, 235–248. doi: 10.1016/j.revpalbo.2004.10.004

CrossRef Full Text | Google Scholar

Bermingham, E., and Avise, J. C. (1986). Molecular zoogeography of freshwater fishes in the southeastern United States. Genetics 113, 939–965.

PubMed Abstract | Google Scholar

Bernatchez, L., and Wilson, C. C. (1998). Comparative phylogeography of nearctic and palearctic freshwater fishes. Mol. Ecol. 7, 431–452. doi: 10.1046/j.1365-294x.1998.00319.x

CrossRef Full Text | Google Scholar

Bishop, P. (1995). Drainage rearrangement by river capture, beheading and diversion. Progr. Phys. Geogr. 19, 449–473. doi: 10.1177/030913339501900402

CrossRef Full Text | Google Scholar

Bruno, M. C., Lizarralde, M., Almirón, A., and Casciotta, J. (2011). Presence of Australoheros facetus (Teleostei: Cichlidae) in the southern Pampean area. Considerations about the ichthyofaunal settlement and distribution in the Southern boundary of the Brazilian Subregion. Ichthyol. Contrib. Peces Criollos. 23, 1–3.

Google Scholar

Bruno, M. C., Mapelli, F. J., Casciotta, J. R., Almirón, A. E., and Lizarralde, M. S. (2016). Phylogeography of Cnesterodon decemmaculatus (Cyprinodontiformes: poeciilidae) in Southern Pampas, Argentina: ancient versus recent patterns in freshwater fishes. Environ. Biol. Fish. 99, 293–307. doi: 10.1007/s10641-016-0474-0

CrossRef Full Text | Google Scholar

Cabrera, A. L., and Willink, A. (1980). Biogeografía de América Latina. Washington, DC: OEA.

Google Scholar

Candela, A. M., Bonini, R. A., and Noriega, J. I. (2012). First continental vertebrates from the marine Parana Formation (Late Miocene, Mesopotamia, Argentina): chronology, biogeography, and paleoenvironments. Geobios 45, 515–526. doi: 10.1016/j.geobios.2012.05.003

CrossRef Full Text | Google Scholar

Casciotta, J., Almirón, A., Cione, A., and Azpelicueta, M. (1999). Brazilian freshwater fish assemblages from southern Pampean area, Argentina. Biogeographica. 75, 67–78.

Google Scholar

Corander, J., Marttinen, P., Sirén, J., and Tang, J. (2008). Enhanced Bayesian modelling in BAPS software for learning genetic structures of populations. BMC Bioinformatics 9:539. doi: 10.1186/1471-2105-9-539

PubMed Abstract | CrossRef Full Text | Google Scholar

Darriba, D., Taboada, G. L., Doallo, R., and Posada, D. (2012). jModelTest 2: more models, new heuristics and parallel computing. Nat. Methods 9:772. doi: 10.1038/nmeth.2109

PubMed Abstract | CrossRef Full Text | Google Scholar

Doyle, J. J., and Doyle, J. L. (1987). A rapid isolation procedure for small quantities of fresh leaf tissue. Phytochem. Bull. 9, 11–15.

Google Scholar

Drummond, A. J., Ashton, B., Cheung, M., Heled, J., Kearse, M., Moir, R., et al. (2007). Geneious v4.8.

Google Scholar

Drummond, A. J., and Rambaut, A. (2007). BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol. 7:214. doi: 10.1186/1471-2148-7-214

PubMed Abstract | CrossRef Full Text | Google Scholar

Excoffier, L., and Lischer, H. E. L. (2010). Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567. doi: 10.1111/j.1755-0998.2010.02847.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Excoffier, L., Smouse, P. E., and Quattro, J. M. (1992). Analysis of molecular variance inferred from metric distances among DNA haplotypes – Application to human mitochondrial-DNA restriction data. Genetics 131, 479–491.

PubMed Abstract | Google Scholar

Felappi, J. F., Vieira, R. C., Fagundes, N. J., and Verrastro, L. V. (2015). So far away, yet so close: strong genetic structure in Homonota uruguayensis (Squamata, Phyllodactylidae), a species with restricted geographic distribution in the Brazilian and uruguayan pampas. PLOS ONE 10:e0118162. doi: 10.1371/journal.pone.0118162

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferreira, M. A., and Suchard, M. A. (2008). Bayesian analysis of elapsed times in continuous-time Markov chains. Can. J. Stat. 36, 355–368. doi: 10.1002/cjs.5550360302

CrossRef Full Text | Google Scholar

Freitas, T. R. O., Fernandes, F. A., Fornel, R., and Roratto, P. A. (2012). An endemic new species of tuco-tuco, genus Ctenomys (Rodentia: Ctenomyidae), with a restricted geographic distribution in Southern Brazil. J. Mammal. 93, 1355–1367. doi: 10.1644/12-MAMM-A-007.1

CrossRef Full Text | Google Scholar

Fu, Y. X. (1997). Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147, 915–925.

PubMed Abstract | Google Scholar

Grassi, D. J., Swarça, A. C., Dergam, J. A., Pastori, M. C., and Fenocchio, A. S. (2017). Cytogenetic characterization of Hoplias malabaricus (Bloch, 1794) from the Ctalamochita River (Córdoba, Argentina): first evidence for southernmost populations of this species complex and comments on its biogeography. Comp. Cytogenet. 11, 15–28. doi: 10.3897/CompCytogen.v11i1.10262

PubMed Abstract | CrossRef Full Text | Google Scholar

Heled, J., and Drummond, A. J. (2010). Bayesian inference of species trees from multilocus data. Mol. Biol. Evol. 27, 570–580. doi: 10.1093/molbev/msp274

PubMed Abstract | CrossRef Full Text | Google Scholar

Herbert, P. D. N., Cywinska, A., Ball, S., and Dewaard, J. R. (2003). Biological identifications through DNA barcodes. Proc. R. Soc. B 270, 313–321. doi: 10.1098/rspb.2002.2218

PubMed Abstract | CrossRef Full Text | Google Scholar

Hrbek, T., and Meyer, A. (2003). Closing of the Tethys Sea and the phylogeny of Eurasian killifishes (Cyprinodontiformes: Cyprinodontidae). J. Evol. Biol. 16, 17–36. doi: 10.1046/j.1420-9101.2003.00475.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, C. P., and Johnson, J. B. (2009). Phylogeography of the livebearer Xenophallus umbratilis (Teleostei: Poeciliidae): glacial cycles and sea level change predict diversification of a freshwater tropical fish. Mol. Ecol. 18, 1640–1653. doi: 10.1111/j.1365-294X.2009.04129.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ketmaier, V., Bianco, P. G., Cobolli, M., Krivokapic, M., Caniglia, R., and Matthaeis, E. (2004). Molecular phylogeny of two lineages of Leuciscinae cyprinids (Telestes and Scardinius) from the peri-Mediterranean area based on cytochrome b data. Mol. Phylogenet. Evol. 32, 1061–1071. doi: 10.1016/j.ympev.2004.04.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuhner, M. K. (2006). LAMARC 2.0: maximum likelihood and Bayesian estimation of population parameters. Bioinformatics 22, 768–770. doi: 10.1093/bioinformatics/btk051

PubMed Abstract | CrossRef Full Text | Google Scholar

Lawler, J. J., Aukema, J. E., Grant, J. B., Halpern, B. S., Kareiva, P., Nelson, C. R., et al. (2006). Conservation science: a 20-year report card. Front. Ecol. Environ. 4:473–480. doi: 10.1890/1540-9295(2006)4[473:CSAYRC]2.0.CO;2

CrossRef Full Text | Google Scholar

Lévêque, C., Oberdorff, T., Paugy, D., Stiassny, M. L. J., and Tedesco, P. A. (2008). Global diversity of fish (Pisces) in freshwater. Hydrobiologia 595, 545–567. doi: 10.1007/s10750-007-9034-0

CrossRef Full Text | Google Scholar

Li, C., Ortí, G., Zhang, G., and Lu, G. (2007). A practical approach to phylogenomics: the phylogeny of ray-finned fish (Actinopterygii) as a case study. BMC Evol. Biol. 7:44. doi: 10.1186/1471-2148-7-44

PubMed Abstract | CrossRef Full Text | Google Scholar

Loureiro, M., Duarte, A., and Zarucki, M. (2011). A new species of Austrolebias costa (cyprinodontiformes: rivulidae) from northeastern Uruguay, with comments on distribution patterns. Neotrop. Ichthyol. 9, 335–342. doi: 10.1590/S1679-62252011000200010

CrossRef Full Text | Google Scholar

Maltchik, L., Lanés, L. E. K., Keppeler, F. W., Rolon, A. S., and Stenert, C. (2014). Landscape and habitat characteristics associated with fish occurrence and richness in southern Brazil palustrine wetland systems. Environ. Biol. Fish. 97, 297–308. doi: 10.1007/s10641-013-0152-4

CrossRef Full Text | Google Scholar

Martin, L., and Dominguez, J. M. L. (1994). “Geological history of coastal lagoons,” in Coastal Lagoon Processes, ed. B. Kjerfve (Amsterdam: Elsevier), 41–68.

PubMed Abstract | Google Scholar

Overbeck, G. E., Muller, S. C., Fidelis, A., Pfadenhauer, J., Pillar, V. D., Blanco, C. C., et al. (2007). Brazil’s neglected biome: the South Brazilian campos. Perspect. Plant Ecol. Evol. Syst. 9, 101–116. doi: 10.1016/j.ppees.2007.07.005

CrossRef Full Text | Google Scholar

Pallarés, O. R., Beretta, E. J., and Maraschin, G. E. (2005). “The south american campos ecosystem,” in Grasslands of the World, eds J. Suttie, S. G. Reynolds, and C. Batello (Rome: Food and agriculture organization of the United Nations), 171–219.

Google Scholar

Ponce, J. F., Rabassa, J., Coronato, A., and Borromei, A. M. (2011). Palaeogeographical evolution of the Atlantic coast of Pampa and Patagonia from the last glacial maximum to the Middle Holocene. Biol. J. Linn. Soc. 103, 363–379. doi: 10.1111/j.1095-8312.2011.01653.x

CrossRef Full Text | Google Scholar

Quattrocchio, M. E., Borromei, A. M., Deschamps, C. M., Grill, S. C., and Zavala, C. A. (2008). Landscape evolution and climate changes in the Late Pleistocene–Holocene, southern Pampa (Argentina): evidence from palynology, mammals and sedimentology. Quat. Int. 181, 123–138. doi: 10.1016/j.quaint.2007.02.018

CrossRef Full Text | Google Scholar

Reznick, D. N., and Bryga, H. (1987). Life-history evolution in guppies (Poecilia reticulata): 1. Phenotypic and genetic changes in an introduction experiment. Evolution 41, 1370–1385. doi: 10.1111/j.1558-5646.1987.tb02474.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rozas, J., Sánchez-DelBarrio, J. C., Messeguer, X., and Rozas, R. (2003). DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics 19, 2496–2497. doi: 10.1093/bioinformatics/btg359

CrossRef Full Text | Google Scholar

Saadi, A., Machette, M. N., Haller, K. M., Dart, R. L., Bradley, L., and Souza, A. M. P. D. (2002). Map and Database of Quaternary Faults and Lineaments in Brazil. U.S. Geological Survey. Open-File Report 02–230. Available at http://pubs.usgs.gov/of/2002/ofr-02-230

Google Scholar

Sorenson, M. D. (2003). Avian mtDNA Primers. Available at: http://people.bu.edu/msoren/Bird.mt.Primers.pdf

Google Scholar

Speranza, P. R., Seijo, G., Grela, I. A., and Neffa, V. G. S. (2007). Chloroplast DNA variation in the Turnera sidoides L. complex (Turneraceae): biogeographical implications. J. Biogeogr. 34, 427–436. doi: 10.1111/j.1365-2699.2006.01622.x

CrossRef Full Text | Google Scholar

Stephens, M., Smith, N. J., and Donnelly, P. (2001). A new statistical method for haplotype reconstruction from population data. Am. J. Hum. Genet. 68, 978–989. doi: 10.1086/319501

PubMed Abstract | CrossRef Full Text | Google Scholar

Swartz, E. R., Skelton, P. H., and Bloomer, P. (2007). Sea-level changes, river capture and the evolution of populations of the Eastern Cape and fiery redfins (Pseudobarbus afer and Pseudobarbus phlegethon, Cyprinidae) across multiple river systems in South Africa. J. Biogeogr. 34, 2086–2099. doi: 10.1111/j.1365-2699.2007.01768.x

CrossRef Full Text | Google Scholar

Swartz, E. R., Skelton, P. H., and Bloomer, P. (2009). Phylogeny and biogeography of the genus Pseudobarbus (Cyprinidae): shedding light on the drainage history of rivers associated with the Cape Floristic Region. Mol. Phylogenet. Evol. 51, 75–84. doi: 10.1016/j.ympev.2008.10.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Tajima, F. (1983). Evolutionary relationship of DNA sequences in finite populations. Genetics 105, 437–460.

Google Scholar

Thomaz, A. T., Malabarba, L. R., Bonatto, S. L., and Knowles, L. L. (2015). Testing the effect of palaeodrainages versus habitat stability on genetic divergence in riverine systems: study of a Neotropical fish of the Brazilian coastal Atlantic Forest. J. Biogeogr. 42, 2389–2401. doi: 10.1111/jbi.12597

CrossRef Full Text | Google Scholar

Tonni, E. P., and Cione, A. L. (1997). Did the Argentine Pampean ecosystem exist in the Pleistocene. Curr. Res. Pleistocene 14, 131–133.

Google Scholar

Trenti, P. S., Gómez, S. E., and Ferriz, R. A. (1999). Capacidad de natación en tres peces pampásicos. Aprona 38, 2–9.

Google Scholar

Turchetto, C., Fagundes, N. J. R., Segatto, A. L. A., Kuhlemeier, C., Neffa, V. G. S., Speranza, P. R., et al. (2014). Diversification in the South American Pampas: the genetic and morphological variation of the widespread Petunia axillaris complex (Solanaceae). Mol. Ecol. 23, 374–389. doi: 10.1111/mec.12632

PubMed Abstract | CrossRef Full Text | Google Scholar

Turchetto-Zolet, A. C., Pinheiro, F., Salgueiro, F., and Palma-Silva, C. (2013). Phylogeographical patterns shed light on evolutionary process in South America. Mol. Ecol. 22, 1193–1213. doi: 10.1111/mec.12164

PubMed Abstract | CrossRef Full Text | Google Scholar

Villwock, J. A., and Tomazelli, L. J. (1995). Geologia costeira do Rio Grande do Sul. Notas Técnicas 8, 1–45.

Google Scholar

Ward, R. D., Woodwark, M., and Skibinski, D. O. F. (1994). A comparison of genetic diversity levels in marine, freshwater, and anadromous fishes. J. Fish. Biol. 44, 213–232. doi: 10.1111/j.1095-8649.1994.tb01200.x

CrossRef Full Text | Google Scholar

Zemlak, T. S., Habit, E. M., Walde, S. J., Battini, M. A., Adams, E. D., and Ruzzante, D. E. (2008). Across the southern Andes on fin: glacial refugia, drainage reversals and a secondary contact zone revealed by the phylogeographical signal of Galaxias platei in Patagonia. Mol. Ecol. 17, 5049–5061. doi: 10.1111/j.1365-294X.2008.03987.x

PubMed Abstract | CrossRef Full Text | Google Scholar