Phylogeographic pattern of Jenynsia multidentata ( Cyprinodontiformes : Anablepidae ) in the southern boundary of the Brazilian Subregion , Argentina

The Atlantic drainage basins are located in the southern Pampean area, in the southernmost part of the Brazilian subregion. Tectonic and paleoclimatic phenomena, occurring during the Quaternary, have currently isolated these basins from the present hydrographic system. Their ichthyofaunal composition is similar to that of hydric systems located further northward. Jenynsia multidentata is a species with wide distribution in the Atlantic drainage basins, as well as in those Northern basins. Here we investigated the phylogeographic pattern of Jenynsia multidentata and analyzed its relationships with the paleoclimatic and geologic history of the region. The analysis of the population’s assemblage showed moderate genetic diversity, demographic equilibrium and marked genetic structure in the populations that occupy the extremes of the distributional range. The populations in the center of the range did not present genetic diversity, sharing a single haplotype. These results allow inferring that the presence of Jenynsia multidentata in the area results from historical demographic processes that are independent but complementary. In turn, these processes could arise from paleoclimatic changes occurred during the Quaternary.


Introduction
The Neotropical freshwater ichthyofauna is the most diverse in the world (Reis, 2003).The actual number of species could reach 8,000 according to Schaefer (1998), Vari & Malabarba (1998) and Reis et al. (2003).The reasons for such a marked diversity are likely to be both historical and ecological, a result of millions of years of evolution from the breakup of Gondwana to the present (Ribeiro, 2006).Much of the diversification of Neotropical freshwater fishes took place in the dynamically changing rivers and watersheds of South America during the Late Cretaceous and Cenozoic.These continental waters were at once agents and the products of landscape evolution (Lundberg et al., 1998).
The Neotropical region is divided into the Austral and Brazilian subregions (Ringuelet, 1975;Arratia et al., 1983).The southern boundary of the Brazilian subregion is represented by the rivers and streams of the Atlantic coastal drainage in the southern Pampean area, Argentina (Casciotta et al., 1999) (Fig. 1).The courses in this drainage basin flow southwards following a parallel pattern and are isolated from each other.Moreover, all these water courses, mostly relatively small, are also separated from the ones located to the north (Almirón et al., 1997;Casciotta et al., 1999).
Most of the streams in the southern Pampean area run south from Tandilia, Ventania and the Llanura Interserrana Bonaerense (Casciotta et al., 1999).When the whole positive area of Tandilia, the Llanura Interserrana Bonaerense and Ventania are considered, the drainage pattern appears to be radial.The hydrogeographic pattern of the region is the result of a combination of tectonics, climate and sea level changes since the Miocene (Casciotta et al., 1999).Due to eustatic and/ or neotectonic changes occurred during the late Pleistocene and Holocene, the current river beds could have been tributaries of different water courses in the past (Frenguelli, 1956;Tonni & Cione, 1997).In fact, the Sauce Grande River has been proposed to have been part of the Colorado River basin during the Quaternary (Casciotta et al., 1999;Ponce et al., 2011).The southern Pampean area has been characterized by alternating arid and cold/humid and warm conditions during the Quaternary and even since the Miocene (Aguirre et al., 1999;Aramayo et al., 2002;Quattrocchio et al., 2008).Paleoclimatic reconstructions suggest an arid-semiarid habitat for the southern Pampean region during much of the Quaternary (Quattrocchio et al., 1988(Quattrocchio et al., , 1993;;Tonni et al., 1999;Nabel et al., 2000) (Fig. 1).
The fish fauna present in the Atlantic coastal drainage in southern Buenos Aires province, Argentina is very poor in terms of species richness.The number of recorded species is no higher than 14, and it is particularly interesting that the ichthyofauna of these basins is similar to that of hydric systems located further north, even though there is at present no connection between these basins.
According to Cione & Barla (1997) and Casciotta et al. (1999) the presence of Brazilian species in these basins may be explained by the following hypotheses: 1.These fishes were present in the area before the last glacial period and resisted unfavorable climates; 2. These fishes entered the area during periods of marine regression, through basins that are currently submerged in the ocean platform.These basins may have been connected with others situated farther north.
The fishes of the genus Jenynsia are small viviparous cyprinodontiforms.This genus comprises thirteen species distributed in highland drainages of southeastern Brazil, lowlands in the La Plata basin, coastal Atlantic drainages of southern and southeastern Brazil, Uruguay, western and northeastern Argentina, and southeastern Bolivia (Ghedotti, 1998(Ghedotti, , 2003;;Lucinda et al., 2006).
Jenynsia multidentata is a species with widespread distribution in Argentina (Liotta, 2005).In addition, its fossil record in the study area dates back to the Middle Pleistocene (Bogan et al., 2009).Because of these characteristics, it is a good model for phylogeographic and historical biogeographic studies, as it can provide information about the arrangement of the basins in this region before the tectonic and paleoclimatic processes determined their isolation from northern basins.
In this work we investigate the phylogeographic pattern and historical demography of Jenynsia multidentata.In addition, we analyze the relationship between this pattern and the geomorphologic and paleoclimatic history of the region, in order to test the hypotheses proposed about the ichthyofaunal population presence in the research area.

Material and Methods
Sample collections.Tissue samples were obtained from 99 Jenynsia multidentata specimens from ten localities that span the entire study area (Fig. 1).Each locality was considered as a population.Once the tissue was removed for DNA extraction, the specimens were fixed in 10% formaldehyde and kept in 70% ethyl alcohol.The tissue samples were deposited in the tissue collection of the Centro Regional de Estudios Genómicos, Universidad Nacional de La Plata, Argentina (CREG-EM tissue bank).The well-preserved specimens were deposited in the ichthyological collection of Museo Argentino de Ciencias Naturales Bernardino Rivadavia, Argentina (MACN-Ict).Voucher identification number is indicated in Table 1.
DNA extraction, amplification and sequencing.DNA was extracted following the protocol of Aljanabi et al. (2007), which consists of protein precipitation in ClNa and subsequent DNA precipitation with isopropanol.
The entire control region of mitochondrial DNA was amplified using polymerase chain reaction (PCR), using the primers K 5' AGCTCAGCGCCAGAGCGCCGGTCTTGTAAA 3' and G 5' CGTCGGATCCCATCTTCAGTGTTATGCTT 3' (Lee et al., 1995).The amplification reaction was done in a total volume of 25 µl with a final concentration of 0.25 µ/l Taq DNA polymerase, 1.5 µl of Cl 2 Mg 3mM, 0.4 µl of dinucleotids 50 mM, 0.25 µl 10mM of each primer and 1µl of DNA as template.The reaction was made under the following conditions: initial DNA denaturalization at 95ºC for four minutes, followed by 35 cycles of denaturalization at 95ºC for 30 seconds, annealing at 57ºC for 30 seconds and extension at 72ºC during 45 seconds, followed by a final extension period of five minutes.Negative controls were performed for all samples to verify absence of contamination.The PCR products were purified following the alcohol purification protocol (Sambrook, et al., 1989).Amplicons were sequenced in a capillary sequencer model ABI 3100 (Macrogen Inc., korea).
Chromatograms were edited using the software Proseq (Filatov, 2002) and aligned using Clustal W (Thompson et al., 1994) with default parameters.Genetic variation.The software ARLEQUIN 3.5 (Excoffier et al., 2010) andDNAsp 5.10 (Librado et al., 2009) were used to estimate genetic diversity on the basis of two different parameters: haplotype diversity (h), defined as the probability that two randomly chosen haplotypes are different in the sample (Nei, 1987), and nucleotide diversity (π) defined as the probability that two randomly chosen nucleotide sites are different (Nei, 1987).For the first approach, a Maximum Parsimony (MP) analysis was performed as implemented in PAUP * 4.0b10 (Swofford, 2002).A heuristic search with 1,000 random stepwise additions and tree bisection and reconnection (TBR) branch swapping was made.On the other hand, a Maximum likelihood (ML) was conducted in RAxML (Stamatakis et al., 2006).Node support was assessed by 1,000 bootstrap replicates using the fast bootstrapping algorithm implemented in the RAxML web-servers (Stamatakis et al., 2008).

Genealogic relationships among haplotypes. To evaluate the phylogenetic relationships among mitochondrial DNA
For the second approach, an analysis based on genetic distances was made using the Neighbor-Joining method (NJ) implemented in PAUP * 4.0b10.Because of, TPM3uf+G is not a model available in PAUP, we select the GTR model, for calculating genetic distances.
Both MP and NJ methods, a 1,000 bootstrap replicates (Felsenstein, 1985) were conducted to assessed node support values by calculating the 70% Majority Rule Consensus Tree.
Sequences from the Cyprinodotiformes Xiphophorus montezumae and Poecilia latipinna (GenBank acc.number: DQ445680.1 and HM567257.1 respectively) were used as outgroup to root the trees.
Because traditional methods of phylogenetic analysis have not been designed to be used at intraspecific level (Posada & Crandall, 2001), evolutionary relationships between haplotype variants were obtained by constructing a haplotype network.This was performed using the software Network 4.5.1 (http://www.fluxus-engineering.com), applying the Median-Joining algorithm (Bandelt et al., 1999).(Fu, 1997) and Tajima's D (Tajima, 1989), and the R 2 statistic (Ramos-Onzins & Rozas, 2002).Fu's Fs and Tajima's D tests evaluate departure from neutrality (mutation-genetic drift equilibrium) as expected under a demographic expansion model.Significant negative values are expected in populations that have undergone recent expansion.The R 2 statistic is based on the difference between the number of singleton mutations and the mean number of nucleotide differences between sequences in a sample.The populations that have undergone large expansion are expected to exhibit low R 2 values.The significance of these tests was assessed using 10,000 replicates of coalescent simulations in DNAsp (Librado et al., 2009).Additionally, a mismatch distribution analysis was performed, which are the distribution of pairwise differences among haplotypes (Slatkin & Hudson, 1991;Rogers & Harpending, 1992).These distributions are ragged and erratic in samples from populations that have been stationary for a long time, whereas they are smooth and usually unimodal in populations that have been growing for a long time or that have experienced a single burst of population growth in the past (Harpending, 1994).To test the validity of the population expansion model of the observed mismatch distribution with respect to the expected one, a goodness of fit test was made using the sum of square deviations (SSD) and a Harpending's Raggedness index (r) (Harpendig, 1994).These analyses were performed using the software ARLEQUIN 3.5 (Excoffier & Lischer, 2010).

Historical demography. Historical demography was assessed by means of neutrality tests including Fu's Fs
Population structure.To estimate the genetic distance between populations the pairwise Φ ST 's among localities was calculated (Reynolds et al., 1983;Slatkin, 1995).
A description of the genetic structure was obtained by means of Analysis of Molecular Variance (AMOVA) (Excoffier et al., 1992) using the software ARLEQUIN 3.5 (Excoffier et al., 2010).For this analysis, three different settings were carried out.The first analysis included three population groups, established on the basis of the geographic location of the watercourse headwaters.These groups are the East or Tandilia (Chapadmalal, Chocorí, and El Moro streams, and Quequén Grande River), the Central or Llanura Interserrana Bonaerense (Cortaderas and Claromecó streams), and the West or Ventania (Quequén Salado River and Sauce Grande, Saladillo, and Sauce Chico streams).The second analysis, included the same three populations groups, however here were considered the lower sections of the watercourses.These groups are East or Tandilia group (Chapadmalal, Chocorí, and El Moro streams), the Central or Llanura Interserrana Bonaerense (Quequén Grande and Quequén Salado Rivers, Cortaderas and Claromecó streams), and the West or Ventania group (Sauce Grande River, Saladillo, and Sauce Chico streams).The third analysis included all populations as a unique group.Statistical significance was assessed using 10,000 permutations.A Mantel test (Mantel, 1967) was performed to test for Isolation by Distance pattern (Wright, 1943;Slatkin, 1993) by testing for correlation between geographic and genetic distances using 1,000 permutations.

Genetic variation.
A 771 base pair (bp) fragment from the entire control region was analyzed for 99 individuals distributed into ten populations (Genbank acc number: KC 485466-485474).It is worth to note that the obtained fragment was initially 880-900 bp long; due to the presence of adenine homopolymers, whose presence results in sequence reading errors, these sites were discarded at alignment.Twenty-two polymorphic sites were identified, of which 16 were transitions and six transversions.Likewise, five indels were detected, to define a total of nine haplotypes.These haplotypes were defined in the same manner whether or not indels were taken into account, and therefore are considered as polymorphic sites for the different analyses.The most widely distributed haplotype was found in eight populations, and it was the only one for six of them.The remaining populations presented between two and four haplotypes (Table 2).
Genealogic relationships among haplotypes.Both the evolutionary method (MP and ML) and the one based on genetic distances (NJ) showed trees with similar topologies (Figs.2-4).Three main clades (MP, ML) or clusters (NJ) were recovered with high bootstrap support values.The first clade/ cluster comprises three of the nine haplotypes found (H1, H3, H9) including the one with widest distribution (H1).The second clade/cluster includes haplotypes that occur in the East or Tandilia group (H2, H4, H5) and the third clade/cluster include haplotypes that occur in the Sauce Grande River (H6, H7, H8) which belongs to the West or Ventania group.
The haplotype network shows the same terminal groups (Fig. 5) recovered in the phylogenetic and distance trees, differentiated by several mutational steps, with high geographic concordance.The ancestral haplotype in the network is occupied by three hypothetical, closely related ancestral haplotypes.The occurrence of intermediate or hypothetical haplotypes may be explained by lack of sampling, or these may be haplotypes that have disappeared before the present.
Historical demography.The mismatch distribution for the set of samples showed an erratic curve (Fig. 6).Although Harpending's Raggedness index (r) was high though not significant (0.24 P = 0.97), the high and significant SSD value (0.28 P = 0.0001) demonstrates that the data do not fit a model of population expansion.The R 2 index was high but not significant (0.11 P = 0.81).Regarding the neutrality tests, both Tajima's D and Fu's Fs were positive and non-significant (0.71 P = 0.81 and 5.23 P = 0.94 respectively).The values for each population are provided in Table 3.   Population structure.The genetic differentiation between population pairs, measured as pairwise Φ ST 's, showed that the highest significant values are those between the East basins (Tandilia group), with respect to the remaining populations.It is worth to note that in the Tandilia group the population from El Moro stream, showed the highest differentiation with values between 0.37 and 1.The values for each population are presented in figure 7. The value of the fixation index Φ ST for all populations was 0.56 (P<0.001).These values indicate a high degree of population structure.
The AMOVA considering the headwater courses grouping indicated that most of the variance was explained between populations within groups (50.11%), while 43.14% of the variance was distributed within populations, and 6.75 % between groups.The fixation indices Φ ST, Φ CT and Φ SC were 0.56 (P<0.001),0.067 (P = 0.22) and 0.53 (P<0.001)respectively.The high and significant value of Φ SC shown differentiation at population level within groups.The variance among groups defined by the Φ CT was not significant.
The analysis considering the lower sections of watercourses, shown that most of the variance was explained among groups (45.52 %), while the 37.82 % among populations between groups, and the 16.66 % within populations.The fixation indices Φ ST, Φ CT and Φ SC were 0.83 (P<0.001),0.45 (P = 0.14) and 0.69 (P<0.001)respectively.By considering this grouping, the fixation indices support evidence of strong genetic structure among groups.
The AMOVA considering a unique group showed an 80.92 % of the variance among populations and a 19.08 % into populations.The Φ ST was 0.80 (P< 0.0001).
The Mantel test yielded a coefficient of correlation r= 0.037 (P= 0.86) that evidences that the populations do not fit a model of Isolation by Distance (Fig. 8).

Discussion
Although the overall analysis of the results of this work supports the affirmation that the populations of Jenynsia  Grande and Saladillo in the West), whereas the remaining populations are represented by a single haplotype and thus it is not possible to estimate genetic diversity values for them.
Unquestionably, the settlement and persistence of the species in this area are part of a complex scenario that cannot be explained only by a single demographic event.Such demographic events were associated with major climatic and geological changes that has been affected the study area during the Late Pleistocene and Holocene.
The decrease of sea level during the last glacial maximum (22,000 calibrated years BP, Ponce et al., 2011) produced a great eastward expansion of the coastal line, prompting changes in the distribution of rivers and the integration of the drainage network.In this sense, the Sauce Grande River could have been part of the Colorado River basin during the Pleistocene (Ponce et al., 2011).These   multidentata are under demographic equilibrium and show strong population differentiation, these results should be interpreted with caution.It is evident that not all the populations show the same pattern.Genetic diversity values obtained in this work appear to be low to moderate for the global samples.Similar pattern was observed in others cyprinodontiforms fishes, such as the family Poeciliidae (Johnson, 2001;Gutiérrez-Rodríguez et al., 2007).In addition, these values were recovered from four populations that occur in the extremes of the geographical distribution (Chapadmalal and Chocorí in the East, and, Sauce authors found fossil evidences of Percichthys sp. for the middle Pleistocene in the southern Pampean area and note that the entrance path for its settlement was through these paleoconnections between different basins while the sea level was below the current costal line, probably during the Pleistocene.This suggests that Jenynsia multidentata became settled in the area at a similar time, achieving wide distribution in the region thanks to the at least partial integration of the drainage network.Indeed, fossil evidence in the area demonstrates the presence of Jenynsia sp. at the locality Centinela del Mar during the middle-late Pleistocene (Bogan et al., 2009); likewise, scales from Cyprinodontiformes have been reported for the Holocene from Napostá stream (Quattrocchio et al., 1998).
The stratigraphic analysis of Quaternary alluvial terraces associated with the main water courses that drain the mountainous area indicate that their activity has been highly discontinuous throughout their geologic history.Thus, although these rivers are the major agent for valley excavation, their fill sediments indicate that they would have remained dry during long periods of their evolution (Zavala et al., 2005).These adverse climatic and geologic conditions could have led to local extinctions, and the survival of some populations in those water courses that remained active during those climate fluctuations.
Given such a scenario, it is possible that the populations of Jenynsia multidentata that are currently structured could have inhabited environments that remained active throughout their evolution.These populations could have achieved population differentiation without or with only restricted gene flow, in a partially disintegrated drainage network resulting from the deactivation of some water beds.
These findings are consistent with idea that, the fishes were present in the area before last glacial period and resisted unfavorable climates (hypothesis 1, see introduction).
For the Pleistocene, the area has been described as having wide flood plains in an arid-semiarid environment (Zavala & Navarro, 1993;Quattrocchio et al., 1993Quattrocchio et al., , 1998;;Aramayo et al., 2002).Around 9,000 calibrated years BP, with the beginning of the Holocene and improved climate conditions, these wide flood plains were occupied by interconnected lakes (Aramayo et al., 2002) and anastomosing river systems (Quattrocchio 1998).These interconnected systems could have acted as the dispersal route by means of which a few populations, colonizing most of the area.This situation is reflected in the lack of fit to an Isolation-by-Distance model, expected in populations that have recently colonized the area that they occupy (Slatkin, 1993).Thus, we find populations that possibly remained in the area under adverse climate conditions in environmental refuges, such as the oldest water courses.
Later on, and due to more favorable climate conditions, they may have undergone demographic expansion with subsequent colonization of the area, as a consequence of the formation of a modern hydrographic system.Evidence for this fact is the presence of a single haplotype widely distributed in the area.This situation could be supporting the second hypothesis about the species occurrence in the area.In fact, the age of some portions of the Sauce Grande River has been dated as Late Pleistocene-Holocene, while other water courses such as the Sauce Chico stream are 4,400 years old (Prieto, 1996;Quattrocchio et al., 2008).
To sum up, the presence of Jenynsia multidentata in the area is the result of independent but complementary demographic processes together with the geologic evolution of the area.In this sense we found support for hypotheses 1and 2 that could account for the presence of the species in the southern Pampean area

Table 2 .
Genetic diversity of Jenynsia multidentata in the Southern Pampean area.N = number of individuals.S = number of polymorphic sites.NH = number of haplotypes.h±SD = haplotype diversity ± standard deviation.π ± SD = nucleotide diversity ± standard deviation.-Not estimated.

Fig. 2 .
Fig. 2. Phylogenetic strict consensus tree obtained by Maximun Parsimony based on mitochondrial DNA control region of Jenynsia multidentata.Numbers under the nodes represents the Bootstrap values.The length branches are proportional to mutations per site.H1-H9: Haplotypes.Black vertical lines next to sample localities are proportional to the individual numbers in each locality.Shading patterns indicates localities.Black: Chapadmalal.Black doted: Chocorí.Cross: El Moro.Dark grey: Quequén Grande.Horizontal lines: Cortaderas.Light grey: Claromecó.Grey: Quequén Salado.Diagonal lines: Sauce Grande.Vertical lines: Saladillo.Grey doted: Sauce Chico.

Fig. 3 .
Fig. 3. Phylogenetic consensus tree obtained by Maximum Likelihood based on DNA control region of Jenynsia multidentata.Numbers under the nodes represents the Bootstrap values.The length branches are proportional to mutations per site.H1-H9: Haplotypes.Black vertical lines next to sample localities are proportional to the individual numbers in each locality.Shading patterns indicates localities.Black: Chapadmalal.Black doted: Chocorí.Cross: El Moro.Dark grey: Quequén Grande.Horizontal lines: Cortaderas.Light grey: Claromecó.Grey: Quequén Salado.Diagonal lines: Sauce Grande.Vertical lines: Saladillo.Grey doted: Sauce Chico.

Fig. 4 .
Fig. 4. Distance tree obtained by Neighbor-Joining method based on mitochondrial DNA control region of Jenynsia multidentata.Numbers under the nodes represents the Bootstrap values.The length branches are proportional to mutations per site.H1-H9: Haplotypes.Black vertical lines next to sample localities are proportional to the individual numbers in each locality.Shading patterns indicates localities.Black: Chapadmalal.Black doted: Chocorí.Cross: El Moro.Dark grey: Quequén Grande.Horizontal lines: Cortaderas.Light grey: Claromecó.Grey: Quequén Salado.Diagonal lines: Sauce Grande.Vertical lines: Saladillo.Grey doted: Sauce Chico.

Fig. 6 .
Fig. 6.Mismatch distribution of mitochondrial DNA control region for Jenynsia multidentata.Thick line: observed distribution.Fine line: expected distribution under a sudden population expansion model.The dashed line represents confidence interval at 95%.

Fig. 7 .
Fig. 7. Population relationships described by Φ ST computed between pairs of populations inferred from mitochondrial DNA control region of Jenynsia multidentata.