Quaternary refugia and secondary contact in the southern boundary of the Brazilian subregion : comparative phylogeography of freshwater fish

Freshwater fish are an ideal model to illustrate how climate-tectonic changes affect the distribution of genetic variation. Freshwater bodies are extensively affected by environmental changes, with streams even changing their courses in the most extreme cases. Fortunately, this situation is reflected in the genetic composition of populations and may currently be inferred from the study of mitochondrial DNA molecular markers. Here we analyze and compare the phylogeographic patterns of the species Corydoras paleatus and Jenynsia multidentata at the southern limit of the Brazilian subregion. These basins are isolated in the current hydrogeographic pattern due to geologic and paleoclimatic changes. Our results support a concurrent pattern for both species. Some lineages have persisted in the area under adverse climate conditions, possibly in environmental refuges, while other lineages may have colonized the area later by means of paleodrainages connections. In addition, the presence of independent, greatly diverging lineages, even within the same watercourses, suggests secondary contact between these lineages. This work represents a first approach to understand how geologic and paleoclimatic changes have affected the distribution of genetic variation in the Southern Pampean Area.


Introduction
The distribution of genetic variation in freshwater organisms is affected by the configuration and the geologic evolution of the basins where they occur.Freshwater fishes show strong genetic structure as a result of their confinement to given hydrologic systems (Faulks et al. 2010;loxterman & keeley 2012).The diverse factors which influence basin morphology and dynamics may be ma-jor barriers that limit the dispersal capacity of freshwater organisms.When such barriers arise, they restrict gene flow, thus increasing the differentiation of populations.In contrast, when no barriers are present, gene flow tends to homogenize the populations (HutcHinson & templeton 1999;allendorF & luikart 2007;loxterman & keeley 2012).The use of DNA markers in combination with paleoclimatic reconstruction provides new information on the evolution of species during the climatic oscillations of the Quaternary.Such analyses can identify genetic subdivision and its geographical pattern across a species range, and provide attenuated signals of past demographic and geographic changes (Hewitt 2004).
The Neotropical region is divided into the Austral and Brazilian subregions according to fish distribution (casciotta et al. 1999).The richest ichthyofauna occurs in the Brazilian subregion, extends to the drainages south of Buenos Aires province, Argentina.The Southern Pampean Area, in the Buenos Aires province (Fig. 1A) is the southernmost limit of the Bra zilian subregion.The area is an extensive plain grassland whose general monotony is interrupted by two northwest-southeast-trending sierras ranges.These ranges are Ventania (Sierras Australes) and Tandilia (Sierras Se pten trionales).The plain between both ranges is known as Llanura Interserrana Bonaerense (Bonarean Intersierran Plain) (casciotta et al. 1999).Currently, the water courses within the Southern Pampean Area run parallel, and are therefore unconnected to each other; and they are also to the basins that limit this area.The fish fauna of these basins is one of the last assemblies of the Brazilian subregion; despite this, knowledge of the ichthyofauna of this area and of the possible causes that led to its isolation is still fragmentary.The hydrogeographic pattern of this region has been subjected to the Quaternary climate cycles, which have influenced the morphology of the drainage network (tonni & cione 1997;aramayo et al. 2002;QuattroccHio et al. 2008).Paleoclimatic reconstructions suggest an arid-semiarid climate for this area during much of the Quaternary, associated in part with glacial and interglacial periods (tonni et al. 1999;QuattroccHio et al. 2008).This probably affected the flow of many water courses, with the smaller ones even disappearing during the driest periods (nagle & simons 2012).The geologic and paleoclimatic changes that resulted in the current isolation of the studied basins probably also influenced the distribution of the fish fauna.
The general goal of this work is to analyze and compare the phylogeographic pattern of two teleostean species that are co-distributed in the basins of the Southern Pampean Area: Jenynsia multidentata and Corydoras paleatus.The sequences for J. multidentata are those of Bruno et al. (2013) with the addition of new individuals.Both species are widely distributed and have documented fossil records within the area.In addition, we examine the relationship between these patterns and the geomorphological and paleoclimatic history of the region.Thus, this approach aims to infer if such geological and paleoclimatic conditions have influenced the distribution of gene lineages in the ichthyofauna of the Southern Pampean Area.

Study area and sample collections
Specimens were collected in the localities El Durazno, Chocorí, El Moro, Cortaderas streams, Quequén Grande river, Claromecó stream, Quequén Salado and Sauce Grande rivers, Sauce Chico and Saladillo streams (Fig. 1A).For the purpose of these analyses, each locality was considered as a population.
Part of the pectoral fin of each specimen was removed and used for the molecular analyses.Tissue samples were deposited in the tissue bank at the Molecular Ecology Laboratory (Laboratorio de Ecología Molecular, CREG, UNLP).

Molecular analyses
Total DNA was extracted from the preserved tissue samples following the protocol outlined by aljanaBi & martínez (1997).The entire control region of the mitochondrial DNA was amplified by the Polymerase Chain Reaction (PCR), from the specimens of C. paleatus using the primers FLR-l (5'-aactcccaaagctaggattc -3') and FLR-H (5'-gctcgtggaactttctagg -3') which were specifically designed for this study using the software Gene Runner v 3.01 (spruyt & BuQuiccHio 1994).The amplification reaction was carried out in a total volume of 25 µl, with a final concentrations of 0.25 µl of Taq DNA polymerase, 1.5 µl of Cl 2 Mg 3 mM, 0.4 µl of dNTP's 50 mM, 0.25 µl of each primer 10 mM and1 µl of DNA template.The reaction was made under the following conditions: initial denaturation at 95 °C during four minutes, followed by 30 cycles of denaturation at 95 °C during 45 seconds, annealing at 59 °C during 45 seconds and extension at 72 °C during 1.15 minutes, followed by a final extension of 5 minutes.
In the case of J. multidentata it was possible to amplify eleven individuals, which were added to the preexisting ones, using the primers and under the conditions described in Bruno et al. (2013).Negative controls were performed in all cases to verify the absence of contamination.
The amplification products were purified by lithium chloride precipitation and absolute ethanol precipitation (samBrook et al. 1989).The amplicons were sequenced in a capillary sequencer ABI 3100 (Macrogen Inc., Korea).Chromatograms were edited using the software Proseq (Filatov, 2002) and aligned using Clustal W (tHompson et al. 1994).

Genealogic relationships among haplotypes
To evaluate phylogenetic relationships among the haplotypes of the studied species, we used analyses based both on evolutionary methods (Maximum Parsimony and Bayesian Inference) and on genetic distances (Neighbour Joining).The software JModelTest (posada & Buck ley 2004; posada 2008) with Akaike's information criterion (akaike, 1974) was used to select the nucleotide substitution model that best fit the data.For J. multidentata the resulting model was TPM1uf + G.For C. paleatus, the resulting nucleotide substitution model was HKY + G.The parameters obtained for the nucleotide substitution models were used to optimize the data in the phylogenetic reconstruction.
A Maximum Parsimony (MP) analysis was implemented in PAUP * 4.0b10 (swoFFord 2002), performing a heuristic search with 100 random stepwise addition and tree bisection and reconnection (tBr) branch swapping.Secondly, a Bayesian Inference (BI) analysis was performed using the software Mr. Bayes (ronQuist & HuelsenBeck 2003).Four independent runs of 3 × 10 6 generations were made for the two species.The first 40,000 trees were discarded as burn-in for J. multidentata, and the first 30,000 for C. paleatus.The Neighbour-Joining (NJ) algorithm was implemented in PAUP * 4.0b10 (swoFFord 2002).Because of TPM1uf + G is not a model available in PAUP, the GTR model was selected in J. multidentata analyses.Node support for the MP and NJ analyses was tested using 1,000 bootstrap replicates (Felsenstein 1985).Sequences of Poecilia latipinna (GenBank accession number DQ445680.1)and Corydoras rabauti (GenBank accession number 29570739.1)were used as outgroup in the J. multidentata and C. paleatus analyses respectively.
Alternatively the evolutionary relationships between haplotype were analyzed by constructing a haplotype network for each species, applying the Median Joining algorithm (Bandelt et al. 1999) implemented in the software Network 4.5.1 (http://www.fluxus-engineering.com).

Historical demography
Historical demography was assessed by means of Fu's (1997) and Tajima's ( 1989) neutrality tests.These analy-ses were made both globally for the pooled samples of each species, and for each locality.Their significance was tested using 1,000 coalescent simulations in Arlequin v 3.5 (excoFFier & liscHer 2010).
To evaluate for a sudden demographic expansion model a Mismatch distribution analysis (rogers & Harpending 1992) were conducted for the different populations and for the pooled samples using the software Arlequin v 3.5 (excoFFier & liscHer 2010).To test the validity of the population expansion model of the observed vs. the expected mismatch distribution we used a goodness of fit test based on square sum deviation (SSD) and Harpending's Raggedness index (rg) (Harpending 1994).These analyses were also carried out using the software Arlequin v 3.5 (excoFFier & liscHer 2010).

Population structure
A description of the genetic structure was obtained from an Analysis of Molecular Variance (AMOVA) (excoFFier et al. 1992) using Arlequin v 3.5 (excoFFier & liscHer 2010).For this analysis, three hierarchical population groupings were established on the basis of the geological province where each water course was located: Eastern or Tandilia group, Central or Bonaerean Intersierran Plain group, and Western or Ventania group.The statistical significance was assessed using 1,000 permutations.In order to evaluate the existence of migration-genetic drift equilibrium at regional scale, we used the approach proposed by HutcHinson & templeton (1999).Under equilibrium conditions, genetic distances are expected to increase with the increase of geographic distance between populations, in a pattern known as isolation by distance (wrigHt 1943).The absence or weakness of one or both of these relationships would indicate lack of regional equilibrium.The relationship between geographic and genetic distances (F st ) was evaluated by means of a Mantel test (mantel 1967).Statistical significance was assessed using 1,000 permutations.

Genetic variation
A 773 base-pair (bp) fragment from the control region was analyzed, using 110 J. multidentata individuals (Table 1).Thirty-seven polymorphic sites were identified, 22 of which were transitions and eight, transversions.Likewise, seven indels were found, defining a total of 12 haplotypes (GenBank acc number KM086733-KM086735).Total haplotype diversity (h) was 0.54 +/-0.05, while total nucleotide diversity (p) was 0.0096 +/-0.005,indicating low to moderate genetic diversity for the populations.Values for each population are show in Table 1.
A 923-base pair (bp) fragment of the control region was obtained for C. paleatus from 72 individuals.Thirty-six polymorphic sites were observed, including 27 transitions, two transversions and ten indels, to define a total of 14 haplotypes (GenBank acc number KM077147 -KM077160).Total haplotype diversity (Table 2) was moderate (0.580 ± 0.066) and so was nucleotide diversity (0.0048 ± 0.0026).Table 2 presents the values for each population.

Genealogical relationships among haplotypes
The phylogenetic relationships among J. multidentata haplotypes (Fig. 2) showed similar patterns according to evolutionary and distance methods.In both cases, three major clades were found.One include the widely distributed haplotype as well as haplotypes from the localities Chapadmalal (East) and Saladillo (West).A second clade restricted to the East (or Tandilia region) of the geographical range.The third clade is formed by haplotypes that occur in the West locality (Sauce Grande).The main difference between the MP and BI methods concerned haplotype 7, which was independent from these clades according to MP, while it had a basal position within Clade 1 in the BI tree.In any case, this latter relationship had low support as indicated by the posterior probabil-ity of the node (data not shown).The haplotype network was highly concordant with the trees obtained previously.The same terminal groups described in the trees were recovered (Fig. 1B).Four major lineages differentiated by several mutational steps can be observed.Of these, one occurs exclusively at Sauce Grande river (haplotypes 8, 9, 10 and 11).Another lineage is formed by a single individual from Quequén Salado river (haplotype 7).A third group was defined by individuals from the populations Chapadmalal, Chocorí, and El Moro (haplotypes 2, 4, 5, 6), which are part of the Tandilia region.Lastly, a widely distributed lineage was recovered, which was observed in all populations with the exception of Chocorí and El Moro.This latter lineage is formed by haplotype 1, and two haplotypes are differentiated from it by one mutational step: one belonging to the Saladillo population and another occurring in both Saladillo and Chapadmalal.Even though signs of local demographic expansion can be observed within each lineage (except for haplotype 7, locality Quequén Salado), the population as a whole shows clear population structure.
In the case of C. paleatus, the MP tree (Fig. 3) showed one clade which includes one group formed by the localities El Moro, Claromecó y Quequén Salado, and the widely haplotype (6) and two Chocorí haplotypes occupying basal positions in this clade.Basal to this group and forming a polytomy are the haplotypes from the locality Sauce Grande.Haplotypes 3 and 4 were independently separate.The BI tree showed similar topology, but in this case the first clade shows a split so that the haplotypes from Sauce Grande form an independent clade, although with low support value.Haplotypes 3 and 4 are also differentiated as an independent clade, again with low support (data not shown).The NJ topology was highly concordant with the one obtained through MP.
The haplotype network for C. paleatus (Fig. 1C) shows marked population structure, with haplotypes that are unique to some populations such as Claromecó, Quequén Salado and Sauce Grande.According to the divergence observed in terms of mutational steps, three haplogroups or lineages may be differentiated.One group is formed by a widely distributed haplotype (H 6), related with the populations from Chocorí, Sauce Grande

Historical demography
The global mismatch distribution for J. multidentata showed an erratic curve (Fig. 4a).This result could sug-gest constant population size, or alternatively a population subdivision for a long period.The goodness of fi t tests yielded high values, although not signifi cant in the case of Harpending's Raggedness index (rg) (Table 1).Such values confi rm that the data do not fi t a sudden demographic expansion model.The neutrality tests were positive and not statistically signifi cant, supporting the results of the mismatch distribution analysis (Table 1).
In the case of C. paleatus the mismatch distribution analysis for the global populations also showed an erratic curve that does not fi t the expected sudden expansion  model (Fig. 4b).This result could suggest constant population size, or a population that has been subdivided for a long period.The neutrality tests were negative in both cases, although only statistically significant for Tajima's D (Table 2).

Population structure
The genetic differentiation between J. multidentata populations, showed the highest and significant F st values in the El Moro and Chocorí populations when compared to the populations from Quequén Grande, Cortaderas and Claromecó.Likewise, high F st values were observed in the comparisons with the Sauce Chico and Saladillo populations (Table 3).In contrast, the lowest values were found in the locality Quequén Salado vs. Chapadmalal and Saladillo (Table 3).The F st fixation index for the pooled populations was 0.80 (p < 0.001), indicating high degree of population structure.
On the other hand, the AMOVA indicated that most of the variance is explained among groups, while the rest of it is distributed among populations within groups, and lastly, within populations (Table 4).
The Mantel test yielded a correlation coefficient r= 0.28 (P = 0.076).Although some localities were shown to be greatly divergent, this result from the Mantel test suggests that equilibrium between flow and gene drift has not yet been achieved (Fig. 5a).
Among the C. paleatus populations, the highest significant F st values were found for the populations Sauce Grande vs. Quequén Grande, Claromecó and Cortaderas, as well as for Claromecó vs. El Durazno, Quequén Grande and Cortaderas.The lowest values corresponded to El Durazno vs. Quequén Grande and Cortaderas, as well as Cortaderas vs. Quequén Grande.Values for the remaining comparisons were medium to low (Table 5).
On the other hand, the AMOVA (Table 4) indicated that most of the variance is accounted for within populations.This result shows the marked differentiation within the Corydoras paleatus populations, given by highly divergent haplotypes within a single population (e.g.Claromecó and Sauce Grande).The remaining variance was distributed among groups and lastly, among populations within groups (Table 4).
The Mantel test (Fig. 5b) resulted in a correlation coefficient r = 0.43 (P = 0.0038) which indicates low fit to a pattern of isolation by distance.

Discussion
The analysis of genetic diversity of C. paleatus and J. mul tidentata yielded moderate values, both for haplotype and nucleotide diversity.According to avise (2000) and grant & Bowen (1998), high values for both estimators suggest a demographically stable population with large effective population size maintained through time.On the contrary, low values for both estimators could suggest a severe population bottleneck or possibly selective sweeps, ultimately resulting in a strong reduction of population size.The moderate values found for both species could suggest that the respective populations are on their way to demographic equilibrium.The patterns shown by the historical demography analyses support a scenario of demographic stability.
The results of the analysis of genealogic relationships were similar for C. paleatus and J. multidentata.Both taxa comprised, on the one hand, diverging haplogroups that were in some cases confined to small populations, and on the other, a widely distributed haplotype, with the presence of both lineages in the same locality in some

A B
cases, suggesting the occurrence of secondary contact between the different lineages.
The lack of fit between mutation-genetic drift equilibrium in J. multidentata could be explained by strong effect of genetic drift.According to HutcHinson & templeton (1999), the expected patterns in regions that have not reached equilibrium may be affected both by the time elapsed since colonization of the region, and by the degree of dispersal within said region.If environmental conditions change so that the population that occupies the colonized region becomes fragmented into small isolated populations, genetic drift will have relatively greater influence than gene flow (Case III, HutcHinson & templeton 1999).
In contrast, for C. paleatus there was low although significant fit between genetic and geographic distances (Case IV, HutcHinson & templeton 1999).Given these results, it is probable that the most diverging lineages fit an isolation-by-distance pattern, but the presence of a widely distributed haplotype hinders a more defined manifestation of such pattern.
The integration of the information available so far indicates a similar phylogeographic pattern for both species: both possess haplogroups that are greatly divergent among them.In both species one of these haplogroups is widely distributed, while the remaining ones are restricted to a few watercourses.
The geological and paleoclimatic changes that have affected the Southern Pampean Area could account for the distribution of genetic variation in the fish species that inhabit the area.For freshwater fishes, a pattern of wide colonization in an area can be explained in terms of basin connectivity.On the other hand, for population structure to occur, enough time is necessary so that the individuals may accumulate differences through mutation and in absence of gene flow, which entails a certain degree of geographic discontinuity (avise, 2000; allendorF & luikart 2007).The existence of the phylogeographic patterns found in this study provide an approximation to a geologic and paleoclimatic scenario with alternating connection and disconnection of the basins in this region.Such phenomena have been characteristic of much of the Quaternary in the Southern Pampean Area.
During the middle Pleistocene, the southeastern Buenos Aires province was characterized by arid conditions with marked eolic influence (aramayo et al. 2002;QuattroccHio et al. 1993QuattroccHio et al. , 2008)).The sea level drop during the last glacial maximum (approx.22.000 years AP, ponce et al. 2011) produced great eastwards expansion of the coastal line.In turn, this expansion favored changes in the distribution of rivers and the integration of the drainage network.By the late Pleistocene (16.000 -12.000AP), during dry periods, the main courses were discontinuous and accessory drainage networks, both permanent and semipermanent, were established in association with shallow lakes, and their extension underwent marked seasonal fluctuations (zavala et al. 2005).These adverse geologic and climatic conditions could have caused, on the one hand, local extinctions, and on the other, the survival of some populations in those courses that remained active during climate fluctuations.Assuming this paleoclimatic scenario, it is possible that the highly genetically differentiated populations of J. multidentata and C. paleatus may have occupied environments that remained active during their evolution.Such populations could have achieved population differentiation in absence of or with restricted gene  et al. 2009).Nevertheless, although they have persisted to the present, it is possible that the populations of these lineages underwent drastic size decreases due to habitat fragmentation brought about by the effects of the abovementioned paleoclimatic conditions on the drainage network.This could account for the current population structure.Furthermore, some environments could represent environmental refuges, such as the Quequén Salado river where both species present lineages that are quite divergent with respect to the others.In addition, both species presented exclusive or unique haplotypes in Sauce Grande river.The star-shaped topology observed in the haplotype networks for this locality suggests demographic expansion.Similarly, both species showed the highest values of genetic distance measured as pairwise F st .It is evident that the characteristics of Sauce Grande river have resulted in the population differentiation of these species.More benign climatic conditions have been suggested to occur around 9,000 years BP, during the Holocene sequence (QuattroccHio & Borromei 1998;zavala & QuattroccHio 2001;aramayo et al. 2002).Due to the higher temperature and consequent increase of precipi-tations, medium-sized interconnected lakes were formed in the area (aramayo et al. 2002).Given this scenario, it is possible that some haplotypes, such as the widely distributed one, entered the system after the LGM, due to the favorable connections of the drainage network.This situation may explain the presence of at least one widely distributed haplotype in both species.
Thus, we find lineages that have remained in the area under adverse climatic conditions, in potential environmental refuges, together with lineages that may have colonized the area later due to the more benign climate conditions, as the result of the formation of the modern hydrographic system (Bruno et al. 2013), showing evidence of secondary contact in some localities.
The geologic and paleoclimatic history of the region has largely influenced the distribution of genetic variation in the studied species, promoting population differentiation in some cases.Some water courses, such as the Sauce Grande river, host a considerable genetic diversity, which should be taken into account when formulating conservation plans in the area.Future studies of the ichthyofauna that occurs in basins outside the Southern Pampean area would contribute to enhance the knowledge of factors at a regional scale, that could have shaped the distribution and diversity of the fish fauna in the southern limits of the Brazilian subregion.

Fig. 3 .
Fig. 3. Phylogenetic strict consensus tree obtained by Maximun Parsimony based on mitochondrial DNA control region of Corydoras paleatus.Numbers on the nodes represents the Bootstrap values.The length branches are proportional to mutations per site.H1-H14: Haplotypes.Each square corresponds to one individual.

Fig. 2 .
Fig. 2. Phylogenetic strict consensus tree obtained by Maximun Parsimony based on mitochondrial DNA control region of Jenynsia multidentata.Numbers on the nodes represents the Bootstrap values.The length branches are proportional to mutations per site.H1 -H12: Haplotypes.Each square corresponds to one individual.

Fig. 4 .
Fig. 4. Mismatch distribution of mitochondrial DNA control region for a: Jenynsia multidentata and b Corydoras paleatus.Thick line: observed distribution.Fine line: expected distribution under a sudden population expansion model.
and Claromecó, as well as with those from Claromecó, Quequén Salado and El Moro.On the other hand, two haplotypes are widely divergent and connected to the network by intermediate haplotypes, which occur in the populations Quequén Salado and Claromecó.

Table 4 .
AMOVA for Corydoras paleatus and Jenynsia multidentata considering three hierarchical population groupings: Eastern or Tandilia group, Central or Bonaerean Intersierran Plain group, and Western or Ventania group.
(Bruno et al. 2013) network that was partially disintegrated due to the deactivation of given water courses(Bruno et al. 2013).In fact, both genera occur in the middle Pleistocene fossil record of the study area (descHamps 2005; Bogan