Meta‐analyzing the likely cross‐species responses to climate change

Abstract Ecological Niche Models (ENMs) have different performances in predicting potential geographic distributions. Here we meta‐analyzed the likely effects of climate change on the potential geographic distribution of 1,205 bird species from the Neotropical region, modeled using eight ENMs and three Atmosphere‐Ocean General Circulation Models (AOGCM). We considered the variability in ENMs performance to estimate a weighted mean difference between potential geographic distributions for baseline and future climates. On average, potential future ranges were projected to be from 25.7% to 44.5% smaller than current potential ranges across species. However, we found that 0.2% to 18.3% of the total variance in range shifts occurred “within species” (i.e., owing to the use of different modeling techniques and climate models) and 81.7% to 99.8% remained between species (i.e., it could be explained by ecological correlates). Using meta‐analytical techniques akin to regression, we also showed that potential range shifts are barely predicted by bird biological traits. We demonstrated that one can combine and reduce species‐specific effects with high uncertainty in ENMs and also explore potential causes of climate change effect on species using meta‐analytical tools. We also highlight that the search for powerful correlates of climate change‐induced range shifts can be a promising line of investigation.


| 11137
ORTEGA ET Al. controlled experiments (Walker et al., 2006). Empirical evidence of the effects of global warming on biodiversity is so overwhelming that the number of synthesis of these impacts is increasing steadily (e.g., Ainsworth, Rosenberg, & Wang, 2007 and references therein; Halupka & Halupka, 2017).
These studies estimated potential range shifts by using one or several methods for niche modeling associated with different Atmosphere-Ocean General Circulation Models (AOGCM) and Greenhouse Gas Emission Scenarios (GES). If one is interested in summarizing these results, in principle, the projected range shift of each species could be used as an unit of information in a meta-analysis.
Different issues arise when applying formal meta-analytic techniques to summarize niche modeling results, but critical is the uncertainty associated with range shift estimates. For example, Diniz-Filho et al. (2009) showed that the estimated shifts in species' geographic ranges varied broadly, mainly due to the use of different Ecological Niche Models (ENMs hereafter). Other sources of uncertainty include the choices of AOGCM, GES, parameterization methods and rules to transform continuous outputs of models into presence and absence estimates (Nenzén & Araújo, 2011). In other words, different model projections are obtained when different combinations of ENMs, GES, and AOGCMs are used. Furthermore, the characteristics of the modeled organisms (e.g., life history and dispersal ability) can also influence model outcomes (Buisson, Thuiller, Casajus, Lek, & Grenouillet, 2010;Lawler et al., 2010;Thuiller, Guéguen, Renaud, Karger, & Zimmermann, 2019).
Here, we propose that the uncertainty of species' range shifts can be used as weights in a meta-analysis. We exemplify our new approach using the estimated biogeographical range shifts of Neotropical birds to future climate change. Further, in a cross-species analysis, we asked whether the weighted mean projected shift has a negative (loss) or a positive (gain of geographic range) signal.
Finally, we tested the predictive ability of species-specific biological traits as correlates of range shifts.

| Data
We exemplify our analytical framework with the dataset originally used by Diniz-Filho et al. (2009). This dataset includes the extents of occurrence (range filling) for 1,205 Neotropical bird species. The data were downloaded from the BirdLife (former Nature Serve; http://www.dataz one.birdl ife.org) and resampled to a grid of 1° latitude × 1° longitude.
We used four bioclimatic variables in our ENMs (mean annual rainfall and variability, average temperature of the warmest and coldest months) because they represent the major drivers of species diversity in large spatial scales (see Hawkings, Porter, & Diniz-Filho, 2003). We obtained data on bioclimatic variables from the World Climate Research Program's (WCRP) Coupled Model Intercomparison Project phase 3 (CMIP3) multimodel dataset (Meehl et al., 2007;https ://www.esgf-node.llnl.gov).
We used a subset of the species (n = 1,205) modeled by Diniz-Filho et al. (2009) because we kept only those for which the following variables were available (see Appendix S1): altitude midpoint, body size, IUCN categories of extinction risk, clutch size, and migratory behavior (migrants or non-migrants). A previous study indicated that these traits are important in predicting bird extinction risk (Machado & Loyola, 2013). Altitude midpoint was defined as the mean between maximum and minimum altitude within ranges. Data on body size, clutch size, and migratory behavior came from different sources (see Machado & Loyola, 2013). IUCN categories of extinction risk were transformed to a discrete scale attributing zero, 1, 2, 3, and 4 to least concern (LC), near-threatened species (NT), vulnerable species (VU), endangered species (EN), and critically endangered species (CR), respectively (Cardillo et al., 2004;Purvis, Gittleman, Cowlishaw, & Mace, 2000).
We built a majority rule consensus phylogeny (Bryant, 2003) among the species with 10,000 random phylogenetic trees with "Hackett constraint" for the backbone topology from Jetz, Thomas, Joy, Hartmann, and Mooers (2012). For matching, we used the function match.phylo.comm of "picante" (Kembel et al., 2010) and the function "consensus.edges" of "phytools" R package (Revell, 2012) to build the consensus phylogeny.

| Statistical analyses
Our analyses were divided in three steps. First, for each species, we fitted eight ENMs (assuming unlimited dispersal): BIOCLIM, Euclidean distance, Random Forest, Generalized Additive Model, Generalized Linear Model, Gower distance, Multivariate Adaptive Regression Splines and Maxent (see Rangel & Loyola, 2012). Using these models, we estimated the potential geographic distributions for baseline and future climates, considering three AOGCM (CCSM3, CSIRO-Mk3.0 (CSIRO), and UKMO-HadCM3 (UKMO); see Diniz-Filho et al., 2009 for details), for the scenario A1 (IPCC, 2000). In sum, for each species, we generated 24 values of potential geographic range size (eight ENMs projected into three AOGCMs), both for baseline and future climates (generating a total of 48 projections). The eight ENMs listed above were originally selected by Diniz-Filho et al. (2009) to study the variation in the results derived from the use of different modeling strategies, from simple (e.g., BLIOCLIM) to more computer-intensive methods (Random Forest).
Because a key aspect of our approach also consists in measuring the variation in range shift and, for comparative purposes, we use the same ENMs.
Second, for each species i, within a given AOGCM and for each ENM j, we estimated the proportional difference between potential geographic distributions for future and baseline climates (D ij ): where y ij is the future potential geographic distribution (number of 1° × 1° grid cells) for species i according to the ENM j; x ij is the baseline potential geographic distribution (number of 1° × 1° grid cells) for the same species i and ENM j. A negative D ij value suggests that species i will lose a proportion in its projected range of occurrence according to ENM j. These proportional differences were then averaged across the eight ENMs, generating D i . The D i was taken as a measure of effect size in our study (difference in geographic distribution; henceforth "effect size" for simplicity and to adhere to the standard nomenclature for meta-analysis). The uncertainty in the estimates was evaluated by the variance (V i ) over the eight values of D ij for each species within a given AOGCM.
Third, the weighted mean range shift (M), summarizing the results across the 1,205 bird species, was estimated assuming a random-effects model, where the weight (W i ) assigned to each species i was given by the reciprocal of the V i plus the between-species variance (T 2 ). In a typical meta-analysis, T 2 is an estimate of the betweenstudies variance in effect sizes (see Borenstein, Hedges, Higgins, & Rothstein, 2009;Borenstein, Higgins, Hedges, & Rothstein, 2017).
In our study, this statistic reflects the variability among the species in terms of their likely responses to climate change (difference between current and future potential geographic distribution). Thus, the larger the variance (V i ) associated to a D i , which, in our study, is due to the use of different ENMs, the smaller the weight of a species in estimating M. The larger the between-species variance (T 2 ), the larger our uncertainty over M.
Our random-effects model consisted of a multilevel model with bird phylogeny as a random effect (Lajeunesse, 2009;Nakagawa & Santos, 2012) and was estimated using the "rma.mv" function in the "metafor" R package (Viechtbauer, 2010). We estimated a phylogenetic variance-covariance matrix describing the expected variance (in the diagonal) and the covariances (in the off-diagonals) of a given trait following a Brownian-motion process (Lajeunesse, 2009), using the "vcv" function of "ape" package (Paradis, Claude, & Strimmer, 2004). We used this phylogenetic variance-covariance matrix to describe the correlation among species' responses in the random-effects multilevel model.
We estimated a T 2 statistic of variability between species (T 2 s ) and due to the phylogenetic random effects (T 2 p ) with the multi-level model (Nakagawa & Santos, 2012). We also estimated a statistic called I 2 , which indicates the proportion of variability in the effect sizes (i.e., range shifts) that comes from true differences between species and, therefore, can likely be explained by species-level variables. We further decomposed I 2 into two components: betweenspecies heterogeneity (I 2 s ) and phylogenetic-level heterogeneity (I 2 p ; Nakagawa & Santos, 2012).
To illustrate the possibility of using meta-analytic tools to explore the reasons for the variability in effect sizes (i.e., the among-species variation in range shifts), we evaluated whether the magnitude of the effects sizes (D i ) was related to the following set of variables: altitude midpoint, body size, species' IUCN categories of extinction risk, clutch size, and migratory behavior. For this, we applied random-effects multilevel meta-regression models, while controlling for phylogenetic effects (Lajeunesse, 2009;Nakagawa & Santos, 2012).
For each meta-regression model, we estimated a quantity analogous to the coefficient of determination (pseudo-R 2 ) following Borenstein et al. (2009). All analyses were performed in the R environment (R Core Team, 2018) using the packages cited above (Appendix S2).
F I G U R E 1 Variation of effects sizes among different atmospheric-ocean global circulation models projections. The vertical solid line indicates effect size equal to zero, and the dashed line indicates the weighted effect size for each projection.
Horizontal lines indicate 95% confidence intervals (95% CI) of each effect size. The size of each circle indicates the weight of each effect size for the weighted effect size calculation. Negative ES: effect sizes which the upper limit of 95% CI does not include zero; Positive ES: effect sizes which the lower limit of 95% CI does not include zero; Non-significant ES: effect sizes which 95% CI includes zero. For simplicity, we omit six, five, and eight very imprecise effects sizes in "a", "b," and "c," respectively There was significant heterogeneity in all AOGCM, and the highest values of the heterogeneity (T 2 ) was due to between-species differences (Table 1). According to the coefficient I 2 t , from 81.3% to 99.8% of the observed variance can be attributable to real differences among species in terms of geographic range shift (D i ). From 65.9% to 74.3% of the observed variance in D i can be due to real differences between species, and from 7.0% to 34.0% can be due to phylogenetic effects (Table 1).

| RE SULTS
The high variability in the differences between potential distributions for baseline and future climates suggests that there is much scope for the study of correlates of range loss/gain due to climate change.
IUCN categories of extinction risk, clutch size, and altitude midpoint (only for UKMO) were the main predictors of differences in species' range shifts ( Table 2)

| D ISCUSS I ON
We examined the potential effects of projected climate change on the geographic range sizes of Neotropical birds. Matching results obtained by many studies (Parmesan & Yohe, 2003;Walther et al., 2002; see Parmesan, 2006 for a comprehensive review), which measured phenological (e.g., advancement of spring events) and distributional changes (e.g., poleward shifts in geographical ranges, range contractions or expansions), we projected that effects of climate change on the range sizes of Neotropical birds can be negative. The negative sign of this forecasted pattern seems to be coherent (sensu Parmesan, 2006;Pecl et al., 2017) for the scale of the Neotropics.
Further, we showed that more than 80% of the observed variance in species-range size change can be attributable to real differences between species and, therefore, can potentially be explained by species-levels (explanatory) variables. However, finding correlates of range shift proved to be difficult.
Even under the unrealistic assumption of unlimited dispersal (which favors the hypothesis of no effect of climate change on the biota), we found that 509 (for CCSM3 and CSIRO) to 628 species (for UKMO) from all 1,205 bird species assessed are projected to experience range contraction. Further, species that were projected to lose less or even to increase their range of occurrence tended to present the most imprecise estimates. Our estimates of the relative number of bird species suffering more than 50% range contractions (4.4%, 3.7% and 27.8% of all bird species considered, respectively for CCSM3, CSIRO, and UKMO) are also comparable to those made by Jetz et al. (2007), who found that from 4.5% to 20.6% of the species analyzed are likely to lose more than 50% of their ranges.
However, an important message of our study is that we should not attribute the same weight for different species when assessing cross-species differences in potential distributions. This is so because, in the context of ensemble forecasting (Araújo & New, 2007),  2012;Sales et al., 2017;Thuiller, 2004). Our results considered these differences because we down-weighted species with variable range shifts estimates due to the use of different ENMs. The reason behind weighting is that species with low uncertainties in terms of range shifts will give estimates of the effects of climate change on geographic ranges that are more robust to the choices of ENMs. Also, weighting the estimates of range shift accordingly will likely increase the likelihood of correctly estimating the potential effects of climate changes. If we did not weight the D i values, we would observe a mean difference in potential distribution of 6.38% (95% CI = ±12.44%), 6.28% (±12.49%) and −4.63% (±16.88%) under CCSM3, CSIRO, and UKMO projections, respectively. Thus, without taking uncertainty into account, our estimates of the effect of climate change on geographic range size would be much smaller and, probably, downward biased. In general, we believe that taking uncertainty into account, while estimating the effects of climate change on species distributions, will be an important step in rebutting climate change denialism (e.g., Boussalis & Coan, 2016).
Our meta-analytic approach provides an alternative method to quantify the relative variation in estimates of range shifts within (i.e., due to the use of different ENMs and AOGCMs) and among species.
Our results are in line with those reported previously by showing that a substantial part of the variation in range shifts (from 25.7% up to 34.1%) can be attributed to the use of different ENMs (for a recent study, see Thuiller et al., 2019). Further, the geographic range losses, as given by the weighted mean effect sizes (M), were estimated to be much higher (ca. 57%) under UKMO than under CSIRO. However, our findings also highlight that most of the variation in range shifts, for a given AOGCM, can be attributable to real differences between species and that a small part of this variation is phylogenetically structured (except under UKMO). Taken together, our results suggest that the main question (or source of uncertainty) is not whether many Neotropical bird species will lose a large proportion of their climatic space (as it has also been projected for northern-boreal land bird; see Virkkala, Heikkinen, Leikola, & Luoto, 2008), but "which species will be under most threat" (Sinclair, White, & Newell, 2010).
The search for correlates of range shifts under contemporary climate change has been an active research line in ecology (e.g., MacLean & Beissinger, 2017;Williams & Blois, 2018). The method we used to decompose the variation in range shifts among species also allowed us to infer that phylogenetic relatedness is a poor predictor of range shifts (but see Comte, Murienne, & Grenouillet, 2014 for a contrasting result with stream fishes). We found that bird species with larger clutch size are expected to lose a smaller range of occurrence than species with smaller clutch size. This pattern was consistent across all AOGCM and with studies showing that species with higher fecundity may be buffered against a decrease in their overall range of occurrence (e.g., Amano & Yamamura, 2007).
Bird species with higher threat of extinction (higher IUCN categories of extinction risk values) may lose a larger range of occurrence projected for CGCM and CSIRO. For UKMO, these species were projected to lose a smaller range of occurrence than species classified with a smaller extinction risk. Altitudinal midpoint correlated F I G U R E 3 Relationship between raw mean difference (D, %) and species midpoint of altitude (a), IUCN categories of extinction risk (b) and clutch size (c) for atmospheric-ocean global circulation models projections of UKMO. Circle size indicates the weight of each effect size to meta-regression parameter estimates. Fitted lines represent partial effects of each moderator variable on D positively with differences in range of occurrence. However, studies have shown that mountain bird species will lose a greater range of occurrence than other species (e.g., La Sorte & Jetz, 2010;Freeman, Scholer, Ruiz-Gutierez, & Fitzpatrick, 2018; see also a meta-analysis by Scridel et al., 2018 argue that our estimates of potential geographic range loss would be lower after allowing for adaptation. However, recent studies indicate that the effect of adaptation (evolutionary rescue) may be limited (Cotto et al., 2017;Diniz-Filho et al., 2019). Valladares et al. (2014) showed that the forecasted effects of climate change on species range shifts are expected to even increase when phenotypic plasticity is incorporated into ENMs. Our estimates of variation in range shifts are likely to be conservative because a recent study by Thuiller et al. (2019) showed that ENMs differ greatly in predicting future distributions even when only ENMs with high predictive accuracies are considered in the analyses.
Despite the uncertainties associated with ENMs, which were considered in the meta-analytic approach implemented here, our results suggest that the effects of human-induced global warming on the geographic distributions of Neotropical birds are concerning. We believe that a promising line of investigation would consist in finding other species traits with high capability to predict likely ranges shifts (as forecasted by ENMs). For instance, an additional bird trait not considered here is the extent of breeding seasons. In a recent meta-analysis, Halupka and Halupka (2017), for example, found that climate change may extend breeding seasons for multibrooded species and reduce it for single-breeding species from the northern hemisphere. Species reproducing for a longer period may increase their population size and consequently the range that this species occurs. Finding strong correlates of range shift would help policymakers to focus not only where the effects of climate changes are likely to be stronger (e.g., Jetz et al., 2007;Lawler et al., 2010;Lawler et al., 2009), but also to list the species (as predicted by their traits) that are at higher risk.

ACK N OWLED G M ENTS
Our work on distribution modeling has been supported by various

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
All the data and script used are available at the Dryad Digital Repository (https ://doi.org/10.5061/dryad.qv5p3r8).