Mismatches between demographic niches and geographic distributions are strongest in poorly dispersed and highly persistent plant species.

The ecological niche of a species describes the variation in population growth rates along environmental gradients that drives geographic range dynamics. Niches are thus central for understanding and forecasting species' geographic distributions. However, theory predicts that migration limitation, source-sink dynamics, and time-lagged local extinction can cause mismatches between niches and geographic distributions. It is still unclear how relevant these niche-distribution mismatches are for biodiversity dynamics and how they depend on species life-history traits. This is mainly due to a lack of the comprehensive, range-wide demographic data needed to directly infer ecological niches for multiple species. Here we quantify niches from extensive demographic measurements along environmental gradients across the geographic ranges of 26 plant species (Proteaceae; South Africa). We then test whether life history explains variation in species' niches and niche-distribution mismatches. Niches are generally wider for species with high seed dispersal or persistence abilities. Life-history traits also explain the considerable interspecific variation in niche-distribution mismatches: poorer dispersers are absent from larger parts of their potential geographic ranges, whereas species with higher persistence ability more frequently occupy environments outside their ecological niche. Our study thus identifies major demographic and functional determinants of species' niches and geographic distributions. It highlights that the inference of ecological niches from geographical distributions is most problematic for poorly dispersed and highly persistent species. We conclude that the direct quantification of ecological niches from demographic responses to environmental variation is a crucial step toward a better predictive understanding of biodiversity dynamics under environmental change.

In 1957, George Evelyn Hutchinson introduced his seminal concept of a species’ ecological niche (1). The Hutchinsonian niche is defined as the set of environmental conditions for which demographic rates result in a positive intrinsic population growth rate and thus permit a species to form self-perpetuating populations (2, 3). This niche concept has received much attention as a theoretical foundation for explaining the geographic distributions of species and forecasting range shifts under environmental change. However, this use of the Hutchinsonian niche concept has been critically revisited in recent years (4, 5). Theoretical models of range dynamics predict that the geographic distribution of a species does not necessarily match the geographic projection of the species’ niche. This is because geographic distributions are structured by extinction and colonization events that arise from a dynamic interplay of spatial variation in demographic rates, local population persistence, and dispersal (6). Specifically, migration limitation can prevent a species from colonizing parts of its potentially suitable range, source–sink dynamics can sustain populations by immigration even if local population growth rates are negative, and time-lagged local extinction can cause species to occur in locations that became unsuitable due to environmental change (4–6). Strong mismatches between niches and geographic distributions can severely bias niche estimates and forecasts of species distribution models (7, 8), which are widely used in global change–biodiversity assessments (9). However, the extent of these niche–distribution mismatches is poorly known, mainly because of a lack of the comprehensive, range-wide demographic data needed to directly infer ecological niches (5, 10–12). Demography-based niches were so far only quantified for a few single species (13–15), which precluded comparative analyses. More than 60 y after Hutchinson introduced his niche concept, it is thus unclear how relevant niche–distribution mismatches are in the real world and how they depend on the dispersal and persistence ability of species.

Here we quantify niches from extensive demographic data for 26 closely related plant species to analyze how niche sizes, geographic range sizes, and, finally, niche–distribution mismatches depend on species’ life-history traits. Our study species are shrubs of the Proteaceae family endemic to the Cape Floristic Region, a global biodiversity hot spot (16). All study species are serotinous: they store their seeds over multiple years in a canopy seedbank until fire triggers seed release, wind-driven seed dispersal, and subsequent establishment of new recruits (17). Fire is also the predominant cause of mortality of established adults (18). This fire-linked life cycle allows the efficient measurement of key demographic rates in single visits to each population (19) and enabled us to collect data that are informative of variation in population growth rates across each species’ geographic range.

Results and Discussion

We analyzed a total of 3,617 population-level records of fecundity, recruitment, and adult fire survival (). The analyses used hierarchical demographic response models that describe for each key demographic rate the species-specific response curves to variation in climatic-edaphic conditions (minimum winter temperature, maximum summer temperature, summer aridity, and soil fertility), fire return intervals, and intraspecific density (Fig. 1 and see , for all study species). The fitted response curves explained much variation in long-term fecundity (Nagelkerke’s R2N 0.29 to 0.91 across species, mean = 0.58; ), recruitment (R2N = 0.05 to 0.84, mean = 0.42), and adult fire survival (R2N = 0.41 to 0.83, mean = 0.62). The response curves of all three demographic rates were then integrated in models of density-dependent local population dynamics that predict variation of intrinsic (low-density) population growth rates (r0) along environmental gradients. This delimits each species’ niche as a hypervolume in environmental space for which r0 is positive (Fig. 1). Within species, different demographic rates often respond similarly to the same environmental variable. However, we also found some opposing responses, in particular between fecundity and recruitment rates (for an example, see responses to aridity in Fig. 1). Such opposing responses indicate demographic compensation that can broaden species’ environmental niches and buffer effects of environmental change (15, 20, 21).

Fig. 1.

The Hutchinsonian niche and geographic distribution of Protea longifolia. (A) Responses of key demographic rates and the resulting annual intrinsic population growth rate (r0) to variation in minimum winter temperature (Tmin), maximum summer temperature (Tmax), indices of summer aridity and soil fertility, and fire return interval. (B) A projection of the niche hypervolume into a 3D environmental subspace (gray) delimits the conditions for which r0 > 0. The marginal 2D heat maps show the predicted r0 when all other niche axes are set to their respective optima. (C) Geographic projection of r0 across the Fynbos biome (colored areas) in comparison to the natural geographic range (dashed line) and to populations established outside the natural range (crosses). (D) Enlarged map showing presence records of natural populations (open circles) and demographic sampling sites (green circles). Model predictions in all subplots are the medians of the respective Bayesian posterior distributions.

Relationship Between Life-History Trait Effects and Niche Sizes.

We first examined how the size of the estimated niches is related to the dispersal and persistence abilities of species. Niche sizes were quantified separately for the 4D environmental niche defined by spatially varying long-term averages of climatic-edaphic conditions and for the disturbance niche defined by fire return intervals that also show strong temporal variation in any given location (22). Dispersal ability was quantified from species-specific parameterizations of a trait-based, mechanistic model of wind-driven seed dispersal (23). Persistence ability was characterized by resprouting as a key functional trait: some of the study species (resprouters, n = 7) possess fire-protected meristems from which individuals can resprout and are thus more likely to survive fire than individuals of species lacking this trait (nonsprouters, n = 19) (18, 19). Since resprouter populations do not exclusively rely on successful reproduction in each fire cycle, they are expected to be less vulnerable to short fire return intervals that prevent the buildup of canopy seed banks (18). We indeed found that resprouting ability had a clear positive effect on adult fire survival rates and disturbance niche size (Fig. 2 and ), whereas dispersal ability had no effect. In contrast, environmental niche size showed a strong positive relationship with dispersal ability (Fig. 2). This finding is consistent with a scenario of correlational selection on niche size and dispersal, where narrower environmental niches select for lower dispersal distances and vice versa (24, 25).

Fig. 2.

Life-history trait effects on niche sizes. (A) Effect of persistence ability on disturbance niche size. (B) Effect of dispersal ability on environmental niche size (points: posterior means; bars: posterior SDs). The line shows the estimated linear regression (posterior means, 90% credibility interval as shaded areas, slope = 1.00, P = 0.006).

Mismatches Between Demographic Niches and Geographic Distributions.

For each species, we projected niches from environmental space into geographic space (Fig. 1 and see , for all study species) and then compared this potential geographic range (the region of demographic suitability where predicted r0 > 0) to independent and extensive distribution records (17). Several species showed a remarkably strong agreement between potential ranges and observed geographic distributions (area under the receiver operating characteristic curve [AUC] values >0.8 for 10 of the 26 study species). However, there also was a high variation in this agreement across species (AUC 0.55 to 0.97, mean = 0.77; ), indicating interspecific variation in mismatches between demographic suitability and geographic distributions. Since processes that can generate these mismatches are expected to act on different spatial scales (26), we further analyzed the relationship between demographic suitability and species occurrence separately at large and small spatial scales.

On large spatial scales, dispersal limitation can cause an incomplete filling of potential ranges since species are unable to reach suitable areas (4–6). We thus tested for a positive relationship between dispersal ability and range filling. Range filling was measured as the proportion of the potential range that is covered by a species’ geographic range (the alpha-convex hull encompassing all natural occurrences; Fig. 1). We found range filling to vary widely across species (0.08 to 0.95, mean = 0.47; ) and to be unrelated to species’ persistence ability (). However, as expected, range filling strongly increased with species’ dispersal ability (Fig. 3; for additional analyses showing that this finding is robust to potential bias in predicted r0 and that species distribution models estimate higher degrees of range filling and a weaker, nonsignificant, relationship between range filling and dispersal ability, see ). Hence, good dispersers not only have larger environmental niches (Fig. 2) and thus tend to have larger potential geographic ranges, but they also fill more of their potential ranges, so that both factors add up to explain the larger geographic ranges of good dispersers ().

Fig. 3.

Life-history trait effects on the mismatch between niches and geographic distributions. (A) Effect of dispersal ability on range filling (points: posterior means; bars: posterior SDs). The line shows the estimated linear regression (posterior means: 90% credibility interval as shaded areas, slope = 0.58, P = 0.022). (B) Relationship between demographic suitability (predicted r0) and occupancy within the range. Points show the mean occupancy in sites that were binned according to deciles of predicted r0 (i.e., 10 points per species). The lines show average predictions of this relationship for species with different persistence ability (posterior means: 90% credibility interval as shaded areas). (C) Variation in species’ mean occupancy of sites within their ranges that are predicted to be unsuitable (r0 < 0) or suitable (r0 > 0) among species with different persistence ability.

Absence from suitable areas could of course also indicate that occurrence is limited by environmental factors not considered in our analyses. Ideally, transplant experiments could be used to test model predictions of suitable areas outside the range (27). Such large-scale transplant experiments do, however, pose substantial ethical and logistic problems (28). Instead, we made use of the fact that transplantation by humans (notably flower producers) caused most of our study species to form naturalized populations in natural ecosystems outside their native geographic range (Fig. 1) (17). When evaluating our model extrapolations for these naturalized populations, we found that the predicted r0 was positive for an average of 80% populations per species (). This quasi-experimental evidence suggests that the demographic niche models capture key factors limiting the geographic distributions of our study species.

On small spatial scales, niche–distribution mismatches can arise from source–sink effects and from time-delayed extinction (4–6). To assess the match between demographic suitability and a species’ occurrence within its geographic range, we regressed spatial variation in occupancy (on a 1′ grid, c. 1.55 × 1.85 km2) against the locally predicted r0. Occupancy generally increased with predicted r0 (Fig. 3). However, the strength of this relationship varied strongly among species (R2N: 0.01 to 0.55, mean = 0.18; ) and significantly differed between species of different persistence ability, where variation in occupancy was better explained for nonsprouters than for resprouters (ANOVA, F1,24 = 4.74, P = 0.039). For resprouters, small-scale niche–distribution mismatches are greater mainly because populations more frequently occur in unsuitable sites (Fig. 3). This can be explained by populations of more persistent species being less vulnerable to adverse conditions in both temporally fluctuating (29) and directionally changing environments (30). Dispersal ability had no positive effect on the occupancy of unsuitable sites. Hence, we found no indication of source–sink effects at the spatial scale of our analysis.

Conclusions

In summary, our comparative analysis of demography-based niches indicates that key life-history traits shape the geographic distributions of plant species by affecting not only niche sizes but also niche–distribution mismatches. Specifically, mismatches between niches and geographic distributions arise because poorly dispersed species are absent from suitable sites beyond their range limits and because species with high persistence ability are present in sites that are unsuitable under current, average environmental conditions. Importantly, this identifies poorly dispersing and highly persistent species as cases where static, correlative species distribution models are more likely to fail. For such species, range forecasts require dynamic species distribution models that incorporate demographic niche estimates (7, 8, 11). From a theoretical perspective, the quantification of spatial variation in species’ intrinsic population growth rates is also the first step toward understanding effects of biotic interactions on range dynamics, large-scale species coexistence (31, 32), and niche evolution (4). Linking properties of demography-based niches to a wider spectrum of functional traits can furthermore elucidate patterns and dynamics of functional biodiversity (33, 34). A demographic quantification of ecological niches, particularly in well-studied model systems, thus holds great promise for better integrating ecological theory and empirical biogeography. This is urgently needed to advance our predictive understanding of biodiversity dynamics under environmental change (10, 12).

Materials and Methods

Study Species and Demographic Data.

We studied 26 species of the Proteaceae family, specifically of the genera Protea (16 species) and Leucadendron (10 species), that are endemic to the Cape Floristic Region (17). These species were chosen to represent variation in geographic distributions as well as variation in dispersal and resprouting ability. For each species we obtained data on between-population variation in key demographic rates across the entire life cycle, namely, the total fecundity of adult plants since the last fire (size of individual canopy seed banks), per capita postfire seedling recruitment (ratio between postfire recruits and prefire adults), and adult fire survival. The latter two rates were measured on recently burned sites (<3 y after fire), where burned prefire adults were still identifiable (35, 36). For each study species, demographic sampling sites were selected to cover both major environmental gradients across the species’ global geographic distribution and variation in population densities (). The final dataset comprised 3,617 population-level records from an average of 99 (median = 85) study sites per species (). For details on the demographic data collection, see ref. 19.

Study Region and Environmental Variables.

Our study area was defined on a regular grid with a spatial resolution of 1′ × 1′ (c. 1.55 km × 1.85 km) and included all grid cells of the Cape Floristic Region in which >5% of the area is covered by Fynbos vegetation (37). Climatic and edaphic variables that are expected to be main determinants of the performance and survival of serotinous Proteaceae were extracted from the South African Atlas of Climatology and Agrohydrology (38). We included January maximum daily temperature (Tmax), July minimum daily temperature (Tmin), and a January aridity index (AI) calculated as the ratio between the mean values of precipitation (P) and temperature (T): AI = P/(T + 10 °C) (39). Climatic variables are averages over the years 1950 to 2000. As an edaphic variable we used a soil fertility index that combines soil texture and base status and ranges from 0 to 10 (38). These climatic and edaphic variables showed no strong collinearity and were also not strongly correlated with population densities of the study species (). Information on the fire return interval was obtained from both observational records and model predictions. For the demographic sampling sites, information on the fire history (time since the last fire and length of the previous fire interval) was inferred from a combination of measured plant ages, historical records and Moderate Resolution Imaging Spectroradiometer (MODIS) satellite observations (19, 40–42). For predictions of population growth rates across the study region (see below), we used probability distributions of fire return intervals predicted from a climate-driven model of postfire ecosystem recovery (22).

Demographic Response Model.

We used a hierarchical Bayesian modeling approach for estimating the species-specific responses of key demographic rates (fecundity, per-seed establishment, and adult fire survival) to environmental covariates (). The model considers effects of climatic and edaphic conditions, variable fire return intervals, and intraspecific density dependence at both the adult and the seedling stage. Below, we describe the submodels for variation in each demographic rate.

Fecundity.

The recorded size of the canopy seed bank (Seed.count) of plant j in population i is described by an overdispersed Poisson distribution,

where the expected value of mean fecundity μ.fec is determined by limiting effects of postfire stand age (Age), environmental covariates (E), and population density (D):

Effects of stand age on fecundity arise from the time of maturation until the first flowering and cone production, increasing accumulation of standing cones on growing plants, cone loss, and possibly senescence of aged individuals:

where Mi is a binary random variable (0, 1) indicating maturity. The probability of population-level maturity is calculated from a Weibull distribution for the age (t.mat) of first cone production:

The species-specific time to reproductive maturity (t.mat) was constrained to be at least 3 y for nonsprouters (17). The effects of the environmental covariates k = 1…K are described by Gaussian demographic response functions:

where opt.fec denotes the optimal conditions and sig.fec measures the width of the response curve. Effects of population density D on fecundity are described as (43)

Establishment.

The number of observed recruits on a recently burned site depends on the fecundity of prefire adults, a density-dependent per-seed establishment rate, and also how long after fire the site was sampled. The establishment of new recruits from seeds is modeled as a binomial process where the number of recruits (#Recruits) in population i depends on the total number of available seeds (#Seeds) in the canopy seed bank at the time of the last fire and the per-seed establishment rate π.est:

Since #Seeds is unknown for recently burned sites where recruitment was recorded, it is modeled as a latent state variable,

where #Parents denotes the number of prefire seed sources (only females for dioecious Leucadendron species) and Fec depends on environmental covariates (E), the postfire stand age (Age), and the adult population density (D) at the time of the previous fire as described in the fecundity submodel. Establishment rate π.est is affected by environmental covariates (E) and by the densities of seeds SD = #Seeds/Area and fire-surviving adults AD = #Adults/Area:

As for fecundity, the effects of different environmental covariates k = 1…K are described by Gaussian demographic response functions:

Density effects on establishment result from the density of seeds (SD) as well as from the density of fire-surviving adults (AD), with different strengths (γ.est.SD resp. γ.est.AD) for each of these density effects:

Density-dependent mortality of recruits (self-thinning) is a continuous process, which for Fynbos Proteaceae generally occurs within the first 3 y after a fire (44). This is described by weighting the density effects with a factor ci that depends on the postfire stand age (pf.Age) at the time of sampling:

Thereby the model accounts for the fact that more seedlings can be observed if a site is surveyed just shortly after germination (19).

Survival.

Adult fire survival is modeled as a binomial process for the proportion of survivors among all prefire adults (#All.Adults):

Similar as for fecundity, the effects of different environmental covariates k = 1…K and of prefire stand age (Age) are described by Gaussian response functions and effects of population density D by a negative-exponential function:

Since adult fire survival rates of nonsprouters are generally low with little intraspecific variation (19), we modeled them as species-specific constants and considered effects of covariates only for the survival rates of resprouters.

Bayesian parameter estimation.

Parameters of the model were estimated independently for each study species in a Bayesian framework (see for details and model code). An overview of parameter prior distributions is given in . In particular, we chose a wide normal prior (mean = 0, variance = 104) for the location of the environmental optima (opt.feck, opt. estk, and opt.survk) and an inverse negative exponential prior (rate = 1) for the environmental response strength (sig.feck, sig.est2k, or sig.surv2k). These prior distributions allow the model to describe, over the range of occurring environmental conditions, unimodal relationships as well as monotonic positive, monotonic negative, or practically flat demographic response curves ().

Model evaluation.

For each species we assessed the model fit separately for each observed demographic variable (fecundity, recruit:parent ratio, and adult fire survival) by calculating Nagelkerke’s general RN2 (45) relative to null models in which demographic rates (π.est, π.surv, and μ.fec) are species-specific constants. The explained variance in each demographic variable for each study species is shown in .

Prediction of Niche Sizes and Geographic Variation in r0.

From the estimated demographic response curves we characterized a species’ niche by predicting the expected intrinsic growth rate of small populations (r0) across environmental gradients. The demographic response model predicts fire survival rate π.surv, fecundity μ.fec, and establishment rate π.est as functions of environmental covariates E, fire interval T, and population density at different stages. Based on these demographic rates the expected population size N after a fire interval of length T can be calculated as the sum of fire survivors and new recruits:

For the dioecious Leucadendron species the parameter p.fem specifies the proportion of female individuals in a population and accounts for the fact that fecundity rate μ.fec was defined per female. The niche of a species is defined as the set of environmental conditions for which the intrinsic growth rate of small populations (r0) is positive. To calculate r0 we set all density variables to zero and first calculated the rate of change in population size per fire interval,

Hence λ0 is the expected population size after the next fire (number of established offspring plus the possibly surviving adult) if a single adult occurs on a site without competition from conspecifics. The intrinsic growth rate r0 was then calculated on an annual basis as

To quantify species’ niches, r0 was predicted on a 5D grid spanned by niche axes according to the four climatic-edaphic covariates in E (Tmax, Tmin, AI, and soil fertility) and the fire return interval T (log-transformed). Since the demographic data analyzed in the demographic response models were measured in natural communities and thus incorporate effects of interspecific biotic interactions, the predicted r0 represents the postinteractive (or realized) niche (2). For commensurability of the different niche axes, these were confined and scaled to the respective range of values that occur throughout the Fynbos biome for each variable (so that each axis ranges from zero to one), and each niche axis was regularly sampled with a resolution of 0.01 (1010 grid points). The niche size was then quantified separately for fire return interval (disturbance niche size) and for the four climatic-edaphic variables (environmental niche size). The disturbance niche size was determined as the range of fire return intervals for which a positive r0 is predicted when the climatic-edaphic variables are set to their optimal values. Likewise, the environmental niche size was determined as the hypervolume in the 4D climatic-edaphic subspace for which the predicted r0 is positive when setting the fire return interval to its species-specific optimum.

To geographically project r0 across the study region, r0 was not predicted for a single fixed fire return interval but integrated as a weighted geometric mean over the site-specific probability distribution of fire return intervals (22). In all cases, values of r0 were predicted from each posterior sample of the demographic response model, yielding full posterior distributions for each predicted value.

Species Distribution Data, Geographic Ranges, and Range Filling.

The Protea Atlas Project was an extensive citizen science project that used a standardized protocol to collect complete species lists from Proteaceae communities across the Cape Floristic Region (17). For our study region, the Protea Atlas database contains 54,642 sampling locations with a total of 126,690 recorded presences of the 26 study species (). We aggregated these occurrence data to the proportion of presence records (occupancy) among the sampled communities in each 1′ × 1′ grid cell. As an overall assessment of how well species occurrence is predicted by demographic suitability, we calculated the AUC (46) for a binary classification of grid cell presence–absence by the predicted r0. We furthermore described the large-scale geographic range of each species by an alpha-convex hull over all presence records, using the alphahull package in R with parameter α = 0.5 (47, 48). The size of each species’ geographic range was calculated as the overlap of this alpha-convex hull with the Fynbos biome (study region). The degree of range filling was calculated as the proportion of potentially suitable area (cells with predicted r0 > 0) that lies within the geographic range. For comparison, we also applied two alternative approaches for quantifying species’ range filling, where we 1) tested whether potential species-specific biases in predicted r0 could affect the analysis of interspecific variation in range filling and 2) quantified range filling based on species distribution models that were directly fitted to the presence–absence data of each species (see for details and , for results).

Statistical Analyses of Relationships Between Life-History Traits, Demographic Rates, Niche Characteristics, and Geographic Ranges.

We performed separate regression analyses to test for effects of key life-history traits on species’ demographic rates (maximum fecundity max.fec, maximum establishment rate max.est, and maximum adult fire survival rate max.surv), the two measures of niche size (disturbance niche size and environmental niche size) as well as on species’ geographic range sizes and range filling. The traits considered as explanatory variables were persistence ability (0 = nonsprouter, 1 = resprouter) and long-distance dispersal ability. The species-specific relative long-distance dispersal ability was derived from a trait-based mechanistic model of primary and secondary wind dispersal (23) and measured as the number of neighboring cells on a 1′ × 1′ rectangular grid that can be reached by dispersal from a source cell with a probability of at least 10−4. In each regression model we also accounted for phylogenetic dependence. A molecular phylogeny of our study species was obtained by pruning a phylogenetic tree of 291 Proteaceae species (; see ref. 34 for numerical data). This phylogenetic tree was constructed from a supermatrix combining molecular markers for Leucadendron, Protea, and related species of the Proteaceae family (49–51). As quantitative measure of the degree of phylogenetic dependence the model estimates Pagel’s λ (52). For each response variable, we estimated a full model that included effects of persistence ability, dispersal ability, and their interaction (see for details on Bayesian parameter estimation). All simplified models nested in the full model were then compared by the deviance information criteria (DIC) (53), and we report parameter estimates for the DIC-minimal models ().

Statistical Analysis of the Relationship Between Occupancy and Demographic Suitability.

We used a binomial nonlinear regression model to analyze the relationship between each species’ occupancy in the 1′ × 1′ grid cells within its geographic range and habitat suitability. The model predicts the occurrence probability ψs,i of species s in grid cell i in dependence of the respective predicted intrinsic population growth rate r0s,i as ψs,i = as/(1 + exp(–bs(r0s,i – cs))), where as, bs, and cs are species-specific regression parameters (see for details on Bayesian parameter estimation).

Data and Code Availability.

Demographic data were partly used under license from CapeNature for the current study and are available from the authors upon reasonable request and with permission of CapeNature. Documented JAGS code for the statistical analyses is provided in . Additionally used R code and data generated during the analyses are available from the corresponding author upon request.