A dynamic continental moisture gradient drove Amazonian bird diversification.

The Amazon is the primary source of Neotropical diversity and a nexus for discussions on processes that drive biotic diversification. Biogeographers have focused on the roles of rivers and Pleistocene climate change in explaining high rates of speciation. We combine phylogeographic and niche-based paleodistributional projections for 23 upland terra firme forest bird lineages from across the Amazon to derive a new model of regional biological diversification. We found that climate-driven refugial dynamics interact with dynamic riverine barriers to produce a dominant pattern: Older lineages in the wetter western and northern parts of the Amazon gave rise to lineages in the drier southern and eastern parts. This climate/drainage basin evolution interaction links landscape dynamics with biotic diversification and explains the east-west diversity gradients across the Amazon.

INTRODUCTION

A deep understanding of biological diversification in the Amazon drainage basin (hereafter AB) is crucial to inferring processes generating such a rich biological diversity, as well as understanding future stability under global change (). A plethora of competing and overlapping hypotheses has been postulated (–). The relative roles of drainage evolution and rain forest expansion and contraction due to climate oscillations over the past 10 million years (Ma) have been the focus of intense debate (, , –).

Wide Amazonian rivers are considered to be notable barriers to gene flow for many organisms, as they mark range limits of many species (), pairs of sister lineages (), and areas of endemism (Fig. 1) (). However, these patterns contrast with evidence of gene flow across rivers via river meanders () and in headwater regions (). Other studies have revealed incongruences between modern river courses and genetic transitions between populations (–) or highly contrasting divergence times in different clades across particular rivers (). Amazonian river barriers are more dynamic than previously thought, as paleoriver barriers can explain incongruences between locations of modern rivers and genetic transitions (), and paleomorphological studies provide evidence of river flow inversions, drainage catchment capture events, and palaeochannels across the AB (, ).

Fig. 1

Amazon region.

Summary of the area of inquiry of this study, with major interfluvia (or endemism areas) associated with different geological provinces depicted in red (Guiana Shield), gray (Amazonian foreland basins), and green (Brazilian Shield) colors.

Alternatively, rivers may simply be zones of secondary contact between populations that differentiated in isolation in forest remnants (“refugia”) driven by climatic fluctuations (). This “refugium” hypothesis has been challenged by mismatches between the timing of splitting events, as measured from molecular phylogenies, and the extensive changes in precipitation regimes over tropical South America in response to glacial boundary conditions that would have led to refugium formation (), as well as by the absence of demographic signatures that would be expected under refugial scenarios (, ). Extensive palynological () and isotopic () evidence makes clear that the AB was not broadly covered by savannah vegetation during glacial maxima, as had been originally proposed (, ). Instead, forest composition was most likely subject to marked changes driven by a steep gradient in precipitation during glacial periods characterized by a more stable humid climate in the west in comparison to a very unstable and drier climate in the east (, , ). Furthermore, recent niche-based reconstructions of the geographic distributions of currently extant species in the Pleistocene under glacial and interglacial conditions highlight that refugia may be the result of inimical climate conditions, changes in forest structure and composition, or increased climate seasonality (, ).

Last, recent data have shown that pure vicariance models are not sufficient to explain diversification in Amazonia, as species-specific dispersal and environmental tolerance attributes can account for different responses to both refugial and riverine barriers (, ). More broadly, these and other recent phylogeographic studies have challenged the view that refugia and riverine models alone can account for the diversification of the Amazonian biota, posing serious hurdles on testing of alternative biogeographic scenarios (). While many patterns of diversification have emerged for Amazonian organisms over the past decades (), their ultimate and proximate drivers remain largely unknown, with current most debated hypotheses falling short from providing a basin-wide general mechanism of biotic evolution.

Here, we marshal comprehensive new collections of biotic material from across the AB to present phylogeographic data and paleodistributional projections that permit a detailed cross-species analysis. Using independent multilocus phylogeographic analyses and ecological niche model–based paleodistributional projections to glacial climate conditions, we analyze evolutionary patterns of lineage origination, dispersal, and extinction in the context of dynamics of suitable areas for each lineage. We present a spatiotemporal model of lineage formation and dispersal across the AB, based on multiple, independent suites of evidence, adding resolution to a “necessarily complex model” of diversification ().

RESULTS

This synthesis encompasses 23 codistributed, upland, humid, terra firme forest bird species, or species complexes (comprising geographically replacing allopatric or parapatric taxa) that are endemic or nearly endemic to the AB (table S1), 157 distinct lineages (species, subspecies, and populations), 1050 individuals sampled and sequenced, and 6527 primary occurrence records (table S1 and fig. S1).

Diversification and ancestral area reconstructions

A dominant pattern in our analyses was that the youngest splits corresponded to the diversification of the southeastern populations [Xingu and Belem interfluvia () located on the Brazilian shield] (Fig. 1). Because of the complexity of our data, we did multiple arrangements of our results, and all highlight this common main pattern (Fig. 2A, fig. S2, and table S2). Eleven of 19 species or species complexes present in southeastern Amazonia shared this pattern—Psophia spp., Galbula albirostris/Galbula cyanicollis, Malacoptila rufa, Thamnomanes caesius, Hylophylax naevius, Phlegopsis spp., Willisornis poecilonotus/Willisornis vidua, Dendrocincla merula, Dendrocolaptes certhia, Sclerurus caudacutus, and Sclerurus rufigularis. Divergence times associated with these southeastern lineages fell solidly in the middle and late Pleistocene [= 0.5 Ma ago (±0.3, SD)] (fig. S2). Moreover, in these 11 and 6 additional groups also occurring in southeastern Amazon (Megascops watsonii/Megascops usta, Cymbilaimus lineatus, Thamnophilus aethiops, Hypocnemis cantator complex, Xiphorhynchus guttatus, and Sclerurus macconnelli/Sclerurus peruvianus), a “counterclockwise” pattern of cladogenesis was evidenced (Fig. 2A, fig. S2, and table S2), in which splits associated with ancestral lineages took place primarily on the Guiana Shield (Imeri and Guiana interfluvia) or in western Amazonia (Napo, Inambari, and Rondonia interfluvia) before 2 Ma ago, with subsequent splits originating distinct lineages across northern and western-central Amazonia [= 1.08 Ma ago (±0.63, SD)], and into the Brazilian Shield in southeastern Amazonia [= 0.79 Ma ago (±0.55, SD)] (fig. S2). The oldest splits across the Brazilian Shield in southeastern Amazonia date to around = 1.67 Ma ago (±0.75, SD) (fig. S2). For all these taxa, ancestral area reconstructions based on geological and endemism areas support a pattern of origination or earlier vicariance in western and central Amazonian interfluvia, followed by subsequent dispersal, colonization, and vicariance in the southeastern interfluvia (Xingu and Belem; Fig. 2A and table S2).

Fig. 2

Summary spatiotemporal diversification pattern of Amazonian terra-firme bird species.

(A) Species tree and ancestral area reconstruction scenario exemplifying a counterclockwise pattern of cladogenesis observed from our analyses. Splitting dates and ancestral area reconstructions for each node were estimated on the basis of coalescent multilocus species trees obtained for each of the 21 species and species complexes showing this pattern of diversification (see table S2.3 for more details). (B) Lineage through time plots divided by nonpasserine and passerine birds, as well as one estimated for all lineages combined, based on our multilocus species tree reconstructions. Time is informed in million years.

The other two groups (approximately 10% of the lineages occurring in southeastern Amazonia; Rhegmatorhina spp. and Phoenicircus carnifex/Phoenicircus nigricollis) showed a distinct “clockwise” differentiation pattern, with ancestral area estimates not ruling out an origin or early presence and vicariance in southeastern Amazonia, followed by subsequent splits in the central and western parts of the AB (fig. S2 and table S2.1). The remaining four lineages analyzed (those that do not occur in southeastern Amazonia) also diversified extensively in northern and western Amazonia in the past 2 Ma, but have no extant populations in southeastern Amazonia, which reinforces the overall trend in which more than 90% of all 23 lineages sampled began to diversify outside of or away from southeastern Amazonia (fig. S2 and table S2.1). An overall variable rate of diversification was estimated for all lineages sampled throughout Amazonia (Fig. 2B and table S3), with a remarkable rate increase starting in the Plio-Pleistocene boundary, i.e., approximately 2.3 Ma ago.

Rivers as barriers: Disparate temporal splits across riverine barriers

The modern courses of main Amazonian rivers coincide with most (i.e., 75%) of the splitting events documented in our phylogenetic analyses (Fig. 3A). Times of diversification associated with putative riverine barriers fall within the Plio-Pleistocene, but there is substantial variation among taxa associated with the same putative barrier (Fig. 3B). For instance, the Solimões (upper Amazon) River could be related to multiple splitting events spanning the past 4 Ma (depicted in cyan in Fig. 3B).

Fig. 3

Modern barriers separating lineages of 23 Amazonian upland terra firme birds.

(A) Percentage of modern barrier types separating lineages of the species and species complexes sampled. (B) Detailed counts for diversification events through time coincident with the presence of modern Amazonian rivers (top graph), within interfluvial regions and across the Andes (bottom graph). Dotted bars represent SDs for time estimates obtained with *BEAST analyses. Note that older lineages were often putatively separated by western rivers.

Refugia: Indistinct demographic signatures between refugia and non-refugia

Present-day and paleodistributional estimates from group-specific ecological niche reconstructions show two regions of consistently reduced suitability across species under glacial conditions (putative non-refugia): (i) a northwest-southeast swath across the Guiana Shield, corresponding to the Jaú and Imeri interfluvia; and (ii) much of the southeastern interfluvia (Xingu and Belem; Fig. 4 and fig. S3). In contrast, the Guiana, Napo, and Rondonia interfluvia showed greater stability (putative refugia; Fig. 4 and fig. S3).

Fig. 4

Climate suitability predicted from past and current species distribution models.

Percentage of reduction in climate suitability during the Last Glacial Maximum (LGM), across the Amazon, averaged from all the 23 terra firme bird ecological niche models considering the (A) Community Climate System Model (CCSM) and the (B) Model for Interdisciplinary Research On Climate (MIROC) climatic conditions.

Demographic signatures of population expansion (77% of groups) are more frequent than stability (22% of groups), not only after the last glacial period but also earlier, throughout the Pleistocene (Fig. 5, fig. S4, and table S1). These demographic events were geographically scattered, i.e., neither expansion events were always coincident with putative non-refugia areas nor was stability mostly found within putative refugia. In other words, no consistent demographic pattern could be related to either putative or nonputative refugia (Fig. 5, fig. S4, and table S1).

Fig. 5

Demographic trends during the past 300 thousand years in putative refugia and non-refugia.

Top row: Putative refugia [(A) Guiana, (B) Napo, and (C) Rondonia]. Bottom row: putative non-refugia [(D) Jau/Imeri and (E) Xingu/Belem). Black dotted lines represent the expected demographic trends under refugia and non-refugia during glacial (0.16 to 0.21 Ma) and nonglacial periods. Straight lines depict median demographic trends inferred from coalescent-based Bayesian analyses for each sampled lineage within a given interfluve, and pie charts depict the percentage of each type of demographic event: black, stability; dark gray, expansion; and light gray, decline (see fig. S4 for more detailed information).

DISCUSSION

Understanding palaeoclimate changes and its influence on Amazonian habitats and biota (), identifying past and current river dynamics (), and describing former and present species diversity across all the AB (, ) are still unsettled, prolific lines of research. Consequently, pure model-based inference approaches might still be limited when studying Amazonian organism (), as the alternative hypotheses to be tested possibly will be insufficiently built given the complex, controversial, and still largely unknown biogeographic history of its biodiversity (, ). Instead, model-building approaches, based on multiple lines of evidence, are still needed to ground the biogeographical history of the AB (, , ). Although many previous treatments of this topic have suggested mechanisms that would apply to the entire AB (, ), our large-scale, multidimensional analyses reinforce that former explanations are overly simple. Rather, different interfluvia may have experienced markedly different climatic histories, and different groups with distinct ancestral distributions were affected by climatic fluctuations differently. The scenario posited here acknowledges these diverse processes () and incorporates the action of each in the diversification of Amazonian biotas.

First, we support that the pure vicariance riverine-barrier model of diversification is insufficient to explain Amazonian bird diversification. As demonstrated earlier (), the lack of temporal congruence among splitting events associated with single riverine barriers challenges a strict vicariant scenario, as envisioned originally (, ). Furthermore, we support that current riverine barriers may or may not have been the initial drivers of diversification, since river dynamics might explain temporal inconsistencies between splitting events and river formation. River course changes are well documented at least for the Madeira and Negro rivers (, ) and were previously inferred for the Japura and Tapajos rivers (, ). Similarly, paleoriver dynamics that do not correspond to current river configurations could also be related to diversification events in the Guiana, Imeri, Inambari, and Rondonia interfluvia (Fig. 3B) (, ). Hence, rather than a common and single sequence of vicariant events, our results support temporally dynamic rivers as having provided semipermeable dispersal barriers through time, which nonetheless played an important role in the generation of differentiation and maintenance of established genetic cohesion among populations (–, ). An overall rate increase of lineage diversification estimated to around 2.3 Ma ago (Fig. 2B) matches the proposed times for an observed peak in deposition of sediments of Andean origin in both western Amazonia () and the Amazon River mouth area (). Together, these biological and geological data support that the consolidation of the modern AB was a major driver of biotic diversification, despite the rivers’ relatively fast spatiotemporal dynamics ().

Second, our present-day and paleodistributional estimates from group-specific ecological niche reconstructions suggest that changes in climatic suitability associated with glacial conditions were not distributed evenly across the AB (Fig. 4 and fig. S3). However, the stable regions that we identify cannot be viewed as general refugia (). We observe demographic trends neither temporally coincident among each other nor geographically coincident with putative refugia [where one would expect demographic stability ()] or non-refugia [in which demographic expansions would be more frequent following glacial maxima ()] (Fig. 5 and fig. S4). Idiosyncratic responses to similar events might mask common processes acting on the local biota (), but these dissimilar demographic trends are enough evidence to refute the refugia hypothesis, at least, as originally suggested (). Conversely to former views of a homogenous Amazon forest (), the composite nature of this biome is currently largely accepted (our results) (, , ). This heterogeneity in the environment should be accounted as an enhancer to the species-specific evolutionary histories of Amazonian biota (). Furthermore, the contrasts in climate stability that we found were spatially consistent with the diversification patterns herein described, separating western Amazonia and the Guiana Shield and southwestern and southeastern Amazonia. Hence, we note close correspondence between climate-driven distribution dynamics and phylogeographic patterns in terra firme bird lineages. Previous studies have indicated that savannah-associated bird and plant taxa did not necessarily expand broadly during glacial periods, particularly at the Last Glacial Maximum (LGM) (), as had been posited originally (). As such, Pleistocene climate fluctuations likely did not produce massive forest-savannah turnover but rather changes in forest structure and composition (). Similar events are expected to have happened during other glacial periods (). These changes in habitat composition might have affected forest bird populations differently, depending on their ecological requirements and dependency on closed canopy rainforest ().

Last, our data corroborate that Amazon terra firme bird species are a young fauna (), with the oldest diversification events dated at less than 6 Ma ago (fig. S2 and table S2.1). Vicariance seemingly related to river formation has likely been driven largely by changing climatic suitability, with synergistic influences of drainage evolution and climate oscillations generating current patterns. Summarizing our evidence, a general scenario for the diversification of Amazon terra firme faunas is as follows (Fig. 6). During late Pliocene, ancestral taxa occupied areas more likely to have maintained more stable upland rain forested habitats in northern and western Amazonia (Figs. 2A and 6A, fig. S2, and table S2) (). The first splits separated lineages situated on opposite banks of the Amazon, Negro, and Madeira rivers (Figs. 2A and 6B, fig. S2, and table S2) (), which were among the first established following the entrenchment of Amazon River Basin during the Plio-Pleistocene border (, ). Rivers subdivided and structured populations as the AB was being formed (Figs. 3 and 6, B to D and F), but cyclical changes in climatic conditions were essential in reinforcing isolation events, sometimes resulting in speciation, others in population reductions and extirpations (our results; Fig. 6, C to F) (, , ). Diversification events for species within the Tepuis (mountains with sharp cliffs on the Guiana Shield) corroborate the relevant role of climate oscillation, and habitat composition changes within northeastern Amazon, throughout the late Pliocene and Pleistocene (Fig. 6, B and C) (). Recolonization, and river dynamics during the middle Pleistocene, then originated undifferentiated populations in adjacent southeastern areas (Fig. 6, D to F), or populations that have differentiated only recently.

Fig. 6

Model of regional biological diversification for Amazonian upland terra firme birds.

(A) By the late Pliocene, most ancestral lineages inhabited northern and western Amazonia. (B) During the Plio-Pleistocene boundary, most of the earliest splits occurred across the Guiana Shield, the western Amazonian foreland basins, and the westernmost part of the Brazilian Shield and involved lineages currently separated by the Amazon-Solimões and Negro rivers. Colonization of central southern Amazonia started during this time (represented by the blue arrows). (C) During the mid-Pleistocene, diversification events were most frequently observed in the western Amazonian foreland basins and the western part of the Brazilian Shield and involved lineages currently separated by the Solimões, Madeira, and Tapajos rivers, with inferred colonization toward southeastern Amazonia. (D) During the late Pleistocene, the center of diversification shifted to the Brazilian Shield, with lineages continuing to spread eastward and most splits taking place across the Madeira, Tapajos, and Xingu drainages. (E) However, because of river course dynamics particularly in the Madeira and Tapajos drainages, the process of diversification was reversed westward for some lineages. (F) Most recently, climate-driven retraction followed by recolonization allowed modern diversification events in southeasternmost Amazonia, with most splits involving lineages situated on opposite banks of the Xingu and Tocantins rivers.

All this resulted in older and richer faunas in western Amazonia, as well as less diverse and younger faunas in eastern sectors of the AB (). Our results corroborate that southeastern interfluvia, in particular, have not been continuously suitable for Amazonian terra firme forest bird lineages and so have either been (re)colonized recently or not at all (Fig. 2 and fig. S2), resulting in less diverse avifaunas in that region (). More generally, the rich biological diversity of the AB is the result of interacting suites of factors: climatic fluctuations both providing episodic isolation and causing local extinctions of populations, and rivers serving as dynamic and semipermeable barriers that limit dispersal and constrain gene flow and recolonization of areas.

MATERIALS AND METHODS

Study taxa and input data

Monophyletic avian species complexes, and species, endemic or with their distributions in largest part within the Amazon region, known to be strongly tied to upland terra firme humid forest, were included in this study on the basis of (i) availability of gene sequence data for at least one mitochondrial marker and at least one nuclear marker, and (ii) availability of sufficient data from across the geographic distribution of the complex/species in the AB [i.e., genetic and occurrence data from all main Amazonian geological provinces and major interfluvia (, , ) occupied by the group].

For each species or species complex analyzed, a concatenated phylogeny was used initially as a guide tree to further test lineage boundaries based on multilocus coalescent methods. Essentially, highly supported, reciprocally monophyletic groups (i.e., subspecies and populations within species or species within complex of species that tend to replace each other geographically) were assumed as hypothesized species (in the case of complexes of species) or independent evolutionary units (for species datasets) whose limits were then tested with multilocus coalescent methods (see below). The gene sequence data were drawn from both published and unpublished studies developed by the Aleixo laboratory group (see table S1 for a complete list of references, GenBank accession numbers, species/complex of species, lineages, and sample sizes).

Occurrence records were drawn from data associated with specimens in the ornithological collection of the Museu Paraense Emílio Goeldi (MPEG), as well as from scientific collections providing access to data via VertNet.org. Sample sizes ranged from 71 to 346 unique occurrence localities per species or complex of species (full dataset available at http://hdl.handle.net/1808/24034).

Climate data for the present day (1950–2000) were drawn from the WorldClim climate archive (). All occurrences were checked carefully for consistency and correct placement with respect to interfluvia and river positions. In particular, to represent overall tendencies and variation in temperature and precipitation, we explored the “bioclimatic” coverages in that data archive. Our concern was that we might end up calibrating niche models in an overly dimensional environmental space, so we plotted 1000 random points across the region of analysis (see next paragraph) and measured pairwise Pearson product-moment correlation coefficients for all variables. We eliminated one of each pair of variables that showed high (r ≥ 0.80) correlations (), except for the maximum temperature of the warmest month and the minimum temperature of the coldest month, wishing to retain both because previous analyses had indicated that these two variables (albeit correlated) are highly informative for AB bird distributions (). Hence, we retained the following for analysis: annual mean temperature, mean diurnal temperature range, isothermality, temperature seasonality, maximum temperature of the warmest month; minimum temperature of the coldest month, annual precipitation, precipitation of the driest month, and precipitation seasonality. All analyses were developed at a spatial resolution of 2.5′ (~5 km at the Equator), reflecting the approximate spatial accuracy of the georeferencing.

To summarize glacial climates, we used LGM (about 20 to 21 thousand years ago) climate datasets developed by R. J. Hijmans () to parallel the WorldClim current climate data, both at the same spatial resolution. Climate model projections are available only for LGM and not for any of the previous glacial maxima; but in climatic terms, if not in terms of accumulating geological effects, the different glacial maxima were relatively similar in terms of the conditions that they presented (). LGM data were drawn from the outputs of general circulation model (GCM) simulations from the Community Climate System Model (CCSM) () and the Model for Interdisciplinary Research On Climate (MIROC) (), downloaded from the website of the Paleoclimate Modelling Intercomparison Project 2 (www.pmip2.cnrs-gif.fr/). The GCM data had a native spatial resolution of 2.8° (approximately 300 km at the Equator), which were downscaled via the following procedure. First, the difference between the GCM output for historical and recent conditions was calculated. These differences were then interpolated to a 2.5′-resolution grid using the spline function in ArcInfo (ESRI, Redlands, CA) with the tension option. The interpolated differences were then added to the high-resolution current climate datasets from WorldClim. Last, established routines (www.worldclim.org/bioclim-aml) were used for generating the so-called bioclimatic datasets from raw monthly temperature and precipitation data. This procedure has the dual advantage of producing palaeoclimatic datasets at resolutions relevant to the spatial scale of analysis and of calibrating the simulated climate change data to the actual observed (present-day) climate data.

We required an approximation of the area that has been accessible to the species over relevant periods of time, which is the appropriate area for model calibration (). For lineage-wide analyses (i.e., covering most or all of the AB), we buffered the limits of the Amazon () by 400 km, in light of the broad extent of interdigitation between Amazon forest and Cerrado and Llanos vegetation types to the south and north, respectively. We delimited this area further to remove parts of the Andes presenting elevations ≥2500 m. All analyses were thus developed within this area but further delimited by the interfluvia (Fig. 1), in which each of the species/species complexes is known to occur.

Because broad spatial autocorrelation in environmental variables can lead to pseudoreplication of environmental signals and exacerbate overfitting in niche models (), we assessed spatial lags (i.e., a radius of nonindependence of points in environmental terms) across the AB. We calculated autocorrelation lags for each climate dimension separately, based on 12 bins of distances. Over the nine variables, lag distances ranged 0.8° to 4.2°, but annual mean temperature and annual precipitation were both below 1° (115 km at the Equator), such that we constrained occurrence in which that no points were closer to any other than 1° and developed five replicate subsamplings for each species.

Ecological niche modeling

We used Maxent () to calibrate ecological niche models for each species or lineage. In view of the difficult calibration challenges, we set the convergence threshold to 10−6. We generated 10 bootstrapped replicate models (50% subsampling) and extracted the median output across the 10 replicates for each. We then extracted the pixel-by-pixel median of the five replicate subsamplings, each representing the median of the 10 bootstrap replicates. Given emerging realization of contrasts between model information content and model transfer ability (), we did not use Akaike Information Criterion (AIC)–based model selection approaches.

Once present-day and the two LGM projections of niche models were in hand, we established appropriate thresholds for separating prediction of presence versus absence on the basis of the least training presence threshold approach (), but modified to take into account the expected error parameter E (). That is, instead of setting a threshold at the value that includes 100% of the training occurrence data, we used the threshold that includes (100–E)% of the training occurrence data, which is thus a more restricted area but takes into account the possible presence of noise in the occurrence data. We explored E values of 0, 5, and 10% to provide a variety of views of confidence in model projections. Thresholds were established on the basis of the present-day models and then applied to LGM maps.

To test general ideas regarding the stability of potential distributional areas through time, we subtracted the present-day data layer from the LGM projection for each species. We averaged these difference maps across all species. After removing all ecoregions () not corresponding to AB biomes, we also removed all areas termed either varzea (seasonally flooded forest) or seasonal forests. We then used zonal statistics approaches (ArcGIS 10.3) to calculate changes in suitable areas within each of the interfluvia across the region.

Phylogeographic analyses and time estimates

Multilocus coalescent Bayesian species trees were reconstructed with BEAST 1.8 () from the gene sequences for each species and species complexes from all Amazonian main geological provinces and interfluvia (, , ) occupied by the group (table S1). Coalescent times were estimated under a relaxed-clock approach with an uncorrelated log-normal distribution and calibrated using the mutation rate estimated for cytochrome b sequences [0.01105 substitutions/nucleotide/Ma ()], and a speciation Yule process as tree prior (). PartitionFinder 1.1.1 () was used to select the best-fitting model of nucleotide substitution. We ran at least two independent chains with a length of at least 108, and the initial 10% were excluded as burn-in. We verified effective sample sizes and convergence of the chains with Tracer 1.5 (). Trees were combined with LogCombiner (), and tree topology was assessed using TreeAnnotator () and FigTree 1.4 (http://tree.bio.ed.ac.uk/software/figtree/). For some species trees, an outgroup was used to achieve proper convergence of the chains: for the Malacoptila fusca/Malacoptila semicincta analysis, sequences from M. rufa were included, as well as a sequence of Hylophylax for the Phlegopsis nigromaculata/Phlegopsis erythroptera complex.

The possible common temporal-spatial splits across rivers, assumed by the riverine hypothesis (, ), were assessed correlating ingroup nodes with putative riverine barriers, except for Psophia, whose confidence intervals for split time estimates were quite wide (table S2.1). Graphs were obtained in R () using ggplot2 package ().

BEAST 1.8 () was also used to construct multilocus extended Bayesian skyline plots (EBSP) () for all major monophyletic lineages with enough sampling (n ≥ 3) and genetic variability (table S1). This method uses a Markov Chain Monte Carlo procedure to estimate the effective population size through time and also the credibility intervals considering both phylogenetic and coalescent uncertainty (). Because of the use of coalescent theory, small sample sizes might be used, particularly in multilocus analyses, such as this (). A strict clock model was used, but calibration, choice of the best-fitting mutation model, and convergence were assessed as aforementioned. In a few cases, an uncorrelated log-normal molecular clock was preferred, since better convergence was reached (as reported in table S1 and fig. S4). Analyses for most of the species trees and EBSPs were run in the CIPRES Science Gateway (). Common temporal-spatial demographic stability in refugia, and instability in non-refugia, as assumed by the refugia hypothesis (), was verified through plotting demographic trends for each putative area. Plots were constructed in R ().

Ancestral areas reconstructions

Two models of ancestral area reconstruction were tested using BioGeoBEARS implemented in R 3.1 (, ), one considering distributions among nine areas of endemism/interfluvia, following Borges and Silva (), and another considering four geological areas (, ) (Fig. 1). These two models were chosen because distribution of most of the reciprocally monophyletic lineages analyzed can be allocated to areas of endemism previously described for the Amazon region (table S1) (). However, in some cases, more than one lineage is present within the same area (e.g., H. naevius/Hylophylax naevioides complex) or a single lineage occurs in distinct areas (e.g., X. guttatus), occasionally with a distribution more coincident with the geological division of the AB (, ) (e.g., P. carnifex/P. nigricollis). Yet, because estimates based on the interfluvia include a much higher number of possible ancestral areas, they might be more informative (as also demonstrated for other BioGeoBEARS estimates) (, ).

For each model of areas, we performed a total of six different analyses including the Dispersal-Extinction Cladogenesis (DEC) model, a likelihood version of the Dispersal-Vicariance Analysis (DIVALIKE), and a version of the Bayesian inference of historical biogeography for discrete areas (BAYAREALIKE), as well as “+J” versions of these three models, which include founder-event speciation, an important process left out of most inference methods (). These biogeographic processes were implemented in a maximum likelihood (DEC, DEC+J, DIVALIKE, and DIVALIKE+J) or Bayesian (BAYAREALIKE and BAYAREALIKE+J) framework, as free parameters estimated from the data. Our multilocus species trees were used to infer the ancestral area probabilities to each lineage, which were computed for each node and subsequently plotted on the majority-rule chronograms. Outgroups were excluded from analyses to avoid insufficient representation of its diversity in our sampling, which could mislead biogeographic estimates. Last, we compared the six different models for statistical fit via comparison of Log-Likelihood (ln L) and AIC values. Dispersal (d), extinction (e), founder (J), ln L, and AIC values were directly obtained from BioGeoBEARS (results are summarized in table S2).

Lineage diversification analyses

We used the “APE,” “GEIGER,” and “LASER” packages (–) in R () to obtain the lineage diversification patterns for each species group. We fitted a total of eight diversification models implementing two constant rates (pure birth and birth-death) and three variable rates (exponential and logistic density-dependent and two-rate pure birth) using LASER 2.4 (). The SPVAR (time-varying speciation), EXVAR (time-varying extinction), and BOTHVAR (time-varying speciation and extinction) models that allow differential extinction and speciation rates were measured following Rabosky and Lovette (). The models fit were ranked by comparison of values of ln L, AIC, ΔAIC, and Akaike weight (ωi). We generated lineage-through-time (LTT) plots from our multilocus species trees to visualize the tempo of diversification within each species group. Because we generally used relatively small species trees with less than 10 lineages for each species/species complex to infer LTT plots, one concern is that the number of branches in each tree were not sufficient to detect patterns deviating from a constant rate of diversification. Despite this caveat, 15 of the 23 groups analyzed showed convergent patterns of rate constancy, which is notable considering the total number of lineages analyzed as well as their disparate phylogenetic affinities. Nonetheless, to minimize this potential bias, we further evaluated different models of diversification by combining all species and lineages sampled into a single LTT, following the same approach of the individual analyses (table S3).