Using ISI Web of Knowledge we conducted a systematic literature search of all relevant articles—published between 1990 and 2015—that reported an observed change in climate in the study area; indicated that birds and/or mammals have undergone a change (for example, in distribution, population size, phenology, behaviour, genotype, phenotype) attributable to climate in the past 100 years; and/or suggested that populations of a species were not affected by recent climate change. For each study and each species considered (70 studies and 120 species for mammals, 66 studies and 569 species for birds), we identified the type of impact experienced.
A negative response was assigned to a species if all (at least one) or >50% of its populations (if the species had both negative and no responses in different portions of its range) were reported to have undergone declines in population size, geographic range size, survival or reproductive rate, and body mass, thus reducing the risk of false attributions. These responses were confidently attributable to recent climate change by the authors of the studies, for instance due to the fact that the most significant change in environmental and biotic conditions reported in the area in which the population of the species was impacted was related to climatic variables. Although we acknowledge that some of the studies may have been more rigorous than others, with such variation in the methods used and the effect size themselves it would have been difficult to adjudicate the level of confidence around the claimed relationship, although we believe that evaluating the strength of attribution is a priority for future work.
A positive response was assigned if the majority of the populations of a species experienced geographic range expansions, increase in population size, survival rate and/or reproductive rate, body mass, and/or changes in phenology. An unchanged response was attributed if no response was observed despite the recorded change in climate. Finally, species that exhibited a combination of the negative and positive (not necessarily in the same proportion) responses in different parts of their range were classified as mixed.
To identify the relationships between the observed response of mammals and birds to climate change and a set of intrinsic and spatial variables (see Supplementary Methods for description of these predictors and a priori hypotheses), we performed a multinomial logistic regression using the ‘nnet’ package in R. This model uses maximum likelihood estimation to evaluate the probability of the different possible outcomes of a categorical dependent variable with more than two classes. To reduce the overdispersion in models and avoid collinearity, we performed Spearman’s correlation tests between the predictors and removed those that were highly correlated (R2 > 0.75) and led to the minimum loss in model performance.
We included taxonomic order as fixed variable of our models, for a total of 11 orders of mammals and 22 of birds. By including taxonomy as a fixed effect, we aimed to control for the non-independence of observed responses across species, and for the latent variables that may affect the responses to climate change that are phylogenetically conserved. We did not include taxonomic family or genus because it resulted in strong underdispersion, as observed data on the response to climate change (which we used as a base for our predictions on threatened species) were often available only for the populations of one species per family/genus. Since we are not aware of frequentist methods to implement phylogenetically corrected models with a multinomial distribution, and concerned that phylogenetic non-independence in the species in our data set could nevertheless be important, we tested for the existence of phylogenetic signal in the residuals of our models. We used phylogenetic trees for mammals and birds34, 35 to estimate Pagel’s lambda, assuming a star-shaped phylogeny and the actual phylogeny (Brownian motion models). We tested whether the value of lambda differed significantly from 0 (no phylogenetic signal) and 1 (trait distribution matches a Brownian model of evolution), by computing the likelihood ratio, and then comparing it to a Chi-squared distribution with one degree of freedom. If the test is significant there is phylogenetic signal in the residuals. However, we found lambda values of 6.73e-05, 5.56e-04 and 2.68e-04, and p-values of 0.51, 0.47 and 0.62 for mammals, birds in breeding ranges and birds in non-breeding ranges, respectively. Therefore we conclude that there is no phylogenetic signal in the residuals of the models and a phylogenetically informed model is not justified.
We performed a model selection using the Akaike Information Criterion (AIC) to find the set of predictors to include in the final model that minimize the Kullback–Leibler distance between the model and the observed values. We applied logarithmic and quadratic transformations to the predictors and included variable interactions in the models, but most of them did not lead to a decrease in AIC or increase in model performance calculated by using the Area Under the Curve (AUC). Finally, to test our models for overdispersion, we calculated the sum of squared Pearson residuals and compared it to the residual degrees of freedom by using a Chi-squared test. P-values close to 1 indicate that the probability of the model being overdispersed approaches 0 (Supplementary Table 13).
On the basis of the relationship between the observed response of species and our independent variables found with the best multinomial models, we predicted the probabilities of the four classes of response to climate change by using the function predict in R. For predictions we considered all threatened birds (1,272 species, as listed on the 2014 IUCN Red List) and terrestrial non-volant mammals (873 species) with available data. We excluded sea mammals from our analysis because the environmental variables that influence the persistence of marine and terrestrial species are different, and most of the variables important for marine species (for example, sea temperature, salinity) were not available for the study period. Chiroptera could not be considered in this study because of the paucity of data available on their life history.
Our model is at the species level, but our data (observed responses to climate change) is at the population level. Because the spatial extent of the study area was not available for the vast majority of studies, we were forced to average the annual temperature change experienced by the species across all of its range. However, the average climatic change might not be representative of the change experienced by the populations we used to train the model, especially with species with large range size. By resampling the response category assigned to each species from the multinomial distribution 100 times and deriving coefficient intervals and mean values of the richness of species with negative responses, we tried to reduce the uncertainty around our predictions. In addition, to identify the taxonomic orders for which our predictions were most reliable, we used a Kolmogorov–Smirnov nonparametric test which quantifies the distance between the empirical continuous distribution functions of two samples, and the null hypothesis is that the samples are drawn from the same distribution. By comparing the distribution of the same numeric trait in both the observed and the predicted sample, if the p-value of the test is above the α threshold, that is, 0.05, we can assume that threatened species in the considered taxonomic order are well represented in the sample of observed data. This means that, for this order, our predictions are more robust.
The authors declare that (the/all other) data supporting the findings of this study are available within the article and its Supplementary Information files.