Herbaceous perennial plants with short generation time have stronger responses to climate anomalies than those with longer generation time
Demographic dataTo address our hypotheses, we used matrix population models (MPMs) or integral projection models (IPMs) from the COMPADRE Plant Matrix Database (v. 5.0.156) and the PADRINO IPM Database57, which we amended with a systematic literature search. First, we selected density-independent models from COMPADRE and PADRINO which described the transition of a population from 1 year to the next. Among these, we selected studies with at least six annual transition matrices, to balance the needs of adequate yearly temporal replicates and sufficient sample size for a quantitative synthesis. This yielded data from 48 species and 144 populations.We then performed a systematic literature search for studies linking climate drivers to structured population models in the form of either MPMs or IPMs. We performed this search on ISI Web of Science for studies published between 1997 and 2017. We used a Boolean expression containing keywords related to plant form, structured demographic models, and environmental drivers (Supplementary Methods). We only considered studies linking macro-climatic drivers to natural populations (e.g., transplant experiments and studies focused on local climatic factors such as soil moisture, light due to treefall gaps, etc. were excluded). Finally, we used the same criteria used to filter studies in COMPARE and PARDINO, by selecting studies with at least six, density-independent, annual projection models. This search brought two additional species, belonging to three additional populations, which we entered in the COMPADRE database.One of the studies we excluded from the literature search because it contained density-dependent IPMs, also provided raw data with high temporal replication (14–32 years of sampling) for 12 species from 15 populations58. Therefore, we re-analyzed these freely available data to produce density-independent MPMs that were directly comparable to the other studies in our dataset (Supplementary Methods).The resulting dataset consisted of 46 studies, 62 species, 162 populations, and a total of 3761 MPMs and 52 IPMs (Supplementary Data 1). The analyzed plant populations were tracked for a mean of 16 (median of 12) annual transitions. To our knowledge, this is the largest open-access dataset of long-term structured population projection models. However, this dataset is taxonomically and geographically biased. Specifically, among our 62 species, this dataset contains 54 herbaceous perennials (11 of which graminoids), and eight woody species: five shrubs, two trees, and one woody succulent (Opuntia imbricata). Moreover, almost all of these studies were conducted in North America and Europe (Supplementary Fig. 1), in temperate biomes that are cold, dry, or both cold and dry (Supplementary Fig. 1, inset). Our geographic and taxonomic bias reflects the rarity of long-term plant demographic data in general. This dearth of long-term demographic data is particularly evident in the tropics. The ForestGEO network59 is an exception to this rule, but to date, no matrix population models or integral projection models using these data have been published.We used the MPMs and IPMs in this dataset to calculate the response variable of our analyses: the yearly asymptotic population growth rate (λ). This measure is one of the most widely used summary statistics in population ecology60, as it integrates the response of multiple interacting vital rates. Specifically, λ reflects the population growth rate that a population would attain if its vital rates remained constant through time61. This metric therefore distills the effect of underlying vital rates on population dynamics, free of other confounding factors (e.g., transient dynamics arising from population structure62). We calculated λ of each MPM or IPM with standard methods61,63. Because our MPMs and IPMs described the demography of a population transitioning from one year to the next, our λ values were comparable in time units. Finally, we identified and categorized any non-climatic driver associated with these MPMs and IPMs. Data associated with 21 of our 62 species explicitly quantified a non-climatic driver (e.g., grazing, neighbor competition), for a total of 60 of our 162 populations. Of the datasets associated with these species, 19 included discrete drivers, and only three included a continuous driver.Climatic dataTo test the effect of temporal climatic variation on demography, we gathered global climatic data. We downloaded 1 km2 gridded monthly values for maximum temperature, minimum temperature, and total precipitation between 1901 and 2016 from CHELSAcruts64, which combines the CRU TS 4.0165, and CHELSA66 datasets. Gridded climatic data are especially suited to estimate annual climatic means45. These datasets include values from 1901 to 2016, which are necessary to cover the temporal extent of all 162 plant populations considered in our analysis. For our temperature analyses, we calculated the mean monthly temperature as the mean of the minimum and maximum monthly temperatures. We used monthly values to calculate the time series of mean annual temperature and total annual precipitation at each site. We then used this dataset to calculate our annual anomalies for each census year, defined as the 12 months preceding a population census. Our annual anomalies are standardized z-scores. For example, if X is a vector of 40 yearly precipitation or temperature values, E() calculates the mean, and σ() calculates the standard deviation, we compute annual anomalies as A = [X − E(X)]/σ(X). Therefore, an anomaly of one refers to a year where precipitation or temperature was one standard deviation above the 40-year mean. In other words, anomalies represent how infrequent annual climatic conditions are at a site. Specifically, if we assume that A values are normally distributed, values exceeding one and two should occur every 6 and 44 years, respectively. We used 40-year means because the minimum number of years suggested to calculate climate averages is 3067.Z-scores are commonly used in global studies on vegetation responses to climate8,68, and they reflect the null hypothesis that species are adapted to the climatic variation at their respective sites. Across our populations, the standard deviations of annual precipitation and temperature anomalies change by 300% and 60%, respectively (Supplementary Fig. 2). Thus, a z-score of one refers to a precipitation anomaly of 50 or 160 mm and to a temperature anomaly of 0.5 or 0.8 °C. Our null hypothesis posits that species are adapted to these conditions, regardless of the absolute magnitude of the standard deviation in annual climatic anomalies. If this null hypothesis were true, each species would respond similarly to z-scores. Z-scores are more easily interpreted when calculated on normally distributed variables. We found our temperature and precipitation z-scores were highly skewed (skewness above 1) only in, respectively, 2 (for temperature) and three (for precipitation) of our 162 populations. We concluded that this degree of skewness should not bias our z-scores substantially.To test how the response of plant populations to climate changes based on biome we used two proxies of water and temperature limitation. For each study population, we computed a proxy for water limitation, water availability index (WAI), and temperature limitation using mean annual temperature. To compute these metrics, we downloaded data at 1 km2 resolution for mean annual potential evapotranspiration, mean annual precipitation, and mean annual temperature referred to the 1970–2000 period. We obtained potential evapotranspiration data from the CGIAR-CSI consortium (http://www.cgiar-csi.org/). This dataset calculates potential evapotranspiration using the Hargreaves method69. We obtained mean annual precipitation and mean annual temperature from Worldclim70. Here, we used WorldClim rather than CHELSA climatic data because the CGIAR-CSI potential evapotranspiration data were computed from the former. We calculated the WAI values at each of our sites by subtracting mean annual potential evapotranspiration from the mean annual precipitation. Such proxy is a coarse measure of plant water availability that ignores information such as soil characteristics and plant rooting depth. However, WAI is useful to compare water availability among disparate environments, so that it is often employed in global analyses68,71. As our proxy of temperature limitation, we use mean annual temperature. While growing degree days would be a more mechanistic measure of temperature limitation48, this requires daily weather data. However, we could not find a global, downscaled, daily gridded weather dataset to calculate this metric.The overall effect of climate on plant population growth rateTo test H1, we estimated the overall effect sizes of responses to anomalies in temperature, precipitation, and their interaction with a linear mixed-effect model.$${mathrm{log}}left( lambda right) = alpha + beta P + eta T + theta P{mathrm{x}}T + varepsilon$$
(1)
where log(λ) is the log of the asymptotic population growth rate of plant population P is precipitation, T is temperature. We included random population effects on the intercept and the slopes to account for the nonindependence of measurements within populations. We then compared the mean absolute effect size of precipitation, temperature, and their interaction. This final model did not include a quadratic term of temperature and precipitation because these additional terms led to convergence issues. This likely occurred because single data sets did not include enough years of data.Population-level effect of climate on plant population growth ratesTo test our remaining three hypotheses, we carried out meta-regressions where the response variable was the slope (henceforth “effect size”) of climatic anomalies on the population growth rate for each of our populations. Before carrying out our meta-regression, we first estimated the effect size of our two climatic anomalies on the population growth rate of each population separately. We initially fit population-level and meta-regression simultaneously, in a hierarchical Bayesian framework. However, these Bayesian models shrunk the uncertainty of the noisiest population–level relationships, resulting in unrealistically strong meta-regressions. We, therefore, chose to fit population models separately, resulting in more conservative results.For each population, we fit multiple regressions with an autoregressive error term, and we evaluated the potential for nonlinear effects in the datasets longer than 14 years. We fit multiple regressions because temperature and precipitation anomalies were negatively correlated, so that fitting separate models for temperature and precipitation would yield biased results72. We fit an autoregressive error term because density dependence and autocorrelated climate anomalies can produce autocorrelated plant population growth rates. The form of our baseline model was$${rm{log}}(lambda )_y = alpha + beta _pP_y + beta _tT_y + varepsilon _y$$
(2)
$$varepsilon _y = rho varepsilon _{y – 1} + eta _y$$
(3)
The model in Eq. 2 is a linear regression relating each log(λ) data point observed in year y, to the corresponding precipitation (P) and temperature (T) anomalies observed in year y, via the intercept α, the effect sizes, β, and an error term, εy, which depends on white noise, ηy, and on the correlation with the error term of the previous year, ρ. When multiple spatial replicates per each population were available each year, we estimated the ρ autocorrelation value separately for each replicate. This happened in the few cases when a study contained contiguous populations, with no ecologically meaningful (e.g., habitat) differences.We compared the baseline model in Eqs. 2 and 3 to models including a quadratic climatic effect and non-climatic covariates. We estimated quadratic climatic effects only for time series longer than 14 years. We choose this threshold because when using a model selection approach to select a quadratic or linear regression model, the recommended minimum sample size is between 8 and 25 data points73. We fit models including a quadratic effect of temperature, precipitation, or both (Supplementary Table 1).Finally, we also tested whether non-climatic covariates could bias the effects of climate on log(λ) estimated in our analysis. Such bias, either upwards or downwards, could result in the case non-climatic co-variates interacted with climate. For example, harvest can have multiplicative, rather than additive effects on the climate responses of forest understory herbs74. We tested for an interaction between a covariate and climate anomaly in 17 of the 21 studies that included a non-climatic covariate. In the remaining three studies, discrete covariates corresponded with the single populations. Because Eqs. 2 and 3 is fit on separate populations, it implicitly accounted for these covariates. For the 17 studies above, we fit a linear effect of the non-climatic covariate and its interaction with one of the two linear climatic anomalies. Thus, including the linear model in Eqs. 2 and 3, the nonlinear models, and the covariate interaction models, we tested up to six alternative models for each one of our populations (Supplementary Table 1). We selected the best model according to the Akaike Information Criterion corrected for small sample sizes (AICc75). We carried out these and subsequent analyses in R version 3.6.176.In the populations for which AICc selected one of the model alternatives to the baseline in Eqs. 2 and 3, we calculated the effect size of climate by adding the effect of the new terms to the linear climatic terms. For example, when a quadratic precipitation model was selected, we calculated the effect size of precipitation as β = βp + βp2. For models including an interaction between temperature and a non-climatic covariate, we evaluated the effect of the interaction at the mean value of the covariate. Therefore, we calculated the effect size as β = βt + βxE(Ci) for continuous covariates. For categorical variables, we calculated the effect size as βp + βx0.5: that is, we calculated the mean effect size between the two categories. We quantified the standard error of the resulting effect sizes by adding the standard errors of the two terms.The effect of biome on the response of plants to climateWe used a simulation procedure to run two meta-regressions to test for the correlation between the effect size of climate drivers on λ, and our measures of water or temperature limitation. These meta-regressions accounted for the uncertainty, measured as the standard error, in the effect sizes of climate drivers. We represented the effect of biome using a proxy of water (WAI) and temperature (mean annual temperature) limitation. For each of our 162 populations, the response data of this analysis were the effect sizes (βp or βt values) estimated by Eqs. 2 and 3 or their modifications in case a quadratic or non-climatic covariate model were selected. In these meta-regressions, the weight of each effect size was inversely proportional to its standard error. To test H2 and H3 on how water and temperature limitation should affect the response of populations to climate, we used linear meta-regressions. These two hypotheses tested both the sign and magnitude of the effect of climate. Therefore, we used the effect sizes as a response variable which could take negative or positive values. As predictors, we used population-specific WAI (H2, only for effect sizes quantifying the effect of precipitation), and mean annual temperature (H3, only for effect sizes quantifying the effect of temperature). The null hypothesis of these meta-regressions is that plant species are adapted to the climatic variation at their respective sites. Such an adaptation implies that a precipitation z-score of one should produce effects on log(λ) of similar magnitude and sign across different climates. This should happen across average climatic values that are connected to substantially different absolute climatic anomalies (Supplementary Fig. 2). On the other hand, our hypotheses posit that at low WAI and MAT values, species are more responsive to z-scores than expected under the null hypothesis.We performed these two meta-regressions by exploiting the standard error of each effect size. We simulated 1000 separate datasets where each effect size was independently drawn from a normal distribution whose mean was the estimated β value, and the standard deviation was the standard error of this β. These simulated datasets accounted for the uncertainty in the β values. We fit 1000 linear models, extracting for each its slope, βmeta. Each one of these slopes had in turn an uncertainty, quantified by its standard error, σmeta. For each βmeta, we then drew 1000 values from a normal distribution with mean βmeta and standard deviation σmeta. We used the resulting 1 × 106 values to estimate the confidence intervals of βmeta. This procedure assumes that the distribution of βmeta values is normally distributed. We performed one-tailed hypothesis tests, considering meta-regression slopes significant when over 95% of simulated values were below zero.The effect of generation time on the response of plants to climateTo test H4 on how the generation time of a species should mediate its responses to climate, we used a gamma meta-regression. We fitted gamma meta-regressions because our response variables were the absolute effect sizes of precipitation and temperature anomalies, |β|, which are bounded between 0 and infinity. To test H4, we therefore fit gamma meta-regressions with a log link, using |β| values as response variable and generation time (T) as predictor. We calculated T directly from the MPMs and IPMs (Supplementary Methods). We log-transformed T to improve model fit. We carried out these meta-regressions using the same simulation procedure described for testing H2 and H3. We also carried out one-tailed hypothesis tests, by verifying whether 95% of βmeta values were below zero.The effect of plant types on estimates of climate effectsWe verified whether certain plant types could bias our results by subdividing our species as graminoids, herbaceous perennials, ferns, woody species (shrubs and trees), and succulents. We ran ANOVA tests to verify whether the effect sizes of precipitation and temperature anomalies differed between plant types. We then tested for significant differences in pairwise contrasts between plants types by running Tukey’s honestly significant difference tests. We carried out these tests on the average effects of climate, without accounting for differences in parameter uncertainty. If Tukey’s test identified significant differences among plant types, we ran additional tests of H2–H4 excluding the plant type, or plant types, whose response to climate differed.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
