Research Article

Modeling the ecology and evolution of biodiversity: Biogeographical cradles, museums, and graves

See allHide authors and affiliations

Science  20 Jul 2018:
Vol. 361, Issue 6399, eaar5452
DOI: 10.1126/science.aar5452

Simulating South American biodiversity

The emergence, distribution, and extinction of species are driven by interacting factors—spatial, temporal, physical, and biotic. Rangel et al. simulated the past 800,000 years of evolution in South America, incorporating these factors into a spatially explicit dynamic model to explore the geographical generation of diversity. Their simulations, based on a paleoclimate model on a 5° latitude-longitude scale, result in shifting maps of speciation, persistence, and extinction (or cradles, museums, and graves). The simulations culminate in a striking resemblance to contemporary distribution patterns across the continent for birds, mammals, and plants—despite having no target patterns and no empirical data parameterizing them.

Science, this issue p. eaar5452

Structured Abstract


Individual processes that shape geographical patterns of biodiversity are increasingly understood, but their complex interactions on broad spatial and temporal scales remain beyond the reach of analytical models and traditional experiments. To meet this challenge, we built a spatially explicit, mechanistic model that simulates the history of life on the South American continent, driven by modeled climates of the past 800,000 years. Operating at the level of geographical ranges of populations, our simulations implemented adaptation, geographical range shifts, range fragmentation, speciation, long-distance dispersal, competition between species, and extinction. Only four parameters were required to control these processes (dispersal distance, evolutionary rate, time for speciation, and intensity of competition). To assess the effects of topographic heterogeneity, we experimentally smoothed the climate maps in some treatments.


The simulations had no target patterns. Instead, the study took a fundamental approach, relying on the realism of the modeled ecological and evolutionary processes, theoretical derivations of parameter values, and the climatic and topographic drivers to produce meaningful biogeographical patterns. The model encompassed only the Late Quaternary (last 800,000 years), with its repeated glacial-interglacial cycles, beginning at a time when South America was already populated with a rich biota, comprising many distinct lineages. Nonetheless, current consensus holds that the contemporary flora and vertebrate fauna of South America include numerous lineages that have undergone rapid diversification during the Quaternary, particularly in the Andes. In our model, over the course of each simulation, a complete phylogeny emerged from a single founding species. On the basis of the full historical records for each species range, at each 500-year interval, we recorded spatial and temporal patterns of speciation (“cradles”), persistence (“museums”), extinction (“graves”), and species richness.


Simulated historical patterns of species richness, as recorded by maps of the richness of persistent (museum) species, proved quite successful in capturing the broad features of maps of contemporary species richness for birds, mammals, and plants. Factorial experiments varying parameter settings and initial conditions revealed the relative impact of the evolutionary and ecological processes that we modeled, as expressed in spatial and temporal patterns of cradles, museums, graves, and species richness. These patterns were most sensitive to the geographical location of the founding species and to the rate of evolutionary adaptation. Experimental topographic smoothing confirmed a crucial role for climate heterogeneity in the diversification of clades, especially in the Andes. Analyses of temporal patterns of speciation (cradles) and extinction (graves) emerging from the simulations implicated Quaternary glacial-interglacial cycles as drivers of both diversification and extinction on a continental scale.


Our biogeographical simulations were constructed from the bottom up, integrating mechanistic models of key ecological and evolutionary processes, following well-supported, widely accepted explanations for how these processes work in nature. Despite being entirely undirected by any target pattern of real-world species richness and covering only a tiny slice of the past, surprisingly realistic continental and regional patterns of species richness emerged from the model. Our simulations confirm a powerful role for adaptive niche evolution, in the context of diversification and extinction driven by topography and climate.

Observed species richness versus modeled (simulated) richness.

Upper map: Contemporary South American bird richness (2967 species). Lower map: Simulated spatial pattern for the cumulative richness of persistent (museum) species, arising from the model. The map show results averaged over all parameter values for an Atlantic Rainforest founder, excluding the climate-smoothing experimental treatments. Simulated species richness is highly correlated with observed species richness for birds (r2 = 0.6337).


Individual processes shaping geographical patterns of biodiversity are increasingly understood, but their complex interactions on broad spatial and temporal scales remain beyond the reach of analytical models and traditional experiments. To meet this challenge, we built a spatially explicit, mechanistic simulation model implementing adaptation, range shifts, fragmentation, speciation, dispersal, competition, and extinction, driven by modeled climates of the past 800,000 years in South America. Experimental topographic smoothing confirmed the impact of climate heterogeneity on diversification. The simulations identified regions and episodes of speciation (cradles), persistence (museums), and extinction (graves). Although the simulations had no target pattern and were not parameterized with empirical data, emerging richness maps closely resembled contemporary maps for major taxa, confirming powerful roles for evolution and diversification driven by topography and climate.

Despite continually improving documentation of the global distribution of biodiversity and increasing awareness of its vulnerability, we remain confronted by our ignorance of the fundamental ecological and evolutionary processes that have shaped the diversity and complex biogeography of continental biotas (13). Narrative accounts (4) and correlative studies (510) suggest underlying causes, and theoretical models demonstrate possible mechanisms (7, 1117), but spatially and temporally explicit, process-based models (18, 19), founded on a comprehensive suite of well-studied, widely accepted mechanisms, have the greatest potential to assess the complex and sometimes indeterminate interactions among underlying processes (2025). Here, we offer such a comprehensive model, for a simulated biota. We applied it to a fine-scale topographical representation of South America—the most climatically and biologically diverse continent on Earth—driven by a spatially explicit paleoclimate model for the past 800 thousand years (ka), for both temperature and precipitation.

In a changing climate, the geography of species distributions is governed by many interacting environmental and biological processes. These processes include the shifting spatial pattern of environmental variables (16, 26), range shifts (27), dispersal (28), the geographical effects of competition between species (29), niche evolution (30), range fragmentation and re-joining (31, 32), speciation (3335), and extinction (36). Our biogeographical simulation model (Fig. 1) incorporated all these processes at the level of geographical ranges of populations, as realistically as feasible, given the inevitable computational limitations. Our principal objective was to evaluate, experimentally, the relative importance of these mechanisms in a multifactorial framework.

Fig. 1 Simulation model structure.

The processes and parameters implemented in the model, all illustrated here, link climate dynamics and topography to emerging biodiversity patterns. Key entities and patterns (tables S3 and S4) appear in rectangles at the population, species, and assemblage levels. Processes are shown in ovals. Control knobs (table S1) represent the four model parameters: Dmax, maximum dispersal distance; Hmax, maximum niche evolutionary rate; Tmin, minimum time for speciation; and Cmax, maximum intensity of competition allowing coexistence, estimated as a function of phylogenetic distance. Climate change, on a realistic topographical template, drives ecological and evolutionary processes, interacting with each population’s environmental niche to determine range dynamics. Dispersal promotes range shift and range expansion. Interactions between climate change, niche, and geographic distribution may result in adaptive niche evolution, range fragmentation, or extinction. Fragments that remain isolated long enough become new species. Closely related species, in sympatry, may coexist or undergo competitive exclusion. Starting from a single, founding species (and its initial climatic niche), the simulation produces temporal and spatial patterns of biodiversity, including times and places of speciation (cradles), extinction (graves), and persistence (museums). See the Methods section, below, and “Model specification: process sequence” in (95).

Current understanding of South American biogeography

The crucial role of the Andes

The rise of the Andes, beginning 25 million years ago (37), launched a biogeographical experiment unique in Earth’s history (38, 39)—the juxtaposition of a long, trans-tropical mountain chain and a tropical rainforest (40). Throughout their history, the environmental heterogeneity of the Andes is thought to have driven species diversification by (i) providing novel, high-altitude mountain environments; (ii) erecting dispersal barriers that promoted vicariant speciation, both between east and west slopes (41, 42) and between internal valleys and peaks along the mountain chain (43, 44); (iii) offering a north-south, climatically driven, biogeographical corridor; (iv) sheltering species threatened with extinction by reducing regional climate velocity (45, 46); and (v) offering refugia from climatic extremes (4, 4749).

Our model (Fig. 1) encompasses all of these drivers of Andean diversification. Separately from our assessment of the relative importance of these processes, we investigated the role of Andean climate heterogeneity, itself, as a driver of diversification, within the experimental design. To do so, we simulated the biogeographical consequences of gradually smoothing the topography of South America, with the expectation that these diversification processes would be progressively eliminated.

Historical biome dynamics

Although biogeographers unequivocally view the Andes as a driver of species diversification (3841), historical linkages among South American biomes are still under debate. The present-day northeast-southwest Caatinga-Cerrado-Chaco “hot-dry diagonal” poses a dispersal barrier between Amazonian and Atlantic rainforests for vertebrates and plants (5052). However, multiple cases of disjunct distributions across this barrier (5356) support Por’s (57) proposal of an ephemeral connection between the Amazon and Atlantic Rainforest during late Quaternary climate cycles (5860). Seasonally dry tropical forests have also been viewed as important drivers of plant diversity in South America, and they offer a potential explanation for disjunct distributions of woody plants between Atlantic Forest and the Amazon and Andes (50, 61).

Within the Amazon, recent empirical and model-based studies have suggested the existence of a large-scale dipole in hydroclimate dynamics between Western and Eastern Amazonia—a consequence of the regionally discordant effect of glacial cycles on patterns of precipitation (62, 63). Together with smaller-scale, patchy dynamics of forest canopy density (64), recent lineage diversification in the Amazon Basin may have occurred principally as a consequence of sporadic dispersal events and subsequent persistence in isolation (32). Although not at the scale of local forest dynamics, our simulations allow us to assess the degree to which these regional patterns may have been driven by each of the processes we modeled (Fig. 1), in the context of Quaternary climate cycles.

Strategy and scope of the study

Predecessors of our simulation model (22, 23) targeted documented patterns of species richness and range size distributions to guide the exploration of parameter space and to assess the sensitivity of outcomes to individual parameters and their values. The present study takes a more fundamental approach, relying on the realism of the modeled climatic and topographic drivers and modeled ecological and evolutionary processes (Fig. 1) to produce meaningful biogeographical patterns. The simulations had no target pattern and were not parameterized with empirical data. Surprisingly, as we will show, richness maps nonetheless emerged from the simulations that closely resemble contemporary richness maps for South American birds, mammals, and plants, including regional details that mirror conjectures in the biogeographical literature, as outlined above.

Although our paleoclimate model extends further into the past than any three-dimensional atmosphere model previously applied at this temporal resolution, the model nonetheless encompasses only the Late Quaternary (800 ka ago to the present), with its repeated glacial-interglacial cycles, extending as far into the past as the high-precision CO2 record from Antarctic ice core data (65). South America was, of course, already populated with a rich biota comprising many distinct lineages—some quite ancient—at the beginning of this period (66).

Nonetheless, current consensus among biogeographers and paleoecologists is that the contemporary flora (67, 68) and vertebrate fauna (4, 33, 34, 69, 70) of South America include numerous lineages that have undergone rapid diversification during the Quaternary, particularly in the Andes. In contrast, the flora and fauna of the South American tropical lowlands, including the Amazon, are generally considered to be more ancient (33, 40, 70). It is likely, however, that the geographical distributions of most species, whether belonging to an old lineage or a young one, have been shaped by Quaternary climate cycles (71, 72) (Movie 1 and fig. S2).

Movie 1. Spatial and temporal dynamics of South American climates.

The four dynamic maps on the left display minimum and maximum annual precipitation (upper maps) and temperature (lower maps) in South America over the last 800 ka. The colored lines in the corresponding time-series plots (center) indicate, from the top to bottom, (i) maximum, (ii) third-quartile, (iii) median, (iv) first-quartile, and (v) minimum annual precipitation (upper time series) and temperature (lower time series) among map cells. For precipitation, minimum and first-quartile time series overlap. In the dynamic temperature-precipitation climate space (right), each cross corresponds to one grid cell in the map. All cells are illustrated. The width of the cross indicates the annual precipitation seasonality (difference between maximum and minimum), while the height of the cross indicates annual temperature seasonality. The gray scale of individual crosses varies to allow climatically overlapping cells to be visually distinguished.

Cradles, museums, and graves

Although conceptually simple (Fig. 1), the model yielded extraordinarily complex patterns of diversity in space and time. To make sense of the simulations, we examined the history of each simulated species and its contribution to these patterns. Over the course of each simulation, a complete phylogeny emerges from a single founding species. On the basis of this phylogeny and on full historical records of each range and range fragment at each 500-year interval of the modeled paleoclimate data (Movie 1), we analyzed and illustrated spatial and temporal patterns of speciation (“cradles”), persistence (“museums”), extinction (“graves”), and species richness within South America.

Stebbins (73) began a long tradition of referring to locations with unusually high rates of speciation as “cradles” of diversity, and to locations with unusually low rates of extinction as “museums.” Although these terms have previously been applied almost exclusively to broad comparisons between tropical and boreal latitudes (16, 7478), here we follow Fjeldså et al. (4) in downscaling these analogies to the regional level, within South America. In addition, the full evolutionary and biogeographical records that arise from our simulations allow us to define and map a third biogeographical category, “graves”—locations with unusually high extinction rates—and to document not only where, but also when cradles, museums, and graves were most and least active.

As we define them here, “cradles” are about speciation, “museums” about persistence, and “graves” about extinctions. Previously, cradles and museums have generally been viewed as fixed geographical places (4). Because our simulations take place in both space and time, we treat all three patterns as driven dynamically by the processes of speciation, persistence, and extinction. A cradle, museum, or grave has extension and intensity in both space and time, and may move though space and change shape, size, and intensity as time passes.

Lifetime trajectories

Every species in the simulation has a lifetime trajectory, in space and time (Fig. 2). Temporally, each species’ lifetime trajectory extends from its time of origination (a range fragmentation event that leads subsequently to speciation) to one of three endpoints: the point in time when the species splits into two isolated populations (range fragments) that eventually become daughter species, the point in time of its extinction, or the present time—if the species is still alive at the end of the simulation. Each species’ lifetime trajectory is subdivided into three consecutive, distinct, and fully inclusive segments: a speciation trajectory, a persistence trajectory, and, if the species goes extinct during the simulation, an extinction trajectory (Fig. 2). (Species that give rise to daughter species or persist into the present lack an extinction trajectory.) The speciation trajectory covers the period of population isolation between range fragmentation and full genetic isolation, Tmin years later. The extinction trajectory begins when a species starts an inexorable decline (defined statistically) toward extinction. The persistence trajectory comprises the time interval between the consolidation of speciation at Tmin and the beginning of the extinction trajectory.

Fig. 2 Lifetime trajectory of species.

Initially, the species on the left (labeled “ancestral population”) is in a persistence trajectory (thick black line), as a single, viable population. Driven by climate change, the population experiences range fragmentation, yielding two, isolated descendant populations (blue and red dashed lines). These two daughter populations enter speciation trajectories. Once they have remained isolated for at least Tmin years, they are considered independent species (speciation event). Each descendant species then enters its own persistence trajectory (blue and red solid lines). In this example, after a short period of persistence, the red species enters an extinction trajectory (thin dashed red line), as its geographic range continuously contracts in a changing climate, ending in full range collapse (species extinction). The blue species will eventually give rise to two daughter species, undergo extinction, or survive to the end of the simulation.

Occupancy maps and time series

At the end of each simulation, we determined the lifetime trajectory of each species (Fig. 2) and its component segments (speciation, persistence, and extinction) by moving backward in time through the records of the simulation. We recorded the number of time steps that each map cell was occupied by each species for each segment of its lifetime trajectory. We then summed these records for all species—separately for speciation, persistence, and extinction trajectories—to produce five cumulative occupancy maps for the entire simulation: a cradles map, a museums map, a graves map, a net diversification map (cradles minus graves), and a total richness map (fig. S15). Each of these maps is a summation over time. The cumulative total richness map is simply the spatial overlay (summation) of the cradles, museums, and graves maps—a map of total occupancy by all descendants of the founding species, summed over the course of the simulation. Each species in the simulation, whether extinct, an ancestor of other species, or living at the end of the simulation, contributes to these maps. To visualize the temporal pattern behind these cumulative maps, in relation to climate and to each other, we plotted occupancy time series for cradles, museums, graves, net diversification, and total richness, by summing occupancy over all cells for each time step, for each map (fig. S13).

Cumulative occupancy maps provide a deep compilation of historical information on emerging patterns of richness and their dynamics, summed over the time course of a simulation. Occupancy time series, by contrast, represent cradles, museums, graves, or total richness, summed over the entire continent, for each time step of a simulation.


We carried out 10,500 simulations, each spanning the entire 800-ka scope of the paleoclimate time series, to assess sensitivity of the model to its parameters and to a battery of initial conditions. As detailed below in Methods, we treated the four model parameters (Fig. 1) and two sets of initial conditions (founder location and climate smoothing) as factors in a fully realized factorial design, specified in table S5. Movie 2 illustrates the structure and dynamics of these simulations for a small clade, together with a corresponding phylogeny. The temporal and spatial dynamics of species richness, cradles, graves, and net diversification for a larger clade are illustrated in Movie 3.

Movie 2. Demonstration of simulated geographic and evolutionary dynamics for a small clade of Andean origin.

In the temperature versus precipitation climate space diagram (top left), the climatic niche of each extant population is indicated by a rectangle, defined by the population’s maximum and minimum climatic tolerance for temperature and precipitation. As the simulation progresses, and populations become fragmented, the niche of each fragment is represented by its own rectangle. Niches of different populations of the same species share the same color, whereas different species’ niches are shown in different colors. The dynamic map (top right) shows the richness of species at each time step. The phylogeny (bottom) records the events of speciation and extinction that emerge from the interaction of climate dynamics, geographic distribution, and the evolutionary response of species.

Movie 3. Emerging spatial and temporal patterns of species richness, cradles, museums, and net diversification for a rapidly speciating Andean clade.

Spatial patterns of instantaneous (top row) and cumulative (bottom row) total species richness (first column), cradle richness (second column), grave richness (third column), and net diversification (fourth column; the difference between cradle richness and grave richness). Cumulative richness is the sum of instantaneous richness over time, capturing—in a single map—an overview of historical spatial patterns.

Impact of parameters and initial conditions

To assess the impact of parameters and initial conditions on emerging spatial and temporal patterns in biodiversity, we partitioned sources of variation (79, 80) among these patterns, based on matrices of quantitative Bray-Curtis dissimilarities computed between pairs of cumulative occupancy maps for cradles, museums, graves, and total richness generated by each simulation (table S7 and figs. S16 to S19). In these analyses, a parameter or initial condition was judged to be influential if it yielded consistent spatial patterns for any particular parameter setting or initial condition, and different patterns for different settings, regardless of the settings of other parameters or initial conditions. Factors that varied little in their influence on the mean behavior of the simulations, regardless of their settings, were judged less important.

Model parameters and initial conditions varied greatly in their impact on the simulations, but none was as influential, overall, as founder location (Figs. 3 to 5, figs. S16 to S19, and table S7), suggesting an underappreciated impact of historical contingencies in current patterns in species richness (81). The second-most-influential model parameter was the maximum sustainable evolutionary rate realizable by a population (Hmax), which limits the adaptability of niche limits and evolutionary rescue (8285) in the face of changing climates (figs. S16 to S19 and table S7). Low Hmax values indicate that niche traits have low genetic variance, low population growth rates, or both—preventing species from tracking and adapting to changing climates. On evolutionary time scales, this limitation yields a pattern of niche conservatism. Although the balance varied among founders, intermediate levels of adaptive evolution promoted the greatest diversification. If Hmax was too low (strong niche conservatism), species and lineages were subject to extinction; if too high (fast niche evolution), a few species became ubiquitous, and little diversification occurred.

Fig. 3 Simulation results for an Andean founder.

Upper panel: Occupancy time series for speciation (cradle richness, green), extinction (grave richness, red), and mean continental temperature (blue) over the course of the simulation (time moves from left to right). The highest five to seven peaks of speciation (green dashed lines) and extinction (red dashed lines) were marked manually, but time-series cross-correlations were analyzed rigorously (table S6). Precipitation time series appear in fig. S2. Lower panel: Cumulative richness maps for cradles, graves, net diversification (cradles minus graves), and total richness. Each map is a summation over the course of the simulation. The figure shows the average of all parameter values for an Andean founder, excluding the climate-smoothing experimental treatments.

Fig. 4 Simulation results for an Atlantic Rainforest founder.

See the caption for Fig. 3.

Fig. 5 Simulation results for an Amazon founder.

See the caption for Fig. 3.

The effect of spatial and temporal heterogeneity in climate, as assessed by sequential levels of experimental climate-smoothing, ranked third in its capacity to drive variation among simulation outcomes (Fig. 6, figs. S14 to S19, and table S7), as higher levels of heterogeneity promoted faster diversification.

Fig. 6 Pooled simulation results for the Andean, Atlantic Rainforest, and Amazon founders of Figs. 3 to 5.

See the caption for Fig. 3.

Maximum dispersal distance (Dmax) ranked somewhat lower in overall impact. Greater dispersal capacity increased speciation (cradle richness) by promoting occupation of disjunct, yet climatically suitable, regions that initially lay within Dmax but that later became isolated through climate change (figs. S14 to S19 and table S7). Extinction rates (grave richness), by contrast, decreased with greater Dmax, as declining populations were rescued from extinction by dispersal to suitable climates. Thus, net diversification increased with larger Dmax, by its combined effects in increasing speciation and decreasing extinction rates.

The remaining model parameters, minimum time in isolation for speciation (Tmin) and maximum intensity of competition allowing coexistence (Cmax), proved to have surprisingly little effect on the simulations, compared with the other parameters (figs. S14 to S19 and table S7). Regardless of the rate of speciation, similar spatial patterns eventually emerged. We surmise that competitive exclusion at the map cell level, as we modeled it, tends to have only local and probably ephemeral effects, as climates change and species adapt to these changes, with little net impact on large-scale patterns of species richness.

Founder location, cradles, and graves

Figures 3 to 5 show the results for Andean, Atlantic Forest, and Amazon founders (each figure summarizes 375 simulations—all those without experimental topographies), and Fig. 6 combines the results for all three founders. The maps in these figures display cumulative occupancy, summed over the course of the simulations, for cradles, graves, net diversification, and species richness. The corresponding occupancy time-series plots in these figures capture the temporal dimension of speciation and extinction by summing occupancy over all cells in South America for each time step, on the basis of the time-specific occupancy maps of cradles and graves.

In space, Figs. 3 to 5 illustrate the strong influence of founder location. The initial niche of the single founder, in each region, necessarily differed among regions for initial survival, and constraints on niche evolution limited continental-scale convergence in pattern, but these figures share many features of spatial pattern and temporal dynamics, independent of starting location.

With regard to time, the most distinctive feature shared by all simulations is the obvious quasi-periodicity of peaks and valleys of speciation (cradles) and extinction (graves), as shown in the occupancy time-series plots in Figs. 3 to 6. Time-series analysis of log-transformed, detrended data for cradles and graves, in relation to mean annual continental temperature, yielded many significant, time-lagged, cross-correlations (table S6), confirming a subtle but certain role for glacial-interglacial temperature cycles [and thus, for orbital forcing (86)] in driving cycles of speciation and extinction for all founders. Both speciation and extinction tended to peak during glacial terminations, as warming climates returned. Peaks of extinction closely followed peaks of speciation for all three founders—by 18 ka for an Andean founder (Fig. 3), by 20 ka for an Atlantic Forest founder (Fig. 4), and by 24.5 ka for an Amazon founder (Fig. 5). (Time-series analysis for precipitation variables yielded very low correlations, so was not pursued further.)

A notable feature of these simulations—for the Atlantic Forest and Amazonian founders (Figs. 3 and 4)—is the spatial coincidence of cradles and graves, best illustrated by the net-diversification maps, which plot net spatial differences in magnitude between speciation and extinction. In contrast, cradles for the Andean founder are concentrated along the Andean slopes, whereas graves tend to be at lower elevations in the upper Amazon Basin (Fig. 3). Population decline drives both speciation and extinction—speciation through range fragmentation and extinction by range contraction. We conjecture that enviromental heterogeneity (driven by topographic complexity and elevational climate gradients) in the Andes promotes range fragmentation, while at the same time offering climatic refugia from extinction compared with lower elevations (46, 49), thereby focusing cradles at mid-elevations and graves at lower elevations.

Experimental topographies

Our experiments with climate smoothing, as a proxy for decreased topographic heterogeneity, yielded clear-cut evidence for the role of spatial heterogeneity in the location and intensity of cradles of speciation and net diversification, with considerably less effect on graves of extinction (Fig. 7 and fig. S21). Overall, nonsmoothed (realistic) climates promoted three times as much diversification as spatially smoothed climates, with the effect tapering off for smoothing kernels larger than 250 km. The strongest effects were experienced by the Andean founder clade (Fig. 7), where even the smallest possible kernel radius reduced diversification by a factor of 7, whereas the reduction for Atlantic Rainforest founder was only by half.

Fig. 7 The effect of topographic smoothing on rates and cumulative spatial patterns of speciation (cradles), extinction (graves), net diversification (cradles minus graves), and total richness, for Andes, Atlantic Forest, and Amazon founders, pooled.

Upper panel: Occupancy time series for speciation (cradle richness, green), extinction (grave richness, red), and mean continental temperature (blue) over the course of the simulation (time moves from left to right). Black time series are for smoothed topographies. Red and green time series are the same as in Fig. 6. Rates of speciation (cradles) and extinction (graves) were both suppressed by smoothing.

Biogeographical interpretations

Emerging role of the Andes

By mapping the distribution of South American cradles of diversification, in time and space, our simulations offer strong support for the role of the Andes as an episodic “species pump” (38, 69, 87). This phenomenon has been documented not only for endemic Andean clades (43), but also for clades later centered in the Amazon and Atlantic Rainforest (40, 52, 88).

By experimentally smoothing South American topography, we gradually eliminated all of the diversification processes driven by Andean topography. The extent and magnitude of Andean cradles declined drastically, whereas South American graves broadened and intensified (Fig. 7 and fig. S21) and net diversification decreased (Fig. 7 and figs. S14 and S20). With experimentally smoothed topography, the simulations ceased to produce realistic spatial patterns of species richness (fig. S21). These experiments offer compelling evidence that topographic complexity—or more accurately, the fluctuating climatic complexity that constantly mirrors it—drives diversification and biogeographical patterns. This result shows unambiguously not only that climatic complexity promotes diversification, but also that diminished complexity drastically slows diversification.

Emerging regional patterns

Regional biogeographical dynamics in our simulations, on a broader scale, also conform to many expectations from empirical studies. Our simulations (Movie 3) directly support Por’s (57) proposal of an ephemeral connection between the Amazon and Atlantic Rainforest during late Quaternary climate cycles (5860), as well as recent suggestions that Atlantic Rainforest birds may have dispersed through the Cerrado during glacial maxima and through the Chaco during interglacial periods (52, 59).

Our simulations display a pattern of ephemeral, circum-Amazonia “arcs” of seasonally dry climates—sometimes patchy and sometimes continuous—connecting the tropical Andes and Atlantic Forest (50, 61, 89). In our simulations, these episodic biogeographical bridges acted for some species as dispersal corridors, for others as refugia from extinction, but for many others as graves. We did not find evidence, at the scale of our current analysis, for the Amazonian dipole hypothesis (63, 64), although a specific, local investigation might be fruitful.

Emerging patterns of speciation and extinction

Our analyses of temporal patterns of speciation (cradles) and extinction (graves) offer strong support for both generative and erosive effects on biodiversity arising from continental-scale climate change driven by Quaternary glacial-interglacial cycles (Figs. 3 to 6). As warming accelerated following glacial episodes, first speciation, then extinctions spiked. These results are in accord with accumulating evidence for episodes of heightened extinction during Quaternary deglaciations (9092) and support concerns about extinction under rapid anthropogenic climate warming.

In space, cradles and graves closely coincided for Amazonian clades (Fig. 5), suggesting that gradual warming promoted speciation and more rapid or extensive warming in the same region drove extinctions. For an Andean clade (Fig. 3 and Movie 3), cradles were concentrated in the highlands, whereas graves accumulated the adjacent Amazonian lowlands, perhaps suggesting that high climate velocity (93) or climatic divergence (94) in the warming lowlands overwhelmed their capacity to escape to the Andean slopes. This process is also likely to have been a cause of Amazonian extinctions.

Comparison with contemporary patterns of richness

If the processes that we have modeled are realistic representations of the processes that have shaped contemporary biogeography, then a comparison between simulated, historical patterns and empirical, contemporary patterns of species richness should be instructive. Our simulations were not carried out with regard to the richness pattern of any particular real-world taxon, nor were model parameters fitted by targeting such patterns. Indeed, not only the model design, but also the ranges of parameter values were defined strictly on first principles of biogeography, ecology, and evolutionary biology, without regard to any empirical data. Thus, any resemblance between model output and real-world richness patterns must be attributed to (i) the underlying processes built into the model, (ii) the topographic and modeled climatic milieux in which they operated, and (iii) theoretically estimated parameter values that regulated the simulated ecological and evolutionary processes.

Each simulation began with a single founder, 800 ka ago, and unfolded over a geologically brief period of time until the present. Thus, in principle, the present time in each simulation might seem the most appropriate for comparison with the empirical maps of contemporary richness. However, preliminary tests showed that many simulated maps of contemporary richness, especially for Andean founders, were surprisingly species-poor, following massive extinctions in post–last-glacial-maximum (LGM) warming. Figure 3 shows that this Holocene extinction peak is just one of several, each associated with an interglacial warming period for the Andean simulation. In contrast, Amazon founders (Fig. 5) do not show this pattern. We conjecture that the model, as it stands, exaggerates episodes of clade extinction by failing to account for the survival of species under outlier climates in some regions, perhaps supporting a role for microrefugia (4, 49) that lie under the radar of the spatial scale of the model. Moreover, our paleoclimate model exhibits LGM cooling at the high end of the range of full-complexity climate models [see “Comparisons against other paleoclimate models” in (95)], perhaps contributing to the high rates of extinction simulated during deglaciations. Finally, our model does not account for potentially important factors hypothesized to promote speciation, such as subcanopy forest dynamics (64) and the historical origins of the hydrographic network, which are thought to have promoted isolation and vicariance effects in some clades (96).

In contrast to the simulated contemporary map of species richness, cumulative maps (Figs. 3 to 6) provide a richer compilation of historical information on emerging patterns of richness and their dynamics over the time course of the simulations. By aggregating patterns from every phase of the diverse (97) glacial-interglacial cycles and averaging over all parameter values, these maps represent the broader range of possible patterns arising from local and regional processes captured by the model. Thus, cumulative maps of species richness offer a much more representative historical account of model behavior than any single point in time in the simulation, including the present, which represents contemporary climatic conditions—an outlier within the distribution of Quaternary climates.

Cumulative cradle maps record the frequency and location of species origination, grave maps point to regions where species tend to collapse under climate change, and museum maps show where species tend to persist over longer historical periods. Thus, if contemporary patterns of species richness for large clades are representative of deep and persistent historical patterns, we should find stronger correspondence with simulated patterns of museum species richness. For rapidly speciating clades, by contrast, we would expect stronger correspondence with cradle richness, and for clades on the decline, we might expect a better correspondence with grave richness.

Because continental-scale maps of contemporary species richness for South America are scarce and fraught with sampling and data quality problems (98), we used a lower resolution 1° latitude by 1° longitude grid of 1659 cells to develop maps for 2967 species of birds, 1342 species of mammals, and 61,724 species of plants. To compare our simulation outputs with these maps of contemporary richness, we rescaled our hybrid-scale richness maps (figs. S1 and S15 to S19) to the same, uniform 1° by 1° resolution.

Simulated historical patterns of species richness, as recorded by maps of cumulative museum richness, proved very successful in capturing, proportionally, the broad outlines of the empirical richness maps (Fig. 8 and Movie 4) for birds (r2 = 0.6337 for an Atlantic Forest founder, fig. S22), mammals (r2 = 0.6548 for an Atlantic Forest founder, fig. S23), and plants (r2 = 0.4146 for an Andean founder; figs. S24 and S25). Tables S12 to S14 provide more detail and make clear that the maps in Fig. 8 and figs. S22 to S25 represent the best of strong patterns, not one-off accidents.

Fig. 8 Observed species richness versus modeled richness.

Left column of maps: Contemporary spatial patterns for South American bird richness (2967 species, upper map) and mammal richness (1342 species, lower map). Middle column map: The simulated spatial pattern for cumulative museum richness, arising from the model (Fig. 1), averaged over all parameter values for an Atlantic Rainforest founder. Right column of maps: The differences between observed (left maps) and simulated (middle map) richness for birds (upper map) and mammals (lower map). Red indicates that the model underestimates richness, and blue indicates overestimation. Simulated species richness is highly correlated with observed species richness for birds (r2 = 0.6337) and for mammals (r2 = 0.6548). Observed species richness was not targeted in any way by the simulations. A qualitative comparison of modeled richness with South American plants appears in “Contrasting empirical and simulated spatial patterns in species richness” in (95).

Movie 4. Emerging spatial and temporal patterns in museum species richness, averaged for the Andean, Amazonian, and Atlantic Rainforest founders.

Cumulative patterns of cradle, grave, and museum species richness (first three columns), for Andean, Atlantic Rainforest, and Amazonian founders (rows), from model simulations. Static empirical maps (on the right) show contemporary patterns of plant, bird, and mammal species richness. Simulated patterns of cumulative museum richness, over the course of 800 ka, closely resemble current patterns in species richness.

This high level of correspondence between modeled and empirical richness raises an obvious question about time scales. Although several lineages are known to have diversified actively during the time span of our Quaternary simulations, particularly in the Andes (67, 68), most living species in South America are much older than any species in our simulations (33, 40, 70). We do not suggest that our simulations over the geologically brief period of 800 ka might reproduce the actual ranges, evolutionary dynamics, or phylogeny of any living species in the South American biota. Our simulated “species,” however, may just as well be viewed as independent evolutionary units below the level of taxonomic species, nonetheless subject to the same ecological and evolutionary processes. As phylogeographical studies (88) and tools for studying ancient DNA reach farther into the past (91), models such as ours can be expected to take on even greater realism.

A second question is why richness patterns of living species, from the present (outlier) interglacial climate—a single point in time, should correspond so well with cumulative richness patterns from the simulations. Our cumulative richness maps pool both glacial and interglacial distributions, and we know from paleoecological studies that Quaternary temperature cycles (including Holocene warming) shuffled many extant species over elevational (71) and latitudinal (72) gradients, just as in our simulations (Movies 3 and 4). We conjecture, first, that the geographical core of richness patterns may have been more persistent over geological time than generally thought, and, second, that the maps of residuals between simulated and empirical richness (right column in Fig. 8 and figs. S22 to S24) may correspond, at least in part, to regions where past richness differed from present empirical patterns—a topic that merits further research.

What does it all mean?

Our simulations have three strengths. First, they take place in a topographically realistic continental landscape, driven by a paleoclimate model built on well-established principles. Second, the biogeographical simulations were constructed from the bottom up, by integrating mechanistic models of key ecological and evolutionary processes, following well-supported, widely accepted explanations for how these processes work in nature (Fig. 1). Third, despite being entirely undirected by any target pattern of species richness, covering only a tiny slice of the past, and being controlled by only four parameters (two of which turned out not to be very important), surprisingly realistic biogeographical patterns nonetheless emerged from the simulations, not only on a continental scale (Fig. 8), but also on regional scales.

Adaptive niche evolution as a biogeographical force

Phylogenetic niche conservatism (30) is the universal tendency of descendant species to retain the fundamental niche of their ancestors (99101). It may be strong or weak. By modeling adaptation to climate in the trailing edge of a shifting range, our simulation model explicitly regulates the capacity of a species’ climatic niche to respond to climate change by adaptive evolution (evolutionary rescue) (102, 103). When the potential for adaptive evolution is weak (low Hmax, Fig. 1), a pattern of strong niche conservatism emerges. Descendant species accumulate in regions that are climatically similar and geographically close to the original range of the ancestor, with a gradual decline in richness of descendant species with increasing distance and decreasing climatic similarity (104). The distribution of descendant species is constrained by higher extinction rates, as species fail to adapt to changes in trailing-edge conditions. In contrast, when the potential for adaptive evolution is strong (high Hmax, weak niche conservatism), a pattern of niche evolution emerges, with adaptive shifts to novel climates and a broader geographical spread of descendant species, but little diversification, because ranges rarely fragment, as niches adapt to all challenges. Our simulations provide unequivocal support for intermediate levels of adaptive niche evolution as a key factor driving realistic patterns of species richness (figs. S12 to S14, fig. S16, and table S7), confirming the findings of previous simulations (22, 23).

Emerging patterns on a continental scale

In the simulations, the facilitating influence of adaptive niche evolution, acting within the constraints of topography and climate, yielded cumulative patterns of species persistence (museum richness) that correspond well with contemporary richness patterns of birds, mammals, and plants (Fig. 8, figs. S22 to S24, and Movie 4). The inference that contemporary empirical patterns of richness have their origins in the same underlying processes, driven by climatic changes in the same landscapes, seems nearly inescapable.

By revealing the regions and periods of speciation, persistence, and extinction that underlie richness patterns, our results illuminate the role of history in shaping contemporary patterns of species richness on broad spatial scales—an emerging theme in the recent macroecological literature (104106). Observed statistical correlations between contemporary richness patterns and current climate variables (9) should be viewed not as a direct causal link, but rather as a consequence of accumulated historical events driven by geographically structured climate dynamics (86, 107).


Geographical domain

The simulations took place on a gridded map of contemporary South America. The computational demands of spatially and temporally explicit simulations impose a limit on the complexity of simulation models at very high spatial resolutions. Nonetheless, at the large spatial and temporal scales at which we model ecological and evolutionary systems here, topographic heterogeneity, expressed as habitat diversity, is thought to be a key driver of species distributions and evolutionary niche dynamics. Thus, to capture as much of the climate heterogeneity of South America as feasible, while accounting for computational limits imposed by the spatial resolution of the geographic domain, we developed a map grid of “hybrid” spatial scale, in which the 4820 square map cells vary in size in inverse relation to topographical complexity. Cell sizes ranged from 625 km2 in rugged areas of the Andean slopes, where environmental conditions vary greatly within short distances, to 10,000 km2 in flatter regions, such as the Amazon Basin and Patagonia, where large areas have relatively similar environmental conditions (fig. S1). Given the ecological and evolutionary mechanisms implemented in our simulation model (see below), the spatial resolution of our hybrid grid constitutes a balanced trade-off between (i) the computational demand imposed by the number of map cells; (ii) the inherent uncertainty in reconstructing terrestrial paleoclimate dynamics at high spatial resolution; (iii) the organizational level of the mechanisms that drive the evolutionary dynamics of geographical ranges and climatic niches; and (iv) the low resolution of current data on the distribution of real-world species (which we used to evaluate the predictions of the simulation model), and the consequent uncertainty of mapping empirical spatial patterns of biodiversity at higher resolutions.

Paleoclimate simulations

A paleoclimate emulator was developed to overcome the computational challenge of simulating 800,000 years of climate. The emulator was built around the PLASIM intermediate-complexity atmospheric general circulation model, coupled to the ENTS dynamic land surface model and to flux-corrected ocean and sea-ice models, at a 5° latitude-longitude resolution (108). Orbitally forced climates at 500-year intervals were estimated using principal component emulation (109). A transient 800-ka simulation (110) with the faster GENIE-1 model was applied to scale these emulated climates for CO2 and ice sheet climate forcing. Spatial resolution matching the hybrid-scale map was achieved by treating modeled climate variables as anomalies from contemporary climate data, spatially referenced to each map cell.

For each 500-year time interval, from 800 ka to the present (1600 time steps through the Late Quaternary glacial-interglacial cycles), the paleoclimate model assigned to each of the 4280 map cells an estimate of the mean temperature of the warmest and coolest quarters (henceforth minimum and maximum annual temperature) and the mean daily precipitation of the wettest and driest quarters (henceforth minimum and maximum annual precipitation) (Movie 1 and fig. S2).

The temporal resolution of the 500-year interval between time steps is compatible with the macroecological framework used in this study. Assuming a species with a generation time of 5 years, one time step would encompass 100 generations, a reasonable resolution for the population-level biogeographical processes that we are modeling, such as dispersal, competition, range dynamics, and niche evolution. Thus, the full temporal scope of the simulation would encompass ~160,000 generations, well beyond a reasonable time for the emergence of medium-sized lineages (4, 33, 34, 69, 70). Indeed, a finer temporal resolution would probably convert the current model from the population or species level closer to the individual level of organization, requiring a full redesign of the implemented mechanisms (e.g., individual movement, birth-death processes).

Our 500-year resolution is also compatible with currently available knowledge of paleoclimate dynamics and the complexity of our paleoclimate emulator. Indeed, our climate model does not capture shorter–time-scale variability in climate dynamics, because the emulator is built from quasi-equilibrium snapshots, forced by only orbit, CO2, and ice sheets and therefore does not capture variability due to relatively short-term phenomena such as glacial meltwater-driven ocean circulation changes, El Niño–Southern Oscillation, or volcanic eruptions.

To evaluate the reliability of our paleoclimate emulator, we compared spatial patterns of temperature and precipitation variables, at specific time steps, against multimodel predictions carried out by the Paleoclimate Model Inter-comparison Project, phases two and three (PMIP2/PMIP3) [see “Comparisons against other paleoclimate models” in (95)]. We focused the validation of our emulator at three specific moments of the Quaternary: the Last Interglacial (LIG) at ~126.5 ka ago (111, 112), the LGM at 21 ka ago (112, 113), and the mid-Holocene (MH) climate optimum at 6 ka ago (112, 113). The LIG and MH interglacial states, with CO2 and ice sheets similar to those of the present day, provide an opportunity to validate our emulated response to orbital forcing, while our estimates of paleoclimate at the LGM test the emulated response to very different CO2 and ice-sheet forcings. The climate patterns predicted by our model are consistent with multimodel ensemble predictions from complex climate models, suggesting that our emulator can be reliably used in biogeographical simulations, given currently available parallel evidence from independent models.

Biogeographical simulations

In the simulations, geographical space was represented by the gridded map of South America (as detailed above), with each grid cell characterized by its area and geographical position. Each simulation began 800,000 years ago, advancing at 500-year time steps. At each time step, each map cell was characterized by four climatic conditions (minimum and maximum temperature and precipitation), as reconstructed by the paleoclimate simulations (Fig. 1 and Movie 1).

Populations and species

The smallest biological unit explicitly modeled was regarded as a population, characterized as a geographically isolated and continuous species range or range fragment. Thus, the complete range of a species might consist of a single population or of multiple, isolated populations. At each time step, the niche of each population was defined as a two-dimensional region within temperature and precipitation axes, in the same space as that of the modeled paleoclimate of South America.

Climatic niche and geographic distribution

We modeled the evolution of the fundamental climatic niche (114) of populations, which defines the extremes of temperature and precipitation that a population can tolerate at any given time step (Movie 2 and fig. S7). Thus, cells occupied by a population must have climatic conditions within the limits of its fundamental niche. The realized climatic niche is an emergent property of the population, defined by the climatic conditions that it actually experiences across the whole set of occupied cells (e.g., population’s range). However, not every cell with suitable climatic conditions is necessarily occupied by the population (Movie 1). Indeed, the bounds of the population’s evolving fundamental climatic niche may at times extend beyond the limits of its realized niche (115), owing either to dispersal limitation (additional, climatically suitable regions currently exist but cannot be reached), niche conservatism (the fundamental niche has not yet responded to the disappearance of previously suitable climates) (26), or exclusion by competing species (see below).


In each realization of the biogeographical simulation, an evolutionary lineage develops on the grid from a single founding species—initially a single population. The initial geographical range of the founder is determined by its assigned geographical location (a single map cell—an initial condition of the model, table S2), and by a preset environmental niche (see section “Initial conditions” below).


Within each time step, each existing population expands its range to occupy not only adjacent cells, but also disjoint cells with suitable climates, as long as these cells are not separated by cells with unsuitable climate spanning a distance greater than Dmax (a model parameter; Fig. 1 and table S1). All subsequent events are autonomous and deterministic (repeatable), driven by the interaction between climate, topography, and modeled processes.

Evolutionary niche dynamics and evolutionary rescue

The four paleoclimatic conditions in each cell change asynchronously over time and space. Climate dynamics may open opportunities for range expansion by turning an unsuitable cell into a suitable one (a leading-edge cell of a shifting range). In this case, the population simply expands its range to occupy any newly suitable cell, whether contiguous or not, as long as the cell lies within Dmax of the existing range. Climate change, however, may also render a suitable cell unsuitable (a trailing-edge cell of a shifting range), imposing selection pressure on the population in the grid cell. The outcome of trailing-edge selection may be (i) local population extirpation in the trailing edge cell, if the population cannot adapt, or (ii) niche evolution (partial or full adaptation to the new environmental conditions), allowing continued occupation of the trailing edge cell. In nature, selection pressure from climate change, especially in the trailing edge, may cause a gradual adaptive niche shift, bringing niche limits closer to new climatic limits within the current geographic range of the local population (30, 116). This process has been called “evolutionary rescue” (8285), as it promotes species persistence by means of evolutionary changes in niche limits in response to selection pressure imposed by climate change.

We implemented niche evolution as a response to climate change in a simple, quantitative evolutionary genetics framework (117, 118). The evolutionary rate (H) required for sufficient niche adaptation to allow the population to persist in the trailing edge cell c may be estimated by comparing the magnitude of climate change between two consecutive time steps (t and t + 1) in cell c,Hc = [(Ec,t+1Ēt) / σt)] / Δtwhere Hc is the adaptive rate necessary for evolutionary rescue [measured in units of Haldanes (119, 120)], Ec,t+1 is the value of an environmental variable (e.g., maximum annual temperature) in cell c after climate change (time t + 1), and Ēt and σt are the average and standard deviation of the same environmental variable, before climate change (time t), in all cells occupied by the local population within the genetic neighborhood of c (modeled as a circle centered at c, with radius Dmax). Thus, in our model, genetic variation within a species’ range is geographically structured, so that the evolutionary potential of each trailing-edge cell is set, at each time step, by the genetic variation for climate (standard deviation) within the genetic neighborhood of the cell.

In the simulation model, a maximum (critical) evolutionary rate in response to climate change (model parameter Hmax, table S1) is defined for all species, in all trailing-edge cells, uniformly throughout the entire time span of the simulation. Thus, for a trailing-edge cell c, if Hc < Hmax the population is rescued at c by adapting to the new climate, expanding its niche. Conversely, if Hc > Hmax the evolution required for persistence is beyond the maximum evolutionary potential of the population in cell c, and the population is extirpated from the cell.


In classical ecology, species with excessively similar resource requirements cannot coexist in sympatry (121). However, models on broad spatial scales must somehow account for the resources for which species compete, without modeling individual consumers and a myriad of resources and their respective depletion rates. We modeled interspecific competition, without explicitly modeling resources, by implementing the classic assumption that competition is an inverse function of phylogenetic relatedness (122), as measured by the explicit phylogeny generated by our model (123). Assuming that the use of resources by species (e.g., food items, foraging time/strategy) evolves at a constant average rate with variance proportional to time (i.e., a Brownian motion model of trait evolution), the expected intensity of competition between two species declines with phylogenetic distance (PD) between species. Once a pair of sister species achieves a threshold phylogenetic age of Pmin (a model parameter, Fig. 1 and table S1) since divergence, they may coexist in sympatry without competing.

Among the species in each map cell, each species competes against all others from which its phylogenetic distance is less than Pmin. We quantified the intensity of competition between a pair of species as 1 – (PD / Pmin). Thus, the total diffuse competition affecting a particular species in a cell is the summation of the pairwise intensities of competition between that species and all other species present in the cell.

Climatic stress

Assuming that the environmental niche of a population is analogous to a fitness function, individuals occurring in cells with extreme environmental conditions (with respect to the environmental tolerances of the population) have lower fitness than conspecific individuals in climatically more-suitable cells, leading to a lower population density. Conversely, because grid cells with environmental conditions near the center of a population's environmental niche are more suitable for the population, individuals in these cells are assumed to have higher fitness, leading to higher population density. Thus, the cells mapping to the niche center for a species can be considered to offer the most suitable (least stressful) environmental conditions, whereas cells mapping near the niche limits can be considered as the most stressful environmental conditions that nonetheless permit persistence.

We calculated an environmental stress index for each population, in each grid cell, at each time step, as the ratio between (i) the environmental distances between maximum and minimum environmental conditions within the cell and the niche center, and (ii) the maximum environmental scope tolerated by the population. [See “Environmental niche and ecological stress” in (95).] Thus, in a cell with little seasonality and with average climatic conditions similar to those in the niche center of the population, the population has a small environmental stress index. Conversely, a population has a large environmental stress index if the scope of conditions in the cell spans the full range of the climatic tolerance of the population (its niche breadth).

Competitive exclusion

If two coexisting species compete intensely in a particular cell, one of them may be extirpated from the cell. The excluded species is likely to be the competitor under stronger environmental stress, as its population density is likely to be lower. Thus, if the intensity of competition and/or environmental stress is high, the population under greater environmental stress will be excluded from the cell [see “Competitive exclusion” in (95)].

For our simulation model, the index of environmental stress and the index of the intensity of competition were calculated for each population, in cell c, at each time step. These two indexes were then added, for each population, resulting in a single index of competition, Cc, for each population in each cell. All populations occupying a particular cell were then sorted according to the magnitude of this combined index Cc. If the population with the highest competition index Cc had a value greater than the maximum intensity of competition allowing coexistence, parameter Cmax, then that population was eliminated from the cell. The competition index was then recalculated for all remaining populations in the cell c, assuming the absence of the eliminated population, and remaining populations were sorted again. If the population with the highest competition index Cc again had an index greater than Cmax, then that population was also removed. The algorithm iterated until the population with highest competition index in the cell had an index Cc that fell below the threshold Cmax.


The potential geographic distribution of species in our model at any given time step was constrained by available climate, niche limits, dispersal limits, and competitive exclusion. A species became extinct if, at any time step, its entire range was extirpated from all map cells, because of either climate change or competition.

Range fragmentation and coalescing populations

Climate dynamics and competition may cause range fragmentation by imposing barriers of unsuitable climate (Movie 2). When the geographic distribution (range) of an ancestor population became fragmented into independent populations, all smaller populations inherited the environmental niche of the ancestor population (Movie 2). However, owing to founder effects and the spatial structure of genetic variability, smaller populations did not inherit exactly the same niche properties as larger populations. Thus, in our model, in the event of range fragmentation of an ancestor population, the niche limits of the newly isolated, descendant populations were determined by the ancestral population’s niche limits, local environmental conditions, and population size. [See “Environmental niche dynamics of fragmenting and coalescing populations” in (95).] Each population (range fragment) subsequently followed its own evolutionary course.

When the ranges of two populations of the same species were separated by a distance less than Dmax, it was assumed that gene flow was reestablished, therefore coalescing the two populations. Although the niches of the two populations each contributed to the definition of the environmental niche of the newly coalesced population, smaller populations contributed less to the coalescent niche than larger populations. To account for this asymmetry, the contribution of each population was weighted by its range size (total area of occupied cells). Thus, the maximum and minimum tolerance limits of the newly coalesced population, for each niche dimension, were the average of the maximum and minimum tolerance limits of all coalescencing populations, weighted proportionally by their respective range sizes. [See “Environmental niche dynamics of fragmenting and coalescing population” in (95).]


Populations that persisted in isolation beyond a threshold age for speciation Tmin (a model parameter, Fig. 1 and table S1) were assumed to be reproductively isolated and were thus subsequently treated as distinct species (Movie 2). As time passed in the simulation, surviving descendant lineages generated an explicit phylogeny and populated the gridded map, developing patterns of species richness, as species ranges came to overlap following evolutionary divergence (secondary sympatry).

Initial conditions

Certain initial conditions for the model were specified before launching each simulation (table S2). A center of origin (one map cell) was defined for the original founder species, as well as its initial niche (minimum and maximum annual precipitation and temperature tolerated). The historical influence of founder species is believed to have great impact across all scales of spatial and temporal biodiversity patterns (124126). However, because our simulations did not aim to reconstruct any specific real-world lineage, we evaluated spatial patterns in South American biodiversity that emerged from four hypothetical founder lineages, covering the major climatic and geographic zones in South America: high-elevation tropical Andes, lowland Amazonia, lowland Atlantic rainforest, and lowland temperate Patagonia. [See “Experimental design and parameter exploration” in (95).]

Spatial and temporal environmental heterogeneity, particularly in the context of climate change, is widely believed to drive both the extinction and diversification of lineages (81, 127129). In South America, the most extreme climatic heterogeneity is driven by the steep and rugged topography of the Andean mountain chain. Our simulations offer a unique opportunity to assess and quantify the role of topography-driven climatic heterogeneity in ecological and evolutionary modeled mechanisms, as manifested in patterns of cradles, museums, and graves. Thus, we applied a spatial smoothing function to the paleoclimate series, effectively simulating alternative experimental topographies in South America. The climate smoothing factor (an initial condition of each simulation) specifies a smoothing level for all minimum and maximum annual precipitation and temperature maps in the paleoclimate series, thereby generating levels of experimental climatic heterogeneity that defined alternative South American topographies.

Experimental design and model evaluation

To understand the role of the mechanisms implemented in the model (Fig. 1) on emergent patterns of biodiversity, we ran 10,500 distinct simulations, with varying combinations of parameter settings and initial conditions. The factorial design of our simulation experiment consisted in running the model with all possible combinations of parameter values, as listed in “Summary of explored parameter levels and initial conditions” in (95). In our experimental design, we integrated two strategies to define the range of values to be explored for each parameter: (i) a biologically informed definition of the minimum, maximum, and intermediate levels for each parameter, based on the biological interpretation and realism of the implemented process; and (ii) a preliminary experimental evaluation of the feasibility of the simulation, carried out by testing the proposed extreme levels of each parameter. Here, we provide a summary of parameter exploration, but the full conceptual justification may be found in (95).

Parameter exploration

(i) Maximum dispersal distance (Dmax) is a parameter that sets the maximum geographic map distance that a population can disperse across unsuitable climate, over one simulation step of 500 years, to occupy a climatically suitable cell. We specified three intermediate steps between the minimum possible Dmax, given our spatial resolution (150 km), and the maximum Dmax that we considered biologically reasonable, given our temporal resolution (750 km): 200, 350, and 500 km. (ii) Maximum niche evolutionary rate (Hmax) is a parameter that sets the upper limit of potential climatic adaptation of the population in a trailing-edge cell. After a preliminary exploration for a meaningful range of Hmax values, we set five levels, ranging between 0.005 and 0.02 Haldanes, which is, respectively, half and twice the theoretical expectation under natural conditions (117, 118). (iii) Minimum time for speciation (Tmin) is a parameter that regulates the time that a population must remain in genetic isolation before being declared a new species. Although there are no theoretical bounds to Tmin values (except zero), we set three levels for this parameter (17.5, 20, and 22.5 ka—or 3500, 4000, and 5500 generations, assuming a generation time of 5 years), which we considered sufficient for an experimental exploration of meaningful variation in simulated diversification rates. (iv) Maximum intensity of competition allowing coexistence (Cmax) is a parameter that sets the maximum intensity of competition that nonetheless permits coexistence among competing species. We set the minimum experimental value of Cmax to 1.5 units, which in practice specifies that a species under maximum tolerable environmental stress can nonetheless coexist with just one competing species that is phylogenetically close. We gradually explored larger values of Cmax, up to five units, a level at which a species under maximum tolerable environmental stress would nonetheless be capable of coexisting with up to four very closely related species. (v) Minimum phylogenetic divergence for coexistence without competition (Pmin) is a parameter that regulates the phylogenetic distance (PD) beyond which a pair of sister species could no longer compete. Because of the mechanistic association between Pmin and Cmax, we held Pmin at the fixed value of 30,000 generations (150,000 years), to optimize the use of computational resources.

Quantifying the importance of modeled mechanisms

Our experimental design for parameter exploration allowed us to estimate the relative importance of the ecological and evolutionary processes implemented in the simulation model, as they were regulated by parameters and initial conditions. The relative importance of these processes was assessed by quantifying the relative magnitude of divergence among the species richness patterns produced by the model as a consequence of experimental variation of model parameters, each of which regulates one or more of the processes implemented. To quantify the relative influence of initial conditions and parameters, we used a series of analyses of molecular variance (AMOVA) of the simulated spatial patterns in species richness.

Evaluating model performance

The exploration of parameter space was not designed to replicate the real-world diversity pattern of any extant or extinct group of species or lineages. Nonetheless, we evaluated the correspondence between the predictions of our model and contemporary, empirical patterns of species richness of birds, mammals, and plants of South America. To compare our results with published macroecological data at similar spatial resolution, and because of uncertainty in the geographic distribution of real-world species, we created a regular grid of 1659 square cells, each measuring 1° of latitude-longitude. We reprojected the maps of simulated species, from the higher-resolution grid used for simulation, into this lower-resolution grid, and recalculated spatial patterns in total, cradle, museum, and grave species richness. Because we aimed to compare predictions of our model against empirical richness patterns, we included in this analysis only the patterns in species richness emerging from the 1500 simulations that used real-world South American topography, excluding from the analysis all simulations that assumed alternative, experimental South American topographies. We used simple ordinary least-squares regression to estimate the coefficient of determination (r2) of the relationship between empirical maps of species richness (response variable) and simulated maps of species richness variables (predictor variable). [See “Contrasting empirical and simulated spatial patterns in species richness” in (95).]

Supplementary Materials

Supplementary Text

Figs. S1 to S25

Tables S1 to S10

References (130176)

Database S1

References and Notes

  1. See supplementary materials.
Acknowledgments: We thank M. B. Bush, K. J. Feeley, R. Jansson, C. N. H. McMichael, A. S. Melo, and members of the Rangel laboratory for discussions and C. Merrow for help with BIEN plant distributions. Funding: T.F.R. and J.A.F.D.-F. have been continuously supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) (grants PQ309550/2015-7, and PQ301799/2016-4). T.F.R. and R.K.C. were supported by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) (grant SWB134/2012). F.S.C. was supported by CNPq (grant DTI380.376/2017-2) and M.T.P.C. by a CAPES graduate fellowship. This project is also supported by INCT in Ecology, Evolution and Biodiversity Conservation, funded by MCTIC/CNPq (proc. 465610/2014-5) and FAPEG (proc. 201810267000023). C.R. thanks the Danish National Research Foundation for its support of the Center for Macroecology, Evolution, and Climate (grant DNRF96). Author contributions: T.F.R. and R.K.C. jointly conceived and led this study, designed the biogeographical simulation model, created the figures, and led the writing, with input from all authors. T.F.R. programmed the biogeographical model and ran the simulations. N.R.E. and P.B.H. designed, programmed, and ran the paleoclimate model. T.F.R., J.A.F.D-F, and M.T.P.C. performed statistical analyses. All authors contributed conceptually to the design of the study and interpretation of results. Competing interests: The authors declare no conflicts of interest. Data and materials availability: Paleoclimate data are available at the PANGAEA repository ( under the identification code PDI-17504, and source code for the biogeographical simulation model and statistical analysis are available at

Stay Connected to Science


Navigate This Article