Empirical Predictability of Community Responses to Climate Change

of such ecology the assumption that community states and dynamics empirical from and suggest that species communities rarely follow equilibrium dynamics with climate change. This discrepancy between the conceptual foundation of many predictive models and observed community dynamics casts doubts on our ability to successfully predict future community states. Here we used community response diagrams (CRDs) to empirically investigate the occurrence of different classes of disequilibrium responses in plant communities during the Late Quaternary, and bird communities during modern climate warming in North America. We documented a large variability in types of responses including alternate states, suggesting that equilibrium dynamics are not the most common type of response to climate change. Bird responses appeared less predictable to modern climate warming than plant responses to Late Quaternary climate warming. Furthermore, we showed that baseline climate gradients were a strong predictor of disequilibrium states, while ecological factors such as species’ traits had a substantial, but inconsistent effect on the deviation from equilibrium. We conclude that (1) complex temporal community dynamics including stochastic responses, lags, and alternate states are common; (2) assuming equilibrium dynamics to predict biodiversity responses to future climate changes may lead to unsuccessful predictions.

Robust predictions of ecosystem responses to climate change are challenging. To achieve such predictions, ecology has extensively relied on the assumption that community states and dynamics are at equilibrium with climate. However, empirical evidence from Quaternary and contemporary data suggest that species communities rarely follow equilibrium dynamics with climate change. This discrepancy between the conceptual foundation of many predictive models and observed community dynamics casts doubts on our ability to successfully predict future community states. Here we used community response diagrams (CRDs) to empirically investigate the occurrence of different classes of disequilibrium responses in plant communities during the Late Quaternary, and bird communities during modern climate warming in North America. We documented a large variability in types of responses including alternate states, suggesting that equilibrium dynamics are not the most common type of response to climate change. Bird responses appeared less predictable to modern climate warming than plant responses to Late Quaternary climate warming. Furthermore, we showed that baseline climate gradients were a strong predictor of disequilibrium states, while ecological factors such as species' traits had a substantial, but inconsistent effect on the deviation from equilibrium. We conclude that (1) complex temporal community dynamics including stochastic responses, lags, and alternate states are common; (2) assuming equilibrium dynamics to predict biodiversity responses to future climate changes may lead to unsuccessful predictions.

INTRODUCTION
Contemporary climate change impacts the dynamics of biodiversity (Parmesan, 2006;Steinbauer et al., 2018) but our ability to predict these impacts remains limited. Many fields of ecology have historically relied on the concept of equilibrium to study and forecast the responses of biodiversity to climate change. The dynamic equilibrium hypothesis assumes that species distributions and assemblages reflect a climate niche optimum in which species climate niches match the observed climate, and that changes in climate induce changes in community composition or species distribution to stay close to this equilibrium state with climate (Webb, 1986;Prentice et al., 1991). This hypothesis assumes a linear relationship between climate and species climate niches, with limited presence of lags, threshold effects, stochastic variations, and transient dynamics in biodiversity responses to climate changes. These processes probably impair the community responses expected from the equilibrium dynamics hypothesis (Jackson and Overpeck, 2000;Krauss et al., 2010). There is growing evidence that biotic responses observed in nature do not match those expected under the assumption of dynamic equilibrium with climate change (e.g., Devictor et al., 2012;Svenning and Sandel, 2013;Ash et al., 2017).
The dynamic equilibrium hypothesis has provided the conceptual foundation of most anticipatory predictions of biodiversity responses to climate change (i.e., a prediction intended to deduce future states from a given model, sensu Mouquet et al., 2015). The underlying assumption is that species' climate tolerances are constant over time (Pearman et al., 2008;Wiens et al., 2010) and that species ranges reflects climatically suitable areas (Leroux et al., 2013;Stephens et al., 2016). Therefore, species ranges are expected to track climatic change (Parmesan et al., 1999;La Sorte and Jetz, 2012), triggering local turnover in community compositions based on species climate tolerance (Devictor et al., 2008;Gaüzère et al., 2015). This is supported by empirical results at multiple scales and in many taxa. Some of the reported examples from the literature include multi-taxon responses to Younger Dryas climate changes in Switzerland (Birks and Ammann, 2000), woody species responses to late Quaternary climate warming (Jackson and Overpeck, 2000); bird (Tingley et al., 2009) or marine taxa (Pinsky et al., 2013) responses to modern climate change.
However, the relevance of the dynamic equilibrium hypothesis has also been challenged. The assumption of species rangeclimate associations is not strongly supported in all taxa (Jackson and Overpeck, 2000;Beale et al., 2008); species might not be able to keep up with the velocity of climate change (Devictor et al., 2008;Bertrand et al., 2011;Svenning and Sandel, 2013); and changes in available climatic space, habitats or biotic interactions might affect expected responses to climate change (La Sorte and Jetz, 2012;Maiorano et al., 2013;Wisz et al., 2013). In general, ecological systems might intrinsically be unpredictable because of their complexity and the amount of chaotic, neutral, or stochastic processes impairing their dynamics (Petchey et al., 2015). In consequence, climate change only partly explains the dynamics of species and communities since the Last Glacial Maximum (Veloz et al., 2012;Blois et al., 2013) or during modern climate change (Zhu et al., 2012;Ash et al., 2017;Currie and Venne, 2017). While most projected shifts in species distributions or biodiversity (e.g., Thuiller et al., 2011) have relied on equilibrium dynamics, there is no general consensus about the taxonomic, spatial or temporal scales at which this assumption is reasonable.
Delineating the limits of predictability and the presence of non-linear responses is a critical prerequisite for advancing predictive ecology in the Anthropocene (Mouquet et al., 2015).
Contemporary climate change highlights the increasing need to forecast the future state of populations, communities and ecosystems to better inform conservation strategies (Clark et al., 2001;Mouquet et al., 2015;Petchey et al., 2015). While tools for researching and communicating ecological predictability already exist (Petchey et al., 2015;Blonder et al., 2017), there is still a weak empirical understanding of when and why predictability could be possible in communities (Blonder et al., 2018). The goal of this study is to provide an empirical assessment of whether and when anticipatory predictions of community responses to climate change are a reachable goal.
While particular scenarios like equilibrium dynamic responses might be predictable, other type of responses might not. Different response scenarios such as constant-lag dynamics or alternate stable states can lead to a deviation from equilibrium with climate condition (Blonder et al., 2017). As a consequence, no-lag or constant-relationship dynamics, where the community response follows the observed climate with a fixed time delay, are predictable. However, transient dynamics and alternate unstable states are not predictable. While recent theoretical work has identified and defined a broad range of possible scenarios (Blonder et al., 2017), limited empirical work has explored the presence of these different scenarios.
Here we address this gap by documenting the relationship between temperature forcing and responses in community composition within birds and terrestrial plants of North America. We sought to delineate contexts in which predictability is reachable by (i) investigating the response of plant communities since the Last Glacial maximum (−21 Ka-present) and the contemporary responses of bird communities to recent climate changes (1970-2012 C.E.), and (ii) understanding how the predictability of community responses to climate change covaries with climatic gradients, human pressures, and dispersalrelated traits.
We studied community predictability (here defined as the ability to provide anticipatory predictions of a community state from climate observations or projections) via time series analyses of the relationship between environmental forcing and the community response. We use a newly developed framework in which response scenarios are detected by sequentially plotting time series of observed and community-inferred climate values, the community response diagram (CRD) framework (Figure 1, Blonder et al., 2017). CRDs can be used to detect deviations between a climate forcing, such as temperature change, and a corresponding community state response (Figure 1). First, for a given site monitored through time (Figure 1A), the community state is defined as the average niche temperature value of all species composing a community at a given time ( Figure 1B). Community-inferred temperature is similar to a community temperature index (Devictor et al., 2008;Lenoir et al., 2013), a floristic temperature (De Frenne et al., 2013), a coexistence interval (Mosbrugger and Utescher, 1997), or an Ellenberg indicator value (Ellenberg and Mueller-Dombois, 1974). Secondly, these community responses are paired with the observed climate change for a given site/time ( Figure 1C). A CRD is the sequential time series plot of communityinferred temperatures and observed temperatures ( Figure 1D). CRDs can be summarized using three statistics that provide complementary insights into community dynamics ( Figure 1E). First, the absolute deviation,¯ , depicts the time-average absolute deviation of a community's state from its expected state under equilibrium with climate. If¯ = 0 then the observed community follows a 1:1 relationship with climate. Second, the deviation change, d , measures the temporal change of the deviation during the survey. It indicates whether and how much the community gets further or closer to the equilibrium state during the study period. If d = 0 the expected community response to climate varies across the time series. Third, the maximum state number, n, counts the maximum number of community climate states observed for a given temperature over the period. When n > 1 there is more than one community state for a given climate value. If there is only one value of community-inferred temperature corresponding to each value of observed temperature, then n = 1 ( Figure 1D). When n = 1 the community has dynamics that can always be predicted from knowledge of the observed temperature. If n > 1, more than one community-inferred temperature value exists for a given observed temperature value. It is therefore not possible to predict the community's state with knowledge only of the observed temperature. From a CRD, the combination of the three statistics characterizes different theoretical response scenarios and quantifies different aspects of predictability (see Box 1).
We applied the CRD framework to species community data in two contrasting settings, (i) the long-term dynamics of plant communities during the Late Quaternary climate warming (21 Ka-present) and (ii) the contemporary responses of bird communities to recent temperature changes (1971-2012 C.E.) ( Figure S1). We hypothesized that the short-term responses of bird communities to Anthropocene climate change would exhibit a lower predictability than the long-term thermal reshuffling of plant communities during the late Quaternary, for three reasons.
First, short-term variation might be harder to predict than long-term changes because stochasticity is most predominant at fine spatial and temporal scales (Levin, 1992). Second, shortterm resistance to unfavorable conditions (Fordham et al., 2016), delayed effect of climate change via indirect effects (e.g., by changing habitats or resource availability, Gaget et al., 2018) and the time needed for species to track climate changes (Alexander et al., 2018) are all expected to increase mismatch in community responses to climate change. In contrast, changes in average temperatures across longer time periods are expected to be balanced with the regional species pools, due to an extended period in which climate induced extinction and colonization could occur (Webb, 1986;Holm and Svenning, 2014). Third, birds might exhibit lower predictability because they are less sensitive to climate change than plants. They have broader thermal tolerances (Araújo et al., 2013), and the relevance of dynamic equilibrium responses as a conceptual model to explain birds responses to climate change is controversial. For example, their realized niche may not necessarily match the observed climate (Beale et al., 2008). Evidence of breeding bird responses to recent climate change are mixed and appear highly idiosyncratic (Stephens et al., 2016;Currie and Venne, 2017). Climatic data are used to extract observed temperature of the site through time (C), and the community-inferred temperature through time. These time series are combined in a time-implicit parametric plot to build the CRD (D). For each time bin, deviation from the 1:1 line value and the number of community-inferred value existing for each observed temperature value are used to compute three statistics (E) summarizing different aspects of predictability.
We also hypothesized that different aspects of community predictability will be structured by environmental and ecological factors, yielding the following predictions: • More extreme temperature (warm or cold) are linked to stronger disequilibrium in community responses to climate change, with low¯ , negative to null d and n = 1. In the coldest (northern) and warmest (southern) areas of North America, regional species pools might not provide climateadapted species to local communities (Bertrand et al., 2016) and therefore limit their response to temperature change (Blonder et al., 2015) • Complex topography is linked to higher predictability, with low¯ , negative to null d and n = 1. In mountainous terrain, we expect strong variation of temperature within small distances to lower the perceived velocity of large scale climate change (Loarie et al., 2009) and therefore increase community thermal response (Bertrand et al., 2011;Gaüzère et al., 2017). • Human impact is linked to lower predictability, with stronḡ , low d and n > 1. Human impact (i.e., direct exploitation, species introduction, and land conversion) influences community dynamics and also the predictability of community responses to climate change (Maxwell et al., 2016;Bowler and Böhning-Gaese, 2017;Gaget et al., 2018).
Beyond these external factors, we expect biological intrinsic factors to influence community predictability. Assuming that the community mean reflects the local species pool, we predict that: • Species characteristics increasing dispersal ability are linked to higher predictability, with low¯ , negative to null d and n = 1. Communities consisting of species with a high dispersal potential (lower seed mass, migrant species) respond quickly to climate change (Jenni and Kéry, 2003;Svenning and Skov, 2007). • Species characteristics increasing persistence in unfavorable environments are linked to lower predictability, with stronḡ , c. null d and n > 1. Life history traits such as longevity, adult height or bodymass increase resistance and persistence to unfavorable climate conditions, therefore decrease the predictability of community.

Community Data
Plant assemblage composition data across the Late Quaternary and Holocene (21 Ka-present) were compiled from the Neotoma paleoecology database (Goring et al., 2015), relying on the fossil pollen data sets used in Maguire et al. (2016). Sites represent highquality assemblages and were primarily located in eastern North America. We used Blois et al. (2011) selection of sites and revised chronologies. Pollen assemblages obtained from lake sediments can provide a rough proxy for the composition of communities, despite issues on spatial and taxonomical scale integration, species abundance vs. pollen abundance, and detectability of rare taxa (Birks and Seppä, 2004). These issues might influence the quantification of disequilibrium state because rarer taxa with lower dispersal abilities might be undetected. More details on the data processing, spatial and temporal distributions of site are provided in Figure S1. The overall process yielded a presence/absence dataset comprising 425 sites, 103 plant taxa, and 45 time bins (500 yr each) spanning 21 Ka-present.
Bird assemblage composition across the last 50 years were compiled from the North American Breeding Bird Survey (BBS, Sauer et al., 2013, data and protocol at http://www.pwrc.usgs.gov/ bbs/). We used data processed from Barnagaud et al. (2017): firstyear observer effects were removed by excluding the first survey performed by a given observer on a given route. The dataset was restricted to 807 routes monitored at least 8 years and once every 5 years during the 1970-2011 period. Coastal, pelagic and species which accounted together for <1% of the records were not taken into account (87 species). The overall process yielded an abundance dataset comprising 807 routes for up to 435 birds species over 43 years. More details on the data processing, spatial and temporal distributions of the sites are provided in Figure S1.

Temperature time series in North-America
Paleoclimate time-series data were obtained from SynTraCE-21, a set of transient simulations runs using the CCSM3 model (Liu et al., 2009), with decadal averages stored from 22 Ka to the present. These simulations are reasonably congruent with sitebased climate reconstruction (Harrison et al., 2014). Simulations were statistically downscaled to a 0.5 • by 0.5 • grid cell following the CRU TS 3.20 dataset, and then for every 500 years from 21 to 0 Ka, average maximum temperature was calculated based on a 200 year window centered on the 500 year time step (Lorenz et al., 2016).
Contemporary temperature time-series data were extracted from the Parameter Regression of Independent Slope Model (PRISM) (Oregon Climate Service, Corvallis, OR, USA) for the continental United States (Daly et al., 2000). Both climate datasets were created using point meteorological station data, digital elevation models, and other spatial data sets to generate interpolated gridded estimates of monthly, yearly and eventbased climatic parameters. Maps had a spatial resolution of 2.5 arcmin (∼3 km cell size at this latitude), as generated by PRISM (Daly et al., 2000).

Global temperature contemporary distribution
Spatial global temperature data used to estimate climate niches (see below) were extracted from WorldClim (http://worldclim. org/version2) raster at 2.5 min of a degree resolution (Fick and Hijmans, 2017). WorldClim rasters at 2.5 min degree provide the best independent descriptors available at global scale, despite the Frontiers in Ecology and Evolution | www.frontiersin.org 4 November 2018 | Volume 6 | Article 186 Box 1 | Response scenarios between species communities and climatic conditions expected from theory, compiled from Blonder et al. (2017).

No-lag scenario
The climate preferences of species in a given community closely match the observed climate. Changes in climate causes an immediate change in community composition so that the climate preference of the present species matches the new climatic conditions.

Expected statistic values :¯
The community response to climate change has a one-to-one mapping such that the assemblage of species has a unique climatic preference for any observed climate value.
Expected statistic values : The climate preference of a community follows the observed climate with a fixed time delay. When the response to climate change is linear, this scenario reduce to the constant-relationship scenario. In other cases, the lag could follow a periodic function or the future response of the system might depend on its past history, e.g., via memory effects (hysteresis).

Alternate unstable states
A generalization of the constant-lag scenario in which the community shows memory effects and the future response to climate depends on the past state of the system. In this scenario, community dynamics follow the observed climate with a variable delay and magnitude. Under such conditions several species assemblages with similar climate preference combinations could occur together with a single value of observed climate. Thus, the future dynamics of the community cannot be predicted only by knowing the current community state.

Stochastic or chaotic dynamics
Community dynamics are uncorrelated with the observed climate, for example due to stochastic dynamics. As a consequence, many species assemblages with different climate preferences might occur at a given observed climate value.
Expected statistic values :¯ > 0; d → 0; n > 1 uncertainty and error associated with the values (Hijmans et al., 2005). We used the 2.5 min degree resolution spatial resolution (this is about 4.5 km at the equator) in order to erase microclimatic variations and local characteristics of the habitat when inferring species thermal niches.

Climate Niches
We inferred realized plant climate niches for mean maximum temperature for each taxon using independent contemporary occurrence data for each species or genus in the paleocommunity dataset from the BIEN3 database. We choose maximum temperature as a unique climatic niche axis because it is a strong predictor of the species' responses to temperature increase (Jiguet et al., 2007;Lorenz et al., 2016;Maguire et al., 2016). Note that maximum temperature was strongly correlated with mean temperature and minimum temperature in both paleo and contemporary temperature data (correlation coefficients ranging from 0.89 to 0.98). BIEN3 contains more than 30,000,000 geo-referenced vascular plant observations (Enquist et al., 2009) from a much broader geographic scope than represented by the pollen dataset. Contemporary occurrence data were filtered to include only New World records that did not come from cultivated areas. Genus-level distributions were pooled for taxa with multiple names (e.g., "Ostrya/Carpinus"). For each taxa, we estimated the niche maximum temperature value as the mean of the average monthly maximum temperature values ( • C, extracted from the WorldClim raster) over all sites where the taxa was detected. Such realized-niche estimates can be affected by sampling bias, but have the benefit of being estimable from broad-scale and commonly available presence-only data. To ensure that sampling bias does not affect our estimates, we tested the correlation between species' thermal niches estimated from raw data bounded to the 95% confidence interval (i.e., incorporating the density of sampling) with species' thermal niches sampled from a uniform distribution bounded to the 95% confidence interval of real distribution (i.e., without sampling density). A Pearson's product-moment correlation test showed a strong positive correlation between species thermal niches estimated from the different approaches (correlation coefficient ± CI = 0.9967 ± 0.001, t = 120.32, df = 97, p < 0.0001). We concluded that spatial sampling bias is unlikely to affect thermal niche estimates.
We recognize that realized niches for plants may shift over 10 3 -10 4 year timescales (Veloz et al., 2012). To check the correlation between paleo and modern niches estimations, we investigated the relationships between paleo-inferred and contemporary-inferred climate niches ( Figure S2) and showed that they were closely related (Pearson's correlation coefficient = 0.96, t = 39.9, df = 94, P < 0.0001).
We inferred American bird species' thermal niches by clipping their global extent of occurrence maps onto temperature layers from WorldClim. We used an independent dataset combining global extent-of-occurrence maps for 9,886 bird species ( Birds of the World, 2017). These data have been compiled from multiple sources, including specimens, distribution atlases, survey reports, published literature, and expert opinion, and currently represent the most comprehensive assessment of global bird species occurrences. BirdLife International and HBW endeavor to maintain clean, accurate, and up-to-date data at all times. However, such vector maps do not contain within-range heterogeneity in species' presence, and errors on range margins might be present, especially on ill-sampled avifaunas (Herkt et al., 2017). For each taxon, we estimated the niche maximum temperature value as the breeding-season (May, June, July) average monthly maximum temperature values ( • C, extracted from the WorldClim raster) over the breeding distribution range of taxa.

Community Response Diagrams
We used community-inferred temperatures and observed temperatures to build CRDs for the plant and bird datasets, as shown in Figure 1.
First, we computed community-inferred temperature values using the climate niches previously inferred via the species distribution data. For each site monitored at a given time, the community-inferred temperature corresponds to the average climate niche value of all species present in a community. community-inferred temperature was calculated as the average maximum temperature niche value of all species present in a community. We also computed the standard deviation associated with each community-inferred community value.
For birds, we additionally calculated community-inferred temperature based on the abundance-weighted mean of specie's maximum temperature values. This method takes into account the relative abundance of each species in the community and could provide a more sensitive estimation on community-inferred temperature change. Abundance weighted community-inferred temperatures were strongly correlated to estimates based on occurrence only (see Figure S3, Pearson's correlation coefficient = 0.95, t = 540.94, df = 32264, P < 0.0001). We therefore used occurrence-based communityinferred temperature in order to increase the comparability of results between plants and birds.
We then extracted observed maximum temperature for each site and each time bin from paleoclimate and modern time-series data (see Data subsection). In order to reduce the inter-annual variability and temporal stochasticity that might undermine the identification of community and climate dynamics, we smoothed community-inferred and observed temperature times series by using Local Polynomial Regression Fitting (LOESS). LOESS was performed for each time series independently using the loess {stats} R function with an α span of 0.75.
Finally, we built CRDs for each site by sequentially plotting raw and smoothed time series of observed and communityinferred climate values, as shown in Figure 1D. We used CRDs to estimate three complementary statistics depicting predictability of community responses to climate change. For each CRD, the absolute deviation,¯ , was calculated as the absolute value of the average deviation, where the deviation is calculated, for each time value, as the difference between the smoothed value of community-inferred temperature and the smoothed value of observed temperature ( Figure 1D). The deviation change, d , was calculated as the difference of deviation between the last and the first time bin of the CRD, divided by the overall timespan of the CRD. The state number, n, was calculated as the number of smoothed community-inferred temperature values (y-axis in CRD) that correspond to a given single value of observed temperature (x-axis in CRD). n counts the number of times a vertical line on the diagram crosses the communityinferred temperature trajectory. Temporal stochasticity and sampling error tend to inflate the state number n through the detection of more than one community state which are statistically not different. To correct for this false detection of n > 1, we tested for the difference between communityinferred temperature values by comparing the difference between 95% confidence interval associated to each community-inferred value. If the 95% confidence interval were overlapping, we inferred that the community-inferred temperature values could not be differentiated, and reduced n by 1 (further called corrected n). More details, formalization, and simulations of CRDs are provided in Blonder et al. (2017). R scripts and functions used to compute these summary statistics and other descriptive metrics can be found at https://github.com/pgauzere/ Predictability_CRD.

Predictors
For both plants and birds, we gathered a set of five explanatory variables associated to community responses to changes in temperature: two abiotic variables (i.e., baseline temperature and topography), one variable related to human influence on ecosystems (i.e., paleo human density for plants, contemporary human influence index for birds), and two variables related to species characteristics known to influence tracking and resistance processes (i.e., seed mass and plant height for plants, body mass, and migration for birds). The distributions and mapping of predictors variables are shown in Figure S4. Note that none of our pairs of predictors were strongly correlated (maximum correlation coefficient was 0.41 for baseline climate-topography).

Baseline Temperature
To test for the effect of climate on community predictability, we extracted the baseline maximum temperature of each site as the first observed maximum temperature value recorded for a given site. See Temperature Data section.

Topography
To estimate topography, we extracted digital elevation data (DEMs) from the NASA Shuttle Radar Topographic Mission (SRTM). The NASA SRTM provides elevation as 3 arc second (∼90 m) resolution. For each site, we calculated a topography index as the difference between highest and lowest elevation within a 50 km buffer around the site. We tested the influence of the buffer size on our estimate by computing topography on 100 and 150 km buffer and comparing these values to the 50 km buffer.

Human Imprint on Ecosystems
To estimate the human imprint on ecosystems over the Holocene period, we obtained population density data over the Holocene for the period ranging from 21 to 0 Ka from the HYDE 3.2 database (Klein Goldewijk et al., 2011). These data are chosen to provide an approximate proxy for the upper bound of human impacts on communities across the last 21 Ka. We extracted population density value for each site from Neotoma dataset. To estimate the contemporary human imprint on ecosystems, we used the Global Human Influence Index (HII) from the NASA Socioeconomic Data and Applications Center. This index estimates the direct human footprint on ecosystems (Sanderson et al., 2002) as a proxy of recent human pressure on biodiversity. HII incorporates nine global data layers corresponding to population density, human land use, infrastructure, and human access. We extracted HII values for each site from BBS dataset using the raster library.

Species Characteristics
For plants, we computed the community mean seed mass and plant height (Honnay et al., 2002;Matlack, 2005;Normand et al., 2011) using traits values compiled from the TRY database (Kattge et al., 2011). For birds, we used mean body mass and the percent of migrants of each community (Jiguet et al., 2007;Angert et al., 2011) using traits from the Encyclopedia Of Life (http://www. eol.org), the Animal Diversity Web (http://www.animaldiversity. org) and the field guide to North American birds (Sibley, 2014) previously compiled in Barnagaud et al. (2017).

Statistical Analysis
We tested the effect of our set of predictors on each CRD summary statistic by implementing generalized additive models (GAM). We ran six models in which each statistic was the response variable (¯ , d , and n. for both plants and birds) regressed over a set of five predictors (for plants: baseline climate, topography, human density, mean seed mass, means plant height; for birds: baseline climate, topography, human influence index, mean body mass, percent of migrant species). We estimated the linear effect of each predictor, and added a two-dimensional spline based on geographic coordinates in order to account for spatial autocorrelation (Wood, 2006). Because we predicted the effect of baseline climate to occur at coldest and warmest conditions, our models included both linear and quadratic terms for baseline climate. Because state number (n) values take positive integer values, we used Poisson GAMs to analyze these data. Strong outliers corresponding to site with very low number of taxa, or strong variations of number of taxa through time were deleted for the analysis (between one and five points depending on models). P-values reported for parametric and smoothed model terms were based on Wald tests. All analyses were performed using the mgcv library in R statistical software (R Development Core Team, 2013). The code written for data analysis can be accessed at https://github.com/pgauzere/ Predictability_CRD.

RESULTS
We found strong support for disequilibrium dynamics in both plant and bird responses to climate change. We described a wide variation in the type of community dynamics, with CRD depicting potential evidence for all scenarios listed in Box 1. CRD for a few representative communities are shown in Figure 2. All CRDs, along with associated time series and statistics values are shown in Figure S5.
The main component explaining the three summary statistics from the CRD was baseline climate. In the multivariate models, baseline climate was consistently related to variations of CRD statistics (baseline climate has a significant effect on statistics in five over six GAMs performed, Figure 2). Apart from baseline climate, no clear general pattern emerged from other predictors. The geographic two-dimensional splines smooth terms included in the model substantially increased the fit of the models for the absolute deviation and the deviation change but did not increase the fit of the model for maximum state number. Overall, our set of predictors explained from 20.9 to 80.5% of the deviance (apart from n for birds).

Scenarios of Responses
Plant communities generally exhibited lagged monotonic and positive relationship between inferred and observed temperatures during the last 21 Ka. A few communities (2.4%) showed no-lag (or low-lag) dynamics, where relationships between inferred community and observed temperature are linear and close to the 1:1 line. These communities were characterized by¯ <5 • C, −2 < d < 2 • C.Ka −1 , and n = 1 (Figure 2A). Many communities (60.0%) showed approximately constant-relationship dynamics, where increases in temperature lead to monotonic and non-linear change in community-inferred temperature. These communities were characterized by¯ >5 • C, d <0 • C.Ka −1 and n = 1 (Figure 2B). In plant communities, a particular constant-relationship dynamic was commonly observed. It was characterized by communities having strong baseline deviation from climatic equilibrium that did not exhibit substantial variation in community-inferred temperature while observed temperature strongly increased. Such dynamics led to a strong decrease in deviation through time ( Figure 2B). Some communities (23.0%) showed approximately constantlag dynamics with partial or complete loops (Figure 2C) or stochastic dynamics (Figure 2D). These communities were characterized by −2 < d < 2 • C.Ka −1 and n >1 (Figure 2C).

Summary Statistics
Generally, plant communities did show a large deviation between observed and inferred climate (¯ = 8.37 ± 3.4 • C [mean ± sd], Figure 3) and did not follow a 1:1 equilibrium state with climate.¯ values were structured in space, with high values of at northern and southern latitudes and lower¯ (< 5 • C, see Figure 3) at mid latitudes. Across the study sites this deviance from equilibrium was constant through the last 21 Ka (d = −0.54 ± 0.63 • C.Ka −1 ). However, ∼85% of the communities did show negative d values (suggesting that the amount of climate Frontiers in Ecology and Evolution | www.frontiersin.org disequilibrium did decrease at these sites). Low d values did peak at mid latitudes (Figure 3). The maximum state number n ranged from 1 to 4. 24.5% of the sites had more than one state number for a specific observed temperature temperature during the time series. Thus, in a quarter of the communities the inferred temperature value can have multiple values for a single observed temperature. Maximum state number n tended to peak in boreal latitudes.
Frontiers in Ecology and Evolution | www.frontiersin.org 8 November 2018 | Volume 6 | Article 186 , estimated as linear coefficients (± 95% confidence intervals) from generalized additive mixed models. Topography and Human Density were square-root transformed. All predictors were scaled to mean = 0 and SD = 1 prior to modeling to ease comparisons. Point and bar colors indicate the significance level associated to the test (light green: non-significant; light blue: significant at α = 5%; dark blue: significant at α = 1%).

Scenarios of Responses
Bird communities generally exhibited non-directional and stochastic dynamics in climate responses between 1966 and 2011. A few communities (1.4%) showed no-mismatch or low-mismatch dynamics, where relationships between inferred community and observed temperature are linear and close to the 1:1 line. These communities were characterized by¯ <2 • C, −20 < d < 20 • C.Ka −1 , and n = 1 (Figure 2E). A few communities (3.4%) showed approximately constant-relationship dynamics, where increases in temperature lead to monotonic and non-linear change in community-inferred temperature. These communities were characterized by¯ >2 • C, d < −20 • C.Ka −1 and n = 1 (Figure 2F). A few communities (2%) showed approximately Frontiers in Ecology and Evolution | www.frontiersin.org constant-lag or stochastic dynamics with partial or complete loops. These communities were characterized by −20 < d < 20 • C.Ka −1 , and n >1 (Figure 2G). However, many bird communities (c. 60%) exhibited non-directional and stochastic dynamics of response to observed temperatures (strong absolute deviation and low deviation change) with n = 1 (Figure 2H). Such stochastic dynamics associated with low state number are mainly due to the correction applied to the estimation of n, and are hard to identify (see Figure S4).

Summary Statistics
Bird communities showed a consistent deviation from climate equilibrium (¯ = 2.41 ± 1.86 • C, Figure 3).¯ was structured in space, with higher values of¯ in the south-eastern part of North America. On average, equilibrium states did not change in the bird communities between 1966 and 2011 (d =10.06 ± 19.64 • C.Ka −1 ). However, 73.2% of the sites did show positive d values (e.g., an increase in climatic disequilibrium). The lowest negative values of d were present in Western North America, while positive values were generally distributed in south-eastern North America.
Ninety-eight percent of the sites had a maximum state number n = 1, indicating that predictability is reachable for a large part of the communities. However, over this interval there was limited directional change in climate (see CRD), which would also be consistent with slow but stochastic dynamics. The low n value was mainly due to the correction we applied to take into account the strong stochasticity associated with these dynamics (see sections Methods-Community Response Diagrams and Discussion). On average, the uncorrected n value was strong (mean ± sd = 3.2 ± 0.76), with 99% of sites having n > 1 and 88% of sites having n > 2. The max state number n did not exhibit any spatial pattern.
Maximum state number n was not related with any of our predictors (Figure 4B, right panel). The geographic splines smooth terms were not improving the fit of the model (edf = 0, Chi.sq = 0, P = 1). The full model explained 2.3% of the deviance (1% without geographical spline).

DISCUSSION
We explored the limits and the determinants of predictability in community responses to climate change in bird and plant assemblages using CRDs. Currently, anticipatory prediction of biodiversity responses to climate change have considered a limited range of dynamics, relying on predictable relationships between species or community dynamics and climate change. While the no-lag equilibrium hypothesis is the implicit foundation of species distribution modeling (Peterson et al., 2011), only recent extensions of this method has successfully considered constant lag or constant relationship by incorporating dispersal limitation and/or properties of species assemblages (Guisan and Rahbek, 2011;Zurell et al., 2016). We here provided potential evidence for all types of community dynamics (see Box 1), including unpredictable dynamics (e.g., alternate states and stochastic dynamics) which are often not considered in current modeling approaches. Our work suggests that the current understanding of community dynamics in relation to climate change is oversimplified. Among the responses described in our study, equilibrium dynamics were the exception rather than the norm. This result challenges the equilibrium dynamic as the fundamental concept for predictive models of biodiversity response to climate change.
Along with the equilibrium dynamics hypothesis, the spacefor-time substitution approach has often been used to predict the effects of climate change on biodiversity. Although this assumed equivalence may be relevant in situations where equilibrium dynamics prevails (Walker et al., 2010), many studies have emphasized substantial differences between spatial and temporal responses (Johnson and Miyanishi, 2008). Our work suggests that because many community dynamics are diverging from equilibrium, the space-for-time substitution approach should be used with caution to infer future community state. Temporal dynamics might provide fundamentally different insights than spatial patterns (Bonthoux et al., 2013;Bjorkman et al., 2018). While ecology is undergoing a major transformation to leverage and synthesize more spatial datasets (Hampton et al., 2013), timeseries data and analysis are more than ever needed to reach a better understanding and predictability of non-equilibrium biodiversity responses to climate change.
For plants, the most common type of dynamics reported were constant-relationship scenarios. Such dynamics were impaired with strong lagged responses. The current distribution of North American plants are heavily affected by climatic conditions at the Last Glacial Maximum (Ordonez, 2013), and plants are known to show dispersal lags in post-glaciated areas (Normand et al., 2011). This explains the persistence of deviation in responses between observed and inferred climate in North American plant communities across the last 21 Ka. Constant-relationship scenarios responses might be predictable by modeling approaches considering dispersal and community assembly rules (Guisan and Rahbek, 2011). However, almost a quarter of the plant communities exhibited unpredictable dynamics such as constant-lag and alternate stable states.
Bird communities mainly exhibited unpredictable dynamics. Community dynamics appeared stochastic, and often uncorrelated with the observed climate. Despite the fact that maximum temperature change was generally not strong or directional (see Figure S5), the complex disequilibrium observed in bird communities compared to the plant communities is in line with our expectations. Two reasons can explain the stochasticity observed in bird responses to modern climate change. Firstly, climatic determinism of northern hemisphere bird communities is questionable (Beale et al., 2008;Currie and Venne, 2017). For example, broad thermal tolerances (Araújo et al., 2013) and phenological responses (Stenseth et al., 2002;Dunn and Winkler, 2010) might buffer the impact of moderate temperature changes on communities dynamics (Gaüzère et al., 2015). Moreover, bird sensitivity to habitat changes probably influences community-inferred temperatures (Clavero et al., 2011;Barnagaud et al., 2013) and overrides the direct effect of climate warming on community dynamics (Eglington and Pearce-Higgins, 2012;Gaget et al., 2018).
Secondly, short-term variation might be harder to predict than long-term changes because stochastic variability is predominant at fine spatial and temporal scales (Levin, 1992). However, consistent short-term directional changes in bird communityinferred temperature have been reported at continental (Devictor et al., 2012;Princé and Zuckerberg, 2015), national (Devictor et al., 2008), and landscape scales (Gaüzère et al., 2015). Our results showed that local-scale changes in communityinferred temperature were not consistently related to observed temperature change. The scale of effect is defined as the scale at which an environmental attribute has the strongest effect on inferred species-environment relationship. While it is known as a strong determinant of explanatory predictions (Holland et al., 2004;de Knegt et al., 2010), many empirical studies might not be conducted at the appropriate spatial scales (Jackson and Fahrig, 2015). Hence, we can hypothesize that the low predictability exhibited by breeding bird communities might be due to the weak climatic determinism of bird community dynamics at local scale. This suggest that using equilibrium dynamics hypothesis as a conceptual model to predict biodiversity responses to climate change requires caution. We argue that a careful assessment of climate determinism focused on the taxon and the scale of study is a prerequisite for successful anticipatory predictions.
We also showed that some aspects of predictability-absolute deviation from climate equilibrium and deviation changewere structured by environmental or ecological factors, while others-number of alternate states-were not. We expected community predictability to decrease with thermal severity. Our results showed that absolute deviation of plant communities was decreasing with temperature, with a curvilinear relationship showing a plateau on warmest values. Conversely, absolute deviation for birds was increasing with increasing temperature before reaching a plateau. This apparent discrepancy between plants and birds is linked to the distribution of Neotoma and BBS sites. Neotoma sites are distributed in northwestern Nearctic margin and are therefore colder than BBS sites distributed in southern Nearctic margin. Merging the BBS and Neotoma estimates showed a quadratic relationship between absolute deviation and baseline temperature ( Figure S6). As expected, overall absolute deviation increased with coldest and warmest temperature. The consistent effect of baseline climate between taxa and spatial scale suggest a strong regional-scale determinism of predictability, structured by the diversity of realized thermal niches in the regional species pool. In France, Bertrand et al. (2016) already reported a strong effect of baseline temperature on deviation from equilibrium state, in link with the absence of climate-adapted species in the regional pool. Climate severity is expected to have an even stronger effect in North America. The distribution of land masses constrains Neartic species' distribution range in their northern and southern boundaries. These geometric constraints on species distribution, also called "mid-domain effect, " are known to shape latitudinal richness gradient (Colwell and Lees, 2000). While its application to non-spatial domains is scarce (but see Letten et al., 2013), the "thermal mid-domain effect" probably have a strong influence on the species' thermal tolerance present in regional species pools (Brayard et al., 2005;Beaugrand et al., 2013). This suggest that long-term biogeographic history and macro-scale processes have a strong influence on community predictability. Further investigations of the thermal mid-domain effect and its consequence on regional pools should clarify its implication in the predictability of community responses to climate change.
Our analysis showed that plant communities composed of taller plants exhibited lower absolute deviation from climate equilibrium. This result is in line with our predictions. Nolag dynamics and predictable responses are thought to occur when species exhibit low persistence through rapid extinction at trailing range edges (Hampe and Petit, 2005), and/or efficient niche tracking through long-distance dispersal. Conversely, disequilibrium responses are thought to occur when species' responses in these domains are opposite (Svenning and Sandel, 2013). In turn, the importance of these processes is linked to species' dispersal ability and life-history traits. For example, species traits related to weak dispersal ability decrease species' niche tracking (Svenning and Skov, 2007) while persistence processes are enhanced by survival of long-lived individuals (Eriksson, 1996;Holt, 2009;Jackson and Sax, 2010). However, we did not find support for the effect of seed mass on the predictability. While lower seed mass is generally considered as a proxy for longer dispersal distance, empirical evidence is mixed (Thomson et al., 2011). Plant height might even be better predictor for seed dispersal distance (Muller-Landau et al., 2008;Thomson et al., 2011). Because dispersal limitation is expected to be a major driver of climate disequilibrium for plants (Svenning and Skov, 2007, but see Bertrand et al., 2016), the improved dispersal of taller plants supported our result.
For birds we showed that, as predicted, communities composed of larger birds exhibited stronger absolute deviation. Despite the fact that body mass is an integrative species characteristic correlated with many life-history traits, this result confirms that life history traits can influence birds responses through niche tracking and persistence processes (Jiguet et al., 2007). Note that side analyses incorporated the percent of insectivorous birds species as a potential determinant of predictability ( Figure S7). This predictor was removed from our final model to keep consistency between plant and bird analyses.
Nevertheless, the influence of species characteristics was generally weak, and not consistent across taxa and CRD statistics. While trait might be good predictors for population responses to climate change (e.g., Julliard et al., 2003;Jiguet et al., 2007), there is weak support for their effect on species' distribution range shifts (Angert et al., 2011;Tingley et al., 2012;Smith et al., 2013). Different reasons such as the stochastic nature of colonization events, novel species interactions and extrinsic effects of habitat availability and fragmentation might explain these weak effects. Moreover, the properties of species assemblages and assembly rules might be more important for community scale dynamics (Guisan and Rahbek, 2011).
Our set of predictors failed to explain variation in maximum state number. This statistic is a key component of predictability. However, accounting for sampling error challenges a straightforward interpretation of n-values when applied to stochastic dynamics. Without sampling error, stochastic dynamics are expected to cause high n values associated with low predictability. For birds, the uncorrected n values were, indeed, high (99% of the sites having n > 1 and 88% of the sites having n > 2), but the necessary correction starkly reduced this estimate (14% of the sites having n > 1 after correcting for sampling uncertainty).

CONCLUSION
A better understanding of the limits to predictability is a crucial step for predictive modeling and applied ecology (Mouquet et al., 2015). Our study showed that the equilibrium dynamic hypothesis to infer community responses to climate change is only sometimes applicable. In many cases, a straightforward application of the equilibrium dynamic hypothesis to predict biodiversity responses to future climate changes may lead to misleading predictions. Equilibrium dynamics across different taxa and scales should be assumed cautiously. We argue that robust anticipatory predictions will require detailed knowledge of the taxa considered, along with the spatial and temporal scales at which key processes are expected to drive biodiversity response to climate change.

AUTHOR CONTRIBUTIONS
PG and BB designed the study and performed the analysis. PG wrote the first draft of the manuscript. All authors contributed substantially to the conceptual aspect of the study, and to the writing and revisions of the manuscript.