|Home | About | Journals | Submit | Contact Us | Français|
The Neotropics contains half of remaining rainforests and Earth's largest reservoir of amphibian biodiversity. However, determinants of Neotropical biodiversity (i.e., vicariance, dispersals, extinctions, and radiations) earlier than the Quaternary are largely unstudied. Using a novel method of ancestral area reconstruction and relaxed Bayesian clock analyses, we reconstructed the biogeography of the poison frog clade (Dendrobatidae). We rejected an Amazonian center-of-origin in favor of a complex connectivity model expanding over the Neotropics. We inferred 14 dispersals into and 18 out of Amazonia to adjacent regions; the Andes were the major source of dispersals into Amazonia. We found three episodes of lineage dispersal with two interleaved periods of vicariant events between South and Central America. During the late Miocene, Amazonian, and Central American-Chocoan lineages significantly increased their diversity compared to the Andean and Guianan-Venezuelan-Brazilian Shield counterparts. Significant percentage of dendrobatid diversity in Amazonia and Chocó resulted from repeated immigrations, with radiations at <10.0 million years ago (MYA), rather than in situ diversification. In contrast, the Andes, Venezuelan Highlands, and Guiana Shield have undergone extended in situ diversification at near constant rate since the Oligocene. The effects of Miocene paleogeographic events on Neotropical diversification dynamics provided the framework under which Quaternary patterns of endemism evolved.
The Neotropics, which includes South and Central America, contains half of remaining rainforests and the largest reservoir of amphibian diversity. Why there are so many species in certain areas and how such diversity arose before the Quaternary (i.e., more that 1.8 million years ago [MYA]) are largely unstudied. One hypothesis is that the Amazon Basin was the key source of diversity, and species dispersed from there to other areas. Here, we reconstruct a time-calibrated phylogeny and track, in space and time, the distribution of the endemic and species-rich clade of poison frogs (Dendrobatidae) during the Cenozoic (more than 65 MYA) across the continental Neotropics. Our results indicate a far more complex pattern of lineage dispersals and radiations during the past 10 MY. Rather than the Amazon Basin being the center of origin, our results show that the diversity stemmed from repeated dispersals from adjacent areas, especially from the Andes. We also found a recurrent pattern of colonization of Central America from the Chocó at 4–5 MY earlier than the formation of the Panamanian Land Bridge at 1.5 MYA. Thus, the major patterns of dispersals and radiations in the Neotropics were already set by ~5–6 MYA (the Miocene–Pliocene boundary), but the ongoing process of Neotropical radiation is still happening now, especially in the Chocó–Central America region and Amazonian rainforest.
Tropical regions contain more than half of biological diversity on just 7% of the Earth's surface [1,2]. Differences in biodiversity between tropical and temperate regions have been attributed to contrasting speciation and extinction rates . Within the Neotropical realm, the Amazon Basin and the Chocoan region contain half of Earth's remaining rainforests and one of the largest reservoirs of terrestrial biodiversity. However, the impact of pre-Quaternary ecogeographic constraints on Neotropical biodiversity is largely unknown and the mechanisms contributing to species richness are unclear [3,4]. For example, the well-documented high α−diversity (species richness) of the flora and fauna of the Amazon rainforest  is usually attributed to local geoclimatic dynamics that promote monotonic accumulation of lineages [6,7]. However, the lower β−diversity (species turnover relative to distance) within the Amazon Basin is puzzling  and vastly underestimated. Current hypotheses are based on restricted, mostly Quaternary, spatiotemporal scales involving paleogeographic or ecological events (e.g., riverine barriers, Pleistocene climate change) , persistence of conservative niches , and analyses of phylogeography and endemicity . In addition to speciation/extinction processes , major paleogeological events promote diversification, yielding complex phylogenetic patterns of vicariance, dispersal, and secondary sympatry . Using phylogeographic analyses of the endemic and diverse clade of poison frogs (Dendrobatidae), we reconstructed Neotropical biogeography from the Oligocene to the present and revealed a widespread and highly dynamic pattern of multiple dispersals and radiations during the Miocene.
Major geoclimatic events have shaped the Neotropics. The most important include the isolation and reconnection of South America, the uplift of the Andes, the extensive floodbasin system in the Amazonian Miocene, the formation of Orinoco and Amazon drainages, and the dry−wet climate cycles of the Pliocene−Pleistocene (Figure 1). The Panamanian Land Bridge (PLB) between the Chocó and Central America, which formed progressively until the Pliocene , was an important biogeographic catalyst of dispersal and vicariance events at the Miocene−Pliocene boundary (e.g., Alpheus shrimps and freshwater teleost fishes) [12,13]. Similarly, the uplift of the Andes advanced the formation of the Amazon River, converting a widespread, northwest-flowing Miocene floodbasin into the current eastward-running Amazon Basin [14,15]. Two Miocene marine incursions into this wetland system isolated several aquatic taxa as living relicts, including the Amazon River dolphin, lineages of marine-derived teleosts and stingrays, and brackish water mollusks [16,17]. However, controversies exist about the magnitude and duration of these geoclimatic events .
Although well known for its megadiversity, no studies of the Neotropics have examined diversification patterns in highly speciose and widespread lineages over broad temporal and spatial scales. A general explanation that associates rates of speciation with paleogeographic events is lacking. Here, we test two general hypotheses about the spatial configuration of biogeographic areas on the origin of Neotropical diversity (Figure 1). First, under the center-of-origin hypothesis, lineages from the currently most diverse area (i.e., Amazon Basin) dispersed to other areas (Figure 1, HA: SM1). Second, under the stepping-stone hypothesis, paleogeographic events constrained the patterns of lineage diversification in the Neotropics among geographically adjacent areas (Figure 1, HA: SM2). Using a recently developed maximum likelihood (ML) procedure that estimates geographic range evolution, we tested both hypotheses against a null biogeographic model (Figure 1 H0: SM0) using a well-sampled Neotropical clade, the poison frogs (Dendrobatidae). We sampled 223 of the ~353 (264 described and 34–89 undescribed) species, distributed from Central America and Guiana Shield to southeast Brazil and from Andean páramos (4,000 meters above sea level [masl]) to lowland rainforests (<300 masl). However, ~40% of the species diversity remains unsampled (Table S11). Because the true diversity (i.e., described, undescribed, and extinct species) cannot be accurately assessed [19,20], macroevolutionary inference should account for missing diversity. Our goals are to (1) infer how geographic range evolution yielded current species distributions, (2) estimate the general patterns of speciation and extinction under necessarily incomplete taxon sampling, and (3) synthesize these findings with paleogeographic events to explain current patterns of species richness.
We reconstructed a ML phylogeny from an extensive taxon sample (Figure S3A–S3D), and a relaxed-clock Bayesian chronogram from a subsample of the previous . The phylogeny is in general agreement with previous studies [22,23], and the significance of branch support is provided in Figure 2 (for detailed support values see Figure 3A–3D). Although our taxon sampling is the most extensive for poison frogs (i.e., ~63.5% of the diversity in the family), we were concerned that incomplete taxon sampling might cause tree imbalance and introduce bias into divergence time estimation (Table 1) [24,25]. We assessed imbalance by testing the null model that each branch has an equal probability of splitting (Equal Rate Markov or ERM) using three tree imbalance tree indices: Colless IC, Sackin IS, and likelihood shape s (Table S1; see Materials and Methods). Overall, our results showed that the imbalance indices for all trees (i.e., ML, chronogram, and super-region phylogenies) have a tendency to depart from, but do not reject, the equal rates model (ERM). Therefore, our inferences about macroevolutionary events (e.g., ancestral area reconstruction, divergence times, and diversification rates) should not be affected significantly by incomplete taxon sampling.
The patterns of spatial and temporal distribution of poison frogs were inferred using dispersal-extinction-cladogenesis (DEC) modeling . We compared three DEC models (i.e., a null and two alternatives) for ten areas (Figures 1 and S6). The null model (SM0) assumes that spatial structure has no effect on biogeographic patterns of evolution. The alternative models either favor the Amazon Basin as a center-of-origin (SM1), or patterns of dispersal/vicariance that reflect the spatial arrangement or connectivity of biogeographic areas (a stepping-stone model; SM2). Our results strongly support the SM2 model over SM1 (large phylogeny Δln[L] = 140.5, p < 0.001; chronogram Δln[L] = 48.6, p < 0.001) and SM0 (large phylogeny Δln[L] = 17.7, p < 0.001; chronogram Δln[L] = 0.5NS, p > 0.05). The nonsignificant comparison of SM2 and SM0 for the chronogram alone is likely due to its reduced taxon sample.
The chronogram and a summary of the significant biogeographic events with confidence limits (Tables S8, S9A, and S9B) from the stepping-stone model are superimposed on the four major clades (Figure 2). The most recent ancestor of Dendrobatidae was distributed in regions that correspond to the current Venezuelan Highlands and Northern Oriental Andes at 40.9 ± 5.4 million years ago (MYA). This ancestral range split into a Venezuelan Highlands ancestor of clade A and an Andean ancestor (clade B + C + D). The most recent ancestor of each major clade occurred during the Oligocene at 34.9 ± 5.4 MYA (clade A), 30.9 ± 3.9 MYA (clade B), 27.1 ± 3.2 MYA (clade C), and 29.9 ± 4.0 MYA (clade D). We inferred 14 dispersals into and 18 from the Amazon Basin to adjacent areas, including three major radiations and a single lineage extinction within Amazonia. We also found five cross-Andean dispersals, five dispersals from Northern to Central Andes, six dispersals from Northern Andes to Chocó, four dispersals from Chocó to the Andes, and three temporal phases of lineage dispersal with two interleaved periods of vicariant events between the Chocó and Central America (Figure 2 and Table S12).
The diversity of Amazonian poison frogs (>70 species) resulted from 14 separate dispersals into this region, in three phases (Figure 2). First, the two oldest dispersals originated independently from the Guiana Shield (23.8 MYA) and from the developing Andes (21.1 MYA), just before and during the existence of the Amazonian Miocene floodbasin. Second, a single dispersal from the Guiana Shield occurred 15.5 MYA, during a low sea-level period associated with reduction of the Miocene floodbasin system. Third, the 11 remaining dispersals from the Andes took place between (1.6–7.3 MYA) during the formation of the modern Amazon Basin river system. Ancestral area reconstructions using a Bayesian multistate procedure similarly support the recent multiple dispersals to the Amazon Basin (Figure 2). Thus, our results suggest that much of the extant Amazonian biodiversity results from relatively recent immigration of distinct lineages followed by in situ radiation during the last 10 MYA.
At least 18 dispersals from the Amazon Basin to other areas were found in three temporal phases. First, the earliest dispersals to the developing Chocoan lowlands (21.8 MYA) and the Andes (15.2 MYA) occurred during the establishment of the Miocene floodbasin system. Interestingly, for the present Chocoan lineage of Dendrobates pumilio (Figure 2), our results suggest a Miocene overwater dispersal from Chocó to the developing Central America archipelago and the extinction of the Amazonian lineage ancestor at ~19.5 MYA during the formation of the Miocene floodbasin system. Second, dispersals to the Guiana Shield (1), the Venezuela Highlands (1), and the Andes (1) took place after the Miocene floodbasin system receded (8.8–10.8 MYA). Third, the 12 remaining dispersals were very recent (0.7–6.0 MYA), to the Guiana Shield (7), Andes (4), and Brazilian Shield (1). Thus, 16 out of 18 occurred <11 MYA, after the establishment of the current Amazonian geomorphology (Figure 2).
At least four major diversifications occurred within Amazonia: (1) the Allobates trilineatus complex (26 species) is the oldest (14.0–15.1 MYA); the three remaining are more recent (5.4–8.7 MYA), (2) the Amazonian Ameerega (19 species); (3) the Dendrobates ventrimaculatus complex (15 species); and (4) the Allobates femoralis complex (nine species). Moreover, all four lineages entered the Guiana Shield area in the Pliocene (Figure 2).
Species diversity in the Andes (71 species) resulted from a continuous diversification process since the late Eocene (Figure 2). However, several Andean events were contemporaneous with establishment of the Amazon Basin. Five cross-Andean dispersals from Oriental to Occidental Andes (2.0–25.4 MYA), five dispersals from Northern to Central Andes (14.9–30.9 MYA), and six dispersals from the Northern Andes to the Chocoan lowlands (8.3–29.9 MYA) (Figure 2) took place before the establishment of the Andes as a major geographic barrier during the Miocene–Pliocene boundary . We also found five dispersals from the Chocó to North and Central Andes (1.1–6.6 MYA) (Figure 2) that took place mostly in the Pliocene when the Andes were already established as a barrier. Our results indicate that lower montane transition zones between Andean and lowland environments (Chocó and Amazonia) promote diversification, as exemplified by the Amazonian Ameerega  and the Chocoan Epipedobates.
Central American and Chocoan species (>40) also show a complex pattern of diversification at the end of the Miocene. Ten dispersals from Chocó to Central America suggest a pattern of recurrent colonization and isolation in three phases (Figure 2). First, the two oldest dispersals (8.3–12.1 MYA) from the Chocó overlap with a proposed earlier exchange of faunas during the late Miocene . A single vicariant event at 6.8 MYA isolated the first wave of immigrants (i.e., ancestors of Phyllobates and Colostethus 1). Interestingly, the contemporaneous divergence of the Trinidad and Tobago species (Mannophryne trinitatis and M. olmonae) from Venezuelan relatives at 8.3 MYA suggests a global period of high sea level. Second, six Pliocene dispersals from South America (3.2–5.4 MYA), immediately followed a proposed low sea-level period after 6.0 MYA . Six vicariant events in the middle Pliocene (1.1–3.6 MYA) are concomitant with a second high sea-level period (1.5–3.0 MYA) that separated Central America from the Chocó . Third, two dispersals in the late Pleistocene (0.5–2.2 MYA) are contemporaneous with the Great Faunal Interchange at 1.2 MYA . Likewise, the endemic poison frog of Gorgona Island (Epipedobates boulengeri), located 50 km off the Pacific coast, was derived from a Chocoan ancestor 2.4 MYA during the same period. Our results strongly support the repeated dispersal of poison frogs into Central America across the PLB before its final Pliocene closure.
Similar results for the ancestral area reconstruction were obtained by dispersal-vicariance analysis (DIVA) . However, DIVA provided unrealistic ancestral reconstructions for basal nodes (Figure S7), and a large number (i.e., ~16 × 106) of equally parsimonious reconstructions (see Material and Methods). Therefore, DIVA analyses were considered exploratory due to its algorithmic limitations [26,31].
We estimated diversification rates of the chronogram (i.e., dendrobatid family clade) using the adjusted γ statistic to account for incomplete taxon sampling [32,33]. The γ statistic compares the relative position of the nodes in a chronogram to that expected under the pure birth model, and different values of γ characterize whether the diversification rate has increased (γ > 0), decreased (γ < 0), or remained constant (γ = 0) over time . However, we emphasize that the γ statistic is an indirect estimation of the rate of diversification , it should only be applied to monophyletic groups, and our inferences from the γ statistic should be considered relative measures of the diversification. Our result from the chronogram of the adjusted γ statistic is −0.662, and we failed to reject (p = 0.508NS) the null of the pure birth expectation of exponential growth, γ = 0. Simulations have suggested that significant positive values of γ have been associated with two alternatives, either increasing diversification or high extinction through time . The general absence of Tertiary frog fossils from the lowland Neotropics is intriguing , but it does not provide evidence for or against increases in extinction/diversification, as suggested by our estimated γ statistic for the dendrobatid family clade. However, our results of the ancestral area reconstructions strongly suggest that the bulk of recent diversification in poison frogs might be due to rapid radiations in the Amazon Basin and the Central American-Chocó super-regions in the late Miocene. Therefore, our inference of recent dispersals to Amazonia and the recent geological origin of the modern Neotropical rainforest  might weigh in favor of a recent increase in diversification. Additional data from other Neotropical biota might be crucial to validate our inferences.
We further evaluated our claim of a significant increase in diversification in the Amazon Basin and the Central America-Chocó super-regions. We tested for significant changes in diversification rate at the genus-supraspecific level (GSPF chronogram and Figure 3; see Material and Methods) under incomplete taxon sampling using ML approach  with two extreme values of the extinction/speciation ratio (i.e., extinction rate fraction α = μ/λ, α = 0, and α = 0.99) (Table 2). The GSPF chronogram rejected the constant-rate model (all lineages with equal diversification rate) in favor of a variable rate model (at least one lineage has a significant higher or lower diversification rate) (Table 2). Additionally, the GSPF chronogram favored the variable-rate model with diversification rate change in one or more lineages (Figure 3; Tables 3, S13 and S14) over an alternative of retained elevated ancestral diversification rate (Table 2). However, we found possible spurious significant rate increases (i.e., nodes 2 and 9 of Figure 3 and Table 2) that might be dependent on more inclusive lower nodes (i.e., 1 and 5, respectively). This “trickle-down effect” artifact can be explained by a significant rate increase detected in daughter clade being carried over to the adjacent parent clade.
The diversification rate increase within Ameerega is 3.23-fold (α = 0) to 7.55-fold (α = 0.99) higher than the rest of the GSPF chronogram (Figure 3; Tables 3, S13 and S14). Interestingly, Ameerega corresponds the most recent (i.e., 8.7 MYA) widespread radiation of poison frogs in the Amazon Basin after dispersal from the Andes (Figures 2 and and3)3) in the late Miocene. Other significant increases in the diversification rate include two in Amazonia, two in the Chocoan region, and one in the Andes (Figure 3 and Table 3). Significant decreases in the rate of diversification correspond to the rare Guiana Shield endemic Allobates (0.0008-fold reduction) and the mostly Andean endemic Clade B (0.4851-fold reduction) (Figure 3; Tables 3, S13, and S14). Therefore, we found that Amazonian and Central American-Chocoan lineages significantly increased their diversification rate since the late Miocene, while the diversification rate in the Andean and Guianan-Venezuelan-Brazilian Shield lineages have been near constant with a tendency to slow down since the Oligocene. However, we acknowledge that these super-regions (i.e., the Andes and the Guiana-Venezuela-Brazilian Shield) might be undersampled (Table 1) and conclusions about their near-constant rate of diversification need further validation.
The unstable coexistence of lineages within a large community for extended periods of time has been hypothesized as a cause of Neotropical diversity . However, our results suggest that such a model is incomplete; rather, the complex pattern of diversification is strongly intertwined with paleogeographic events. Our inferences about the past history of the poison frogs using ancestral area reconstructions and diversification analyses provide new insights on speciation and extinction patterns in the Neotropics. Three species richness patterns are potential explanations for the extant diversity differences among regions of the Neotropics: (1) high immigration into one region after suitable geoclimatic conditions are established; (2) gradual in situ diversification of old endemic clades, regardless of the geoclimatic conditions, promoting species accumulation; or (3) rapid in situ diversification of endemic clades after favorable geoclimatic conditions are established. We found that all three patterns might apply to different areas depending on historical context.
All extant Amazonian species descended from 14 lineages that dispersed into the Amazon Basin, mostly after the Miocene floodbasin system receded. The recurrent immigrations that originated mostly in the adjacent Andes (species-richness pattern 1), combined with an increased rate of diversification, explain the high α–diversity of Amazonia. Later, from the Miocene-Pliocene boundary to the present, a rapid in situ diversification (pattern 3) gave rise to the extant Amazonian endemic biota. Therefore, most species in Amazonia originated in the last 10 MY. Moreover, lineages immigrating into Amazonia at <8.0 MYA radiated rapidly, resulting in widespread species and young clades (e.g., Ameerega, Allobates trilineatus, A. femoralis, and Dendrobates ventrimaculatus groups).
The diversity in the Chocoan-Central American super-region derived from scattered immigrations (pattern 1) from Andes to the early Chocoan rainforest during the late Miocene. However, starting at the Miocene-Pliocene boundary, significant orogenic events gave rise to the Central American archipelago [11,37] followed by sea level fluctuations , which provided the conditions for repeated dispersal and vicariance events in pre-PLB islands. Evidence of rapid in situ diversification (pattern 3) is supported by the high genetic diversity observed among poison frogs and other lineages especially between Western and Eastern Panamá [12,39,40]. Interestingly, our results might explain the high β–diversity of other endemic clades within the Chocó-Central America super-region  as originating initially from long-distance dispersals between disconnected islands, with diversification later during isolation by high sea levels.
The Andes have undergone extended in situ diversification (pattern 2) since the late Eocene. However, our analyses also provided evidence of decline in the diversification rate since the middle Oligocene, which has important implications for history and conservation of the endemic Andean fauna. First, the Andes uplift at the Miocene–Pliocene boundary caused significant changes in the rate of diversification in the lowland transition zone. We found that several poison frog lineages distributed on one or both sides of the Andes had dispersed repeatedly before the Miocene uplift (i.e., five cross-Andean and five Northern to Central Andes migrations). Paleogeological evidence supports introgression of shallow seas across the northern Andes during the Miocene , suggesting a historical connection between the Amazon Basin and the Chocó. Second, the Pliocene Andean uplift (>2,000 m above sea level)  formed a significant barrier to dispersal, because no other cross-Andean dispersals were found. The uplift also was associated with dramatic ecological changes  and a decrease in diversification rates. These results suggest a role for niche conservatism [41,42], in that some lineages may have gone extinct because of failure to adapt. Alternatively, despite greater sampling effort in the Andes region than in other areas, we failed to find some previously common Andean species (e.g., Hyloxalus jacobuspetersi and the Ecuadorian H. lehmanni). Consequently, it is difficult to separate a natural decrease in diversification rates from the current trend of amphibian species extinctions at high altitudes due to anthropogenic habitat alteration , increased UV radiation , climate change , or pandemic infection . In contrast, the montane transition zones of the Andes and adjacent lowlands (Chocó and Amazonia) have become centers of rapid cladogenesis (pattern 3), and species richness in these transition zones might be underestimated because many Neotropical lineages have been shown to contain several cryptic species . Therefore, dispersals within or across the Andes diminished during the Pliocene, but diversification has intensified in the Andes-lowlands interface.
Although some of the oldest lineages of poison frogs originated in the Guiana Shield and the Venezuelan Highlands (>30 species), our results suggest extended in situ diversification (pattern 2) followed by a decline in the rate of diversification of endemic clades in both areas since the early Miocene. Along the same lines, the Guiana Shield has high poison frog endemism, which is mostly restricted to the summits of the sandstone tepuis , while recent Amazonian poison frog immigrants occupy lowlands adjacent to the tepuis. Our results suggest that the decline of endemic Guianan diversity might be associated with ecological changes in habitat due to the collapse of the ancient tepuis  and repeated dispersals from Amazonian lineages since the Pliocene. However, the diversity of poison frogs in the Guiana Shield is only beginning to be revealed . In contrast, diversification in the Venezuelan region most likely reflects the oldest vicariant event in Dendrobatidae, at 40.9 MYA. The costal ranges of Mérida, Cordillera de la Costa, and Paria peninsula are species rich but their total area is less than 5% of that of the Amazon Basin. No lineage of this endemic fauna has dispersed out to other regions since the early radiation of the poison frog family in the late Eocene. However, Eocene floristic paleoecological reconstruction of the Venezuelan Highlands area showed that it was more diverse than at present , suggesting that the ancestral habitat of the first poison frogs might have been lowland. The depauperate dendrobatid fauna of the Venezuelan llanos and Brazilian Shield plateau is puzzling, but might be related to Holocene aridity .
The recurring dispersals to Amazonia suggests that a large part of dendrobatid diversity results from repeated immigration waves at <10.0 MYA, followed by a rapid in situ diversification after geoclimatic conditions suitable for a rainforest ecosystem were present. The biota of Amazonia was not isolated during the process of diversification, but finely intertwined with the development and export of biodiversity across the entire Neotropical realm. Poison frog diversity in the Chocoan-Central American super-region was significantly associated with formation of the PLB in the Pliocene. Repeated dispersals between disconnected islands followed vicariance by cyclic high sea-level periods, promoted rapid in situ diversification and endemism of poison frog lineages. The extant Andean, Guianan, and Venezuelan Highlands fauna most likely originated after prolongated in situ diversification since the origin of the poison frog clade, but the pace of species formation within these areas has slowed down. Phylogenetic analyses on tropical biota such as birds  and the species-rich genus Inga  as well as models of diversification , have argued that the Amazon might accumulate older lineages; however, the origin of those lineages is not clear. Our results are the first to provide evidence, to our knowledge, of the major involvement of the Andes as a source of diversity of both the Amazon and the Chocó–Central America region. Because 23.5% of endemic Amazonian amphibian species are dendrobatids (i.e., ~70 of 298) , our results may generalize to other Neotropical terrestrial biota with similar distribution. Moreover, these results provide a crucial broad spatiotemporal framework that, coupled with realistic phylogeny-based explanations of the current richness in Neotropics, explains why species occur where they do and how they came to get there. Thus, the major patterns of dispersals and radiations in the Neotropics were already set by the Miocene–Pliocene boundary, but the ongoing process of Neotropical radiation is occurring now, in the Chocó–Central America region and especially in the Amazonian rainforest.
Because there are no fossil poison frogs, a large phylogeny of amphibians was constructed to calibrate the age of the root of the dendrobatid tree. We used a total of 89 terminals including 80 species of anurans (30 families), three species of salamanders (three families), three species of caecilians (three families), and three outgroups (lungfish, human, and chicken) (Figures S1 and S2; Table S2). The amphibian classification partially follows that of Frost et al. . Conflicts with Frost et al.  are indicated as paraphyletic families (e.g., Dicroglossidae 1 and 2).
Molecular data include the mitochondrial rRNA genes (12S and 16S sequences; ~2,400 bp) and the nuclear protein-coding gene RAG-1 (approximately 495 bp). Sequences were retrieved from GenBank (74 terminals) or sequenced (15 terminals) from total genomic DNA (Table S2). The primers and protocols for amplification, purification, and sequencing of PCR products are provided in previous studies [22,55,56]. PCR products were sequenced in both directions and compared to GenBank sequences using BLAST (http://www.ncbi.nlm.nih.gov/BLAST/). By this procedure, we were also able to validate sequences in GenBank and exclude contaminated or mislabeled submissions. GenBank accession numbers are given in Table S2.
A total of 406 individuals for described (137) and undescribed (34–89) species of poison frogs and 12 outgroups (from Hyloidea) were used to estimate the phylogeny (Figure S3A–S3D; Table S3). The estimate of undescribed species (34–89 spp.) corresponds to the estimated minimum and maximum number of new species. Therefore, the described diversity of poison frogs (264 species) plus the diversity discovered by our analysis yields a maximum of 353 (264 known + 89 maximum number of undescribed species), which is a better estimate of the true extant diversity. Of the 127 described poison frog species not sampled (Table S11), we were able to identify their closest relative or species group in 92.1% of the cases (117 species). Thus, we are confidant that we have not missed any crucial lineage and that our conclusions will hold as more data are incorporated. Furthermore, the conservation status of the unsampled species is based on the Global Amphibian Assessment  is as follows (Table S11): 51.3% are data deficient, 16.6% have been described recently (no category), 28.9% are in one of four “threatened” categories, and 3.2% are of least concern.
The classification here partially follows that of Grant et al. . Species placements that conflict with this taxonomy are indicated as paraphyletic genera (e.g., Colostethus 1 and 2). Proposed taxonomic changes and corrections are explained in “Corrections to Taxonomy” (Text S1). Our taxon sample included species from throughout the distribution of dendrobatid frogs, with regions of higher diversity sampled more extensively (Table S3) and the 117 species that were not sampled were assigned to a taxonomic group (Tables S11 and S13) and included in the Materials and Methods section “Speciation and extinction patterns under incomplete taxon sampling.” Molecular data were generated using the same protocols indicated in the first paragraph of Materials and Methods. We included only the mitochondrial rRNA genes (12S and 16S sequences; ~2,400 bp), from which 374 individuals (121 species) have complete sequences. GenBank sequence accession numbers are given in Table S3; because some outgroup sequences will appear in other papers currently under review, their numbers are not listed. Although sequences from other genes are available in GenBank, these were not used because data from the data are highly incomplete; complete data for only 80 species were available. Moreover, simulation studies have indicated that large phylogenetic analyses including many terminals with incomplete sequences might bias branch length estimation, and cause topological inconsistencies due to unrealistic estimates of the rate of evolution .
Contigs were assembled using Sequencher 4.7 . Sequences were initially aligned using ClustalX 1.81  and manually adjusted to minimize informative sites using MacClade 4.08 . The final matrices included 2,688 characters (i.e., 2,192 from mitochondrial genes, 495 from RAG-1) for the Amphibian Tree Matrix (ATM) and 2,380 characters (mitochondrial genes) for the Poison Frog Tree Matrix (PFTM). A total of 228 (ATM) and 113 (PGTM) ambiguously aligned characters were excluded. For the ATM, we divided the data into mitochondrial gene and RAG-1 partitions. For the PFTM, we used 12S rRNA, tRNA-Valine, and 16S rRNA partitions. The best model of nucleotide substitution for each partition was determined using ModelTest 3.7 [62,63]. For the ATM, GTR + Γ + I and TrN + Γ + I were selected for the mitochondrial genes and RAG-1 segments, respectively. For the PFTM, GTR + Γ + I was selected for the 12S and 16S rRNA segments and the GTR + Γ for the tRNA-valine segment.
ATM and PFTM were analyzed with ML methods under a genetic algorithm in GARLI 0.951  and with Bayesian sampling of tree space with MrBayes 3.1.2 [65,66]. For the ATM analyses, the amphibian species were constrained to be monophyletic. For the GARLI analyses, a total of 40 independent runs were used to infer the best tree and 200 nonparametric bootstrap searches were used to estimate support for the nodes. For the MrBayes analyses, tree topology estimation, branch lengths, and Bayesian posterior probabilities (PP) were determined from five independent runs of four incrementally heated chains. Runs were performed for 35–45 × PFTM 106 generations under partitioned models using default settings as priors; the sampling frequency was 1 in 1,000 generations. The convergence of the runs and the optimal burn-in was determined to be 1.238 × 106 (ATM) and 4.282 × 106 (PFTM) generations using MrConverge . This program estimates the point where the likelihood score becomes stationary and the overall precision of the bipartition posterior probability is maximized.
For molecular dating analyses the strict molecular clock model was rejected from both ATM (χ2 = 4,762.6, df = 87, p < 0.001) and PFTM (χ2 = 3822.6, df = 404, p < 0.001) datasets using a likelihood ratio test (LRT) that compared the best unconstrained GARLI trees to those estimated under a strict molecular clock. Therefore, a relaxed molecular clock Bayesian method in MULTIDIVTIME  was used to estimate chronograms for the Amphibian Tree and the Poison Frog Tree.
For chronogram estimation, all taxa in the ATM dataset (Amphibian Chronogram) and a pruned PFTM dataset (Poison Frog Chronogram) were used. The pruned PFTM dataset excluded multiple individuals of the same species to improve computational efficiency. For each analysis (ATM and pruned PFTM), the aligned matrix and the rooted ML topologies without branch lengths were input into MULTIDIVTIME . Branch length estimates under the F84 + Γ model of molecular evolution and variance/covariance matrices were calculated using the BASEML and ESTBRANCHES components of PAML 3.15  and MULTIDIVTIME , respectively. Calibration points (Figures S4 and S5; Tables S4–S6), relaxed-clock model priors, and variance/covariance matrices were then input into MULTIDIVTIME .
For the Amphibian Chronogram, three different sets of time constraints were used to assess the robustness of the dating estimates; these were based on paleogeography, vertebrate fossils, and amphibian fossil records (Table S4). The ingroup tip-to-root distances needed for the estimation of the MULTIDIVTIME rtrate and rtratesd priors  were calculated using TreeStat v1.1 . The relaxed-clock model priors were 344 MYA for the expected age between tip and root (rttm) and 20 MYA for its standard deviation (rttmsd). The rttm prior corresponds to the divergence of Amniota and Amphibia and is based on fossil and molecular analyses [69–72]. The expected molecular evolution rate at the ingroup root node (rtrate) prior and its standard deviation (rtratesd) were estimated at 0.00345 substitutions/site/MY by dividing the median of ingroup tip-to-root distances by the rttm prior as suggested by the MULTIDIVTIME documentation . The priors for the expected value of the Brownian motion constant υ (nu) (brownmean) and its standard deviation (brownsd) were estimated to be 0.058 by setting rttm * brownmean equal to 2.0 (on-line suggestion of Frank Rutschmann ). The bigtime parameter was set to twice the estimated time divergence of Amniota and Amphibia (i.e., 700 MYA). Markov chain (newk, othk, thek) and beta (minab) priors were set to default values. Each MCMC chain was run for 1 × 106 generations with sampling frequency of 1 per 100 generations and burn-in of the first 100,000 generations. All analyses were run twice to ensure convergence of the time estimates. The divergence time estimated for each node of the amphibian chronogram was described by its mean age and 95% confidence interval (CI) (Table S7).
The Poison Frog Chronogram was estimated from the pruned PFTM dataset that included 224 individuals from 157 poison frog species (131 described and 26 undescribed) and 12 outgroups. The fossil record of Tertiary Neotropical frogs is minimal and no fossils of poison frogs have been found . For this reason, we used a three-part strategy to date the Poison Frog Chronogram. First, the expected ages and 95% CIs of the split of the Dendrobatidae from other Hyloidea and the age of the most recent common ancestor of this clade were estimated from the Amphibian Chronogram. Second, a list of geological time constraints (Table S5) was developed based on paleogeological evidence (Table S6). Third, ancestral area reconstruction was inferred and explained in the Materials and Methods section “Ancestral area reconstruction”. To test the overall accuracy of these approaches, three sets of progressively less inclusive time constraints were used (Table S5). The relaxed-clock model priors were 71.4 MYA for the expected age between tip and root (rttm) and 18.6 MYA for its standard deviation (rttmsd). This value, estimated from the Amphibian Chronogram, corresponds to the mean age of the split between Dendrobatidae and its sister clade. The expected value (rtrate) and standard deviation (rtratesd) priors were set to 0.0056 substitutions/site/MY. The priors for the expected value (brownmean) and its standard deviation (brownsd) of the Brownian motion constant υ (nu) were set to 0.2. These priors (rtrate, rtratesd, brownmean, and brownsd) were obtained similarly as in the Amphibian Chronogram (see above). The bigtime parameter was set to three times the estimated age of the Dendrobatidae node from the Amphibian Chronogram (40 × 3 = 120 MYA). Markov chain priors (newk, othk, thek), beta prior (minab), and MCMC chain parameters were the same as for the Amphibian Chronogram. All analyses were run twice to ensure convergence of the time estimates. The divergence time estimated for each node of the Poison Frog Chronogram was described by its mean age and 95% CI (Table S8) and major taxonomic events (Table S9).
We assessed the robustness of the calibrations (Table S5) with three approaches. First, we recalculated the chronogram by using penalized likelihood approach (PL)  implemented in r8s . Because the penalized likelihood method requires at least one fixed node age, the nodes were fixed for root of the poison frog tree at 42.9 MYA or for the major split between clades B + C + D and A at 36.3 MYA. These values correspond to the mean age for each node obtained by averaging the ages from the three estimates of the Amphibian Chronogram (Table S7). We calculated the Poison Frog Chronogram with all constraints and then by removing one calibration constraint at a time (Tables S9 and S10). We use five different smoothing parameters (1, 10, 100, 500, and 1,000), an additive scale for rate penalty, and 20 random starts to estimate the chronogram. The average of all runs without the respective constraint was used as the estimate of node age (Table S9B). Second, we recalculated the chronogram using the same parameters of the relaxed-clock Bayesian method but removing one calibration constraint at a time (Table S9A). Finally, we recalculated the Poison Frog Chronogram to determine the effect of the combination of bigtime and rttmsd priors, which have been shown to have the substantial impact , by changing the estimates of nodes to make them older or younger. We used the same rttmsd prior (18.6 MYA) and two bigtime prior values of 60 and 80 (Table S9A). Based on all tests, the estimates of 36 exemplar nodes that correspond to major phylogenetic events in the family (Table S9A and S9B) and the differences in MYA from the Poison Frog Chronogram node mean are provided (Table S10). We found that all time estimates from each test were within two standard deviations of the mean of poison frog chronogram (Table S10). The time estimates obtained by using the penalized likelihood method mostly were within one standard deviation from those estimated in the Poison Frog Chronogram. In the case of the tests done using the relaxed molecular clock Bayesian method, the removal of constraint A (2.4–15.0 MYA, which was assigned to nodes that correspond to dispersals between South America and Central America) was found to provide older time estimates for all nodes, especially those close to the root. The effect of changing the rttmsd/bigtime priors combination was negligible.
The reconstruction of ancestral areas of the poison frog clade was determined by three methods: a maximum-likelihood inference of geographic range evolution [26,31,77], DIVA , and Bayesian analysis of ancestral areas [78–80]. Ten areas (designated by letters) were delimited on the basis of geological barriers, areas of endemism [81–83], and distribution maps  (Figure 1; Table S3). The Andes were divided into four adjacent regions: Central Occidental Andes (H), Central Oriental Andes (G), North Oriental Andes (E), and North Occidental Andes (F). Northern and Central Andes are divided perpendicular to the southern limit of Carnegie Ridge (parallel 2 °S) in southern Ecuador [84,85]. The Oriental and Occidental Andes are divided along the Interandean valleys that separate both parallel mountain chains. The Guiana Shield (B) and Brazilian Shield (K) regions are located in the eastern shoulder and center of South America, respectively. Both regions are ancient Precambrian plateaus with lowland tropical rainforest, dry forest, and cloud forests from 200–2,800 m. The Venezuelan Highlands (D) region includes the current Caribbean costal cordilleras of Venezuela (Mérida, Cordillera de la Costa, and Paria peninsula) and Trinidad and Tobago Islands. Paleontological and stratigraphical evidence suggests that Venezuelan Highlands region had strong similarities to the western Amazonian Tertiary fossil fauna [86–89]. However, the Venezuelan Highlands biogeographic distinctiveness is evidenced by the Miocene uplift , episodic Miocene floodings, and the formation of the Llanos , the separation from the Guiana Shield region by the current Orinoco river drainage and from the northern Andes by the Táchira depression. The Amazon Basin (C) region includes the river drainage and its extensive lowland tropical rainforest <300 m. The Central America (J) region corresponds to the lowlands and highlands of western side of the PLB to southern Nicaragua (northernmost distribution of poison frogs). The Chocó (I) region includes the eastern side of the PLB and the costal lowland tropical forest below 500 m of the Pacific Coast of Ecuador and Colombia on the western side of Andes. The Magdalena river drainage and the Gorgona Island were also included in the Chocó region based on the paleontological and biotic resemblance to the Chocó [92–94].
Ancestral area reconstructions, time of diversification, and rate of diversification estimates require a fully bifurcate tree [95–98]. We used the best GARLI tree to reconstruct the ancestral areas. Each species in the phylogeny was assigned to one or more regions based on distribution (Table S3). For the first method, we estimated a DEC model using Lagrange package [26,77]. The DEC model is a continuous-time stochastic model for geographic range evolution in discrete areas, with ML parameters estimated for rates of dispersal between areas (range expansion) and local extinction within areas (range contraction). The DEC model considers geographic scenarios of lineage divergence (including scenarios involving within-area speciation), allowing a widespread ancestral range to persist through a cladogenesis event as ancestral states at internal nodes on a phylogeny with observed species ranges at the tips. In all cases, ancestral ranges were assumed to include no more than two areas, the maximum observed for extant species. Moreover, spatial and temporal constraints (e.g., area distances, dates of geological origin) may be imposed in the DEC model estimation, providing a more accurate estimation of the ancestral areas and hypothesis testing of specific geographic scenarios.
We tested three biogeographic scenarios based on the hypothesized origin of the group (Figure S6). First, the SM0 null model has all areas as equiprobable ancestral ranges and assigns them to be adjacent pairs (i.e., separated by one step in the matrix), plus the individual areas themselves. Second, the SM1 (center-of-origin model) favors an Amazon Basin origin of the poison frog clade and assigns all non-Amazonian areas to be adjacent to the Amazon Basin (i.e., separated by one step in the matrix) and nonadjacent to each other (i.e., separated by two steps). However, the one exception was to make the Chocó and Central America adjacent, because it is not possible to reach Central America without passing through the Chocó. Third, the SM2 model (stepping-stone model) assigns the rates of dispersal between areas to be inversely scaled by their relative distance and connectivity (e.g., the distance between Guiana Shield and Central America is four steps, so the rate of dispersal was constrained to be 0.25). Each analysis estimated the global rate of dispersal and local extinction on the phylogeny with species ranges at the tips, considering all possible range inheritance scenarios (ancestral states) at internal nodes without conditioning on any particular values. Then, the global rates were used to calculate and rank the likelihoods of all ancestral states at every internal node on the phylogeny. The ranking is done by calculating the likelihood at the root of the tree, given the global rates, with one node fixed at one ancestral range scenario, and all other nodes free to vary. The ML scores for all nodes are compared to the overall ML score of each geographic scenario under a “global” method of ML ancestral state reconstruction . The likelihood of each model was optimized against the observed ranges of species and their phylogenetic relationships, with differences greater than 2 log-likelihood units considered significant; the reconstruction with the worst score was rejected . Tests of the three hypotheses (SM0, SM1, and SM2) were repeated in both the complete ML phylogeny and chronogram (reduced number of taxa). All tests were repeated in both large ML phylogeny and chronogram including and excluding Allobates alagoanus whose phylogenetic position at the base of clade A is not well supported (Figure S3A). In all cases similar results were found, so the results excluding A. alagoanus are presented.
For the second method, we performed ancestral area reconstruction by DIVA using DIVA 1.1 . DIVA assumes that speciation occurs as a consequence of vicariance and reconstructs these events at no cost. Therefore, the ancestral geographic distributions are reconstructed by minimizing the number of dispersal or extinction events necessary to explain the actual distribution pattern. Because we assigned all nine regions to favor dispersals equally, no single solution for the large tree was found (~16 × 106 possible reconstructions). Hence, this analysis is considered exploratory due to its limitations. DIVA, which is parsimony-based, only optimizes historical events on cladograms without regard to relative or absolute time (i.e., branch lengths), and is less flexible about geographic lineage divergence scenarios (e.g., widespread ancestors are assumed to undergo vicariance); specific biogeographic hypotheses cannot be tested. In consequence, the strict consensus of the dispersal/vicariance events on >75% of all possible reconstructions was used as the best solution and mapped on the ML poison frog phylogeny (Figure S7).
For the third method, estimates of the values of states at ancestral nodes were derived by point estimates (log-likelihoods) using an ML approach in BAYESMULTISTATE, a component of BAYESTRAITS [79,100]. The reduced chronogram (236 terminals) described in the Materials and Methods section “Determination of divergence times” was used for the analysis. We coded all terminal distributions under two alternatives, Amazonian Basin (C) origin, and non-Amazonian Basin origin (all other areas). We calculated the proportion of likelihood of both alternatives at each node with 10,000 samples under default parameters. The results were mapped onto the chronogram (Figure 2). We also tested whether an Amazonian origin was present at each node by constraining the node to this state using the option “fossilizing” in BAYESMULTISTATE . The likelihoods of the constrained and unconstrained reconstructions were compared, and a difference of ≥2 log-likelihood units was considered significant; the model with the worse score was rejected .
We calculated the tree imbalance of the ML poison frog phylogeny, the reduced-taxon chronogram, and the tree of each super-regions: the Amazon Basin (region C), the Andes (regions E–H), Chocó-Central America (regions I and J), and Guiana Shield-Venezuela-Brazilian Shield (regions B, D, and K). We used conservative imbalance metrics [25,101]: IC , IS , and s [104,105]. All standardized indices and probability of rejecting the null model of each branch having the equal probability of splitting (Equal Markov Rate or ERM) were calculated using functions colless.test, sackin.test, and likelihood.test, implemented in the R package apTreeshape .
An indirect estimate of diversification rates assuming incomplete taxon sampling was explored using the γ statistic  and adjusted for actual phylogenies by excluding the distance between the most recent node to the present (i.e., gn node of a phylogeny of size n), which does not come from the same distribution . The adjusted γ statistic value was obtained from the Poison Frog Chronogram (family level).
The adjusted γ statistic [32,33] was calculated using following functions of the R package LASER : (1) the gn node was excluded from the chronogram (i.e., entire family and major region subtrees) by using truncateTree function; (2) the γ statistic of the truncated tree was calculated using gamStat function; (3) a null-distribution of γ under a pure birth model was obtained by 10,000 Monte Carlo replicates using mccrTest.Rd. This function generates full size trees at the species level and prunes randomly terminals to the actual size of the empirical tree (i.e., simulates incomplete taxon sampling); (4) the empirical γ statistic was adjusted by subtracting the mean value of the simulated null-distribution of γ, which is expected to be 0 ; and (5) the p-values of the adjusted γ statistics were computed from a normal distribution; values outside the ±1.96 standard deviation boundaries are significantly different (alpha = 0.05) from the null pure birth expectation.
We tested for significant changes in diversification rate within the poison frog clade using a ML methodology under the assumption of incomplete taxon sampling . First, we produced a GSPF level chronogram by pruning all but a single lineage per taxonomic group (Figures 3 and S8; Table S13) from the Poison Frog Chronogram, yielding a tree of 78 terminals. The total species richness per taxonomic group was assigned to each terminal based on previous taxonomic and phylogenetic studies [22,23,108–112]. Second, we estimated a constant diversification rate r (i.e., the difference between speciation λ and extinction μ rates) across the phylogeny using a ML estimator that incorporates both known taxonomic diversity and phylogenetic data . We calculated the constant-rate model fit statistics (log likelihood and AIC score) and r using the fitNDR_1rate.Rd function of LASER . Third, we tested for shifts in diversification rate within the poison frog phylogeny by comparing likelihood of the GSPF chronogram under constant and rate-flexible diversification models . Two alternative hypotheses for rate-flexible model may explain the shifts in rate of diversification: (1) an increase within a particular clade (rCL) from the ancestral diversification rate r (flexible-rate model) or (2) a clade-specific decrease (rCL) from the ancestral diversification rate r (rate-decrease model) . We calculated both rate-flexible alternative model fit statistics, r and rCL values using the fitNDR_2rate.Rd function of LASER . The best fitting model was determined using a likelihood ratio test (LRT) between the constant-rate and the flexible-rate models (nested), and by ΔAIC scores between the flexible-rate and rate-decrease models (not nested). All analyses were performed under two extremes of the relative extinction rate (α = μ/λ, α = 0 and α = 0.99) as a fixed parameter to determine the robustness of the results to variation in the extinction fraction .
The phylogram was produced under a genetic algorithm in GARLI 0.951 . We also inferred the tree topology and branch lengths using a Bayesian sampling of tree space with MrBayes 3.1.2 [65,66]. Support values from the nodes were constructed with 200 nonparametric bootstrap replicates (above branch) and Bayesian posterior probabilities (below branch). An asterisk indicates a support value of 100.
(5.20 MB JPG)
(a) Relaxed clock chronogram estimated from the ML tree using time constraints (date from seven fossils and one paleogeographic date constraints as in Table S4); support values for the nodes were constructed with 200 bootstrap replicates (ML) using GARLI; Bayesian posterior probabilities (PP) were estimated using MrBayes; gray bars are 95% CI of the estimated node age. Six major geological events are indicated by red squares: (1) rifting of Pangaea; (2) Gondwana break-up; (3) separation of South America from Africa; (4) Cretaceous–Tertiary mass extinction; (5) initial orogeny of the Andes; (6) formation of the modern Amazon Basin; (7) Stem Dendrobatidae; and (8) Dendrobatidae.
(4.53 MB JPG)
(A–D) The phylogram is the ML methods under a genetic algorithm in GARLI 0.951 . We also inferred the tree topology and branch lengths using a Bayesian sampling of tree space with MrBayes 3.1.2 [65,66]. Support values from the nodes were constructed with 200 nonparametric bootstrap replicates (above branch) and Bayesian posterior probabilities (below branch). (*) indicates that the support value was 100. The coding of the individual terminals includes the generic, species, locality, country code, biogeographic area, and museum, field series, or GenBank accession numbers (see Table S3).
(2.55 MB PDF)
See Table S4.
(4.38 MB JPG)
(A–D) The nodes preceded by “N” correspond to the ML phylogeny and those that were preceded by “#” correspond to the chronogram. See Table S5.
(2.36 MB PDF)
(1) SM0, the null model, assumes no spatial structure, with equal rates of dispersal among all areas and no constraints on ancestral range composition; (2) alternative center of origin model SM1 models the dispersal rate into and out of the Amazonia as twice the rate between all other areas, with widespread ancestral ranges constrained to include the Amazon Basin; and (3) a stepping-stone model SM2 accounts for area adjacency, scaling rates of dispersal between areas inversely by relative distance, and constraining widespread ancestors to spatially adjacent areas. Each hypothesis test was repeated in the large ML phylogeny and in the chronogram. The area codes correspond to Guiana Shield (B), Amazon Basin (C), Venezuela Highlands and Trinidad and Tobago Islands (D), Northern Andean Cordillera Oriental (E), Northern Andean Cordillera Occidental (F), Central Andean Cordillera Oriental (G), Central Andean Cordillera Occidental (H), Chocó, Magdalena Valley, and Gorgona Island (I), Central America (J), and Brazilian Shield (K).
(814 KB JPG)
(7.57 MB JPG)
(5.04 MB JPG)
An asterisk indicates significance at α = 0.05.
(49 KB PDF)
(82 KB PDF)
“Use” identifies which samples were used for the phylogeny (P), chronogram (C), and Langrage (L) analyses. Other columns provide GenBank accession numbers, locality of collection, and biogeographic region of the Neotropics.
(305 KB PDF)
(94 KB PDF)
(107 KB PDF)
(111 KB PDF)
(74 KB PDF)
(124 KB PDF)
(A) corresponds to the node age estimations using MULTIDIVTIME constraint jackknifing and “bigtime” prior variation.
(142 KB PDF)
(72 KB PDF)
(196 KB PDF)
(65 KB PDF)
The species sampled for the analyses are indicated by “*”.
(62 KB PDF)
Node numbers correspond to Figure S8.
(94 KB PDF)
For the loan of tissue samples we thank C. Aguilar (Museo de Historia Nacional San Marcos, Lima); I. de la Riva and I. Rey (Museo Nacional de Ciencias Naturales, Madrid); R. Ibañez (Smithsonian Tropical Research Institute); L. Trueb, W.E. Duellman, and J. Simmons (University of Kansas); S.B. Hedges (Pennsylvania State University); C. Cicero (Museum of Vertebrate Zoology, University of California at Berkeley); and D.L. Dittmann (Louisiana State University Museum of Natural Science). JCS thanks N. Biani, R.F. Guerrero, J. Gómez, J.M. Guayasamin, I. Tapia, E. Tapia, F. Ayala, C. Aguilar, and R. Schulte for field assistance and collection of specimens. JCS thanks the Smithsonian Tropical Research Institute and A. Amézquita (Universidad de los Andes, Bogotá) for field guidance and lab space. LAC thanks E.R. Wild and M. Read for field assistance. JCS and DCC thank N. Biani, J. Brown, B. Caudle, C. Guarnizo, T. LaDuc, A. Lemmon, E. Lemmon, G. Pauly, D. Scantlebury, B. Symula, and D. Zwickl for suggestions about data analysis; B. Moore for insights into diversification analyses and the “trickle down effect”; and J. Thorne for suggestions about MULTIDIVTIME. For research and collection permits we thank the Dirección Nacional de Áreas Protegidas y Vida Silvestre (Panamá; permit SE/A-4–6); Ministerio de Ambiente (Ecuador; permits 0004-IC-FAU-DNBAP/MA, 0006-IC-FAU-DNBAP/MA, 016-IC-FAU-DNBAP/MA, 027-IC-FAU-DNBAPVS/MA); Instituto Nacional de Recursos Naturales-INRENA (Perú; permits 061-2003-INRENA-IFFS-DCB, 002765-AG-INRENA, 003-2005-INRENA-IFFS-DCB, and Convention on International Trade in Endangered Species of Wild Fauna and Flora [CITES] permit number 4326); and the Conselho Nacional de Desenvolvimento Científico e Tecnológico (Brazil; CNPq, Portaria MCT number 170, de 28/09/94) and the Instituto Brasileiro do Meio Ambiente e dos Recursos Naturais Renováveis (Brazil) (IBAMA, number 073/94-DIFAS), under a research agreement between the Oklahoma Museum of Natural History and the Museu Paraense Emilio Goeldi. We also acknowledge the comments from the editor and three anonymous reviewers.
Author contributions. JCS conceived, designed, and performed the experiments. JCS, RR, and DCC analyzed the data. JCS, RR, and DCC contributed analysis tools. JCS, LAC, KS, JPC, RR, and DCC wrote the paper.
Funding. The National Science Foundation provided support to KS (DEB-0134191), JPC (DEB-9200779 and DEB-9505518 to L.J. Vitt and JPC), DCC (EF-0334952), and JCS (DDIG DEB-0710033 under Institutional Animal Care and Use Committee [IACUC] protocol 05111001). RR acknowledges feedback from working groups at the National Evolutionary Synthesis Center (NESCent, NSF EF-0423641). LAC was supported by funds provided by Unidad de Cooperación para el Desarrollo de los Pueblos (UCODEP) (AIDCO/B7–6200/0380/01/TF), and the Dirección General Académica of the Pontificia Universidad Católica del Ecuador. JCS was also supported by University of Texas Ecology, Evolution and Behavior (EEB) Graduate Research Fellowships.
Competing interests. The authors have declared that no competing interests exist.