A stochastic, evolutionary model for range shifts and richness on tropical elevational gradients under Quaternary glacial cycles

Robert K. Colwell, Thiago F. Rangel

Abstract

Quaternary glacial–interglacial cycles repeatedly forced thermal zones up and down the slopes of mountains, at all latitudes. Although no one doubts that these temperature cycles have left their signature on contemporary patterns of geography and phylogeny, the relative roles of ecology and evolution are not well understood, especially for the tropics. To explore key mechanisms and their interactions in the context of chance events, we constructed a geographical range-based, stochastic simulation model that incorporates speciation, anagenetic evolution, niche conservatism, range shifts and extinctions under late Quaternary temperature cycles along tropical elevational gradients. In the model, elevational patterns of species richness arise from the differential survival of founder lineages, consolidated by speciation and the inheritance of thermal niche characteristics. The model yields a surprisingly rich variety of realistic patterns of phylogeny and biogeography, including close matches to a variety of contemporary elevational richness profiles from an elevational transect in Costa Rica. Mountaintop extinctions during interglacials and lowland extinctions at glacial maxima favour mid-elevation lineages, especially under the constraints of niche conservatism. Asymmetry in temperature (greater duration of glacial than of interglacial episodes) and in lateral area (greater land area at low than at high elevations) have opposing effects on lowland extinctions and the elevational pattern of species richness in the model—and perhaps in nature, as well.

1. Introduction

From our viewpoint in the present, with anthropogenic global warming challenging the evolutionary adaptations of many species (Clark et al. 2003; Feeley et al. 2007; Deutsch et al. 2008; Sinervo et al. 2010) as the Earth's temperature nears its highest point in the past million years (Hansen et al. 2006), it is easy to forget that adaptations in contemporary lineages have also been repeatedly tested and shaped by cold during the glacial phases of the past 2.6 Ma (the Quaternary Period). In fact, because of the asymmetry of the peaks and valleys of temperature, the total time elapsed during cold glacial periods is more than twice the time that the Earth has spent in warm interglacials since the advent of dominant 100-thousand-year temperature cycles in the late Pleistocene (Rutherford & D'Hondt 2000; Liu et al. 2008; figure 1). Although continental plate movements and mountain uplift were relatively insubstantial over this geologically short period (Garzione et al. 2008), the freezing and melting of ice masses produced profound changes in shorelines, inland lakes, river courses, continental shelf islands and land bridges. Meanwhile, continental ice-sheets repeatedly altered the poleward limits and alpine glaciers the elevational bounds for terrestrial life, and regional patterns of precipitation shifted, sometimes profoundly, under the influence of temperature-driven changes in global circulation (Marchant et al. 2009). The biogeographical drama played out on this shifting physiographic and climatic template, among essentially modern lineages, continues to be revealed by stoichiometry, by pollen and phytoliths, by plant and animal macrofossils, and by phylogeographic analysis of living populations (e.g. Colinvaux et al. 1996; Schluter & Rambaut 1996; MacFadden et al. 1999; Marchant et al. 2002, 2009; Bush et al. 2004; Hooghiemstra & Van der Hammen 2004; Bush & Flenley 2007; Dick & Heuertz 2008; Blois et al. 2010; Voelker 2010).

Figure 1.

Late Quaternary temperature fluctuations. (a) The EPICA Dome C Antarctic Ice Core 800 kyr dataset (Jouzel et al. 2007). (b) The EPICA series resampled at 1605 kyr equal intervals by linear interpolation, rescaled to the estimated Quaternary range of equatorial, continental temperatures at sea level.

It is now known that the tropics were not spared (Rutherford & D'Hondt 2000; Marchant et al. 2009). Estimates for the scope of terrestrial temperature change in the tropics, at any given elevation, range from 4°C to 10°C between glacial lows and interglacial peaks (Wille et al. 2001), accompanied by complex changes in the amount and seasonality of precipitation (Bush & Silman 2004; Marchant et al. 2009) and the advance and retreat of high-elevation glaciers (Hooghiemstra & Van der Hammen 2004). Although mean annual temperature (MAT), at any given elevation, declines almost linearly with latitude outside the tropics, elevation-specific MAT varies little across tropical latitudes (Ahrens 2006), a pattern that has apparently persisted, on average, throughout the Pleistocene (MacFadden et al. 1999). By contrast, the decline of temperature with elevation (the lapse rate) is similar at all latitudes (Ahrens 2006). Thus, among temperate and boreal taxa, temperature-driven Quaternary range shifts were both latitudinal and elevational (Davis & Shaw 2001), a pattern now being repeated under current warming trends (e.g. Hickling et al. 2006; Moritz et al. 2008). By contrast, Quaternary changes in geographical distributions in the tropics were dominated by elevational range shifts and by responses to changing regional gradients in precipitation (Bush & Flenley 2007; Marchant et al. 2009). Documented upslope elevational shifts in the tropics are already emerging under current global warming trends (e.g. Chen et al. 2009).

Within the time frame of the late Quaternary, contemporary biogeographical patterns of terrestrial species can be viewed as the product of a complex interplay of ecological and evolutionary processes operating largely on pre-existing lineages, constrained by space and climate, all in the context of chance events. In this study, we aim to capture the essence of some of these processes in a simple, stochastic model that simulates both ecological and evolutionary processes affecting the elevational ranges of tropical species through late Quaternary time. Like all models, this one simplifies reality. We model tropical temperatures realistically (albeit approximately), based on an empirical temperature series (figure 1) and estimated lapse rates, but do not attempt to model precipitation (still a daunting challenge for palaeoclimates; Marchant et al. 2009), edaphic factors, habitat specialization, changes in carbon dioxide or changes in sea level. In the model, geographical ranges are represented by their projection on an elevational axis, with their lateral extent modelled as a function of elevational band area. Because temperature is the only environmental variable considered, the evolution of the thermal tolerances of each species, which determine its potential elevational range, is modelled in a one-dimensional Hutchinsonian niche space. We model the tension between thermal niche conservatism, which forces elevational range shifts in response to temperature change, and in situ adaptation to shifting temperature, which promotes range stasis. The model tracks species richness and the phylogenetic structure of assemblages over time, as lineages undergo both speciation and extinction. To assess the ability of the model to produce realistic elevational patterns of species richness, we examine the parameter values that optimize fit to the contrasting elevational patterns of species richness of four contemporary species groups that share a 2900 m elevational gradient in Costa Rica. Throughout, we emphasize the role played by stochastic events in producing enduring patterns.

2. Material and methods

(a) Model design

Our model tracks the geographical distribution of species on a bounded, tropical elevational gradient (the geographical domain) over the course of temperature cycles of the past 800 kyr (thousand years). The model is explicitly constructed upon the duality between the environmental tolerance limits of species in niche space and the projection of those limits on an environmentally heterogeneous geographical domain (Rangel et al. 2007; Colwell & Rangel 2009). In the model, niche space consists of a single environmental axis representing temperature, and the geographical domain (in the simplest version of the model) consists of a single elevational transect somewhere in the wet tropics, from sea level to the top of a mountain massif. Within this spatial domain, only elevation (the geographical z-coordinate) is modelled explicitly, and not distances along the ground (geographical x- and y-coordinates). As an option, we use the areal (x,y) extent of land within elevational bands on the massif to model its lateral extent and shape. The niche space is continuous; the geographical domain is discretized into an arbitrary (large) number of ordered, equal elevational steps. Changes in sea level are not considered, and we have not attempted to model patterns of precipitation or carbon dioxide (Marchant et al. 2009), short-term, abrupt climate shifts (Alley et al. 2003) or climate refuges based on within-elevation spatial heterogeneity (Ashcroft 2010).

The model was implemented in Delphi in the Embarcadero RAD Studio 2010 for Windows environment, as a multi-threaded, compiled application, using visual display libraries from Spatial Analysis in Macroecology (Rangel et al. 2010).

(b) Palaeotemperature data

Palaeotemperatures from 800 kyr before present until 1950 were derived from the EPICA Dome C Antarctic Ice Core dataset (Jouzel et al. 2007), which models palaeotemperature based on oxygen isotope ratios (figure 1a). Because the time-spacing of EPICA samples increases sharply with core depth (age), we resampled the original time series (5800 data points) at 160 constant intervals of 5 kyr, by linear interpolation (figure 1b), for use as the discrete time steps required in the model.

Based on estimates of tropical, terrestrial temperatures at the Last Glacial Maximum (LGM; Flenley 1998; Wille et al. 2001; Bush et al. 2004), the coldest point in the EPICA series, we assumed that sea-level temperature at the equator varied between a minimum of 23°C at the LGM and a maximum of 28°C at the previous (Eemian) interglacial, the warmest point in the EPICA series. For the model, the resampled EPICA series was proportionally rescaled to this range (figure 1b).

The temperature at each elevation on the geographical domain at each point in the resampled, rescaled EPICA series was estimated by assuming a 5°C temperature decline for each 1000 m elevation above sea level, a wet adiabatic lapse rate typical for contemporary tropical climates (e.g. Wille et al. 2001; Bush et al. 2004; Colwell et al. 2008). Although tropical lapse rates may, at times, have been somewhat greater than at present during the late Quaternary and may under some conditions be somewhat larger at high elevations (Wille et al. 2001), these differences are not sufficiently well documented or well understood to include them in our model. The scope of niche space represented on the geographical domain over the course of the time series covers the full range of temperatures, from the coldest, full-glacial mountaintop temperature (which depends upon the height of the modelled massif) to the warmest, interglacial sea-level temperature.

As a control, a flat temperature time series was constructed, the same length as the resampled EPICA series, but representing a constant sea-level temperature of 25.5°C, with corresponding elevation-specific values based on a 5°C lapse rate.

(c) Founders

In the model, before launching the time-based simulation, the geographical domain is first populated with F (a model parameter, table 1) founders, representing pre-existing (Early Pleistocene) lineages that collectively become the ancestors of all subsequent species on the domain. For each potential founder, a niche centre (temperature optimum) is selected at (uniform) random along the thermal niche axis. To allow truncated niches (Colwell & Rangel 2009; Feeley & Silman 2010), the thermal niche axis is modelled to encompass temperatures spanning three times the range of elevational temperatures expressed on the domain at Time 0 in the modelled time series. (The niche axis extends symmetrically both above and below the temperature extremes initially expressed on the domain.)

View this table:
Table 1.

Model parameters and settings.

A thermal niche breadth (in degree-Celsius) for each founder is selected from a normal distribution, with a mean at the niche centre and standard deviation σ (a model parameter, table 1). The niche extends symmetrically to each side of the niche centre. The niche is then projected onto the geographical domain (in its Time 0 climatic state) to determine the initial elevational range of the founder, regardless of whether the niche is fully expressed or, instead, truncated at sea level or mountaintop (Colwell & Rangel 2009; Feeley & Silman 2010). If the modelled niche does not correspond to any portion of the geographical domain, the species is not a successful founder, and a new random niche centre and a new random niche breadth are selected to replace it. This process, which repeats until the specified number of founders is reached, guarantees a stochastically uniform pattern of species richness across the full domain (model 1 of Colwell & Hurtt 1994), avoiding any mid-domain effect (Colwell & Lees 2000) among the founders. Once the domain has been populated with founders and environmental changes begin, the niche centre, niche breadth and distribution of all species on the domain are completely under the control of the model dynamics, and the initial niche breadth parameter σ no longer plays any role.

(d) Anagenetic (within-species) adaptive evolution and niche conservatism

Once the simulation is launched, the temperature at each elevation on the domain changes at each time step in accord with the tropically rescaled EPICA time series (figure 1b). At each step in the simulation, following the change in temperature, each living species is evaluated for discordance between its temperature optimum in niche space (before the environmental temperature change) and the mean temperature of its distribution on the domain after the temperature change.

When niche and distribution are discordant, selection favours an adaptive evolutionary shift in the niche that will decrease or eliminate the difference between the niche optimum and the mean environmental conditions (Smith et al. 1995; Kirkpatrick & Barton 1997; Davis et al. 2005). On the other hand, an upslope or downslope range shift may also bring the geographical distribution into accord with the niche, without any evolution—but only if the range is not already too close to sea level (in a cooling climate) or mountaintop (in a warming climate; Colwell & Rangel 2009). The model assumes that range shifts into contiguous, thermally suitable portions of the domain are always successful. Sometimes, however, following environmental temperature change in a single time step, the new potential range of a species with narrow thermal tolerances maps on the domain outside its range in the previous time step. As detailed in §2e, these ‘range shift gaps’ (Colwell et al. 2008) are assumed to be insurmountable, and extinction ensues.

When ranges track climatic zones with little or no adaptation, the range shifts may be described as evidence for niche conservatism (Martinez-Meyer et al. 2004; Losos 2008). Based on an approach introduced by Rangel et al. (2007), two niche-evolution parameters control how much adaptive anagenetic evolution is permitted (or conservatism enforced) between time steps in the model. If D1 is the distance in niche space (°C) between the old niche optimum and the current mean environmental temperature of a species' realized range, and D2 is the distance between the old niche optimum and the new optimum after selection, then D2 = D1 × Δ. Thus, for perfect niche conservatism, model parameter Δ = 0, whereas for perfect adaptation, Δ = 1, although Δ may take any value between these limits (table 1). A second model parameter, β, regulates evolution or conservatism in niche breadth: W2 = W1 × |N (1, β)|, where W1 is a species' niche breadth before and W2 is its niche breadth after evolutionary change, and N (1, β) is a normal distribution with mean 1 and variance β. Thus, for perfect niche conservatism, parameter β = 0, whereas values of β > 0 allow stochastic evolutionary change in niche breadth, which may prove adaptive (table 1).

Once these evolutionary adjustments in each species' niche have been made, the niche is freshly projected onto the domain to establish the realized geographical range. If little or no evolution has been permitted, the new range will have shifted to track the corresponding upslope or downslope shift of climatic thermal zones, to the degree possible given the bounds of the spatial domain.

(e) Extinction

The model incorporates two kinds of extinction. After each step of temperature change and subsequent niche evolution (as described above), each newly re-projected range is considered for stochastic extinction, based on its elevational range size as a proxy for population size, which is commonly considered the primary risk factor for such extinctions (Purvis et al. 2000). The relationship between range size and probability of extinction is governed by model parameter α (table 1) by evaluating Pext = −ln(r/v)(1/α), where Pext is the probability of extinction, r is the range size and v is the full extent of the domain (figure 2; Rangel et al. 2007). Variables r and v are measured either in metres of elevation or, as an option, as the areal extent of land within elevational bands bounded by the species' elevational limits (r) and the sum of all elevational areal bands on the massif (v). To explore the effect of lateral area on model outcomes, we used a right cone with base radius equal to its height (approximating a stratovolcanic cone), as well as an empirical mountain gradient in Costa Rica, as described in a later section.

Figure 2.

The extinction function. Based on a single parameter, α, the function scales probability of stochastic extinction, per species per time step, to range size r and geographical domain size v (Rangel et al. 2007).

In contrast with stochastic extinctions, deterministic extinctions occur when temperature change during a single time step of the model leaves a species without any portion of the domain that corresponds with the scope of its thermal niche (Williams et al. 2007), leaving no population to evolve. Three kinds of deterministic extinctions occur in the model: mountaintop extinctions during warming phases, sea-level extinctions during cooling phases and range-shift gap extinctions (Colwell et al. 2008) that befall species at mid-elevations (described in §2d). Narrow ranges (initially as a result of small σ, but ranges may become stochastically narrow as well) predispose species to deterministic distinctions as well as stochastic ones. Note that the number of time steps (160 in our model) chosen for the EPICA resampling affects the rate of deterministic extinctions as well, by controlling the mean shift in climate in an average step—perhaps, a reasonable proxy for the rate of environmental change in relation to biological response.

(f) Speciation: cladogenetic adaptive evolution and niche conservatism

Speciation in the model is conceived as allopatric or parapatric, with the niche and distribution of a child species diverging, stochastically, from its parent species, which itself remains initially unchanged. The child species is envisioned as inhabiting a neighbouring but not overlapping region of the massif, perhaps an adjacent watershed (Jansson & Dynesius 2002; Fjeldså & Rahbek 2006). A parent species is chosen at random from among species that have survived extinction in the present time step, with probability γ (a model parameter, table 1). The degree to which the thermal niche, and therefore the elevational distribution, of the child species resembles its parent's niche and distribution after cladogenesis is governed just as for anagenetic evolution and range shifts between time steps. The child's niche evolves towards the optimum for its realized range to the degree permitted by the niche centre conservatism/evolution parameter Δ and expands or contracts stochastically to the degree allowed by the niche breadth conservatism/evolution parameter β.

As an option, a ceiling of K species may be imposed on the model, imposing a degree of diffuse interspecific competition and forcing a degree of lineage selection on the model (and as a practical matter, preventing the runaway accumulation of species). If a ceiling has been set, once it is reached, speciation occurs only to the degree that extinction reduces the total number of living species below K.

(g) Model outputs and behaviour

A model run (realization) begins at EPICA Time 0 (800 kyr ago) and ends at Time 160 (1950), or earlier if all lineages go extinct. The programme maintains a species ‘vital statistics’ table, recording for each species its time of origination (0 for founders), time and cause (deterministic or stochastic) of extinction, parent species (if not a founder), number of descendant species, and initial and final values for niche size and breadth (°C) and range size (in m elevation). Biogeographic data, consisting of the presence or absence of each species at each time step and at each elevation, yield a time versus elevation plot of species richness (figures 3d and 4), as well as time-based trends in patterns of range size and elevational location (figure 5).

Figure 3.

Visualization of results for a single realization of the model driven by the EPICA temperature series. (a) Cladograms for the 10 founders and their descendant species (each clade a different colour). Vertical location of clades in the figure is arbitrary. (b) Extinctions for each 20 kyr interval for the clades in (a). Note that more extinctions occur during cold phases than during interglacials. (c) The resampled EPICA temperature series (from figure 1b). (d) Time-elevation map of species richness for the clades in (a). Cool colours indicate lower richness, hot colours higher richness (0–100 species). (e) Time-elevation map of the phylogenetic composition of local assemblages. Hot colours indicate significant (p < 0.025, red) or nearly significant (p < 0.05, orange) phylogenetic clustering, whereas dark blue indicates phlyogenetic mixing (‘overdispersion’), although no assemblages quite reach p > 0.975 in this example. The ‘average taxonomic distinctness’ (Warwick & Clarke 1998) of assemblages was computed for each elevation and time, based on the cladogram in (a). Significance was assessed by a permutation test (500 permutations) for each cell, based on the pool of all living species at that time interval. See table 1 for parameters and settings for this replicate.

Figure 4.

The signature of chance in time-elevation maps of species richness. All nine examples share exactly the same values of all settings and parameters except for the extinction parameter, α. The three examples within each column, however, are true replicates, sharing the same value of α as well as all other parameters. See table 1 for other parameters and settings.

Figure 5.

Elevational range shifts as a function of niche conservatism/evolution. (a) The three time-elevation maps show representative richness patterns for three values of Δ, the niche centre evolution parameter (table 1). (b) The small graph below each richness map plots, for that example, the mean of the elevational range midpoints for all living species as a function of sea-level temperature, through the full course of the EPICA glacial cycles (figure 1). Strong niche conservatism (Δ = 0.1) forces ranges to shift upslope during warming phases and downslope during cooling phases, whereas rapid niche evolution (Δ = 0.9) permits ranges to remain relatively fixed on the elevation gradient regardless of temperature. In (c), the slopes of ordinary least-squares regressions (the lines in (b)) for 10 replicates at each of seven values of parameter Δ (spanning the range between 0.01 and 1.0) are plotted as a function of Δ, demonstrating a linear decrease in elevational range shift as niche conservatism weakens. See table 1 for other parameters and settings.

From parent–child relationships and times of origination and speciation, a fully dichotomous phylogeny is recorded (in Newick notation) for each founder's descendants, with time-based branch lengths (figure 3a). On the basis of phylogeny and the modelled biogeographic data, the phylogenetic structure of assemblages at each time step is assessed for each elevation (figure 3e) based on the ‘average taxonomic distinctness’ (Warwick & Clarke 1998) among species at each elevation at each time step (see figure caption for details). This pattern can be compared with the corresponding pattern of richness (figure 3d) as an indication of lineage selection and lineage mixing.

We combined visualization and statistical tools to explore the behaviour or the model and its parameters (table 1). Simple statistical procedures are documented in context in §3 or in figure captions. The more complex methods we used for pattern matching are explained in §2h.

(h) Target pattern matching and parameter assessment

To assess the potential of the model to reproduce empirical richness patterns on a tropical elevational gradient, we used four contrasting elevational records of species richness as ‘target patterns’ and systematically analysed the model parameter space to maximize the correspondence between observed and modelled patterns. The empirical datasets represented four ecologically and phylogenetically distinct groups of organisms on the same gradient, the Barva Transect in Costa Rica (tables S1–S4 in Colwell et al. 2008): 568 species of epiphytes (Cardelús et al. 2006), 82 species of understorey Rubiaceae (Gilman 2007), 739 species of geometrid moths (Brehm et al. 2007) and 495 species of ants (Longino et al. 2002).

For all Barva Transect model runs, a maximum elevation of 2900 m was set, to match the actual scope of the transect, with the domain divided into 140 steps of 20 m elevation each, with a lapse rate of 5°C per 1000 m elevation (electronic supplementary material and methods in Colwell et al. 2008). To assess the effect of lateral area on extinction and patterns of richness, we applied an estimate of the land area in each 100 m elevational band on the Atlantic slope of Costa Rica (Gilman 2007) to the extinction function (figure 2).

For each of the four target richness patterns, a parameter space comprising nearly 50 000 combinations of parameter values was considered and formally explored, with a set of 10 replicate simulations (a ‘replicate set’) per parameter value combination (table 2). The results for different replicates within a replicate set vary among replicates only as result of the stochastic elements of the model. We did not attempt to match the absolute number of species in the target patterns, only the shape of the distribution (figure 7). Because it is sensitive to all aspects of the shape of distributions, but not to their absolute magnitude, we minimized the Kolmogorov–Smirnov (K–S) statistic (Sokal & Rohlf 1995) as a criterion of fit.

View this table:
Table 2.

Parameter space and model selection for Barva Transect richness simulation. (Parameter values that appeared in the most frequent parameter sets among the best-fitting 1% of replicates are indicated by x or X. Parameter sets for the examples in figure 7 are marked X. Entries in the ANOVA column indicate factor importance averaged across groups. Preliminary explorations having shown that the probability of speciation, as long as not zero or extremely high, matters little to the outcome, γ was fixed at 0.06 for all parameter sets. Model settings as for figure 3 (see table 1)).

Separately for each target group of organisms, a multi-way analysis of variance (ANOVA) was carried out on the full set of simulation runs, with K–S values as the response variable and 10 runs under identical parameter settings as replicates, to determine which parameters were most important in distinguishing good fit from poor for each target pattern, given the parameter combinations considered. Each parameter (table 1) was treated as a factor and its settings as levels. Standardized F-statistics were compared among parameters and target groups, but no significance tests were appropriate or performed for the ANOVA.

To identify the parameter value combinations that most accurately reproduced each target pattern, the model outcomes for all replicates (for all parameter values) were first sorted by the K–S statistic, for each target pattern separately. The best-fitting (smallest K–S) 1 per cent were selected for further analysis. Next, to help eliminate chance successes, this subset was searched for the most frequently represented parameter combinations (those closest to 10 replicates per combination) and the parameter settings for these sets were compared among the four target groups. No significance tests were appropriate or performed for the K–S statistics.

3. Results

(a) Visualization of patterns over time and elevation

Each run (realization) of the model simulates the evolutionary, ecological and biogeographic fate of one or more founders and all descendant species over the 800 kyr of the EPICA temperature series (figure 1), under some particular fixed value of each of the model parameters and settings (table 1). Another run using exactly the same values of the parameters is considered a ‘replicate’ of this biogeographical experiment. We refer to groups of runs that share the same parameter values and model settings as ‘replicate sets’. Only by running the model repeatedly under the same parameter settings can we distinguish between the roles of chance (stochasticity) and mechanism. Figure 3 illustrates, in detail, the results for a single replicate.

A key visualization tool for this study is the time versus elevation map, which plots the magnitude of richness for each discrete time step (160 for the EPICA series) and each elevational step (we used 100–140). For example, figure 3d plots local species richness, at each time and elevation, as a function of time (x-axis), and elevation (y-axis) under EPICA temperature cycles (figure 3c), for a single realization of the model. Species richness, which is simply the number of overlapping species ranges at a particular time and elevation, is indexed by colour according to the scale below the time-elevation map (figure 3d).

(b) Stochastic variation in pattern among replicates

The stochastic elements of the model (table 1) are limited to (i) initial location of founder niche centres (uniform random), (ii) the magnitude of initial niche breadths (normal distribution tuned by parameter σ), (iii) the niche breadths of child species (normally distributed differences from parent niche breadth, tuned by parameter β), (iv) the identity of species chosen to speciate (uniform random with probability γ), and (v) the identity of species chosen to undergo stochastic extinction (weighted by range size and tuned by parameter α). All other processes and algorithms in the model are fully deterministic, given the prior settings (table 1) of the model. Nevertheless, temporal and spatial (elevational) variation in the pattern of species richness, range size, range location and local phylogenetic structure among replicates sharing the same parameter values proved to be remarkably rich.

Figure 4 illustrates the striking stochastic variation, among replicates within replicate sets, in species richness patterns as a function of time and elevation, for three levels of the extinction parameter α. The three replicates within each column share all parameter values. These examples show how the stochastic events of history become consolidated by speciation and the inheritance of niche characteristics.

(c) An emergent pattern: mid-elevation richness peaks

Despite the variation in spatial and temporal features of richness and ranges among replicates and between parameter settings, consistent patterns do emerge. The most striking is the tendency for richness to reach its maximum near mid-elevations, or at least to decline towards the boundaries of the domain, a pattern frequently seen along elevational transects in nature (Rahbek 2005).

When two or more founders initiate the simulation and a species ceiling is imposed, under most parameter settings subsequent biogeographic and phylogenetic history shows a pattern of domination by one or a very few clades, while others become extinct. The surviving clade or clades tend to be those that either began or (if parameters permit) have shifted towards mid-elevations, away from the extinction perils of the bounds of the domain. For example, in figure 3 (which is typical in this regard), of 10 founders, only two produced lineages that survived to the end of the simulation (figure 3a). Both started out at moderate elevations: 750 m and 1680 m on a 2900 m domain. The phylogenetic composition of local assemblages routinely shows the increasing dominance of the surviving clades as time passes and other clades become extinct. In the example, figure 3e plots the significance level of phylogenetic clustering or mixing in the assemblage at each elevation at each time. Following an early period of coexistence with other clades (see the cladogram in figure 3a), the two surviving clades come to dominate different elevations (red areas of figure 3e), briefly mixing at intermediate elevations (blue areas), eventually broadly overlapping (increasing yellow and green areas) to produce the final richness pattern. A formal, quantitative analysis of such repeated patterns is beyond the scope of the present study, but they strongly suggest lineage-level selection for biogeographic location on the domain.

Moreover, although the EPICA temperature fluctuations certainly strengthen this effect, we found that stochastic extinctions near the domain boundaries, where ranges tend to be smaller, can produce mid-elevation richness peaks even with a flat temperature series replacing EPICA (maintaining the usual lapse-based elevational temperature gradient). Rangel & Diniz-Filho (2005; further discussed by Colwell & Rangel 2009) showed a similar result for a related, niche-based stochastic model. Previous stochastic models that incorporate speciation within a bounded domain have all yielded mid-domain peaks of richness (Bokma et al. 2001; Davies et al. 2005; Arita & Vazquez-Dominguez 2008; Connolly 2009). The differential extinction of small populations near domain boundaries was first proposed by Colwell & Hurtt (1994) as a biological explanation for the mid-domain effect (the tendency for larger ranges to overlap towards the middle of a bounded, spatial domain owing to geometric constraints on range location).

(d) The effect of niche conservatism versus niche evolution on range shifts

In the model, in the complete absence of anagenetic (within-species) niche evolution (Δ = 0 and β = 0; table 1), species ranges are forced to shift downslope during cooling phases and upslope during warming phases in lockstep response to the constantly changing temperatures of the EPICA glacial–interglacial cycles. At the other extreme, if a species is capable of rapid evolutionary response to directional selection imposed by discordance between its niche optimum and the current mean temperature within its elevational range (Δ = 1 and β ≥ 0), its range centre remains fixed on the gradient. Figure 5 shows that the model yields a linear decrease in the magnitude of elevational range shift as niche conservatism weakens, for values of Δ spanning the range between 0.01 and 1.0, confirming that the model behaves as expected.

(e) The effects of niche conservatism, temperature asymmetry and lateral area on extinction

Extinctions are more likely with strong niche conservatism, for two reasons. First, if little or no adaptive evolution is possible, temperature change may leave a species without any region of the domain that corresponds to its thermal niche (deterministic extinction). Second, species with ranges abutting either boundary of the domain (sea level or mountaintop) may express only a portion of their fundamental niches, guaranteeing discordance between niche optimum and mean environmental conditions in the expressed range. If selection cannot correct this discordance, the range remains small, and is thus susceptible to stochastic extinction. At extremes of the glacial–interglacial cycles, both lowland and mountaintop species in the model are therefore more vulnerable to both kinds of extinction than species with mid-elevation ranges, especially so when niche conservatism is strong. However, two opposing, well-known asymmetries complicate this simple picture.

The first asymmetry lies in the long-term pattern of temperatures during the late Quaternary. For the resampled EPICA series (figure 1b), the temperature lies at or below the middle of its long-term range for 73 per cent of the 5 kyr intervals. At the extremes, the asymmetry is even more striking: 33 per cent of the time intervals lie within the coldest 20 per cent of the EPICA temperature range—the lingering glacial episodes, whereas only 6 per cent of intervals reach the warmest 20 per cent of the temperature range—the relatively brief spikes of the interglacials. Thus, even if extinctions are random across time intervals, far more extinctions would be expected to occur during glacial episodes than during interglacials. The pattern of extinction shown in figure 3b illustrates this asymmetry, which would be expected, over time, to impose stronger selection for cold tolerance among lowland lineages than for heat tolerance among mountaintop lineages.

The second asymmetry arises from the shape of mountain massifs, which generally (but not always) offer decreasing land surface area with increasing elevation, yielding smaller populations and thus greater vulnerability to stochastic extinction at higher elevations. We modelled the effect of lateral area on extinction by treating the domain as a cone, reckoning species ranges as the surface area of an elevational band, as a proxy for population size. The net result of the opposing effects on extinction of these two asymmetries is thus a quantitative matter, not predictable from first principles. Under the conditions of the model, figure 6 shows that increasing extinction pressure lowers the elevation of the centre of mass for species richness, evidently overcoming the contrary effect of greater numbers of extinctions in the lowlands during glacial maxima. Figure 4 illustrates a variety of examples from these simulations.

Figure 6.

The effect of lateral surface area on the elevational pattern of richness. The elevational centre of mass for species richness is plotted for 10 replicate simulations for each of four extinction intensities (parameter α, table 1). The extinction function (figure 2) was modelled for a hypothetical mountain shaped as a right cone, 5000 m high with a 5000 m base radius. Species elevational ranges were transformed to the corresponding area (r in the extinction function) of a band on the surface of the cone, as a proxy for relative population size, with the area of the full domain (v in the function) equal to the entire surface of the cone. With increasing extinction pressure, lower-elevation lineages are increasingly favoured by reason of their larger population sizes. The slope of the regression line, which passes the test for linearity, is significant at p = 0.003 (Sokal & Rohlf 1995, p. 476). See table 1 for other parameters and settings.

(f) Target pattern matching for Barva Transect plants, moths and ants

The empirical patterns of species richness with elevation for four groups of organisms on the Barva Transect (blue profiles in the graphs on the right of figure 7) vary greatly in shape and in the elevation at which species richness peaks, but all four richness patterns can be reasonably well reproduced from the model (red profiles in figure 7). Not surprisingly, given the role of chance in the richness patterns produced by the model (as illustrated in figure 4 for a different analysis), even within a particular replicate set, some replicates fit better than others. Nevertheless, patterns of success and failure do emerge, and the range of parameter settings that yield the best fit differs among the four groups for some parameters (table 2), but not for others. We endeavoured to minimize the hazard of finding one-off fits by considering only parameter sets that repeatedly produced a good fit to the empirical patterns (based on the K–S statistic). Moreover, the examples illustrated in figure 7 were not selected from the realizations of the model used to define ‘best sets’ of parameter values, but were produced subsequently by independent runs based on one or more of those sets. In effect, the procedure can be viewed as an informal test of the hypothesis that the empirical profiles could have been drawn from the distribution of richness patterns produced by the model.

Figure 7.

Modelling target richness patterns for four groups of species on the Barva Transect, Costa Rica. The time-elevation maps illustrate the modelled temporal history of richness for each simulation, on a colour scale from 0 (blue) to 100 (red) species. The blue curves on the right show the smoothed empirical profile of richness with elevation for four groups of species, and the red curves show the modelled richness profile at the end (right edge of the richness map) of the simulation, on a relative scale. See the text for methods and table 2 for parameter settings for each group. (a) geometrid moths, (b) epiphytes, (c) Rubiaceae and (d) ants.

As table 2 shows, a surprisingly high level of the niche centre evolution parameter, Δ (0.7 or 0.9, where 1.0 represents full adaptive evolution at each time step), proved to be crucial for a good fit of modelled to empirical richness for all groups. By contrast, a good fit called for relatively low levels of the niche breadth evolution parameter, β, although a somewhat higher level of this parameter proved optimal for the broad, centred distribution of the moths. High rates of extinction proved to be optimal for all groups, perhaps serving to enhance the rate of evolution of niche centres under the constraint of the species ceiling (K). Groups differed in the optimal level of founder number (few for the low-elevation ants, many for the broad, convex richness pattern of the epiphytes) and lateral area effect (which tends to produce lower elevation richness peaks, figure 6). The levels for speciation rate (γ, results not shown) and founder niche breadth (σ) proved to matter little, as long as not zero.

4. Discussion

Integrating determinism with chance and ecology with evolution, we have attempted to model a brief moment in the history of the Earth for a simplified biogeographical gradient, based on a single environmental factor. By simplifying both niche space and geographical space, we were able to model and visualize environmental effects on species distributions on a geological time scale (figure 3), taking ‘general simulation models’ in macroecology in a new direction (Gotelli et al. 2009). Although we incorporated lateral extent into the extinction function by taking elevational band area into account, a more explicit accounting for ‘x, y space’ (where elevation is the z dimension) would clearly be desirable. As for the environmental side of the model, the most significant limitation of our model is its lack of accounting for patterns of precipitation. While no species can escape its limits of thermal tolerance, those limits are of course often modified by water availability, particularly for plants.

To us, the most surprising result was not that the constraints and mechanisms built into the model succeeded in producing reasonable patterns of phylogeny and macroecology, but the remarkably rich variation among replicate realizations under identical parameters, as speciation consolidates the differential survival of lineages and inherited ecological tolerances perpetuate their location in space (figure 3). This result, which is unlikely to be changed in a more complete model, suggests that the search for mechanistic explanations for macroecological patterns may inevitably be constrained by the idiosyncratic nature of historical events, perhaps to a greater degree than many of us have appreciated.

Our attempt to reproduce the empirical elevational richness profiles of four groups of organisms on the Barva Transect within the parameter space of the model (figure 7 and table 2) is no more than an exploratory first step towards mechanistic modelling of empirical patterns on tropical elevational gradients. Ideally, parameter-set selection should target multiple characteristics (Grimm et al. 2005), including range-size frequency distributions, diversity and dispersion fields, characteristics of the phylogeny and the phylogenetic structure of local assemblages (Gotelli et al. 2009).

Consistently high rates of niche centre evolution combined with low rates of niche breadth evolution were required for a good fit of modelled to empirical richness profiles, together confining each founder clade to a relatively stable range of elevations, in spite of the repeated temperature cycles. Thus, two founders were more likely to produce a good fit for ants (a single insect family), with their relatively narrow lowland richness peak, whereas 50 founders did a better job for epiphytes (50 plant families; Cardelús et al. 2006), with their broad, convex richness profile. Clearly, more than two lineages of ants and more than 50 lineages of epiphytes have contributed to the contemporary biotas of the Costa Rican Atlantic slope, but the qualitative difference in phlyogenetic, niche and functional diversity may be realistic.

Our results suggest that niche conservatism for temperature optimum over the time scale of the Quaternary might not be a winning strategy, contrary to other views (e.g. Martinez-Meyer et al. 2004; Kozak & Wiens 2010). Experimental and genetic studies of tropical montane species may reveal unexpected levels of intraspecific and intraclade genetic variation within and between elevations (Gilman 2007), just as within-species adaptation to local conditions may play a substantial role in latitudinal range shifts (Davis & Shaw 2001). In fact, although we have consistently referred to the taxonomic units of modelled clades (figure 2), as ‘species’, we recognize that speciation in most groups of organisms has not proceeded at such a rapid pace within the relatively short span of the Quaternary, although there are exceptions (Richardson et al. 2001; Hooghiemstra & Van der Hammen 2004; Kay et al. 2005; Fjeldså & Rahbek 2006). However, the same units may instead be viewed as clades at an infraspecific level, with the caveat that, in the model, lineages do not rejoin.

The opposing effects on lowland extinctions of temporal asymmetry in temperature (greater duration of glacial than of interglacial episodes) and spatial asymmetry in lateral area (greater land area at low than at high elevations) have apparently not previously been recognized. Our results, for a hypothetical conical mountain (of a particular height and slope; figure 6), indicate that area comes to the rescue, yielding richness peaks below mid-elevation, but of course that result may be model-dependent and is certainly dependent on the actual elevational profile of particular mountain massifs. We conjecture that, in general, the steeper and more isolated the mountain, the weaker the area effect in countering lowland extinctions during glacial episodes, forcing the richness profile upslope, compared with patterns on broader mountain massifs. In any case, in terms of selection for traits that allow species to tolerate extreme temperatures, the inescapable inference appears to be that the Quaternary has imposed far more continuous selection for cold-tolerance than for heat-tolerance, at all elevations (Colwell et al. 2008). If, indeed, tropical species are better adapted to glacial cold than to interglacial heat, these opposing asymmetries might help explain the ‘out of the tropics’ (Jablonski et al. 2006) Pleistocene expansion of certain lineages towards higher latitudes. Meanwhile, this inference is not good news for tropical species under the current anthropogenic global warming trend, which comes at an unfortunate point in the orbitally forced cycles of the Quaternary.

Acknowledgements

The comments of R. Burnham, D. Nogués-Bravo, G. Brehm, J. Fjeldså, C. Rahbek and two anonymous reviewers greatly improved the manuscript. We are grateful to G. Brehm, C. Cardelús, A. Gilman and J. Longino for sharing their Barva Transect data. We acknowledge the Danish National Research Foundation for support to the Center for Macroecology, Evolution and Climate, University of Copenhagen (Denmark), the National Science Foundation (USA, DEB-0639979 and DBI-0851245) and the National Center for Ecological Analysis and Synthesis (NCEAS, USA). T.F.R. was supported by the Department of Ecology and Evolutionary Biology, University of Connecticut (USA), the Coordenação de Aperfeiçoamento de Pessoal de Nivel Superior/Fulbright Fellowship (CAPES, Brazil); and the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq, Brazil).

Footnotes

References

View Abstract