Plot and environmental data<\/p>\n
We retrieved 1,903,668 vegetation plots and 40,468,042 observation records from the European Vegetation Archive (EVA; Project 189; accessed 24 August 2023; ref. 24<\/a>), including plots from the ReSurveyEurope initiative53<\/a> and associated LOTVS, GLORIA and forestREplot databases54<\/a>,55<\/a>,56<\/a>; this included single-sampled and resampled plots, which were not treated differently. We selected for plots within Europe, but excluding Ukraine, Russia, Belarus and Turkey. Plots were further filtered to exclude those with location uncertainty exceeding 5\u2009km, noted as experimentally manipulated, lacking date of recording or sampled outside our study period (1981\u20132018) (Supplementary Section 1<\/a>). The lower bound of 1981 ensured that all plots had at least 20\u2009years of preceding climate data (from 1961) for lag analyses. After filtering, 1,138,587 plots remained (Supplementary Fig. 6<\/a>); 91,260 did not contain coordinate uncertainty information but were retained.<\/p>\n A total of 22,083,578 observation records of vascular plants were obtained from these plots. We removed observations with taxonomic identification above the species level (including section, series and subgenus), hybrids, uncertain identities and cultivated individuals. To standardize species identities for modelling, we consolidated infraspecific taxa and taxonomic qualifiers under their corresponding species names, grouping subspecies, varieties, forma and other infraspecific entities that are often not distinguished or misidentified57<\/a>. Although this approach simplifies the data and facilitates modelling, we acknowledge that it may not fully capture the ecological variability within some species complexes. A total of 21,639,385 observations remained. All subsequent analyses were performed at the scale of individual vegetation plots.<\/p>\n All climate data were sourced from the UERRA regional reanalysis for Europe dataset (retrieved 18 October 2023; Copernicus Climate Change Service 2019 (ref. 58<\/a>)). Daily near-surface (2\u2009m above ground) temperature and precipitation records from 1961 to 2018 were obtained from the MESCAN-SURFEX system (0.5\u00b0 by 0.5\u00b0), representing macroclimatic conditions. Although older palaeoclimatic datasets exist, major anthropogenic climate shifts began around the 1960s, making data before this period less relevant for capturing biotic response to recent climate change (last century). Accordingly, we limited our analysis to the post-1960 period to enhance detectability of biologically meaningful lags in response to modern climate trends. Daily climate data were aggregated into 21 annual climate variables for each year: 15 bioclimatic variables (excluding precipitation variables based on temperature-defined periods and vice versa), three growing degree day variables (0\u2009\u00b0C, 5\u2009\u00b0C and 10\u2009\u00b0C thresholds), tree growth season length, snow cover days and frost change frequency (Supplementary Section 1<\/a> and Supplementary Table 2<\/a>). This resulted in annual layers for all 21 variables across 1961\u20132018, providing temporally explicit climate data for dynamic modelling. After excluding highly correlated climate variables (r\u2009\u2265\u20090.7), seven predictors were retained: annual mean temperature, mean diurnal range, isothermality, temperature seasonality, annual precipitation, precipitation seasonality and frost change frequency. These variables were used in both static and dynamic SDMs.<\/p>\n Light indicator values were sourced from a harmonized dataset of Ellenberg-type values at the European scale26<\/a>, using species aggregate values when available (that is, values averaged across closely related or difficult-to-distinguish taxa within a species complex). These values reflect the typical light conditions, or light availability, under which species are found, ranging from deep shade to full sunlight (1\u20139). Disturbance indicator values were sourced from a European database25<\/a>, including disturbance severity, herb-layer disturbance severity, disturbance frequency and herb-layer disturbance frequency. For our final species list, herb-layer disturbance severity was highly correlated with both disturbance frequency (r\u2009=\u20090.70) and herb-layer disturbance frequency (r\u2009=\u20090.76) and thus, frequency indicators were excluded. Disturbance severity represents the expected proportion of biomass removed or destroyed at the whole-community level during a disturbance event; herb-layer disturbance severity reflects the same process but restricted to the herbaceous field layer, and is therefore always equal to or lower than disturbance severity25<\/a>. These disturbance indicators do not record disturbance events or biomass loss directly; instead, they reflect the disturbance conditions characteristic of the habitats in which a species typically occurs, and therefore describe ecological associations of species with disturbance rather than disturbance per se. In total, three disturbance-related predictors were used: light availability, disturbance severity and herb-layer disturbance severity.<\/p>\n Following ref. 26<\/a>, a minimum of five species per plot was required for reliable community-averaged indicator value calculation. To avoid circularity\u2014where presence of a species would influence its own disturbance predictor\u2014we excluded the focal species when calculating plot-level indicators for that species (that is, plot-level indicators were calculated separately for each species)17<\/a>. Thus, plots were required to contain at least six species, ensuring that five remained even after focal species exclusion. This approach kept the set of eligible plots consistent across all species. Applying this criterion excluded 136,580 plots.<\/p>\n To calculate plot-level community-averaged indicator values, we used the unweighted average across species within a plot. When testing weights, we found a high correlation (r\u2009=\u20090.95) between unweighted and weighted values (based on percentage cover of the plot by species), supporting the use of the simpler unweighted approach. We focused on light availability and disturbance regimes because of their strong, known links to disturbance processes25<\/a>.<\/p>\n To assess whether the disturbance and light indicators reflected independently observable forest structural conditions, we compared a subset of plots against remote sensing-derived European tree extent and tree height products59<\/a>. Herb-layer disturbance severity and light availability showed stronger relationships with local-scale remotely sensed forest structure than disturbance severity, with relationships generally weakening at larger spatial buffers (Supplementary Table 3<\/a>).<\/p>\n Species distribution modelling<\/p>\n We developed SDMs to investigate climate-change and disturbance disequilibrium across European plants (Supplementary Section 2<\/a>). Given the focus on detecting disequilibrium and lagged responses, we designed the modelling framework to ensure temporal independence between model training and evaluation. Plots were separated into two temporal bins: a training bin (1981\u20132003) for model training and a testing bin (2004\u20132018) for climate-window selection and model validation. This separation served three purposes. First, it ensured temporal independence between model calibration and evaluation and allowed genuine forward projection across a period of known climate change. Second, it reduced temporal heterogeneity within each bin, limiting noise from short-term variability or uneven sampling that could obscure longer-term trends. Third, it maintained sufficient sample size and spatial coverage in both periods while providing a conservative test for temporal disequilibrium\u2014models trained on earlier data and evaluated on later data perform worse when species\u2013environment relationships have shifted. This two-bin structure therefore provides a robust yet data-efficient framework for quantifying disequilibrium under recent environmental change.<\/p>\n Our modelling framework explicitly tests the assumption of climatic equilibrium inherent to conventional SDMs. Static climate models represent the equilibrium expectation, assuming that occurrences of species reflect current climatic suitability averaged over multidecadal periods. By contrast, dynamic models allow quasi-equilibrium states to shift between temporal bins by incorporating time-varying climate predictors. The difference in predictive performance between these two model types quantifies the degree to which observed occurrences of species deviate from static climatic expectations under ongoing climate change. In this context, occurrence data are interpreted not as perfect reflections of present-day climate suitability but as outcomes of past or lagged climatic conditions, enabling our framework to detect both the strength and temporal depth of climate-change disequilibrium.<\/p>\n Species presence data were spatially thinned with a 10-km buffer to reduce spatial autocorrelation and the over-representation of densely sampled regions; that is, presences occurring within 10\u2009km of another were removed, prioritizing more recent records (Supplementary Section 2<\/a>). Thinning was applied separately for each temporal bin to ensure comparable spatial distributions between training and testing datasets. This procedure substantially reduces clustering bias in the underlying vegetation plot data, improving the spatial independence of observations and preventing inflated model performance due to pseudoreplication. By thinning within bins rather than across them, we also minimize confounding spatial structure with temporal turnover, ensuring that modelled changes in species\u2013environment relationships reflect true temporal dynamics rather than uneven spatial sampling. Species with \u226535 training and \u226510 testing presences were retained, yielding 3,047 species for modelling. Absence data comprised plots without the target species and within 300\u2009km of its presence plots to better represent true absences due to climatic- or disturbance-related factors rather than dispersal limitations. This accessible-area constraint helps minimize the inclusion of absences that were probably dispersal-limited\u2014particularly those arising from environmental changes that predate our study period\u2014while retaining broad climatic and biogeographic representativeness across Europe. Absences were binned and thinned identically to presence plots to maintain consistent spatial independence and sampling density across both temporal periods.<\/p>\n Four SDM types were fitted to the presence\u2013absence plot data of each species, differentiated by their predictor set. (1) Standard SDMs used static climate variables averaged over the temporal bins (1981\u20132003 for training and 2004\u20132018 for projecting and validation). Standard SDMs served as baselines and were deliberately kept simple to avoid embedding disturbance or temporal dynamics within the reference model, ensuring that any improvement in performance directly reflected the added ecological processes rather than differences in model complexity. (2) Dynamic SDMs used 5-year moving-window climate averages, calculated over ten lag intervals spanning 0\u201318\u2009years (in 2-year increments) before the date of recording or year of sampling for a plot; 5-year climate windows spanned \u00b12\u2009years from the selected lag year. For zero-lag windows, climate data after the sampling year were ignored, resulting in a 3-year average. Because sampling years varied across plots, the climate years included in a given lag window also varied between plots. As a result, lag effects were captured at the species or distributional scale (reflecting overall shifts in climatic associations of species) rather than at the level of local population- or plot-specific responses. Only one lag period was considered at a time because testing different lag periods among climate variables would result in an unmanageable number of models (from 10 to 107 dynamic models per species). The best-performing window was selected by evaluating resulting model projections on a subset of presence\u2013absence data from the testing bin (see below); models with different windows were assessed and compared using AUC. (3) Disturbance-informed SDMs used static climate variables (as done for standard SDMs) and disturbance-related predictors (as calculated above, unweighted plot-level community-averaged indicator values of light availability, disturbance severity and herb-layer disturbance severity). Target species were excluded from indicator calculations for presence plots to prevent circularity, so calculations were performed separately for each species17<\/a>. (4) Dynamic, disturbance-informed SDMs combined both lagged climatic predictors and disturbance indicators.<\/p>\n Models were trained using a randomly selected 90% subset of presence\u2013absence data from the training bin (1981\u20132003). We used GAMs, fitted using the mgcv R package (v.1.9.0) with a binomial family, logit link function and smoothing parameter estimation via restricted maximum likelihood; other settings were left as defaults. The basis dimension (k) for smooth terms was limited to four to constrain model flexibility, improving biological plausibility by preventing overly complex response shapes and reducing overfitting to noise by limiting the number of allowable inflection points independent of smoothing parameter optimization60<\/a>. To keep predictor numbers constant across modelling protocols, disturbance-uninformed SDMs used the same disturbance-related predictors but with values that were randomly permutated (for both training and testing data bins).<\/p>\n Testing data (2004\u20132018) were split into two equal subsets: one used to select the best-performing climate window for SDMs using dynamic climate variables; and the other reserved for independent model validation. This separation avoided information leakage between calibration (window selection) and validation steps, ensuring that model evaluation remained unbiased. Using the same data for both window selection and model validation would have artificially inflated model performance and exaggerated the apparent significance of climate-change disequilibrium by overfitting\u2014or, more specifically, over-calibrating\u2014the models to the specific windowed climate conditions. Model performance was assessed using AUC, a popular discriminatory metric that evaluates probabilistic accuracy (or continuous habitat suitability predictions). Model performance was assessed using plots reserved for independent model validation; all downstream analyses relied on AUC scores as assessed based on this independent data. For analyses requiring binary classification of suitability predictions (for example, to assess spatial patterns of overpredictions and underpredictions across plots), we applied the tenth-percentile minimum training presence threshold.<\/p>\n To account for model-inherent errors (false presences or absences arising from statistical noise or sampling limitations), we used a paired-model framework in which each disequilibrium-informed SDM was directly compared with a corresponding non-disequilibrium-informed SDM built using identical data and parameters; only the relevant predictors were different. Because both models share the same structural assumptions and data limitations, model-inherent errors are effectively held constant, allowing changes in performance (\u0394AUC) to reflect ecological mismatches rather than artefacts of model fit. Note, AUC values were logit-transformed before calculating deltas (\u0394logit-AUC), which were used as measures of disequilibrium in downstream analyses.<\/p>\n To quantify the severity of distributional disequilibrium at the species level, we ran ten replicates for each modelling protocol and performed species-specific pairwise t-tests comparing the performance of standard SDMs against dynamic SDMs, disturbance-informed SDMs and dynamic, disturbance-informed SDMs (\u0394logit-AUC). We applied the Benjamini\u2013Hochberg correction to control for false discovery rate for each set of comparisons (n\u2009=\u20093,047) (for example, standard SDMs against dynamic SDMs). Species were classified as exhibiting significant climate change- or disturbance-driven disequilibrium when the relevant model significantly outperformed the baseline.<\/p>\n Between-protocol differences in model performance were tested using an LMM with logit-AUC as the response, modelling protocol as the fixed effect and species as the random intercept (ten replicates per species). Estimated marginal means were then used to obtain pairwise comparisons among modelling protocols with Tukey adjustment for multiple testing.<\/p>\n Relationship between severity and attributes of species<\/p>\n To explore species-level correlates of disequilibrium severity, we modelled \u0394logit-AUC values using LMMs with species attributes as predictors (Supplementary Section 3<\/a>) and species as the random effect. Twenty attributes were selected after removing highly correlated variables (r\u2009\u2265\u20090.7) (Supplementary Table 4<\/a>). These included: three data characteristics (total number of presence plots and median sampling year of training and testing presence plots), two geographic variables (latitude and elevation), three bioclimatic variables (annual mean temperature, mean diurnal range and isothermality), three disturbance-related predictors (community-averaged indicator values of light availability, disturbance severity and herb-layer disturbance severity) and nine plant functional traits (plant height, seed mass, seed number per reproduction unit, stem specific density, leaf dry matter content, specific leaf area, leaf area, leaf nitrogen and leaf nitrogen to phosphorus ratio). We also included a binary categorical predictor indicating if the other disequilibrium type was accounted for in the model; this tests if accounting for the other disequilibrium influenced the detectability of the focal disequilibrium type.<\/p>\n Data characteristics were included to assess potential data-related and temporal influences on disequilibrium detection. The number of presence plots, which was strongly correlated with estimated range size (r\u2009=\u20090.86), tested for possible data-dependence or range-size effects on disequilibrium detection. Median sampling years (training and testing) captured potential temporal biases in the plot data. These variables also allowed evaluation of whether species with more recent observations exhibited greater disequilibrium, which would reflect accelerating impacts of contemporary climate change.<\/p>\n Except for data characteristics and functional traits, predictor values were calculated as the median across presence plots for each species. Elevation data were extracted from EuroDEM 2023 (2\u2009arcsec; retrieved 5 March 2024; www.mapsforeurope.org\/datasets\/euro-dem<\/a>). Community-averaged indicator values were used over original species-level indicators because of broader data coverage (>250 species lacked original values). Where available, original and community-averaged values were highly correlated (light, r\u2009=\u20090.85; disturbance severity, r\u2009=\u20090.94; herb-layer disturbance severity, r\u2009=\u20090.96), supporting their use. Traits data were obtained from a gap-filled dataset compiled by the TRY initiative61<\/a>,62<\/a>. All predictors were standardized before analysis.<\/p>\n Quadratic terms were included for latitude, light availability, disturbance severity and herb-layer disturbance severity to test for potential unimodal relationships. Although we expected greater disequilibrium among species associated with open habitats, Bond\u2019s concept of ecosystem uncertainty16<\/a> implies that disequilibrium should occur at both ends of the habitat spectrum\u2014open and closed vegetation alike. Therefore, we included quadratic terms to allow for nonlinear, potentially symmetric responses across these gradients. These squared terms were uncorrelated with the other 20 predictors and each other, resulting in a final predictor set of 24 attributes.<\/p>\n Although we prefiltered highly correlated predictors, variance inflation factors (VIFs) remained elevated (Supplementary Section 3<\/a> and Supplementary Table 5<\/a>). We therefore removed annual mean temperature (VIF\u2009=\u200931.1). Disturbance severity (VIF\u2009=\u20099.6) and light availability (VIF\u2009=\u20098.6) exceeded the strict threshold of 5, but we retained both because their VIFs remained below the less strict threshold of 10 and excluding them would hinder the interpretability of their quadratic terms. For completeness, we also explored alternative models with both of these variables removed (Supplementary Fig. 4<\/a>). We also tested longitude to contextualize bioclimatic effects (for example, continentality). Longitude was moderately related with isothermality (VIF\u2009=\u200913.8), but excluding isothermality reduced VIFs to acceptable levels (longitude\u2009=\u20091.7) (Supplementary Table 5<\/a>).<\/p>\n Although we used the same climate and disturbance covariates in our SDMs, we first collapsed them into independent species-level summaries (for example, the median value across all occurrence plots) and then used those summaries as predictors in LMMs for disequilibrium severity (\u0394logit-AUC)\u2014a response wholly separate from occurrence probability or habitat suitability\u2014thus avoiding circularity.<\/p>\n Note, we chose not to rely solely on resampled data to analyse climate-change disequilibrium because of the unknown and highly variable spatial bias in the distribution of resurveyed plots (a smaller subset of the EVA plots used); both single-sampled and resampled plot data were used and treated similarly. Hence, we were unable to determine the type of climate-change disequilibrium occurring, that is, between extinction debts and colonization credits.<\/p>\n Spatial distribution of disequilibrium<\/p>\n To identify where disequilibrium occurred along environmental gradients, we directly compared predictions of the test data between paired models\u2014each dynamic and\/or disturbance-informed SDM against its static climate and\/or climate-only baseline. Disequilibrium plots were defined as cases where the baseline model misclassified presence or absence, but the disequilibrium-informed model predicted correctly. This approach isolates locations where accounting for disturbance or temporal climate dynamics resolves mismatches that static, climate-only assumptions cannot, thereby linking model improvement to ecological rather than statistical causes. For each species and selected environmental variables, we calculated median values across all presence plots and all disequilibrium plots. From the previously defined set of species attributes, we assessed patterns for the three disturbance-related predictors, four bioclimatic variables and three geographic variables. Cross-species patterns were visualized by plotting species-level medians (presence versus disequilibrium plots) and fitting GAMs to capture nonlinear trends. Analyses were restricted to species that exhibited significant disequilibrium (as assessed by the pairwise t-test and after adjusting for false discovery rate using the Benjamini\u2013Hochberg) and at least five disequilibrium plots; because climate-change and disturbance-driven disequilibrium varied depending on whether the other was accounted for, significance tests and plot predictions were taken from models that included the other effect when it was significant. For overpredictions, this accounted for 3,015 species for disturbance-driven disequilibrium and 1,578 species for climate-change disequilibrium. For underpredictions, this accounted for 1,711 species for disturbance-driven disequilibrium and 590 species for climate-change disequilibrium. Variables shown in the main text were selected on the basis of statistical significance (severity\u2013attribute relationship), ecological relevance or the emergence of distinct spatial patterns. Although both overpredictions and underpredictions were explored, we focused on the most relevant type for each disequilibrium; full results are presented in Extended Data Figs. 3<\/a> and 4<\/a>.<\/p>\n Analyses of lagged responses<\/p>\n We randomly selected 12 species with significant climate-change disequilibrium (n\u2009=\u20091,434) to illustrate variations in climate-window performances. Two species were selected for each combination growth form and range size. Plant growth forms were obtained from the TRY database and range size was estimated using an alpha-hull approach (Supplementary Section 3<\/a>). Narrow-ranged species were defined as those below the 25th percentile (~175,782\u2009km2) and widespread species as those above the 75th percentile (~1,336,012\u2009km2) of alpha-hull range sizes, based on the distribution across species with significant climate-change disequilibrium. To confirm that observed patterns were not an artefact of this stratified selection, we also examined a fully random set of 24 species (Extended Data Fig. 5<\/a>).<\/p>\n Best-performing climate windows were then grouped into three categories: recent (0\u20138\u2009years before sampling), older (10\u201318\u2009years) and null (static climate average over training and testing periods). For each species, the lag interval corresponding to its best-performing climate window was extracted and modelled as a response variable against species attributes, using the same modelling framework and considerations as for disequilibrium severity.<\/p>\n Reporting summary<\/p>\n