Research article | Open | Published:
Cryptic diversity and deep divergence in an upper Amazonian leaflitter frog, Eleutherodactylus ockendeni
BMC Evolutionary Biologyvolume 7, Article number: 247 (2007)
The forests of the upper Amazon basin harbour some of the world's highest anuran species richness, but to date we have only the sparsest understanding of the distribution of genetic diversity within and among species in this region. To quantify region-wide genealogical patterns and to test for the presence of deep intraspecific divergences that have been documented in some other neotropical anurans, we developed a molecular phylogeny of the wide-spread terrestrial leaflitter frog Eleutherodactylus ockendeni (Leptodactylidae) from 13 localities throughout its range in Ecuador using data from two mitochondrial genes (16S and cyt b; 1246 base pairs). We examined the relation between divergence of mtDNA and the nuclear genome, as sampled by five species-specific microsatellite loci, to evaluate indirectly whether lineages are reproductively isolated where they co-occur. Our extensive phylogeographic survey thus assesses the spatial distribution of E. ockendeni genetic diversity across eastern Ecuador.
We identified three distinct and well-supported clades within the Ecuadorean range of E. ockendeni: an uplands clade spanning north to south, a northeastern and central lowlands clade, and a central and southeastern clade, which is basal. Clades are separated by 12% to 15% net corrected p-distance for cytochrome b, with comparatively low sequence divergence within clades. Clades marginally overlap in some geographic areas (e.g., Napo River basin) but are reproductively isolated, evidenced by diagnostic differences in microsatellite PCR amplification profiles or DNA repeat number and coalescent analyses (in MDIV) best modelled without migration. Using Bayesian (BEAST) and net phylogenetic estimates, the Southeastern Clade diverged from the Upland/Lowland clades in the mid-Miocene or late Oligocene. Lowland and Upland clades speciated more recently, in the early or late Miocene.
Our findings uncover previously unsuspected cryptic species diversity within the common leaflitter frog E. ockendeni, with at least three different species in Ecuador. While these clades are clearly geographically circumscribed, they do not coincide with any existing landscape barriers. Divergences are ancient, from the Miocene, before the most dramatic mountain building in the Ecuadorean Andes. Therefore, this diversity is not a product of Pleistocene refuges. Our research coupled with other studies suggests that species richness in the upper Amazon is drastically underestimated by current inventories based on morphospecies.
Species richness and genetic diversity within species are proposed to co-vary and understanding the details of this relationship is critical to unifying biodiversity theory . Though the forests of the upper Amazon basin are a renowned hotspot for amphibian species richness , thus far there have been few explorations of the phylogeographic and population-level patterns in the region's amphibian taxa [3–9]. Consequently, we have only a preliminary knowledge of the distribution of genetic diversity from merely a handful of amphibian species in the megadiverse upper Amazon. With the aim of augmenting our understanding of spatial patterns of genetic diversity and details of evolutionary history in this famously species-rich area, we assessed the diversity in a widespread upper Amazonian frog species across the upper Amazon of Ecuador.
Historical and environmental characteristics of a region influence multiple levels of diversity, both in its origin and maintenance [1, 10–12]. Numerous regional historical, topographical, and ecological factors of the Andes and Amazon have been suggested as influential in the diversification of species, for example: riverine barriers [13–16]; large uninterrupted area [10, 17]; an Andean "species pump" promoted by complex, isolating montane topography and vegetation ; Pleistocene forest refuges [11, 19–21]; complex biotic interactions [22, 23]; ancient ridges and other palaeogeographic features [4, 24]; ancientness ; and relative youth of the Andes . For amphibians, complex topography may restrict vagility and produce marked population subdivision, and ultimately patterns of adaptive radiation or isolation patterns that are potentially implicated in high local rates of speciation . There are many confounds in the upper Amazon that make exclusive testing of any one of these factors challenging because multiple causes at various historical timescales are surely involved in the origination of diversity [8, 28]. Nonetheless, contributing temporal and geographical patterns of diversity in phylogeographic studies allows us to discern among some of the diversification hypotheses (e.g., ancientness versus youth, orogenesis events versus recent climatic change). Assessing these patterns is critical to disentangling the causes of within and among species diversification and the genesis of anuran biodiversity, especially as this region increasingly suffers deforestation.
Here we use mitochondrial and nuclear DNA to quantify the phylogeography and population structure of a terrestrial upper Amazonian leaflitter frog, Eleutherodactylus ockendeni (Leptodactylidae). We sampled from 13 localities across the species range in megadiverse eastern Ecuador to develop a thorough regional phylogeographic survey. Further, we used phylogenetic and multiple coalescent methods to estimate the depth of divergence among clades and suggest temporal context for the divergence revealed by our analyses.
We sequenced 105 individuals for cytochrome b (cyt b) and 45 individuals for 16S, plus two outgroup taxa (Additional file 1). Ingroup base frequencies are similar to those found in other frogs : cyt b: A = 0.240, C = 0.313, G = 0.145, T = 0.301; 16S: A = 0.299, C = 0.251, G = 0.194, T = 0.255. Ingroup cyt b sequences collapsed into 21 unique haplotypes and 16S into 15 unique haplotypes. For the 16S-cytb data combined (32 individuals; 20 haplotypes and two outgroup taxa; 1248 bp), 227 included characters were parsimony informative.
We found no conflicting phylogenetic signal between the cyt b and 16S data (partition homogeneity test, P = 1.0), justifying the use of the combined data partition in maximum parsimony (MP) analyses. Both MP and Bayesian analyses of the 16S-cytb data produced well-supported trees identical in all major topological details and with three major clades: an Upland Clade (Figure 1 localities Hola Vida, Santa Clara, EBJS, Llanganates, Cando, and Chonta Yacu), a Lowland Clade (Yasuní, Puca Chicta, EBJS, Serena, Auca 14, La Selva, and Cuyabeno); and a Southeastern clade (Kapawi and Auca 14) (Figure 2). The topology of the 16S-cytb trees is congruent with trees resulting from separate analyses of the full 16S and cyt b data (results not shown).
Corrected p-distances between E. ockendeni cyt b haplotypes were high, ranging up to 20% between some haplotypes (Additional file 2). Mean net divergence ± standard error  between the Southeastern Clade and the Upland Clade was 15.45% ± 1.83, between the Southeastern Clade and the Lowland Clade was 15.26% ± 1.90, and between Upland and Lowland Clades was 12.01% ± 1.72. Average sequence divergence was low within the Upland and Lowland clades, at 1.35 and 1.77% respectively, while the Southeastern Clade, which has fewer haplotypes and more geographically distant sampling, showed more intraclade diversity at 5.34%.
For the cyt b data, the identity of haplotypes varied among localities but seven of the 13 sites had only a single haplotype (Table 1). Nucleotide diversities within localities range from 0 to 0.105 ± 0.079. However, high nucleotide diversity within localities EBJS (5) and Auca 14 (7) is an artefact of finding sympatric but genetically distinct clades at those localities. When calculations of locality nucleotide diversity are separated by clade, the values range from 0 to 0.0030 ± 0.003 (Table 1). HKY+G was the best-fit model of nucleotide substitution for cyt b as chosen by hLRT.
The ratio of nonsynonymous to synonymous mutations between all pairs of clades as calculated with a McDonald-Kreitman test does not suggest that selection is responsible for the cyt b divergence among lineages (Upland versus Lowland, P = 1.0; Upland versus Kapawi, P = 1.0; Kapawi versus Lowland, P = 0.3). Overall non-significant value of Tajima's D (0.255; P > 0.1) and within clades (Upland Clade D = -1.707, P > 0.05, Lowland Clade D = -0.834, P > 0.1, insufficient sample size to test Southeastern Clade) also implies neutrality (Table 2).
A survey of five tetranucleotide microsatellite loci across individuals of the Upland and Lowland clades in the Napo River area (localities 3 – 11), including a locale where both mitochondrial clades are present, shows that microsatellite genotypes at all loci are very different between Upland and Lowland clades. This implies complete reproductive isolation. The microsatellite library was developed based on Lowland Clade frogs. Consequently, in individuals from the Upland Clade the microsatellite loci are either non-functional (not amplifying after repeated attempts) or allele sizes are non-overlapping or have very different range of sizes, with the Upland Clade always having the larger mean allele size (Table 3). Sequences of large alleles from a subset of samples from three of the five loci (Eloc-Batman&Robin, Eloc-Laurel&Hardy, and Eloc-Thelma&Louise) demonstrate that very large allele sizes in Upland Clade individuals are due to an increase in the number of microsatellite repeats (unpubl. data).
Estimating divergence and expansion
We used three different methods to estimate divergence: the net divergence method to estimate species divergence time; a Bayesian MCMC method to estimate lineage divergence (time to most recent common ancestor, TMRCA); and, a coalescent MCMC approach to estimate migration and TMRCA. Of these, the net divergence method and the Bayesian MCMC method (as implemented in the program BEAST) gave similar temporal results. Divergence estimates from coalescent analyses (as implemented in the program MDIV) were typically less than half the age of the other two estimates.
By the net divergence method, time of divergence between the Upland and Lowland lineages using our slowest estimated rate of evolution (0.6%/MYR) was 20.02 ± 2.87 mya (early Miocene) while the faster rate (1.0%/MYR) estimates the same split at 12.01 ± 1.72 mya (mid-Miocene) (Table 4). Southeastern Clade split from Upland and Lowland clades approximately 25 ± 3 mya (late Oligocene) at the slower substitution rate and 15 ± 2 mya (mid-Miocene) by the faster rate.
From BEAST, assuming a constant molecular clock and rates of 0.6 and 1.0% substitutions per million years, we estimated the TMRCA of the entire ingroup to be 24.39 mya (late Oligocene) and 14.61 mya (mid-Miocene), respectively. For the TMRCA of the Upland and Lowland clades, the constant clock TMRCA estimates 15.19 mya (mid-Miocene) or 9.11 mya (late Miocene), respectively. The uncorrelated, relaxed clock estimates were not substantially different: assuming mean substitution rates of 0.6 and 1.0%, the TMRCA estimates for the entire ingroup were 27.01 and 15.10 mya, respectively, while those for the Lowland/Upland clade were 15.23 and 9.03 mya, respectively (Table 5).
Coalescent calculations (implemented in the program MDIV) of the Upland and Lowland clade divergence resulted in an average θ of 6.23 ± 0.13 and an average T of 8.78 ± 0.31 (Table 6). Estimates of TMRCA are 7.74 mya assuming a substitution rate of 0.6%/MYR and 4.64 mya assuming the faster substitution rate of 1.0%/MYR. In models of Southeastern Clade versus Upland Clade divergence, θ averaged 8.43 ± 0.27 and T 7.73 ± 0.52, suggesting a TMRCA of 10.55 mya under the slower substitution rate and 6.33 mya, with the faster rate. Models of Southeastern versus Lowland clades identified an average θ of 9.45 ± 0.38 and an average T of 8.98 ± 0.69, or 9.99 mya or 7.77 mya TMRCA, depending on substitution rate. For all clade comparisons M modelled best as 0 suggesting there is no gene flow among clades. The estimated time since population divergence among all clades is less than the TMRCA. MDIV estimates of TMRCA and species divergence are non-overlapping with BEAST and the net divergence method.
Estimating population expansion
We used four different methods to try and identify and estimate the timing of population expansion. Based on mismatch analyses, the sudden expansion model of population growth cannot be rejected for the Upland or the Lowland clades, although the SSD probability was marginal for the Lowland Clade (P = 0.056) (Table 2). Using a mutation rate of 1.0%/MYR and the peak of the mismatch distribution, τ, to estimate the time of population expansion suggests an Upland Clade expansion began in the latter half of the Pleistocene (793 000 YBP), a Lowland Clade expansion (if it occurred) somewhat later (154 000 YBP). The model of sudden population expansion is rejected for the Southeastern Clade. Raggedness indices suggest population expansion (curves are not significantly different than smooth) in all three clades. Fu's F and Tajima's D are not different than would be found under a stable population size. Therefore, we have conflicting evidence from different methods but some suggestion of population expansion particularly in the Upland Clade.
Phylogenetic and phylogeographic relationships
Divergence among the three strongly supported clades in Ecuadorean E. ockendeni is deep and well supported (Figure 2, Additional file 2). In contrast, there is relatively low divergence within the three clades. The Southeastern Clade is basal relative to the Upland and Lowland clades. The timing of the divergences for all three lineages suggests causal events that vastly predate the Pleistocene and the latter phases of orogenesis in the northern Andes. We discuss these below.
Five independent nuclear markers corroborate that the Upland and Lowland clades are evolutionarily distinct. The microsatellite loci were either non-amplifying or composed of twice as many repeats in the Upland as Lowland Clade, even in frogs from the same geographic locality (Table 3). Such differentiation at microsatellite loci between clades even when they are sympatric and the potential for interbreeding exists provides strong indirect support that these are reproductively isolated. Unfortunately, we still know little of the mate recognition system for this complex of species.
The major E. ockendeni clades are geographically restricted and found in the east Andean uplands (approx. 400 to 1000 masl; Upland Clade), the central and northern lowlands (220 to 500 masl; Lowland Clade), and the central and southern lowlands (260 to 240 masl; Southeastern Clade). Our sampling is limited to Ecuador and therefore we cannot infer whether the clades are more broadly distributed. The headwaters of the Napo River appears to be part of a lineage contact zone: at EBJS (locality 5) both the Upland and Lowland Clade can be found; slightly upriver at Serena (3) and Cando (4), the Lowland Clade is found on the south side of the Napo River while the Upland Clade is found immediately on the opposite side of the river. Further east at Auca 14 Road (locality 7), again two lineages can be found sympatrically (Upland and Southeastern clades). At all other geographic localities sampled we found only one of the three clades. Physalaemus petersi was also found to have an upper Napo and lower Napo genetic division and a basal southeastern lineage . Further, an upper and lower Napo basin phylogeographic break has been suggested in Bolitoglossa salamanders . Therefore, evidence from other amphibians may indicate this to be a cryptic geographic break, as has been found in other upper Amazon localities .
In E. ockendeni, the lack of haplotypic diversity at many localities and general lack of haplotype sharing among localities suggests fine-scale restricted gene flow (Table 1), as might be expected from a small terrestrial amphibian . Population structure at finer geographic scales in the upper Amazon has also suggested restricted gene flow in this species .
We have included E. ockendeni from Cuyabeno (locality 2) in northeastern Ecuador in the Lowland Clade; however, they do display a notable genetic divergence from the rest of the Lowland Clade (> 6% corrected p-distance; Additional file 2) in a group that is otherwise characterized by low within-clade diversity. Further sampling from the surrounding geographic area (e.g., northern Peru and southern Colombia) and more molecular markers will either strengthen our inclusion of Cuyabeno E. ockendeni in the Lowland Clade or differentiate it as a separate clade.
Time of divergence
Our divergence estimates suggest that diversification among the "species" E. ockendeni is ancient and predates the Andes at their current height in Ecuador. Divergence between the Southeastern and Upland/Lowland lineages occurred in the late Oligocene or early Miocene and between the Upland and Lowland Clade in the early, mid-, or late Miocene, depending on substitution rate.
Of the three estimates, we consider the Bayesian phylogenetic-coalescent method (Table 5) to be the most accurate because that method accommodates the greatest number of parameters, incorporates molecular evolutionary complexities such as rate heterogeneity, allows for differences in rates among lineages, and allows for tests of the molecular clock . Interestingly, and in support of a long-standing method, the much less complicated net divergence method (Table 4) has yielded very similar species divergence time estimates to the Bayesian method TMRCA. We suggest that the MDIV coalescent method (Table 6) is underestimating time of divergence in this case, perhaps because there is not enough historical information to suit a population coalescent method when species are reciprocally monophyletic and deeply diverged, or because MDIV does not accommodate rate variation among sites.
Some studies have found the rate of molecular evolution in the tropics to be equivalent to temperate areas [34, 35] while others have suggested a faster rate of molecular evolution in the tropics [36–38]. Obviously, the timing of the divergence among E. ockendeni clades will become younger given the same genetic distance if true substitution rate is faster than we have estimated.
Historical population change
Across three different tests, there is weak and somewhat conflicting suggestion of recent population expansion under models of neutral evolution (Table 2). For the Upland Clade, we cannot reject sudden expansion from mismatch distributions or raggedness though Fu's F and Tajima's D do not suggest population expansion. For the Lowland Clade, we can probably reject population expansion since mismatch and raggedness significance values are only barely non-significant (P = 0.056 for mismatch, P = 0.107 for raggedness) and F and D do not suggest expansion. In the Southeastern Clade, sample sizes are very small (two localities, two haplotypes), but mismatch and F do not suggest population expansion while raggedness index does fit the expectation of expansion.
Our results give some support for an Upland Clade expansion in E. ockendeni approximately 800 000 years ago. It is difficult to determine historical population expansion; inconclusive findings from multiple methods are not unexpected because rate heterogeneity, the scale of historical population expansion if there were any, sample size, and number of polymorphic sites are all influential factors . If there has been expansion in the Lowland Clade, this may have occurred 150 000 years ago. In upper Amazonian clades of the terrestrial frog P. petersi, inconclusive results about a population expansion have also been found and, like our results here, some support for an Upper Napo clade expansion .
Geological history and its influence on phylogeographic patterns
The upper Amazon of Ecuador has a history of variations in climate and dramatic geology. The northern Andes are the youngest in the Andean chain: orogenesis occurred in the Neogene (Miocene through Holocene) , with the principle period of upheaval only in the past five million years [40–42].
Throughout the Miocene there was a radical reorganization of drainage patterns in northern South America. In the early Miocene the fluvial system drained to the northwest but by the late Miocene the effects of the rising Eastern Cordillera shifted drainage eastwards, ultimately allowing the Amazon River to reach the Atlantic . Episodic marine incursions as a result of eustatic sea level changes occurred throughout the Miocene [24, 41]. In the Pliocene, Pleistocene, and Holocene, fluvial deposits of many sorts and very large alluvial fans (megafans) were laid throughout the upper Amazon in the wake of Andean tectonism, volcanism, and later glaciation ([41, 43–45] and references therein). These overlays and heterogeneity often leave conflicting evidence and cause geological events east of the Andes to be poorly dated [43, 46].
Climate also has varied over time, with great fluctuations throughout the Tertiary and Quaternary . Alternating wet and dry cycles in the Pleistocene and Holocene may have caused forest contractions and expansions in some areas, though there is disagreement on the extent of this effect in the upper Amazon [21, 47–49]. These forest contracts are the source of the refuge hypothesis of Amazonian diversification .
Clearly there are multiple possible ancient historical geological and climatic influences causing isolation and subsequent expansion in E. ockendeni, ultimately resulting in the speciation we see here. We can infer that diversification between clades of E. ockendeni predates the altitude and shape of the Andes as we know them. Instead, diversification is much older, perhaps precipitated by dramatic changes in the Miocene. For example, mountain-building caused numerous thrust faults to develop and the eastern Subandean Zone fault ( and references therein) approximately matches the geographic location of our Upland/Lowland break between clades of E. ockendeni and may be historically relevant. Geological change in the Miocene has recently been suggested as influential in South American Eleutherodactylus frog diversification in general . Importantly, there is no obvious contemporary barrier between the extant clades of E. ockendeni, such as a major river as would be suggested by the river barrier hypothesis [13–16]). The current lack of elevation difference in the distribution of two of the three clades suggests that elevational gradient differences are not driving the divergence (e.g., under an environmental gradient hypothesis [51, 52]). Contemporary complex topography (as suggested by, for example [27, 53]) does not seem to be relevant to patterns of deep divergence among clades, since divergence far predates existing topography, clades are not isolated by major topographic features, and the distribution of the Upland Clade extends through the headwaters of at least three major river valleys (Agua Rico, Napo, and Pastaza). Further, the genetic diversity in E. ockendeni dates from the Miocene and therefore cannot be attributed to climatic and associated vegetation changes of the past two million years, as has been suggested by proponents of the refuge theory [11, 20, 54]. However, recent population expansion in the Upland Clade may coincide with late Pleistocene climate change.
Cryptic species richness in E. ockendeni
Our findings show that the leaflitter frog species E. ockendeni is three distinct species with apparently extreme morphological conservativism. Furthermore, this area is only a portion of the species' range, which extends from southern Colombia to southern Peru  and Bolivia , so the actual species diversity and richness within "E. ockendeni" is likely much greater than that demonstrated here. When the specimens were collected and catalogued for museum deposition, there was no apparent morphological difference among specimens and there has been no published indication that there would be multiple cryptic species within this species, except for a mention that E. ockendeni from Cuisime (southern Ecuador) are smaller in snout-vent length than those collected in other areas of Ecuador and Peru . It is clear that the divergence among each of the three clades is biologically real and not an artefact of stochastic variance in mtDNA masking recent gene flow  because our coalescent methods model best with no gene flow among lineages (in MDIV M = 0) and five nuclear microsatellite loci are significantly different between the Upland and Lowland clades. Given this new molecular information, a detailed morphological revision of Ecuadorean E. ockendeni is underway (Elmer and Cannatella, unpubl.).
Incidences of sympatric, parapatric, and allopatric cryptic species have been recently discovered in southeast Asian frogs ( and references within). Also, although there are exceedingly few intraspecific molecular phylogenies of neotropical amphibians, those that do exist tend to encounter new species and/or previously unanticipated species diversity ([9, 59, 60] and references within). These findings together suggest that widespread species of amphibians in the tropics have an evolutionarily history that is much more complicated than suggested by morphologies. Consequently, attempts at biological conservation according to current estimates of the number of morphological "species" will drastically underestimate the actual biodiversity in this already species-rich region.
We found deep phylogenetic divisions among clades in this common leaflitter frog, E. ockendeni, suggestive of distinct species. Based on microsatellite genotype profiles for distinct mitochondrial clades and modelling of historical gene flow, we suggest that there is complete reproductive isolation among these clades, even when they are sympatric. These cryptic species have an ancient divergence estimated to have originated in the Miocene. Diversification among these clades coincides approximately with periods of dramatic northern Andean orogenesis and predates the Andes at their current height. Though multiple environmental occurrences surely have been historically influential, Pleistocene climate change refuges as drivers of allopatric speciation are not relevant to extant specific diversification E. ockendeni. Our research strongly suggests that current estimates for the renowned species richness in the upper Amazon of frogs in general and Eleutherodactylus in particular may be a substantial underestimate of the actual phyletic diversity present.
Field and Laboratory Methodology
Eleutherodactylus ockendeni is a small terrestrial leaflitter frog that is relatively abundant at many upper Amazon localities [32, 61, 62]. Its range extends throughout the upper Amazon of Colombia, Peru, Ecuador, Brazil , and Bolivia . We collected from 12 localities across eastern Ecuador (Figure 1) in 2001 and 2003. Inter-locality distances ranged from 1 to 300 km apart (straight-line distance) and included all major river basins in Ecuador. Individuals were euthanized using MS-222 and tissue samples for genetic analyses removed and stored in pure ethanol. Specimens were fixed with 10% formalin and stored in 70% ethanol. Vouchers are deposited at the Museo de Zoología, Pontificia Universidad Católica del Ecuador (QCAZ) (Additional file 1). Samples from the Parque Nacional Yasuní locality were taken from the existing QCAZ collection.
Genomic DNA was extracted using either standard phenol-chloroform methods  or a Qiagen DNEasy kit according to the manufacturer's protocol. Two mtDNA fragments were amplified: 16S rRNA with primers 16Sar-L and 16Sbr-H (ca. 560 bp; numbers 88 and 96 in ) and cyt b with primers MVZ15L and MVZ16H (ca. 790 bp; numbers 141 and 165 in ). Volumes for each PCR reaction were 50 μL and contained: 0.3 μM of forward and reverse primer, PCR enhancing buffer (2.5 mM MgCl2, 10 mM Tris pH 8.4, 50 mM KCl, 0.02 mg BSA, 0.01% gelatin; adapted from ), 0.3 mM dNTP, 0.625 units taq DNA polymerase (Fermentas Life Sciences), and approximately 1 to 3 ng stock DNA. All reaction sets included a negative control. Cycling parameters for 16S were: 94°C for 2 min, 35 cycles of denaturing 94°C for 30 sec, annealing 60°C for 20 sec and extending 72°C for 20 sec, a final extension at 72°C for 5 min followed by extended cooling at 10 or 4°C. Cyt b parameters were: initial denaturing at 92°C for 3 min, 38 cycles denaturing 92°C for 1 min, annealing 51°C for 1 min and extension 72°C for 1 min, a final extension of 72°C for 5 min followed by extended cooling at 4 or 10°C. PCR product was cleaned for sequencing using a Qiaquick Gel Extraction kit according to the manufacturers' instructions for agarose electrophoresis-separated fragments or Pall AcroPrep 96 Filter Plates for PCR products that were not electrophoresed. Samples were capillary sequenced using the BigDye Terminator version 3.1 Cycle Sequencing (Applied Biosystems) chemistry on an Applied Biosystems 3100 or 3730XL Gene Analyzer.
Forty-seven individuals were sequenced for 16S: 23 in both directions and the remainder only in the forward direction. One hundred and seven individuals were sequenced for cyt b: 41 in both directions and the rest only in the forward direction (GenBank accession numbers: 16S EU130581 – EU130626, EU130579, EU130580; cyt b EF581013–EF581063, EU130577, EU130578, EU130627–EU130680). Detailed comparison by eye of forward and reverse sequences in the overlapping regions showed no discordance in the DNA sequence used in subsequent analyses. Thirty-two individuals plus outgroup taxa had sequences for both 16S and cyt b (hereafter, 16S-cytb). In the 16S-cytb fragment, concatenated sequences were trimmed to equal length (1248 bp), except for five sequences that remained shorter (≥ 1197 bp). In the cyt b fragment, all haplotypes were the same length (709 bp) except three (≥ 647 bp). Terminal gaps were coded as missing data.
Sequences were assembled in MACCLADE version 4.07  and aligned in CLUSTAL X version 1.81  using default settings. A hyper-variable loop portion of the 16S alignment could not be aligned with confidence, so nine internal positions were excluded. Alignment of cyt b was not problematic and included no internal gaps. The best-fit models of nucleotide substitution of cytochrome b, each codon position separately, and 16S were estimated using MODELTEST version 3.7  in PAUP* .
Because "intraspecific" diversity in these frogs is so high, we analyzed these data using maximum parsimony (MP) and Bayesian phylogenetic tree methods (rather than network methods ). Before MP analysis, we performed a partition homogeneity test  with 100 replicates and 10 random addition replicates per replicate in PAUP*  to assess whether 16S and cyt b have congruent phylogenetic signal. The 16S-cytb fragment MP analyses were run in PAUP* with outgroups Eleutherodactylus acuminatus and Eleutherodactylus quaquaversus (unistrigatus group ), using a heuristic search strategy with TBR branch-swapping and 1000 random addition replicates. The shortest tree was saved from each replicate. The topology of the 11 most parsimonious trees was tested with 1000 heuristic search bootstrap pseudoreplicates, with 10 random-addition replicates each, and merged into a strict 50% consensus tree. MP analyses were repeated with the full cyt b and 16S data sets separately. For our Bayesian phylogenetic analyses, we used a two partition (16S and cyt b) analysis in MRBAYES version 3.1.2  with likelihood parameters nst = 6, rates = gamma. Two independent MCMC chains were run for 4.8 million generations with reduced temperature difference among chains (t = 0.09) and trees sampled every 100 generations. After assessing for apparent convergence , 25 000 trees were discarded as burn-in and a 50% majority-rule consensus tree was built.
Time of divergence estimates
We used cyt b to estimate divergence times because we had more data for this fragment (individuals and haplotypes), there are more published estimates of rate of evolution available than 16S, and cyt b is less subject to inter-study variation because of alignment-specific exclusion of hyper-variable regions typical of rDNA alignments. Divergence time among lineages or species was estimated in three ways. First, net among clade divergence ( p. 276) ± standard error corrected with a TN model  of nucleotide substitution was calculated in MEGA3 . This method subtracts the intraspecific diversity from overall divergence between two species and is unbiased when lineages are reciprocally monophyletic and ancestral population sizes is equal to the average of the two descendent species [reviewed in ]. Species divergence time was calculated from the net divergence estimate (net divergence × substitutions/site/MYR). The TN model is a more general case of the HKY model , which is implemented in MDIV (below) and was selected as an appropriate model of sequence evolution in MODELTEST.
Second, we estimated lineage divergence time, population divergence, and migration rate between major pairwise clades in MDIV  on CBSU Web Computing Resources. Initial runs were tested under a finite sites (HKY) model of evolution and default priors (M = 10, T = 5) to approximate the posterior distribution of scaled migration rate (M) and time since divergence (T), while allowing MDIV to estimate θ. Once appropriate parameter values to bound a "well-behaved" posterior distribution  were identified (M = 0, T = 30), we ran the MCMC for two million generations with 500 000 generations discarded as burn-in. Convergence was determined by evaluating the consistency of model values for each of the three parameters across five runs, which were then averaged to calculate mean θ and T values ± standard deviation. Time of divergence was calculated as (following ): tdiv = (T θ/2L)/(1/μ) where T (or TMRCA) and θ were estimated by the height of the posterior distribution, L is the sequence length analyzed (709 bp of cyt b), and μ is the mutation rate (here, substitution rate). Substitution rates may be less than mutation rates because neutral mutation rates include population polymorphism that will not eventually be fixed in phylogenetic lineages (reviewed in [82, 83] but see ). Given the high levels of divergence among lineages here, we consider substitution rate more realistic than mutation rate.
As a third approach to calculating the timing of diversification, we estimated time to most recent common ancestor (TMRCA) for various clades using a Bayesian approach with the program BEAST version 1.4.1 . All analyses were performed using an HKY model of nucleotide substitution with gamma distributed rate variation among sites and six rate categories. We ran four separate sets of analyses, first assuming a constant population size and a constant global molecular clock (i.e. no rate differences among lineages) of either 0.6 or 1.0%/MYR, and second using an uncorrelated, relaxed clock again assuming constant population size and mean clock rates of either 0.6 or 1.0%/MYR. Results from two independent runs (10,000,000 generations with the first 1,000,000 discarded as burn-in and parameter values sampled every 1000 generations) for each combination of settings were combined and the effective sample size for parameter estimates and convergence checked using the program Tracer version 1.3 .
The substitution rate of E. ockendeni cyt b is unknown, so for all three methods we used the same two estimates of 0.6 and 1.0 substitutions/site/100 MYR. We inferred our substitution rate from two lines of evidence. First, the mtDNA gene ND2 has been found to be consistent across diverse poikilothermic vertebrate lineages (such as fish, frogs, lizards, and salamanders ) and to have a substitution rate similar to cyt b . A mean ND2 substitution rate of 0.957%/MYR has been suggested in Costa Rican Eleutherodactylus inferred from Eurasian toads . Second, multiple calibrated salamander studies have found cyt b substitution rates in the range of 0.6 to 0.8%/MYR [90–92]. Frogs and salamanders share many characteristics that might be important in determining molecular rates of evolution (e.g., generation time, body size, ectothermy [92–94] and constancy among taxa has been shown in other mtDNA genes , so salamander substitution rate can be considered a reasonable approximation for frog substitution rate.
Population genetics analyses – mtDNA
We used a McDonald-Kreitman approach  to test for selection among clades and calculated Tajima's D  to test for neutrality in DNASP version 4.10.4 . Tajima's D is also used for estimating population expansion.
Cyt b haplotypic (gene) diversity, number of polymorphic sites, and nucleotide diversity  per locality were calculated in ARLEQUIN v. 2.000 . Inter-haplotype TN corrected p-distances were calculated in MEGA3 . Cyt b mismatch distributions and raggedness index (r)  were calculated by clade in ARLEQUIN under a model of sudden expansion using the parametric bootstrap approach (α = 0.05; 1000 bootstraps) . An empirical mismatch distribution that does not deviate from a unimodal distribution of pairwise differences among haplotypes and has a smooth distribution  suggests recent population expansion [101, 102]. The beginning of the population expansion can be estimated from τ, the crest of the mismatch distribution: τ = 2μt , where t is the generation time (unknown, but estimated as 1 year ) and μ is the upper estimate divergence rate (1.0%/MYR × number of bp). As an alternative test of population expansion, in ARLEQUIN we calculated Fu's F  and tested its significant with 1000 bootstrap replications. Significantly negative values of Fu's F suggest population expansion .
Population genetics analyses – microsatellites
Five microsatellite loci (Eloc-Bert&Ernie, Eloc-Beauty&Beast; Eloc-Thelma&Louise; Eloc-Laurel&Hardy and Eloc-Batman&Robin)  were amplified for all individuals in the Napo River area (localities 3 – 11 in Figure 1) for which we also have mtDNA sequences (n = 21 from Upland Clade; n = 76 from Lowland Clade). Samples were amplified and genotyped using published conditions . Amplification and scoring of a subset of samples was repeated to confirm genotypes. Four samples that resulted in large microsatellites in the loci Eloc-Thelma&Louise and Eloc-Batman&Robin were sequenced on an ABI 310 capillary sequencer using the microsatellite primers to determine exact composition of the microsatellite. A homozygous individual for locus Eloc-Laurel&Hardy was amplified by PCR and cloned using pGEM-T Vector System II kit (Promega) and the inserts sequenced in an ABI 3100 sequencer using M13 primers.
Vellend M, Geber MA: Connections between species diversity and genetic diversity. Ecol Lett. 2005, 8: 767-781. 10.1111/j.1461-0248.2005.00775.x.
Duellman WE: Patterns of Distribution of Amphibians: a global perspective. 1999, Baltimore: The Johns Hopkins University Press
Chek AA, Lougheed SC, Bogart JP, Boag PT: Perception and history: molecular phylogeny of a diverse group of neotropical frogs, the 30-chromosome Hyla (Anura: Hylidae). Mol Phylogenet Evol. 2001, 18: 370-385. 10.1006/mpev.2000.0889.
Lougheed SC, Gascon C, Jones DA, Bogart JP, Boag PT: Ridges and rivers: a test of competing hypotheses of Amazonian diversification using a dart-poison frog (Epipedobates femoralis). Proc R Soc B. 1999, 266: 1829-1835. 10.1098/rspb.1999.0853.
Symula R, Schulte R, Summers K: Molecular systematics and phylogeography of Amazonian posion frogs of the genus Dendrobates. Mol Phylogenet Evol. 2003, 26: 452-475. 10.1016/S1055-7903(02)00367-6.
Boul KE, Funk WC, Darst CR, Cannatella DC, Ryan MJ: Sexual selection drives speciation in an Amazonian frog. Proc R Soc B. 2007, 274: 399-406. 10.1098/rspb.2006.3736.
Funk WC, Caldwell JP, Peden CE, Padial JM, de la Riva I, Cannatella DC: Tests of biogeographic hypotheses for diversification in the Amazonian forest frog, Physalaemus petersi. Mol Phylogenet Evol. 2007, 44: 825-837. 10.1016/j.ympev.2007.01.012.
Noonan BP, Wray KP: Neotropical diversification: the effects of a complex history on diversity within the poison frog genus Dendrobates. J Biogeogr. 2006, 33: 1007-1020. 10.1111/j.1365-2699.2006.01483.x.
Lougheed SC, Austin JD, Bogart JP, Boag PT, Chek AA: Multi-character perspectives on the evolution of intraspecific differentiation in a neotropical hylid frog. BMC Evol Biol. 2006, 6: 23-10.1186/1471-2148-6-23.
Darwin C: The Origin of Species (1928 ed.). 1872, London: Dent
Duellman WE: Quaternary climatic-ecological fluctuations in the lowland tropics: frogs and forests. Biological Diversification in the Tropics. Edited by: Prance GT. 1982, New York: Columbia University Press, 389-402.
Avise JC: Phylogeography: the history and formation of species. 2000, Cambridge: Harvard University Press
Wallace AR: On the monkeys of the Amazon. Proc Zool Soc Lond. 1852, 20: 107-110.
Sick H: Rios e enchentes como obstäculo para a avifauna. Atas Simp Biota Amazonica (Zool). 1967, 5: 495-520.
Gascon C, Lougheed SC, Bogart JP: Genetic and morphological differentiation in Vanzolinius discodactylus : a direct test of the riverine barrier hypothesis. Biotropica. 1996, 28: 376-387. 10.2307/2389201.
Gascon C, Lougheed SC, Bogart JP: Patterns of genetic population differentiation in four species of Amazonian frogs: a test of the riverine barrier hypothesis. Biotropica. 1998, 30: 104-119. 10.1111/j.1744-7429.1998.tb00373.x.
Knapp S, Mallet J: Refuting refugia?. Science. 2003, 300: 71-72. 10.1126/science.1083007.
Fjeldså J: Geographical patterns for relict and young species of birds in Africa and South America and implications for conservation priorities. Biodiv Conserv. 1994, 3: 207-226. 10.1007/BF00055939.
Haffer J: Speciation in Amazonian forest birds. Science. 1969, 196: 131-137. 10.1126/science.165.3889.131.
Simpson BB, Haffer J: Speciation patterns in the Amazonian forest biota. Ann Rev Ecol Syst. 1978, 9: 497-518. 10.1146/annurev.es.09.110178.002433.
Haffer J: General aspects of the refuge theory. Biological Diversification in the Tropics. Edited by: Prance GT. 1982, New York: Columbia University Press, 6-24.
Tuomisto H, Ruokolainen K: The role of ecological knowledge in explaining biogeography and biodiversity in Amazonia. Biodiv Conserv. 1997, 6: 347-357. 10.1023/A:1018308623229.
Rohde K: Latitudinal gradients in species diversity: the search for the primary cause. Oikos. 1992, 65: 514-527. 10.2307/3545569.
Räsänen ME, Linna AM, Santos JCR, Negri FR: Late Miocene tidal deposits in the Amazonian foreland basin. Science. 1995, 269: 386-390. 10.1126/science.269.5222.386.
Gaston KJ, Blackburn TM: The tropics as a museum of biological diversity: an analysis of the New World avifauna. Proc R Soc B. 1996, 263: 63-68. 10.1098/rspb.1996.0011.
Hughes C, Eastwood R: Island radiation on a continental scale: Exceptional rates of plant diversification after uplift of the Andes. Proc Natl Acad Sci USA. 2006, 103: 10334-10339. 10.1073/pnas.0601928103.
Wake DB: Diversity of Costa Rican salamanders. Ecology and Evolution in the Tropics: A Herpetological Perspective. Edited by: Donnelly MA, Crother BI, Guyer C, Wake MH, White ME. 2005, Chicago: University of Chicago Press, 65-80.
McKenna DD, Farrell BD: Tropical forests are both evolutionary cradles and museums of leaf beetle diversity. Proc Natl Acad Sci USA. 2006, 103: 10947-10951. 10.1073/pnas.0602712103.
Veith M, Kosuch J, Vences M: Climatic oscillations triggered post-Messinian speciation of Western Palearctic brown frogs. Mol Phylogenet Evol. 2003, 26: 310-327. 10.1016/S1055-7903(02)00324-X.
Nei M: Molecular Evolutionary Genetics. 1987, New York: Columbia University Press
Elmer KR: Genetic diversity across spatial and evolutionary scales in some neotropical amphibians. PhD thesis. 2006, Queen's University, Department of Biology
Elmer KR, Dávila JA, Lougheed SC: Applying new inter-individual approaches to assess fine-scale population genetic diversity in a neotropical frog, Eleutherodactylus ockendeni. Heredity. doi:10.1038/sj.hdy.6801025,
Drummond AJ, Ho SYW, Phillips MJ, Rambaut A: Relaxed phylogenetics and dating with confidence. PLoS Biology. 2006, 4: e88-10.1371/journal.pbio.0040088.
Bromham L, Cardillo M: Testing the link between latitudinal gradient in species richness and rates of molecular evolution. J Evol Biol. 2003, 16: 200-207. 10.1046/j.1420-9101.2003.00526.x.
Crawford AJ: Relative rates of nucleotide substitution in frogs. J Mol Evol. 2003, 57: 636-641. 10.1007/s00239-003-2513-7.
Farias IP, Ortí G, Sampaio I, Schneider H, Meyer A: Mitochondrial DNA phylogeny of the family Cichlidae: Monophyly and fast molecular evolution of the neotropical assemblage. J Mol Evol. 1999, 48: 703-711. 10.1007/PL00006514.
Hoegg S, Vences M, Brinkmann H, Meyer A: Phylogeny and comparative substitution rates of frogs inferred from sequences of three nuclear genes. Mol Biol Evol. 2004, 21: 1188-1200. 10.1093/molbev/msh081.
Wright S, Keeling J, Gillman L: The road from Santa Rosalia: A faster tempo of evolution in tropical climates. Proc Natl Acad Sci USA. 2006, 103: 7718-7722. 10.1073/pnas.0510383103.
Aris-Brosou S, Excoffier L: The impact of population expansion and mutation rate heterogenity on DNA sequence polymorphism. Mol Biol Evol. 1996, 13: 494-504.
Coltorti M, Ollier CD: The significance of high planation surface in the Andes of Ecuador. Uplift, Erosion and Stability: perspectives on long-term landscape development. Edited by: Smith BJ, Whalley WB, Warke PA. 1999, London: The Geological Society of London, 239-253.
Hoorn C, Guerrero J, Sarmiento GA, Lorente MA: Andean tectonics as a cause for changing drainage patterns in Miocene northern South America. Geology. 1995, 23: 237-240. 10.1130/0091-7613(1995)023<0237:ATAACF>2.3.CO;2.
van der Hammen T: Paleoecology of tropical South America. Biological Diversification in the Tropics. Edited by: Prance GT. 1982, New York: Columbia University Press, 60-66.
Räsänen ME, Salo JS, Jungner H, Pittman LR: Evolution of the western Amazon lowland relief: impact of Andean foreland dynamics. Terra Nova. 1990, 2: 320-332. 10.1111/j.1365-3121.1990.tb00084.x.
Christophoul F, Baby P, Soula J-C, Rosero M, Burgos J: Les ensembles fluviatiles néogènes du bassin subandin d'Équateur et implications dynamiques. CR Geoscience. 2002, 334: 1029-1037. 10.1016/S1631-0713(02)01825-4.
Bès de Berc S, Soula JC, Baby P, Souris M, Christophoul F, Rosero J: Geomorphic evidence of active deformation and uplift in a modern continental wedge-top-foredeep transition: Example of the eastern Ecuadorian Andes. Tectonophysics. 2005, 399: 351-380. 10.1016/j.tecto.2004.12.030.
Pratt WT, Duque P, Ponce M: An autochthonous geological model for the eastern Andes of Ecuador. Tectonophysics. 2005, 399: 251-278. 10.1016/j.tecto.2004.12.025.
Räsänen ME, Salo JS, Kalliola RJ: Fluvial perturbance in the western Amazon Basin: regulation by long-term sub-Andean tectonics. Science. 1987, 238: 1398-1401. 10.1126/science.238.4832.1398.
Weng C, Bush MB, Athens JS: Holocene climate change and hydrarch succession in lowland Amazonian Ecuador. Rev Palaeobot Paly. 2002, 120: 73-90. 10.1016/S0034-6667(01)00148-8.
Bush MB, Silman MR, Urrego DH: 48,000 years of climate and forest change in a biodiversity hot spot. Science. 2004, 303: 827-829. 10.1126/science.1090795.
Heinicke MP, Duellman WD, Hedges SB: Major Caribbean and Central American frog faunas originated by ancient oceanic dispersal. Proc Natl Acad Sci USA. 2007, 104: 10092-10097. 10.1073/pnas.0611051104.
Graham CH, Ron SR, Santos JC, Schneider CJ, Moritz C: Integrating phylogenetics and environmental niche models to explore speciation mechanisms in dendrobatid frogs. Evolution Int J Org Evolution. 2004, 58 (8): 1781-1793.
Lynch JD, Duellman WE: The Eleutherodactylus of the Amazonian slopes of the Ecuadorian Andes (Anura: Leptodactylidae). Misc Publ Univ Kansas Mus Nat Hist. 1980, 69: 1-86.
Nevo E, Beiles A: Genetic diversity and ecological heterogeneity in amphibian evolution. Copeia. 1991, 565-592. 10.2307/1446386.
Haffer J: Alternative models of vertebrate speciation in Amazonia: an overview. Biodiv Conserv. 1997, 6: 451-476. 10.1023/A:1018320925954.
Padial JM, Gonzáles L, Reichle S, Aguayo R, de la Riva I: First records of five species of the genus Eleutherodactylus Duméril and Bibron, 1841 (Anura, Leptodactylidae) for Bolivia. Graellsia. 2004, 60: 167-174.
Lynch JD: A taxonomic and distributional synopsis of the Amazonian frogs of the genus Eleutherodactylus. Am Mus Nov. 1980, 2696: 1-24.
Hey J, Machado CA: The study of structured populations – new hope for a difficult and divided science. Nat Rev Genet. 2003, 4: 535-543. 10.1038/nrg1112.
Stuart BL, Inger RF, Voris HK: High level of cryptic species diversity revealed by sympatric lineages of Southeast Asian forest frogs. Biol Lett. 2006, 2: 470-474. 10.1098/rsbl.2006.0505.
Fouquet A, Vences M, Salducci M-D, Meyer A, Marty C, Blanc M, Gilles A: Revealing cryptic diversity using molecular phylogenetics and phylogeography in frogs of the Scinax ruber and Rhinella margaritifer species groups. Mol Phylogenet Evol. 2007, 43: 567-582. 10.1016/j.ympev.2006.12.006.
Ron SR, Santos JC, Cannatella DC: Phylogeny of the tungara frog genus Engystomops (= Physalaemus pustulosus species group; Anura: Leptodactylidae). Mol Phylogenet Evol. 2006, 39: 392-403. 10.1016/j.ympev.2005.11.022.
Pearman PB: Correlates of amphibian diversity in an altered landscape in Amazonian Ecuador. Conserv Biol. 1997, 11: 1211-1225. 10.1046/j.1523-1739.1997.96202.x.
Duellman WE: The biology of an equatorial herpetofauna in Amazonian Ecuador. 1978, Lawrence: University of Kansas
AmphibiaWeb: Information on amphibian biology and conservation. [http://amphibiaweb.org]
Sambrook J, Fritsch EF, Maniatis T: Molecular Cloning: a Laboratory Manual. 1989, Cold Spring Harbour: Cold Spring Harbor Laboratory Press, 3
Goebel AM, Donnelly JM, Atz ME: PCR primers and amplification methods for 12S ribosomal DNA, the control region, cytochrome oxidase I, and cytochrome b in bufonids and other frogs, and an overview of PCR primers which have amplified DNA in amphibians successfully. Mol Phylogenet Evol. 1999, 11: 163-199. 10.1006/mpev.1998.0538.
Hagelberg E: Mitochondrial DNA from ancient bones. Ancient DNA: Recovery and analysis of genetic material from paleontological, archaeological, museum, medical, and forensic specimens. Edited by: Herrman B, Hummel S. 1994, New York: Springer-Verlag Telos
Maddison DR, Maddison WP: MacClade 4: Analysis of phylogeny and character evolution. Version 4.07. 2003, Sunderland: Sinauer Associates
Thompson JD, Gibson TJ, Plewniak F, Jeanmougin F, Higgins DG: The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Res. 1997, 25: 4876-4882. 10.1093/nar/25.24.4876.
Posada D, Crandall K: MODELTEST: testing the model of DNA substitution. Bioinformatics. 1998, 14: 817-818. 10.1093/bioinformatics/14.9.817.
Swofford DL: PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods). version 4. 2003, Sunderland: Sinauer Associates
Templeton AR, Crandall KA, Sing CF: A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics. 1992, 132: 619-633.
Farris JS, Källersjo M, Kluge AG, Bult C: Testing significance of incongruence. Cladistics. 1994, 10: 315-320. 10.1111/j.1096-0031.1994.tb00181.x.
Lynch JD, Duellman WD: Frogs of the genus Eleutherodactylus (Leptodactylidae) in western Ecuador: systematics, ecology, and biogeography. Misc Publ Nat Hist Mus Univ Kansas. 1997, 23: 1-236.
Ronquist F, Huelsenbeck JP: MRBAYES 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003, 19: 1572-1574. 10.1093/bioinformatics/btg180.
AWTY: A system for graphical exploration of MCMC convergence in Bayesian phylogenetic inference Version 0.8. [http://ceb.csit.fsu.edu/awty]
Tamura K, Nei M: Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. Mol Biol Evol. 1993, 10: 512-526.
Kumar S, Tamura K, Nei M: MEGA3: Integrated software for molecular evolutionary genetics. Brief Bioinform. 2004, 5: 150-163. 10.1093/bib/5.2.150.
Arbogast BS, Edwards SV, Wakeley J, Beerli P, Slowinski JB: Estimating divergence times from molecular data on phylogenetic and popuation genetic timescales. Ann Rev Ecol Syst. 2002, 33: 707-740. 10.1146/annurev.ecolsys.33.010802.150500.
Hasewaga M, Kishino H, Yano T: Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J Mol Evol. 1985, 22: 160-174. 10.1007/BF02101694.
Nielsen R, Wakeley JW: Distinguishing migration from isolation: an MCMC approach. Genetics. 2001, 158: 885-896.
Brito P: The influence of Pleistocene glacial refugia on tawny owl genetic diversity and phylogeography in western Europe. Molec Ecol. 2005, 14: 3077-3094. 10.1111/j.1365-294X.2005.02663.x.
Nichols R: Gene trees and species trees are not the same. Trends Ecol Evol. 2001, 16: 358-364. 10.1016/S0169-5347(01)02203-0.
Ho SYW, Phillips MJ, Cooper A, Drummond AJ: Time dependency of molecular rate estimates and systematic overestimation of recent divergence times. Mol Biol Evol. 2005, 22: 1561-1568. 10.1093/molbev/msi145.
Emerson BC: Alarm bells for the molecular clock? No support for Ho et al.'s model of time dependent molecular rate estimates. Syst Biol. 2007, 56: 337-345. 10.1080/10635150701258795.
BEAST v1.4. [http://beast.bio.ed.ac.uk]
TRACER version 1.3. [http://evolve.zoo.ox.ac.uk]
Macey JR, Strasburg JL, Brisson JA, Vredenburg VT, Jennings M, Larson A: Molecular phylogenetics of western North American frogs of the Rana boylii species group. Mol Phylogenet Evol. 2001, 19: 131-143. 10.1006/mpev.2000.0908.
Austin JD, Lougheed SC, Moler P, Boag PT: Phylogenetics, zoogeography, and the role of vicariance and dispersal in the evolution of the Rana catesbeiana (Anura: Ranidae) species group. Biol J Linn Soc. 2003, 80: 601-624. 10.1111/j.1095-8312.2003.00259.x.
Crawford AJ: Huge populations and old species of Costa Rican and Panamanian dirt frogs inferred from mitochondrial and nuclear gene sequences. Molec Ecol. 2003, 12: 2525-2540. 10.1046/j.1365-294X.2003.01910.x.
Tan A-M, Wake DB: MtDNA phylogeography of the California Newt, Taricha torosa (Caudata, Salamandridae). Mol Phylogenet Evol. 1995, 4: 383-394. 10.1006/mpev.1995.1036.
Mueller RL: Evolutionary rates, divergence dates, and the performance of mitochondrial genes in Bayesian phylogenetic analysis. Syst Biol. 2006, 55: 289-300. 10.1080/10635150500541672.
Caccone A, Milinkovitch MC, Sbordoni V, Powell JR: Mitochondrial DNA rates and biogeography in European newts (genus Euproctus). Syst Biol. 1997, 46: 126-144. 10.2307/2413640.
Martin AP, Palumbi S: Body size, metabolic rate, generation time, and the molecular clock. Proc Natl Acad Sci USA. 1993, 90: 4087-4091. 10.1073/pnas.90.9.4087.
Estabrook GF, Smith GR, Dowling TE: Body mass and temperature influece rates of mitochondrial DNA evolution in North American fish. Evolution. 2007, 61: 1176-1187. 10.1111/j.1558-5646.2007.00089.x.
McDonald JH, Kreitman M: Adaptive protein evolution at the Adh locus in Drosophila. Nature. 1991, 351: 652-654. 10.1038/351652a0.
Tajima F: Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989, 123: 585-595.
Rozas J, Sánchez-DelBarrio JC, Messeguer X, Rozas R: DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics. 2003, 19: 2496-2497. 10.1093/bioinformatics/btg359.
Schneider S, Roessli D, Excoffier L: Arlequin ver. 2.000: A software for population genetics data analysis. 2000, Geneva: University of Geneva, Genetics and Biometry Laboratory
Harpending HC: Signature of ancient population growth in a low-resolution mitochondrial DNA mismatch distribution. Hum Biol. 1994, 66: 591-600.
Schneider S, Excoffier L: Estimation of past demographic parameters from the distribution of pairwise differences when the mutation rates vary among sites: Application to human mitochondrial DNA. Genetics. 1999, 152: 1079-1089.
Slatkin M, Hudson RR: Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations. Genetics. 1991, 129: 555-562.
Rogers AR, Harpending H: Population growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol. 1992, 9: 552-569.
Fu Y-X, Li W-H: Statistical test of neutality of mutations. Genetics. 1993, 133: 693-709.
Fu Y-X: Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997, 147: 915-925.
Elmer KR, Dávila JA, Lougheed SC: Isolation of simple and compound polymorphic tetranucleotide microsatellites for the neotropical leaflitter frog Eleutherodactylus ockendeni (Leptodactylidae). Molec Ecol Notes. 2006, 6: 891-893. 10.1111/j.1471-8286.2006.01389.x.
We thank F. Ayala Varela, T. Sugahara, I. G. Tapia, P. Menéndez, S. Padilla, M. Bustamante, G. Vigle, and P. Grefa for assistance collecting frogs in the field. Jatun Sacha Foundation, Arajuno Jungle Lodge, P. Grefa, As. Chonta Yacu, Cabañas Lidia, M. Tapia, La Selva Jungle Lodge, Kapawi Jungle Lodge, Municipio de Arosamena Tola, M. Terry, M. Foley and D. Duncan provided lodging and support in the field. Samples were collected with permits 006-IC-FAU-DBAP/MA, 016-IC-FAU-DNBAP/MA, and 010-IC-FAU-DFP from Ministerio del Ambiente, Ecuador and permission of the Federación Interprovincial de Nacionalidad Achuar del Ecuador at Kapawi locality. L.A. Coloma (QCAZ) identified specimens and advised on field planning. A.J. Baker (Royal Ontario Museum) hosted KRE during the writing of this manuscript. S. Mann designed Figure 1. Some of this work was carried out using the Computational Biology Service Unit resources from Cornell University that is partially funded by Microsoft Corporation. S. O'Reilly also provided computational resources. Funding for field and lab work was provided by NSERC postgraduate scholarships (KRE), a Queen's Doctoral Travel Grant (KRE), and an NSERC Discovery Grant (SCL).
KRE designed the study, collected samples, carried out the molecular studies, performed statistical analyses, and drafted the manuscript. JAD assisted with microsatellite primer design, genotyping, and microsatellite analyses, sequenced microsatellites, and interpreted results. SCL participated in the design of the study, performed statistical analyses, and helped draft the manuscript. All authors read and approved the final manuscript.