Conflicting Evolutionary Patterns Due to Mitochondrial Introgression and Multilocus Phylogeography of the Patagonian Freshwater Crab Aegla neuquensis

Background Multiple loci and population genetic methods were employed to study the phylogeographic history of the Patagonian freshwater crab Aegla neuquensis (Aeglidae: Decopoda). This taxon occurs in two large river systems in the Patagonian Steppe, from the foothills of the Andes Mountains east to the Atlantic Ocean. Methodology/Principal Findings A nuclear phylogeny and multilocus nested clade phylogeographic analysis detected a fragmentation event between the Negro and Chico-Chubut river systems. This event occurred approximately 137 thousand years ago. An isolation-with-migration analysis and maximum-likelihood estimates of gene flow showed asymmetrical exchange of genetic material between these two river systems exclusively in their headwaters. We used information theory to determine the best-fit demographic history between these two river systems under an isolation-with-migration model. The best-fit model suggests that the Negro and the ancestral populations have the same effective population sizes; whereas the Chico-Chubut population is smaller and shows that gene flow from the Chico-Chubut into the Negro is four times higher than in the reverse direction. Much of the Chico-Chubut system appears to have only been recently colonized while the Negro populations appear to have been in place for most of the evolutionary history of this taxon. Conclusions/Significance Due to mitochondrial introgression, three nuclear loci provided different phylogeographic resolution than the three mitochondrial genes for an ancient fragmentation event observed in the nuclear phylogeny. However, the mitochondrial locus provided greater resolution on more recent evolutionary events. Our study, therefore, demonstrates the need to include both nuclear and mitochondrial loci for a more complete understanding of evolutionary histories and associated phylogeographic events. Our results suggest that gene flow between these systems, before and after fragmentation was through periodic paleolakes that formed in the headwaters region. Fragmentation between the Negro and Chico-Chubut systems was driven by the disappearance of these paleolakes during the Patagonian Glaciation.


Introduction
Phylogeographic inference is a powerful tool for understanding an organism's evolutionary history [1,2]. Mitochondrial markers have been the primary data for most of the discipline's history [3]. Mitochondrial loci dominance in the field of phylogeography has been recently challenged because this marker provides only a portion of and possibly misleading picture of the history of the organism under study [4]. Despite these concerns mitochondrial markers are and will continue to be the data of choice for discovering phylogeographic patterns [5]. However, independent markers, mostly in the form of nuclear loci are required to confirm mitochondrial patterns and provide more robust estimates of processes, such as gene flow and demographic changes [6,7]. In this study we employ both mitochondrial and nuclear markers to infer the evolutionary history of a geographically widespread Patagonian organism, the freshwater crab Aegla neuqeunsis (Aeglidae: Decopoda).
Aegla neuquensis is found in two major river systems (Negro and Chico-Chubut) in the Argentinian Steppe (Fig1: [8]). The Negro system is comprised of three rivers, Neuquen, Limay and Negro where as the Chico and Chubut rivers form the Chico-Chubut system. These two river systems are geographically distant from one another along most of their lengths (,500 km). Only in the headwater regions of the Limay and Chubut rivers, in the foothills of the Andes, do these two systems come into close proximity (,2 km) to one another. The headwater region is the only part of the distribution of Aegla neuquensis that was directly impacted by glacial ice during glacial cycles. A paleolake called Lake Cari Laufquen may have formed in the region (,41uS) during the last several glacial cycles (i.e., ,23-25 ka and 0.7 Ma) [9]. Lake Cari Laufquen would have provided suitable habitat for gene flow between systems [9]. Gene flow between these two systems has been proposed in several fish groups [10][11][12][13][14] as well as other decapods [15]. Beyond the limits of glacial ice changes in stream dynamics could have affected the evolutionary history of this taxon. For example, gene flow could also have occurred via the deltoic mosaic that occurred on the continental shelf exposed during previous glacial periods [16,17].
Five major glacial events occurred in the region over the last four million years [18,19,20,21]: A widespread glaciation of the middle-Pliocene (,3.5 Ma); The largest Patagonian glacial advance (1.1 Ma); The coldest Pleistocene glaciation (,0.7 Ma); The last southern Patagonian glaciation that reached the Atlantic coast (180 ka); and lastly the last glacial maximum (LGM), which occurred between 23-25 ka years ago. Each of these periods may have played a role in shaping the phylogeographic history of Aegla neuquensis in general and formed paleolakes near the headwaters specifically. However we predict little impact of glaciation on this taxon because most of the distribution of Aegla neuquensis occurs beyond the limits of glacial ice. In these ice-free regions, we expect to recover stable demographic histories and genetic structuring between rivers and populations. In contrast we predict glacial cycles will have significantly impacted the headwaters region of this taxon. In this region we expect to recover the genetic signature of historical fragmentation when no paleolakes were present and when there were paleolakes in the region gene flow between the Negro and Chico-Chubut systems. Because the signature of previous glacial cycles is likely to be overwritten by later events, we expect these events to date to recent cycles [20,22].
The natural history of Aegla crabs makes them excellent candidates for phylogeographyic study. Aegla species in general occur in freshwater habitats throughout the Patagonian region [21][22][23][24][25][26][27]. These crabs are small (,60 mm) and spawn from February to March. Females are fecund, laying as many as 1500 eggs. Once laid, eggs attach to the females pleopod where they remain until they hatch approximately six to eight months later [24]. Even after hatching, young crabs remain with the female for several days. The absence of a free-floating larval stage and the long-period of attachment to the female suggests that dispersal is restricted [24]. Limited dispersal and widespread distributions make Aegla, including A. neuquensis ideal subjects of phylogeographic study [24] and good indicators on how abiotic events have shaped the freshwater systems where it occurs.
We employ three nuclear and one mitochondrial (three genes) locus and a variety of population genetic approaches to recover phylogeographic patterns and processes in Aegla neuquensis. Comparing patterns and process across multiple markers will allow us to make rigorous inferences [6,28,29] and will provide the most detailed examination of the evolutionary history of a Patagonian freshwater taxon to date and add to our growing understanding of the evolutionary history of the region [16,17,[30][31][32][33].

Sampling and Sequence Data
Mitochondrial and nuclear haplotype diversity was high. We recovered 146 unique haplotypes from 295 individuals in our combined mtDNA data set (Table 1and Fig . 2). Most of the haplotypes (109 of 146, ,75%) were singletons and only nine occurred in more than one location. The complete mtDNA alignment of fragments of the 16 S, COI, and COII genes is 1699 bp long. The COI (659-bp), COII (568-bp) and 16 S (472bp) fragments delineated 82, 79 and 39 haplotypes, with 97, 91 and 41 variable sites, respectively, with an average sequence divergence across genes of 2.7% (0-8.5%). In a 1395-bp concatenated alignment of nDNA fragments, 67 haplotypes were identified from 103 crabs (Fig. 3), with an average sequence divergence across loci of 0.6% (0-1.4%). As in the mtDNA, most of the haplotypes (51 of 67, ,76%) were singletons and only six occur in more than one location. The length of EF1a exon, EF1a intron and ANT intron were 640, 370 and 385 bp respectively defining 13, 21 and 22 haplotypes (phased alleles).

Phylogenetic Relationships
Unexpectedly, the phylogeographic resolution of mitochondrial and concatenated nuclear phylogenies differed significantly (Figs. 2, 3). The mitochondrial phylogeny recovered a non-monophyletic Limay river with Chico-Chubut haplotypes embedded within it. In contrast, the nuclear phylogeny recovered a monophyletic Chico-Chubut system suggesting that this system has been genetically isolated from the Negro system. Similar mtDNA phylogenies were recovered using Maximum likelihood and Bayesian methods (Fig. 2). A McDonald-Kreitman test for selection was not significant using both Fisher's exact test (p = 0.64) and G test (G = 0.221, p = 0.63). Thus, there is no evidence that the mitochondrial data are experiencing natural selection that might obscure phylogeographic signal. The Vaca and Telsen populations clustered with the outgroup taxa and are more genetically distant from the remaining ingroup samples ( Table 2: and see nuclear phylogeny below). These individuals are therefore treated as an unrecognized taxon (new species) and were not included in subsequent analyses. Their distinctive and nonsister phylogenetic position suggests these populations are a distinct species from other aeglids. These populations are currently being described as a new species (Jara in prep). All remaining A. neuquensis samples formed the basal clade. Three main clades are found within this A. neuquensis clade. Two of these clades contain all the Negro and Neuquen individuals. Samples from Covunco, Bocatoma, Chimpay and El Alamito form a clade sister to a clade containing all the Ingeniero Ballester individuals. Samples from Cullin are shared between these two sister clades. The remaining samples from the Limay, Chubut and Chico systems form the remaining clade. Within this clade individuals from the Limay system are basal. Individuals from Chubut and Chico rivers are derived, monophyletic and sister to three haplotypes from the Pichileufu population in the upper reaches of the Limay River. The concatenated nDNA tree recovered three well-supported clades. One clade contained all the individuals from the Limay, Negro, and Neuquen rivers. There is relatively little support for relationships within this clade. As in the mitochondrial tree topology, all Vaca and Telsen are monophyletic and genetically distant from the remaining samples. All the Chico and Chubut haplotypes form the third clade and are sister, but with no support, to the Vaca and Telsen clade.
The mitochondrial TMRCA for all A. neuquensis (minus Vaca and Telsen samples as they are a distinct species) was 235 ka (149-322 ka; Fig. 2). This range overlaps two glacial cycles, the last Patagonian glaciation (,180 ka) and the previous cycle (,270 ka). The mitochondrial TMRCA for the clade containing nearly all the Negro and Neuquensis river individuals was 219 ka (140-300 ka; Fig. 2). The youngest mitochondrial TMRCA is 29 ka (17-41 ka: Chico-Chubut clade, Fig. 2), which is consistent with the LGM (,20 ka). See EBSP results for the TMRCA from the combined data sets.

Networks
Individual haplotype networks from both genomes (Fig. 4,5,6,7,8,9) show relationships that are consistent with those recovered in the phylogenetic analyses (Figs. 2, 3); that is, haplotypes from the Neuquen, Negro and Limay rivers are more closely related to one another than to those found in the Chubut-Chico systems. The ANT intron network is an exception because it has one widespread haplotype (#15) that is shared between these two systems and a Limay river haplotype (#20) that is more closely related to Chico-Chubut than it is to the other Negro system haplotypes ( Fig. 9). Overall, networks group geographically close haplotypes. Nuclear networks further illustrate the relatively low haplotype diversity observed within the Chico-Chubut river system.

Multilocus Nested Clade Phylogeographic Analysis
While NCPA found multiple significant inferences (Table 3), only three were detected across two or more loci (Table 4). Only these three will be considered for further discussion. Fragmentation (either past-fragmentation [PF] or allopatric-fragmentation [AF]) between the Limay and Chubut rivers was inferred in all four loci. As predicted this fragmentation event occurred between   Table 1) and number of individuals are shown after that haplotype. doi:10.1371/journal.pone.0037105.g002 the headwaters of the Limay and Chubut rivers near the base of the Andes (Fig. 1). Restricted gene flow (RGF) was detected in the Chubut (three loci: ANT Intron, EF1 Exon, EF1 Intron) as well as in the Neuquen/Negro River (two loci: mtDNA and ANT Intron) systems.

Demographic Model Selection
The highest ranked demographic model using Akaike Information Criterion (AIC) was ABADE (Tables 5 and 6: four additional models could not be rejected using chi-square with or without Bonferroni correction). This model had unequal rates of gene flow  Table 1) and number of individuals are shown after that haplotype. Color codes indicate the six major rivers. doi:10.1371/journal.pone.0037105.g003 between the Chico-Chubut and Negro systems with gene flow (where M = m/m) from the Chico-Chubut into the Negro system (0.0404) greater than the reverse (0.0118). The estimated ancestral (qa) population (where h = 4 Nem) size was the same as the Negro system (8.3701), while the Chico-Chubut has the smallest population size (1.8135). The time of population divergence between these two river systems was 137 ka. The geometric mean of the mutation rate across loci was 4.008839610 25 /combined loci/generation. The 95% posterior probability distribution (not shown) for time (t) does not include zero, supporting the hypothesis of historical isolation (i.e., fragmentation) with post isolation gene flow. Each of the three IMa runs converged on nearly identical parameter estimates and had relatively high ESS values (.50: values for t tended to be lower). Thus, we conclude that these analyses were well sampled and converging on correct values and that we were justified in combining the three runs for the model selection procedure.

Gene Flow
Asymmetrical gene flow was observed in both mtDNA and nDNA with a greater magnitude of downstream flow detected within the Neuquen and Negro river system (Fig. 10). Similar patterns between mtDNA and nDNA suggests that there is no sex biased gene flow. Both mtDNA and nDNA suggest that the majority of populations in the upper reaches of the Limay River are exchanging genes. In the Chubut and Chico rivers, only mtDNA possessed sufficient variation to provide reliable estimates. In the Chubut system mitochondrial gene flow is occurring both upstream and downstream from the Los Altares population. Conversely, in the Chico River all gene flow is following an

Extended Bayesian Skyline Plots
Contrasting demographic histories between the two major river systems were estimated using EBSP. We cannot reject a constant demographic model for the Negro system ( Fig. 11a inset). However, the 95% HPD includes one demographic change. Inspection of the EBSP suggests that this event could have occurred over the last 65 ka (Fig. 6a). In contrast, we can reject a constant demographic model for the Chico-Chubut system because the mode number of changes is one (Fig. 11b inset). The EBSP shows an explosive increase over the last 2-3 thousand years (Fig. 11b). These plots show that the Negro system's TMRCA is much older than that of the Chico-Chubut, 300 ka and 48 ka, respectively. The effective population size for Chico-Chubut was much lower (near zero) than the Negro's for most of the history of this population. Each of the three EBSP runs conducted for each system converged on similar estimates for our parameters of interest and these parameters had high (.200) ESS values. See the appendix for the BEAST xml files and plots of TMRCA for each locus.

Population Genetic Analysis
Population genetic diversity as indicated by haplotype diversity (Hd) was generally high except those of Bayas, Covunco and Chubut, but nucleotide diversity (h p ) was relatively low, especially in the populations from the Chico-Chubut River (Table 1). These observations indicate a higher genetic diversity and deeper genetic divergence in populations from the Negro system than those from the Chico-Chubut system, results consistent with the phylogenies and networks. Stronger genetic divergence was observed between A. neuquensis populations (F ST = 0.774, P = 0; AMOVA test), and pairwise F ST values were large and significant (P,0.05) among all populations, excluding those between Bocatoma and Chimpay, Bocatoma and Cullin, Bocatoma and Pichileufu, Piedra del Aguila and La Rinconada, and Los Altares and Ameghino (Table 7). Genetic differentiations between the Neuquen-Negro River and Chico-Chubut River populations (F ST = 0.70820.988) were the largest, followed by those between the Neuquen-Negro River and Limay River (0.503-0.986), while those between the Limay River and Chico-Chubut River (0.566-0.927) were relatively small. Total genetic diversity as measured by theta was much greater in Negro system versus that found in the Chico-Chubut system ( Table 8).
Tajima's D tests were significant (P,0.05) and negative for Bayas, Sarmiento (not shown) and Chico-Chubut River (Table 8). Fu's Fs was significant and negative for Rinconada, Los Altares, El Blanco, Mayo, Sarmiento (not shown) and Negro and Chico-Chubut systems (Table 8). Significant and negative values of Tajima's D and Fu's Fs suggest that these populations and rivers have experienced recent demographic growth.

Fragmentation Event
We find considerable support for our primary hypothesis of gene flow and fragmentation between the Negro and Chico-Chubut systems. These two systems are reciprocally monophyletic in our nuclear phylogeny. Evolutionary relationships primarily from genetic networks suggest that this fragmentation event and the between systems gene flow occurred between populations within the headwaters region of these two large river systems. Process inference from NCPA and an isolation-with-migration analysis provide corroboration for our primary hypothesis and additional details on the demographic history of this event. For example, multilocus NCPA analysis inferred this fragmentation event in all four loci. The highest-ranking demographic isolationwith-migration model suggests that the ancestral (pre-fragmentation) population was equal in size to the current size of the Negro population and that both are larger than the Chico-Chubut population. These relative population sizes are consistent with the patterns of haplotype diversity. The isolation-with-migration model suggests that gene flow between systems is asymmetrical, with up to four times as many individuals moving into the Chico-Chubut system from the Negro system as occurred in the opposite direction.
The fragmentation event between the Negro and Chico-Chubut systems occurred roughly 137 ka. Of the five major climatic events outlined by Ruzzante et al. [17,19], only the Patagonian glaciation (,180 ka) is consistent with our estimate. We posit that a paleolake formed in the headwaters region of these two systems during the Patagonian glaciation [16,18]. The existence of a lake in this region is supported by the detection of limited gene flow between populations in the headwaters. Our isolation-with-migration analysis confirms post-isolation gene flow. Flooding of the region during the Patagonian glaciation would have provided opportunities for gene flow between systems. These two systems were then isolated again following the retreat of glaciers in the region. We find no evidence that these two systems came into contact in the deltoic mosaic that would have been present on the exposed continental shelf during the last glacial maximum [16].
Our estimate of the timing of this fragmentation event could be much older than 137 ka because mitochondrial introgression (see below) occurred after the two systems and the three nuclear loci diverged. Therefore our estimate of 137 ka could be biased when we combined these two genomes into a single analysis. Observation of reciprocal monophyly in nuclear loci over a relatively recent period further suggests that our divergence estimate is biased by mitochondrial introgression, however, the extremely low historical effective population size (a population bottleneck) observed in the Chico-Chubut system suggests that monophyly could have occurred over short evolutionary time.
Samples from the Vaca and Telsen Rivers are morphologically similar to Aegla neuquensis but are in fact more closely related to other Aegla taxa. These rivers are geographically isolated from other systems and have not been connected to rivers that contain other A. neuquensis. This evolutionarily independent group was treated as a new taxon in our analyses and will be formally described in a forthcoming work.

Post Fragmentation History
Post fragmentation histories of these two systems were different. First, the Negro system appears to have been demographically stable and has maintained a greater amount of genetic diversity throughout its history. In contrast, the Chico-Chubut system has relatively lower levels of diversity even after a period of demographic expansion over the last several thousand years. In addition, the most recent common ancestor is younger in the Chico-Chubut system than in the Negro system.
Gene flow recovered using MIGRATE suggests very different patterns in these two systems and within individual rivers. Downstream gene flow in both mitochondrial and nuclear markers is extensive within most of the Negro system, with very little movement upstream. Conversely, gene flow is bidirectional between populations in the headwaters of the Limay River. Gene flow in the Chico River is exclusively upstream. Bidirectional flow is seen in the Chubut River, with the Los Altares population serving as the primary source of emigrants within this river. Both systems exhibited restricted gene flow. Differences in gene flow appear consistent with the overall patterns of genetic diversity and demographic history, in the following manner. First, the stable demographic history and larger genetic diversity of the Negro system is characteristic of populations at equilibrium and therefore not colonizing new habitats [34]. These observations would suggest that both females and males are being carried downstream within this system and that very few individuals are emigrating out of the Negro River upstream into the Limay or Neuquen Rivers. In Figure 5. Haplotype networks with nested clade design for mtDNA (COI & COII) sequence data comprised mostly of individuals from the Neuquen River. Haplotype circle size is proportional to its frequency. Each haplotypes is numbered. Haplotypes that exceed twenty individuals are labeled with the total number of individuals of that haplotype. Arrow indicates the root haplotype. Black empty circles represent extinct (or unsampled) haplotypes. Color codes indicate the six major rivers. doi:10.1371/journal.pone.0037105.g005 stark contrast, the large amounts of upstream movement detected in the Chico-Chubut system especially in the Chico River is indicative of recent colonization, an inference consistent with the observation of low genetic diversity and recent demographic expansion. The best demographic model detected some gene flow between systems after they were fragmented. This gene flow would have occurred within the last 137 ka. Only during the last glacial maximum (23-25 ka) would there have been the opportunity for individuals to move between systems via another paleolake. Our best model suggests that the majority of individuals moved south from the headwaters of the Neuquen River into the Chico River, though some appear to have dispersed in the reverse direction.

Comparative Patterns and Processes Between Cis and Trans Andean Aegla Taxa
Patagonian freshwater systems are unique and geographically diverse. Andean orogeny (16-23 million years ago [Ma]) played a major role in shaping these systems primarily due to the mountain range's relative proximity to Pacific and Atlantic coasts ( Fig. 1:   Figure 7. Haplotype network with nested clade design for mtDNA (16 S) sequence data comprised of all samples. Haplotype circle size is proportional to its frequency. Each haplotypes is numbered. Haplotypes that exceed twenty individuals are labeled with the total number of individuals of that haplotype. Arrow indicates the root haplotype. Black empty circles represent extinct (or unsampled) haplotypes. Color codes indicate the six major rivers. doi:10.1371/journal.pone.0037105.g007 [35]). Pacific systems are short in length (,200 km) and fastflowing due to steep elevational gradients [36]. In comparison, Atlantic drainages are longer slow-moving rivers that are geographically distant from one another. Glacial cycles likewise affected Pacific and Atlantic river systems [37]. The impact of glaciation cycles on the two regions was remarkably different [36,38]. The LGM covered a significant portion of Pacific drainages, but only a small portion-primarily in the headwatersof Atlantic drainages were covered. These differences in geology and glacial impacts suggest that freshwater organisms in these two regions experienced different evolutionary histories [8].
Glacial cycles and geography played fundamental but different roles in the evolutionary history of A. neuquensis (this paper) and A. alacalufi [24]. Overall phylogenetic, NCPA, and demographic Figure 8. Haplotype networks with nested clade design for the nuclear EF1a intron and exon. Haplotype circle size is proportional to its frequency. Each haplotypes is numbered. Haplotypes that exceed twenty individuals are labeled with the total number of individuals of that haplotype. Arrow indicates the root haplotype. Black empty circles represent extinct (or unsampled) haplotypes. Color codes indicate the six major rivers. doi:10.1371/journal.pone.0037105.g008 Figure 9. Haplotype networks with nested clade design for the nuclear ANT intron data. Haplotype circle size is proportional to its frequency. Each haplotypes is numbered. Haplotypes that exceed twenty individuals are labeled with the total number of individuals of that haplotype. Arrow indicates the root haplotype. Black empty circles represent extinct (or unsampled) haplotypes. Color codes indicate the six major rivers. doi:10.1371/journal.pone.0037105.g009 Asterisks denote inferences in common across two or more loci (see Table 4). doi:10.1371/journal.pone.0037105.t003 Table 4. Multilocus NCPA [28,29] inferences shared across more than one loci. analyses suggest that glacial cycles were key drivers in shaping the evolutionary history of A. alacalufi [24]. For example, six clades with deeper divergences and more genetic structuring were detected in non-glaciated island populations. A single-shallow clade was recovered in the continental-glaciated populations.
Populations of A. alacalufi persisted in the relatively ice-free western islands, whereas much of the mainland distribution was recently colonized. Xu et al. [24] proposed a novel colonization route where individuals followed glacial waters flowing from higher to lower elevations. Remarkably, this relatively complex history was formed in a relatively short time period over the last two glacial cycles. Glacial cycles played a very different, but nonetheless important, role in shaping the evolutionary history of A. neuquensis. In this taxon, only a small portion of the total distribution was directly impacted by glacial cycles [9,18]. Regardless, this limited exposure to the direct effects of glacial cycles dramatically shaped the evolutionary history of this taxon by fragmenting it into two (three including the Telsen+Vaca clade) lineages. Later, during the lastglacial maximum gene flow occurred when another paleolake formed in the headwaters of the two systems. Demographic expansion was experienced after the last-glacial maximum in the Chico-Chubut system.

A Case of Mitochondrial Introgression
Mitochondrial loci have long been the standard marker for phylogenetic and phylogeographic analyses because they are more likely to recover recent evolutionary history [3]. In this study, we used both nDNA and mtDNA to provide independent lines of evidence on the evolutionary history of Aegla neuquensis. Nearly all phylogeographic studies that have employed both nuclear and mitochondrial data have observed greater resolution in mitochondrial loci, including relatively deeper levels of genetic divergence in the mitochondrial versus the nuclear loci [3]. In this study of Aegla neuquensis, the opposite pattern was observed, though within rivers  AIC calculated using -2(log) +2k, (where k is the number of free parameters in the model). Models in bold as in table 5. Based on AIC score the model (ABADE) provides the best fit to our data. This model had equal population sizes between the Negro and ancestral populations and unequal gene flow between Negro and Chubut systems ( Table 5). Estimates of the time of divergence (years before present [to the nearest thousand]) between river systems were obtained by dividing t by the geometric mean of the mutation rate across loci ( = 4.008839610 25 : see methods for further details). doi:10.1371/journal.pone.0037105.t006 Figure 10. Gene flow estimates for mtDNA and a combined nDNA analysis using MIGRATE [71]. A thick arrow indicates gene flow that is at least two-times greater in one direction than the opposite direction. Color codes indicate the six major rivers. doi:10.1371/journal.pone.0037105.g010 the mitochondrial locus was more informative. Multiple analyses of independent nuclear loci suggest that populations of Aegla neuquensis in the Chico-Chubut system are evolutionarily independent from those in the Limay-Negro system. Phylogenetic analysis of our mitochondrial locus recovered a non-monophyletic Limay system with Chico-Chubut haplotypes embedded within it. Additionally, and again contrary to theoretical expectations, genetic diversity was greater (between Limay-Negro and Chico-Chubut) in the nuclear loci than in the mitochondrial locus. Why are there conflicting signals between genomes? Several alternative explanations have been proposed to explain conflicting signals between genomes [39,40]. First, because of the stochastic nature of the evolutionary process one might detect greater resolution in a nuclear locus over mitochondrial simply by chance [3]. We reject this hypothesis because it is unlikely that three independently sorting nuclear loci would, by chance, recover a stronger signal than a more rapidly evolving mitochondrial locus. Another possible explanation is that the mitochondrial genome is under selection [39]. However, we did not detect selection in the mitochondrial locus. Extremely high female, relative to male, dispersal could explain the differences in loci. But to our knowledge no data exist to suggest a sex-based dispersal bias in this or any Aegla taxon and the breeding system of Aegla suggests that dispersal would be much greater in males than in females. Furthermore, our gene flow analysis suggests similar local patterns of gene flow for both nuclear and mitochondrial data sets.
The most likely explanation for the contradictory signal between genomes is interspecific mtDNA introgression with restricted nuclear gene flow. This phenomenon of ''mitochondrial introgression'' is observed when the mitochondrial genome of one species occurs against a predominant nuclear background of another species [40,41,42]. During hybridization, each species' mitochondrial genome was exposed to the genetic backgrounds of the alternate species, and one species' mtDNA variant (here individuals from the Limay River) may have relatively higher fitness than individuals from the Chico-Chubut system, causing the displacement of mtDNA in the Chico-Chubut by samples from the Limay River [40]. For the nuclear genomes involved in hybridization, selection against nuclear gene introgression in backcross generations could maintain their divergence and distinctness. Because of recombination of nuclear genes, backcross offspring could have various degrees of nuclear genes from the two parental species. If the genes of a parental species have co-adapted for a long time, it is conceivable that backcross offspring having more genes from a single species are favored over those having a mixture [43], i.e., reduced fitness in backcross generations caused by disruption of co-adapted parental gene complexes would be expected [41,44,45].

Conclusions
A detailed understanding of the patterns and processes that shaped the evolutionary history of Aegla neuquensis was possible because we employed multiple analytical methods and loci from both nuclear and mitochondrial genomes [46]. Multiple loci allowed us to detect a recent fragmentation event between two of the major river systems of the region. We suggest that this fragmentation event was driven by the disappearance of a paleolake that had formed at the headwaters of these two river systems. By comparing phylogenetic patterns between mitochondrial and nuclear loci, we identified a rare case of mitochondrial introgression. Our interpretation of the phylogeographic history would have been misleading without the addition of nuclear loci. Future study of aquatic taxa in these two systems will determine the taxonomic breadth of the impact of this fragmentation event.

Sampling and Sequence Data
A total of 295 specimens were collected during 2006-09 from 21 locations spanning the entire distribution of A. neuquensis, an area approximately 3.4610 5 km 2 ( Fig. 1 and Table 1). Specimens were caught by dipnet or hand. Gill or muscle tissues were removed for DNA analysis. Individuals were preserved in ethanol and deposited in the crustacean collections of the Monte L. Bean Life Science Museum at Brigham Young University and the arthropod collections of Universidad Austral de Chile. As invertebrates, no specific permits were required. Locations were not privately owned nor protected and this species is not endangered.

Phylogenetic Analysis
Phylogenetic analyses for the mitochondrial and nuclear data sets were conducted separately using maximum-likelihood (ML) and Bayesian methods. The mitochondrial and nuclear data consisted of haplotypes of concatenated sequences from A. neuquensis and sequences of outgroup species (Aegla ringueleti, A. septentrionalis, A. humahuaca, A. jujuyana, A. sanlorenzo, A. inercalata, A. scamosa, A. platensis, A. uruguayana, [23]). The best-fit models of sequence evolution were determined using Akaike information criterion (AIC) as implemented in JModeltest [50]. ML analyses were performed using RAxML 7.0 [51], with 1000 bootstrap runs and searching for the best-scoring ML tree. Individual a-shape parameters, GTR-rates, and base frequencies were estimated and optimized for each partition. Bayesian analyses were performed in MrBayes 3.2 [52], using four independent runs with random start trees. We ran 10 million generations using four chains, with sampling frequency of one per 1000 for each chain. Results were visualized and checked with Tracer 1.5 [53] and a burn-in of 20% discarded. Divergence times (time to the most recent common ancestor -TMRCA) of major nodes were estimated with BEAST 1.6.1 [54]. Phylogenetic analyses of the mitochondrial data  suggested that they might be under selection (see results). We used a McDonald-Kreitman test [55] as implemented in DnaSP 5.0 [56] to determine if the mitochondrial data were under selection.

Multilocus Nested Clade Phylogeographic Analysis
Multilocus nested clade phylogeographic analysis (NCPA) was used to test for historical and demographic events that are geographically congruent across two or more loci [28,29]. This method has been recently criticized for its high false-positive rates in simulation studies [57,58]. In this study NCPA was used in conjunction with other approaches [59] to examine both patterns and processes of divergence across loci and geography [46]. We consider inferences drawn from multiple approaches as strong evidence. Furthermore only inferences found in two or more loci will be considered [28,29]. Haplotype networks were constructed with TCS [60]. Associations between haplotypes and geography were tested with GEODIS [61]. Ambiguous connections or loops in the networks were resolved following the three criteria listed in Pfenninger and Posada [62]. We calculated distance along rivers instead of direct distance for geographical analysis [63]. To show all relationships in the nuclear loci, we set the connection limit to 40 steps for the nuclear data. We used the default 95% connection limit for the mtDNA. Allelic phase for nDNA loci for this analysis as well as the IMa and MIGRATE analyses (next sections) were determined using PHASE [64,65] with default settings.

Model Selection and Demographic History
Demographic history between Chubut (Chico and Chubut) and Negro (Limay and Neuquen) river systems, using nDNA and mtDNA, was examined using an isolation-with-migration analysis implemented in IMa [66]. This analysis was conducted in successive stages. First, exploratory analyses were performed to determine appropriate prior settings (i.e., posterior distributions were contained within the prior range) and the number of chains that gave consistent results across multiple runs. The priors, chains, and other starting conditions determined from exploratory analyses were: q1 = 5, 2q2 = 5, -qA = 10, 2m1 = 1, 2m2 = 1,t = 18, and -n = 3. Three independent analyses were then run for 72 hours using MCMC simulations to sample genealogies and obtain parameter estimates. Genealogies from these three analyses were then combined into a single data set. Next, we evaluated the relative fit of 16 nested-demographic models [67] using a portion ( = 20,000) of genealogies from the combined dataset using the likelihood option in IMa. Five population parameters are estimated in each model. These are: ancestral theta h A , theta of descendent population one h 1 , theta of descendent population two h 2 , and two measures of gene flow (m 1 and m 2 ) between the two daughter populations. In addition, the time of the split between daughter populations were estimated using IMa. The best models were determined using a chi-square distribution with and without Bonferroni correction. Akaike Information Criterion (AIC) was also employed to rank the models (see [67] for further details).
Time (t) was converted to years before present using a rate of 0.118 substitutions/site/million years that was estimated from the closely related A. alacalufi [24]. Mutation rates for each nuclear locus were obtained by multiplying the ratio of the mutation rate scalars (mtDNA/nDNA) estimated during the IMa runs to the mtDNA mutation rate [68]. The final estimates of divergence times were obtained by dividing t by the geometric mean of the mutation rate across loci per generation [68].

Population Analysis
Population genetic diversities, as indicated by haplotype diversity (Hd), nucleotide diversity, current genetic diversity (h p ) and historical genetic diversity (h W ), were assessed using DnaSP 5.0 [56]. Genetic divergence between populations was estimated using F ST index and population structure was measured with an analysis of molecular variance (AMOVA), as implemented in [69]. Significance levels for multiple tests were adjusted with the sequential correction [70].

Estimating Gene Flow
Patterns of gene flow between populations were determined using MIGRATE 3.2.6 [71,72]. Estimates were obtained by averaging the mean values across three independent analyses. A single analysis consisted of 10 short chains with every 20 steps retained with a total of 1000 steps recorded. Then three long chains were run for 5 million generations, with every 100th steps recorded. Each long chain consisted of four adaptively heated chains determined using default settings. Each analysis was replicated four times with the final estimate averaged across replicates. By design we set the burnin to 10,000 steps. Due to relatively low sample sizes, four pairs of populations were grouped into single populations. These were Chimpay and Bocatoma; Chubut and Leleque; Ameghino and Dolavon; and the Mayo and Sarmiento populations. Due to low nDNA haplotype diversity in the Chico-Chubut system, all samples were grouped into a single population. The migration matrix was designed to allow the exchange of genes only between populations occurring in the same river system. The only exception to this assumption was the allowance of gene flow between populations at the headwaters of the Limay and Chico-Chubut systems.

Extended Bayesian Skyline Plots
Historical demography of the Chico-Chubut and Limay-Neuquen-Negro river systems was examined using Extended Bayesian skyline plots (EBSP) as implemented in BEAST 1.6.1 and visualized in Tracer 1.5 [53]. These analyses employed all four loci. Each analysis was run for 50 million generations with every 1000 generations retained. The first 10 million generations (20%) were removed as burnin. Multiple runs were conducted on each data set to check for consistency. Final parameter estimates were obtained by combining three runs. We concluded the analyses were well sampled when three independent runs converged on similar values and had .100 ESS values for our focal parameters. The first of the three runs was used to generate a skyline plot. Analyses were run with a coalescent tree and relaxed uncorrelated lognormal clock priors. We employed the mtDNA mutation rate of 0.118 substitutions/site/million years [24] and allowed the analysis to estimate the rates, relative to the mtDNA, for each nuclear loci under a lognormal prior. Site models determined using jModeltest [50] for each of the four loci were HKY+I+G (mtDNA), JC+I (EF1 Exon), GTR+I (ANT), and HKY+I+G (EF1 Intron). We provide the Beast xml files (Data File S1 and S2) and TMRCA plots for each loci ( Figure S1). Figure S1 Plots of time of most-recent-common-ancestor for individual loci from three independent runs and from the three runs combined. Analyzed using Extended Bayesian skyline plots (EBSP) as implemented in BEAST 1.6.1 and visualized in Tracer 1.5. Each analysis was run for 50 million generations with every 1000 generations retained. The first 10 million generations (20%) were removed as burnin. Analyses were run with a coalescent tree and relaxed uncorrelated lognormal clock priors. We employed the mtDNA mutation rate of 0.118 substitutions/site/million years [17] and allowed the analysis to estimate the rates, relative to the mtDNA, for each nuclear loci under a lognormal prior. Site models determined using jModeltest [35] for each of the four loci were HKY + I + G (mtDNA), JC + I (EF1 Exon), GTR + I (ANT), and HKY + I + G (EF1 Intron). All ESS values were greater than 500.

Supporting Information
(DOC) Data File S1 Beast input file (.xml) for Extended Bayesian Skyline Plot analysis of the Negro River system. (PDF) Data File S2 Beast input file (.xml) for Extended Bayesian Skyline Plot analysis of the Chico-Chubut River system. (PDF)