Less favourable climates constrain demographic strategies in plants

Abstract Correlative species distribution models are based on the observed relationship between species’ occurrence and macroclimate or other environmental variables. In climates predicted less favourable populations are expected to decline, and in favourable climates they are expected to persist. However, little comparative empirical support exists for a relationship between predicted climate suitability and population performance. We found that the performance of 93 populations of 34 plant species worldwide – as measured by in situ population growth rate, its temporal variation and extinction risk – was not correlated with climate suitability. However, correlations of demographic processes underpinning population performance with climate suitability indicated both resistance and vulnerability pathways of population responses to climate: in less suitable climates, plants experienced greater retrogression (resistance pathway) and greater variability in some demographic rates (vulnerability pathway). While a range of demographic strategies occur within species’ climatic niches, demographic strategies are more constrained in climates predicted to be less suitable.

these, 34 species across 93 populations had presence data adequate for fitting quality Species Distribution Models, which resulted in our final dataset (Figure S1.1 and Table   S1.1). www.compadre-db.org (Salguero-Gómez et al. 2015). We imposed the following species and study selection criteria: (1) field observations were carried out for at least three years (two matrix transitions), 2) matrices were under control conditions (experimental manipulations were excluded) (3) individual matrices were reported for each year and population separately so that we were able to compute stochastic demographic metrics and measures of temporal variability, (4) the study was conducted within the species' native geographic range (5) the geographic location of populations was reported with a coordinate accuracy of at least 10 km (i.e., at least one digit) to ensure matching with grid size of environmental variables, (6) matrices were complete, with all vital rates, including fecundity, recorded, (7) matrices were quality checked i.e., survival was <1.09 in "survival issue" column in COMPADRE (stage-specific survival probabilities cannot be much greater than 1, see Salguero-Gómez et al. 2015) and matrices were not suspected of seed stage error i.e., fecundity was reported as the fate of organisms over a time interval and not simple seed count, (8) matrices were irreducible and ergodic (see Stott et al. 2010 for details) (9) species belonged to either the "tree" or "herbaceous perennial" growth types (other growth forms were excluded due to low sample size after all other selection criteria). Precision of population coordinates was crosschecked and verified on Google Earth for all populations where it was possible based on the original paper. Subsequently another seven species (eight populations) with study length longer than 25 years and matrix dimension higher than matrix dimension of other populations.
Species occurrence data selection. Overall, the incomplete and sometimes poor quality data of global biodiversity repositories limited the number of species selected from COMPADRE for which we could fit quality distribution models. To overcome the reliability issue of single-source occurrence data for producing realistic distribution maps (Duputié et al. 2014), we compiled the dataset from two major global and continental data repositories: GBIF (Global Biodiversity Information Facility, accessed on 28 November 2014) and BIEN (Botanical Information and Ecology Network, www.bien3.org 1 , accessed on 16 November 2014), from local and regional herbaria and digitized species distribution maps from atlases (Meusel and Jäger 1992) (Table S1.1 below). There are a number of limitations to using occurrence records from big data repositories in species distribution modeling (Newbold 2010), therefore we carefully cleaned the GBIF, BIEN and regional herbaria data and selected locality records as reliably as possible, including only occurrences reported to at least one digit (~10 km) accuracy, with a coordinate precision of 10 km or no reported precision, and deleting cultivated specimens, non-native occurrences and duplicated records. For species outside Eurasia, we approximated based on independent descriptions whether the occurrence dataset covered roughly the whole species area of distribution, and if gaps spanning > ca. 75% of the range were identified, we assumed that the probability of the existing dataset representing the species' ecological requirements was low and the 1 http://bien.nceas.ucsb.edu/bien/; http://vegpath.org/BIEN3/wiki/conditions_of_use; species was omitted. In Eurasian taxa, where global datasets were incomplete we used digitized species distribution maps from atlases to ensure we covered species' full geographic range (Meusel and Jäger 1992). We extracted all presence points reported as such and another set of presence points from a large number of random points within the polygon of the distributional range (see maps below for a visual comparison of occurrence data completeness of the pooled GBIF and BIEN datasets and the Meusel-Jäger dataset for species with available maps in our dataset).
To reduce spatial autocorrelation due to sampling bias, occurrences closer to each other than 10 km were removed using a disaggregation index in the ecospat R package (Broennimann et al. 2015). Various taxa were checked for taxonomic consistency across the different data sources and infrataxa with poorly known geographic distribution were not modeled (e.g., Anthyllis vulneraria ssp. alpicola).
Taxa for which the final number of occurrences was lower than 30 and/or were suspected of not covering the range of environmental conditions across the area of distribution were not used because of increased risk of biased predictions (Wisz et al. 2008). The minimum number of occurrence points available per species was 56 (Table   S1.1), ensuring sufficient modeling power. We did not use the COMPADRE population locations to build SDMs and we deleted any occurrence points that overlapped with COMPADRE populations. Appendix S1.2. Calculation of population growth rates, time to quasi-extinction,

transient population dynamics and demographic processes (matrix transitions)
Mean and variation of population growth rate (λ). Matrix projection models (MPMs) were used to calculate five different measures of population growth rate. Two mean deterministic measures were calculated for each population: the arithmetic mean of individual population growth rate estimates over time (λ arith ) and the geometric mean of individual population growth rate estimates over time (λ geom ). In contrast to deterministic population performance metrics stochastic metrics include the effect of observed temporal variability in estimates. We calculated stochastic population growth rate (λ iid ) for each population by simulation. We started from the stable stage distribution i.e., a theoretical state in which the proportion of individuals in each developmental stage remains constant over time and projected population growth over 50,000 time intervals by randomly choosing one of the individual matrix models recorded for a population. λ iid was the arithmetic mean of log [N (t+1) /N (t) ] of each successive interval over all time steps (Morris and Doak 2002). To account for differences in the effect of the temporal variation in individual matrix entries on population growth, we calculated Tuljapurkar's analytic approximation for population growth rate, λ Tulja , which estimates the stochastic population growth rate from matrices and covariances among matrix entries (Tuljapurkar 1982). However because this method assumes that the variation among individual matrices is low, results need cautious interpretation (Morris and Doak 2002).
Both the simulation and Tuljapurkar's analytic approximation assume lack of temporal autocorrelation between matrix elements and population growth rates (i.e., identical and independently distributed environments). We estimated a third, sequential growth rate (λ seq ) in temporally autocorrelated environments. For this metric, each species' population, the population matrix models were used in the same sequence that they were used from data from the field, and we report results only for the latter. To obtain robust estimates, extinction probability curves were averaged over 10 separate simulations, with each simulation based on 5,000 iterations.
Transient population dynamics. Disturbances like fire, extreme weather, disease epidemics, pest outbreaks or human activity can change a population's age-, size-or stage structure if some stages are affected more than others. When a population is not at stable state, it can exhibit rates of growth that are faster (population 'amplification'), or slower (population 'attenuation') than predicted under stable conditions. These phenomena are due to a relative over-or under-representation of reproductive individuals compared to the stable structure. Transient bounds (Stott et al. 2011) describe the limits of these dynamics: the largest or smallest population size achievable following a disturbance, relative to size of the population in the absence of disturbance (i.e. growing at stable asymptotic rate). We calculated four bounds: the upper and lower bounds on reactivity, which respectively describe the largest and smallest possible abundance relative to stable state in the first time step of projection immediately following a disturbance and the upper and lower bounds on inertia, which respectively describe the largest and smallest ratios of population abundance relative to stable state that a population may settle to ad infinitum as a result of disturbance. We determined the range of transient dynamics which can occur in the first time step and ad infinitum respectively using log(upper bound) -log(lower bound) for each of reactivity and inertia.
The indices were calculated from the element-by-element arithmetic mean MPM for each population.
Demographic processes (matrix transitions). Matrix population models compiled in COMPADRE are constructed from field measurements and represent the proportion of individuals that transition across different size-, stage or age categories in a standard format across species (Salguero-Gómez et al. 2015). From each projection matrix we extracted the following basic demographic processes: fecundity, progression, retrogression and stasis (Silvertown et al. 1993). Fecundity represents reproductive effort and/or recruitment given that plants have survived, progression is when plants survive and grow, retrogression is when plants survive but lose some modules or shrink in size, stasis is when plants survive but neither grow or shrink. Clonality was not reported separately for either species. We calculated a weighted average of the corresponding matrix entries in the transition matrix. Weights corresponded to relative abundances of each stage at stationary equilibrium, retrieved from the right eigenvector of the projection matrix (the stable stage distribution). For example, to calculate mean progression in Figure S1.2, we calculated the sum of the entries in the corresponding columns in green, which were multiplied by the corresponding element of the stable stage distribution.

Study extent
In comparative distribution modeling studies, it is important to choose a large enough area, covering the niche of the modeled species that has been accessible to the species over relevant time periods (Barve et al. 2011 Table S1.1. for details) and a weight was applied in the calibration of the models so that presences and pseudo-absences have equal (0.5) prevalence. Models were calibrated with 70% of the data and evaluated with the remaining 30% of the data. The procedure was replicated 25 times per species with random training and evaluation datasets, adding up to 100 models per species (4 techniques x 25 runs each). We compared four different evaluation metrics: Roc AUC (Fielding and Bell 1997), TSS (Allouche et al. 2006), KAPPA (Allouche et al. 2006) and the "presence-only" evaluator Boyce index (Boyce et al. 2002, Hirzel et al. 2006     Aster amellus GBIF_BIEN