Abstract

Climate warming is shifting biological communities, with warmth-demanding species being favoured at the expense of cold-adapted species in a process referred to as thermophilization 1 , 2 , 3 , 4 . Because biodiversity responses often lag behind climate warming, climatic debts are accumulating in many ecosystems across the world 5 , 6 , 7 . Although we might expect that thermophilization and climatic debts will vary among habitats, standardized quantification across ecosystems is lacking. Here we analysed multidecadal data from 6,067 resurveyed vegetation plots over 12–78 years in forests, grasslands and on alpine summits across Europe. We demonstrate that forest understory and grassland plant communities experienced positive thermophilization, although not significantly different from zero. By contrast, alpine summit vegetation showed much stronger (up to five times) and significant thermophilization. Thermophilization was driven largely by increases in warmth-demanding species in grasslands, by declines in cold-adapted species on alpine summits and by both processes in forests. Significant climatic debts have accumulated in forests and alpine summits, but less so in grasslands, with debts positively correlated with macroclimate temperature changes. Our findings uncover divergent thermophilization trajectories and increasing climatic debts across ecosystems. Moreover, we highlight the mechanisms that enable some communities to track climate change more closely than others and provide a basis for projecting future shifts in plant communities under accelerating climate warming.

Main

Global climate change is having substantial impacts on biodiversity 8 , 9 , 10 . Range shifts are being observed in many species and ecosystem types as species track their thermal requirements 11 , 12 , 13 , 14 . As a result, directional shifts in community composition at the expense of cold-adapted species and in favour of warmth-demanding species (referred to as thermophilization) are increasingly evident in terrestrial and aquatic plant and animal communities 1 , 2 , 3 , 4 , 15 , 16 . For example, significant thermophilization of vegetation has been documented across alpine summits and grasslands in Europe 3 , 17 , and in plant communities of temperate and tropical forests despite the buffering effects of forest canopies 1 , 18 , 19 , 20 . However, in many biological communities, thermophilization proceeds more slowly than contemporary warming 13 , 18 , 21 , 22 , which results in a climatic debt in community composition 2 , 5 , 6 , 15 , 23 .

Climatic debts reflect the inherent life-history and ecophysiological characteristics of species, such as limited dispersal and establishment capacity and long lifespan 7 , 24 , 25 . However, they may also arise from the buffering effects of microclimatic processes or interactions with other global change drivers 1 , 18 , 23 , 26 . Apparent debts may also be shaped by adaptive shifts in climatic niches of species, although such processes are generally slow 7 . Microclimates in non-forest ecosystems, such as alpine summits and grasslands, are mainly shaped by topography and vegetation feedbacks on snow cover, solar radiation, air mixing and evapotranspiration 27 , 28 , 29 , 30 . By contrast, in forests, the tree canopy structure largely governs microclimatic conditions and thermal buffering 1 , 18 , although topography is also important 31 . As such, the magnitude of thermophilization might appreciably differ among ecosystems. However, comparing studies of thermophilization rates is difficult, as they often use disparate methodologies to quantify thermal niches and climatic debts (for example, as a result of different sources, resolutions, seasons and base periods of climate data). Thermophilization rates in response to climate change may differ across ecosystems 4 . However, we lack consistent large-scale, multitaxa and long-term comparisons of thermophilization rates and climatic debt magnitudes among vegetation types based on a consistent methodology. Here to address this knowledge gap, we use vegetation resurveys with a consistent methodology to answer two key questions: (1) how thermophilization rates and climatic debts differ among different ecosystem types; and (2) how these differences are related to variation in local climate warming across contrasting ecosystems.

We report cross-biome thermophilization for Europe’s temperate forests, grasslands and alpine summits. These ecosystems represent widespread and ecologically distinct vegetation types with contrasting structures and microclimates relevant to responses of species to climate change. We compiled plant occurrence data for 1,919 species from 6,067 resurveyed vegetation plots (Fig. 1). In each plot, we quantified the thermophilization rate of the plant community 1 , 18 using a consistent methodology (Methods). In total, 4,372, 1,209 and 486 plots from forests, grasslands and alpine summits, respectively, were analysed, with sampling intervals (baseline to resurveys) of 12–78 years (median of 42 years), 16–78 (median of 35 years) and 12–15 (median of 14 years), respectively (Extended Data Table 1 and Extended Data Fig. 1a). We used these data to estimate rates of thermophilization and climatic debt (°C per decade) for each plot using an established method that relies on the thermal affinities of species inferred from their geographical distributions relative to long-term climate data 1 , 18 , 32 . To further disentangle the mechanisms that underlie thermophilization, we separately calculated changes in the 5th and 95th percentiles of the floristic temperature distribution to represent declines in cold-adapted species compared with increases in warmth-demanding species, respectively 1 .

Forest cover (green) and grassland cover (yellow) are based on CORINE Land Cover 2018 data 50 . Elevation data (grey) are from the NASA Jet Propulsion Laboratory 51 .

First, we estimated thermophilization rates using plant abundance data (excluding rare species with range sizes fewer than 1,000 grid cells at 20 × 20 km resolution; Methods), which mainly reflect changes in species dominance and community reordering. We found positive, although not significantly different from zero, rates of overall thermophilization in forests (0.03 °C per decade, 95% credible interval (CI) of −0.02–0.07) and grasslands (0.03 °C per decade, 95% CI of −0.06–0.11). By contrast, the overall thermophilization rate on alpine summits was much higher and significantly different from zero at 0.12 °C per decade (95% CI of 0.05–0.19) (Fig. 2a). When estimated using plant presence/absence data, which are more sensitive to species turnover, including rare or transient species, similar results were found. The thermophilization rate on alpine summits was five times higher than in forests and grasslands (see the section ‘Sensitivity analysis’ in the Methods and Extended Data Fig. 2a). When assessed using abundance data and excluding rare species, thermophilization in grasslands was driven largely by an increase in warmth-demanding species (0.08 °C per decade, 95% CI of −0.03–0.19). In forests, thermophilization reflected both declines in cold-adapted species (0.04, 95% CI of −0.01–0.09 °C per decade) and increases in warmth-demanding species (0.03, 95% CI of −0.03–0.09 °C per decade). On alpine summits, however, we found a significant decrease in cold-adapted species (0.26 °C per decade, 95% CI of 0.18–0.34) and a decrease in warmth-demanding species (−0.16 °C per decade, 95% CI of −0.25 to −0.06). These trends did not show marked changes when rare species and/or when vegetation presence/absence data were included in the assessment (see the section ‘Sensitivity analysis’ in the Methods and Extended Data Fig. 2b).

a, Vegetation thermophilization rates quantified as declines in cold-adapted species (5th percentile shift), overall (median shift) and increases in warmth-demanding species (95th percentile shift). b, Climatic debts (°C per decade) calculated according to climate change rates and vegetation thermophilization rates of declines in cold-adapted species, overall and increases in warmth-demanding species. Thermophilization rates and climatic debts were estimated using plant abundance data that excluded rare species. Circles represent means (n = 4,372 plots, 1,209 plots and 457 plots for forests, grasslands and alpine summits, respectively) with 80% (thick line) and 95% (thin line) CIs and posterior distributions obtained from Bayesian mixed-effects models.

Our second analysis involved calculations of climatic debts by subtracting thermophilization rates from local rates of climate change 2 , 24 . These climatic debts were mostly positive (ranging from 0.13 to 0.34 °C per decade across ecosystem types; Fig. 2b and Extended Data Fig. 2), which indicated that thermophilization lags behind macroclimate changes. The magnitude of these lags did not strongly depend on whether they were estimated from declines in cold-adapted species, increases in warmth-demanding species or overall thermophilization rates. Climatic debts were generally highest in forests, followed by alpine summits and were lowest in grasslands. Thermophilization rates were weakly correlated with rates of increase in temperature regardless of the ecosystem type. Conversely, climatic debts strongly covaried with rates of temperature change across all ecosystem types (Fig. 3).

a–c, Relationships between the growing season mean temperature (T mean ) change rate (°C per decade) and the thermophilization rate (°C per decade) quantified as declines in cold-adapted species (5th percentile shift) (a), overall (median shift) (b) and increases in warmth-demanding species (95th percentile shift) (c). d–f, Relationships between the growing season T mean change rate (°C per decade) and climatic debt (°C per decade) calculated according to the thermophilization rate quantified as declines in cold-adapted species (d), overall (e) and increases in warmth-demanding species (f). Thermophilization rates and climatic debts were estimated using plant abundance data that excluded rare species. Each circle represents a plot (n = 4,372, 1,209 and 457 for forests, grasslands and alpine summits, respectively). Lines are fits from Bayesian mixed-effect models, and shaded areas denote 95% CIs.

The overall rates of thermophilization on alpine summits were generally higher (up to five times) than those in forests and grasslands. The divergence in thermophilization rates among ecosystem types may be partially attributed to differences in plant lifespans, with communities composed of shorter-lived plants potentially responding more rapidly to climate warming 33 , 34 . However, currently available lifespan datasets (for example, from the TRY Plant Trait Database 35 ) are still limited in taxonomic coverage and resolution. This limitation is especially relevant for ecosystems such as alpine summits and forests, where many plant species are perennials and have multidecadal to centuries-long lifespans, which prevented us from quantitatively testing this hypothesis. The high thermophilization observed on alpine summits may predominantly result from topographical compression of thermal belts, which creates strong spatial gradients in microclimate and species composition such that even small shifts in species distribution can lead to marked community-level changes in thermal affinity 36 . In the alpine summit dataset, however, we did not detect significant relationships between thermophilization rates and variables that reflect rugged topography, such as the terrain ruggedness index (TRI) or slope (Extended Data Fig. 3a). This result is probably due to the fact that summit-level coordinates in our database and mid-scale topographical metrics did not capture the fine-scale microhabitats that more strongly shape community responses, especially in alpine environments. Nevertheless, thermophilization was lowest on north-facing slopes (Extended Data Fig. 3b). This finding supports the hypothesis that microhabitats are generated by topographical gradients and emphasize the role of slope exposition to sun radiation in shaping thermophilization patterns of treeless alpine vegetation.

Our results further suggest that thermophilization of the understory plants in forests is driven by both declines in cold-adapted species and increases in warmth-demanding species, which is expected with warming. By contrast, only increases in warmth-demanding species contributed to overall thermophilization in grasslands. Unlike forests, grasslands lack tree canopy buffering, and grassland plant communities experience hotter maximum daily temperatures than forest understoreys. Together with higher species turnover and recruitment rates, this factor could have favoured warmth-demanding taxa more as temperatures rise 37 , which results in thermophilization being driven mainly by increases in warmth-demanding species. For alpine summits, however, our results suggest a simultaneous decline in both cold-adapted and warmth-demanding species, which could potentially be linked to recurring severe drought during the survey periods 38 , 39 , although it may reflect additionally distinct ecological pressures. Cold-adapted species may decline or be lost because of substantial climate warming and habitat contraction 40 , 41 . Conversely, the decline in warmth-demanding species may result from reduced persistence after colonization, which may be due to environmental stress or biotic interactions 42 . Despite declines in both types of species, net thermophilization can still result from stronger declines in cold-adapted species than warmth-demanding species (Fig. 2 and Extended Data Fig. 2) and from compositional shifts that favour warmth-demanding species through both topography-facilitated colonization and differential persistence 43 . By partitioning thermophilization into the effects of gains, losses and abundance shifts of persisting species, we found that abundance shifts, along with losses, were the primary drivers of thermophilization on alpine summits (see the section ‘Sensitivity analysis’ in the Methods and Extended Data Fig. 4a). Our results also confirmed that thermophilization is jointly driven by losses in cold-adapted species and gains in warmth-demanding species in forests and largely by gains in warmth-demanding species in grasslands.

Substantial climatic debts existed in all the types of ecosystems studied here. Climatic debts in forests were substantial and generally higher than those in alpine summits and grasslands. This difference is probably because tree canopies provide microclimatic buffering of climate warming, which is not captured by the macroclimate data used to assess global climate warming 23 , 44 (and in our study). Moreover, forest plant life-history traits (for example, long lifespan) probably impede their ability to rapidly adapt to changing conditions, thereby delaying their response to and increasing their climatic debt 6 , 45 . The lower climatic debt in grasslands indicates that shifts in grassland plant communities are more closely tracking climate warming 37 , especially for warmth-demanding species that showed limited climatic debts. We also quantified the thermal suitability of the background climate for plant communities and assessed the extent to which species at the plot level were exposed to climatic conditions beyond their thermal tolerance ranges (see the section ‘Sensitivity analysis’ in the Methods for details). We found significant declines in the probability density of thermal suitability in forests and marginal declines in grasslands, but not on alpine summits (Extended Data Fig. 5a). The proportion of mismatched species significantly increased in forests but not in grasslands or alpine summits (Extended Data Fig. 5b). These results show that thermal mismatch progresses unevenly across ecosystems and that comparing the difference between changes in climate and floristic temperature may have overestimated climatic debts in grasslands and alpine summits. Notably, disturbance levels among ecosystems could also contribute to these differences. Although we were unable to assess this parameter with the current data, it represents a significant direction for future research. Together, these results highlight the ecosystem-specific responses of vegetation to climate warming.

Our results also showed that climatic debts increase with the amount of temperature change, which suggests that rapid climate change increases the risk of species reshuffling and ecosystem disruptions 5 , 7 . This outcome highlights the urgency of proactive strategies to mitigate and adapt to future climate change. As physical and biological drivers of microclimate are substantially different in forests, grasslands and alpine summits 46 , future research should quantify how climatic debt scales with different rates of warming in each of these contrasting ecosystems. In addition to climate warming, other globally changing factors, such as droughts, nitrogen and sulfur deposition, land-use change and elevated CO 2 , may interact with or modulate thermophilization patterns 14 , 47 . Part of the thermophilization signal we detected might be due to niche axes that are correlated with, but not causally linked to, temperature. Multifactor assessments, including other global environmental changes, are therefore needed to more precisely predict the responses of plant communities and inform more integrative conservation strategies.

In summary, we uncovered that thermophilization rates vary significantly among major terrestrial ecosystem types—forests, grasslands and alpine summits—in Europe. Thermophilization is driven largely by increases in warmth-demanding species in grasslands, by decreases in cold-adapted species and increases in warmth-demanding species in forests, and by decreases in cold-adapted species and abundance shifts of persisting species on alpine summits. Although time lags between climate warming and changes in the composition of plant communities have previously been noted 6 , 7 , 12 , 18 , 48 , 49 , we demonstrated that climate change has generated contrasting shifts in plant communities among different ecosystem types. These contrasting patterns suggest that there is asynchronous reorganization among ecosystems, which reflect differences in exposure, biotic resistance and response lags. Notably, ecosystems with accumulating climatic debts are likely to face increased risks of ecological mismatches, functional instability and local species extirpations. Given rapid rates of climate change, climate lags will probably persist over many decades even if the climatic debt must ultimately be paid off. The differing thermophilization rates among ecosystems also highlight the importance of comparing several ecosystems when assessing the impacts of global change on biodiversity, as such quantification can enable forecasting of the duration of biodiversity changes that will continue even if climate warming stabilizes. Our findings provide a basis for projecting community change with accelerating climate change, which will facilitate the development of targeted management strategies, such as assisted migration in forests, adaptive grazing practices in grasslands and conservation of cold-adapted species, to enhance ecosystem resilience and mitigate biodiversity loss under ongoing climate change.

Methods

Databases and vegetation resurveys

We used data from three databases, and each database is a collation of vegetation resurvey data for a certain ecosystem type in Europe (Fig. 1 and Extended Data Table 1). Data for forests were obtained from the forestREplot database (v.2.3; www.forestreplot.ugent.be) with 4,372 plots from 76 sampling sites. The forestREplot is a resurveyed vegetation plot network that focuses on all vascular plants in the understory community (<1 m height, including tree seedlings) across temperate forests in Europe 52 , as they represent on average the most important plant biodiversity compartment in forests 53 . All plots are located in natural (never cleared for other land uses than forests according to the oldest available sources such as historical maps, typically at least 150–300 years 54 ) or semi-natural forests according to a published definition 55 . Grasslands were represented by 1,209 resurveyed vegetation plots from 22 sites that are situated in semi-natural grasslands, which is defined as unploughed mown or grazed areas with limited fertilization 56 , from the GRACE database 57 . Data for alpine summits were from the Global Observation Research Initiative in Alpine environments (GLORIA, ), with 486 plots in natural vegetation with a complete list of vascular plant species from 64 sites 58 . In each site, plant communities were surveyed across permanent (that is, plots were permanently marked in the field) or quasi-permanent (that is, plots were not permanently marked but have known coordinates or clear indications on physical maps for later accurate relocation) 59 plots at two time points, namely the baseline survey and resurvey. Between the two surveys, there was no major change in land use such as clearcutting, replanting, conversion to agricultural land or urbanization, among others. Moreover, at least at the plot level, the same methodology was used in both surveys, such as the number and size of plots and cover-abundance scale (for more details, see refs. 52 , 57 , 58 ). The time spans between the baseline survey and resurvey ranged between 12 and 78 years with a median of 42 years. Specific lengths of the study period per ecosystem type are provided in Extended Data Table 1. The median number of plots per site was 42.5 and 39 in forest and grassland surveys, respectively, and the median size of plots was 200 and 25 m 2 in forest and grassland surveys, respectively 60 . Alpine summits were always resurveyed in 8 spatial sections that together covered the entire area from the highest summit point to the contour line 10 m in elevation below this point, and the median summit area was around 300 m 2 . A complete list of vascular plant species was recorded for all plots in each ecosystem type. Species abundance was recorded as the percentage cover in forest and grassland plots, and as an ordinal scale from one (very rare) to eight (dominant) in summit plots 58 . The ordinal abundance scores were standardized for each plot and survey time to a continuous 0–100 scale via min–max transformation and was subsequently converted to relative abundances, thereby rendering them directly comparable to percentage cover estimates. Vegetation surveys were generally carried out during the peak of the vegetation season (May–August), and plots with a year difference between the baseline survey and resurvey less than 10 years were excluded in our study.

Temperature data

To determine the rate of temperature change between the baseline survey and resurvey, we used monthly mean temperature data to obtain mean growing season (April–September) temperatures, given that this period broadly corresponds to the peak of plant physiological activity, although on alpine summits, this period is shorter (June–August), thereby better reflecting the climatic condition that drives species turnover and community dynamics 1 , 3 , 61 , 62 . We used the data of CRU TS (v.4.07), which provides updated gridded climate datasets at a spatial resolution of 0.5° (around 50 km at the equator) since 1901, based on monthly observations at meteorological stations across the land areas of the world. This dataset was selected because of its temporal depth and consistency, which enabled robust estimation of temperature change over the full period between the baseline survey and resurvey (from 1931 to 2020). For each plot across the ecosystem types, we calculated the climatic temperature change rate (expressed in °C per decade) between the baseline survey and the resurvey as the difference between the mean growing season temperature during the 5 years preceding the baseline survey and the 5 years preceding the resurvey, divided by the time span between the baseline survey and the resurvey 1 , 3 , 18 . The resulting distributions of mean growing season temperature change rates are presented in Extended Data Fig. 1b.

Thermophilization

To calculate thermophilization rates, we followed established and widely used methods in the literature based on thermal niche and distribution range data of species 1 , 3 , 6 , 15 , 18 , 20 , 63 . First, we extracted gridded long-term (1970–2000) climatic data on monthly mean temperatures of the growing season (April–September) from WorldClim (v.2) as a raster file with 10 arcmin resolution (about 20 × 20 km) 64 . Then, to obtain distribution data of the species, we estimated species range sizes as the area of occupancy (AOO) 65 using point occurrence records in the Global Biodiversity Information Facility (GBIF, www.gbif.org) 66 . We retained records from 1930 onwards to ensure full temporal coverage of the distributions of species across the baseline survey and resurvey intervals, the earliest of which began in 1931 (Extended Data Table 1). To minimize spatial sampling bias, occurrences were aggregated to 20 × 20 km resolution, and the AOO was calculated on the basis of the number of unique occupied grid cells per species, thereby not only avoiding overrepresentation of well-sampled areas but also aligning with the resolution of the climate data used. Records with missing or erroneous coordinates, high spatial uncertainty (>10 km) or flagged as outliers were removed using standard data-cleaning procedures. After filtering, a total of 1,919 species (85.8%) out of the 2,237 species identified to the species level included in our study across ecosystem types were matched and retained. This GBIF-derived range sizes for plant species in Europe has been found to correlate strongly with expert-drawn range maps and can therefore provide good relative estimates of the distributions of species 60 , 67 . Next, to estimate the realized thermal niches of species, defined as the climatic niche space that is genuinely occupied by a species 66 , we set the study area extent to all European land between 15° W and 60° E and 25° N and 75° N, with grid cells of 20 × 20 km. For each species, we randomly sampled 1,000 grid cells, with replacement, within the range size represented by the AOO 18 , 32 . To generate consistent data series for all species (1,000 values), we used random sampling with replacement instead of a random sample or all grid cells, which allowed species with small distribution ranges (that is, <1,000 grid cells) to be included under the same sampling scheme and ensured comparability of climatic niche estimates across species. However, because rare species (defined here as species with range sizes smaller than 1,000 grid cells, which accounted for 36.4% of the total species pool; Extended Data Fig. 1c) are less frequently encountered, sampling efforts may not consistently detect them and therefore lead to potential underestimation or overestimation of their actual contribution to thermophilization of vegetation communities. We therefore first estimated the thermophilization rate and climatic debt only with species with AOO ≥ 1,000 grid cells and then reran the analyses with all the species included in our database to check whether the exclusion of rare species biased our results (see the section ‘Sensitivity analysis’ below).

We calculated the floristic temperature (T plant , °C) for each plot by using the thermal niches of all species present during each survey (baseline survey and resurvey). To account for the variability and uncertainty in the thermal preferences and niche widths of species, we generated a distribution of plot-level floristic temperatures at each survey by resampling 1,000 times from the thermal niches of species (that is, from the 1,000 grid cells), with weights defined by species abundance, following a previously described method 68 . Then, we estimated the overall thermophilization rate for each plot as the difference between the median of the floristic temperature distribution of the resurvey minus the baseline survey (ΔT plant ), divided by the time span (Δt, in decades) between the two survey periods. Given that thermophilization represents a slow ecological response to climate warming, we therefore used decadal rates (°C per decade) in this study to better reflect long-term trends and reduce sensitivity to interannual variability. A positive value therefore reflects an overall temporal increase of floristic temperature in °C per decade.

To further evaluate the contribution of declines in cold-adapted species and increases in warmth-demanding species to the overall thermophilization rate, we quantified shifts in the left (5th percentile) and right (95th percentile) tails of the plot-level distribution of floristic temperatures per survey period, respectively 1 . A positive value in the left versus right tail shift indicates thermophilization due to a decline in relatively cold-adapted species versus an increase of the frequency of more warmth-demanding species, respectively. As we were specifically interested in this per-plot thermophilization, it is important to note that this calculation of shifts in the cold-adapted versus warmth-demanding species is a per-plot relative quantification based on the left and right tail of floristic temperature distribution shifts per plot 1 . As a result, a cold-adapted species in a plot in a warm site can be a warmth-demanding species in a plot in a colder site. Finally, as a sensitivity analysis, we also converted all vegetation data into presence/absence data and recalculated thermophilization rates and climatic debt (see the section ‘Sensitivity analyses’ below).

Climatic debt

Climatic debt is defined as the degree to which the rates of compositional changes in vegetation are lagging behind the rates of climate warming. The climatic debt per plot between the baseline survey and resurvey was calculated as (ΔT clim /Δt) – (ΔT plant /Δt) (following refs. 2 , 24 ), where ΔT clim and ΔT plant are changes in the actual mean climate temperature of the growing season and floristic temperature (the median, 5th and 95th percentile of the distribution of floristic temperatures) between the baseline survey and resurvey, respectively. A large climatic debt implies a large mismatch between the dynamics of plant communities and climate change 7 , which indicates that a relatively large number of species are occurring under suboptimal thermal conditions. Positive or negative climatic debts indicate that the response rates of plant communities are slower or faster, respectively, than climate warming rates.

Statistical analyses

We used Bayesian mixed-effects models with the brms package (v.2.22.0) 69 in R (v.4.4.2) for all the statistical analyses and the 95% CI to determine statistical clarity 70 . We adopted a Gaussian regression model with an identity link to evaluate whether thermophilization rates and climatic debts (°C per decade) were significantly different from zero in each type of ecosystem by fitting ecosystem type as a fixed-effect factor and ‘site’ as a random-effect term. The suitability of Gaussian assumption was evaluated using quantile–quantile plots of the response variables, which indicated approximate normality (Extended Data Fig. 6). Also, to account for the potential impacts of regression to the mean effects 71 , we included the mean climate temperature of the growing season at the baseline survey (T clim baseline ) as a covariate, given the detected significant linear relationship between T clim baseline and temperature change rate (Extended Data Fig. 7a). This approach helped separate baseline effects from true temporal trends. We also fitted plot size as a fixed-effects factor in the model, which was specified as the following formula (model 1):

Thermophilization rate or climatic debt ~ ecosystem type + plot size + T clim baseline + (1|site).

Then, to assess the relationship between actual climatic temperature change rate and the thermophilization rate or climatic debt based on vegetation surveys across ecosystem types, we fitted the interaction of actual climatic temperature change rate and ecosystem type as fixed-effects factors (model 2):

Thermophilization rate or climatic debt ~ ΔT clim × ecosystem type + plot size + T clim baseline + (1|site).

For each model, we fitted 4 chains of 8,000 iterations, excluding the first 2,000 iterations as the warm-up in each chain. We used the Gelman–Rubin statistic (hatR < 1.01) to ensure model convergence, and all the models in our analyses converged. Data visualization was performed using the ggh4x package (v.0.3.1).

Sensitivity analyses

Thermophilization was primarily estimated using abundance-weighted data, which emphasize shifts in community structure by giving more weight to dominant species. This approach provides insights into changes in ecological dominance and potential functional consequences. However, it may also underrepresent the role of rare or newly colonizing species, which can be important indicators of climate-driven species turnover and bias the results because of methodological and observer differences in abundance measurement. To address these limitations and to assess the robustness of our findings, we therefore converted all abundance data into robust presence/absence data following previous studies 1 , 18 and reran the analyses for estimations of thermophilization rate and climatic debt. This complementary approach captured changes in species composition regardless of abundance, which enabled a more comprehensive assessment of community thermal shifts. The rationale for performing this analysis is that the time intervals between the baseline surveys and resurveys were mostly sufficiently long (all >10 years) to allow for species replacements due to colonization and extinction. We found similar results of thermophilization rate and climatic debt generated from vegetation abundance-weighted and presence/absence data (Extended Data Fig. 2).

Rare species with small range sizes are often ecological outliers with distinct habitat requirements and thermal tolerances, and their presence may therefore disproportionately influence the calculation of average thermal preferences of plant community 72 . Also, rare species are less likely to be consistently detected across sampling events owing to their low abundance, and this variability in detection can introduce noise in the data, which makes it difficult to accurately assess trends in species composition over time 73 . We therefore excluded rare species with AOO < 1,000 grid cells in our analyses. However, rare species themselves are important components of the vegetation community and can be important indicators of climate-driven species turnover. Therefore, to test whether the exclusion of rare species in our analyses may affect estimations of thermophilization rate and climatic debt, we reran such analyses by including all the species in our dataset. The results showed that the inclusion of rare species did not affect the overall patterns of thermophilization rates and climatic debts across all ecosystems (Extended Data Fig. 2b), which indicated that the observed patterns are not an artefact of sampling bias towards rare taxa and reflect consistent community-level responses to climate warming.

In our study, we classified cold-adapted species and warmth-demanding species on the basis of the lower and upper 5% percentiles of the floristic temperature distribution. Although this approach focused on species at the thermal distribution margins, it may also result in relatively small groups that are sensitive to outliers or sampling variation. To evaluate whether our results were disproportionately influenced by these few extreme cases, we conducted a sensitivity analysis using the lower and upper 20% percentiles of the floristic temperature distribution to define cold-adapted and warmth-demanding species, respectively. The results were consistent with our original findings, showing similar patterns in the direction of thermophilization responses across ecosystems (Extended Data Fig. 7b). However, the magnitude of the effects was slightly attenuated, as expected, owing to the inclusion of species closer to the median of floristic temperature distribution. These results confirm that our conclusions are robust to the choice of quantile thresholds used to define thermal adaptation.

As a supplementary analysis for the method used to quantify the contributions of shifts in cold-adapted versus warmth-demanding species to thermophilization, we further partitioned the observed changes in community-weighted, at the plot level, mean temperature (ΔT CWM , °C) into three additive components, namely, species gains, species losses and abundance shifts of persisting species 40 , 74 . This framework enabled us to attribute community-level floristic temperature change into different forms of species turnover and reassembly. A gains effect was defined as the average thermal niche of species that were newly recorded in the resurvey but absent from the baseline survey minus the baseline T CWM . The losses effect was calculated as the baseline T CWM minus the mean thermal niche of species that disappeared between the baseline survey and resurvey. The gains effect describes the directional impact of the arrival of novel species, and the losses effect reflects the impact of species loss on thermal restructuring. The abundance shifts effect was quantified as the ΔT CWM between the two surveys using only the set of species identified in both baseline surveys and resurveys, weighted by their relative abundances in each survey, which captures compositional shifts among persisting species that contribute to warming or cooling, independent of species turnover. For each plot, we estimated gains and losses effects using the median of species-specific thermal niche distribution as the respective mean thermal difference relative to the baseline T CWM . The abundance shift effect was estimated as the change in T CWM over time based only on the subset of shared species in both surveys. The sum of the three effects produced a total estimate (estimated ΔT CWM ), which was then compared with the observed ΔT CWM to assess the explanatory power of mechanism. We also estimated the relative contribution of each mechanism by dividing the absolute value of its effect by the sum of the absolute values of all three effects, and the contributions were then expressed as percentages to facilitate comparison among mechanisms. We identified the dominant mechanism for each plot as the one with the largest percentage contribution. The results revealed significant linear relationships between estimated and observed ΔT CWM across ecosystems, with R 2 values ranging from 0.23 to 0.53 (Extended Data Fig. 4b), which indicated that the three mechanisms collectively provide a reasonable explanation for the observed thermophilization patterns. Building on this, the results further suggested that gains and losses were the dominant contributors to thermophilization of understory vegetation in forests, losses played the primary role in grasslands, whereas losses and abundance shifts of persisting species were the major contributor on alpine summits (Extended Data Fig. 4a). These results highlight the fact that the relative importance of gains, losses and abundance shifts in driving thermophilization is ecosystem-dependent and reflect distinct community dynamics and responses to climate warming across ecosystems.

Despite the widespread use of the climatic debt metric by comparing ΔT clim and ΔT plant , it could neglect the variation in the thermal niche space in a community because one can obtain a climate debt but the background observed climate of the resurvey remains in the thermal niche space of all species in a community. Therefore, to test the robustness of the climatic debt results in our study, we further quantified the thermal suitability of background climate for plant communities and assessed the extent to which species in each community, at the plot level, were exposed to background climatic conditions beyond their thermal tolerance ranges. Specifically, for thermal niche suitability, we estimated the probability density of floristic temperature distribution (that is, T plant ) at the observed climatic temperature (that is, T clim ) for each plot and survey time. T plant distribution values were generated using the 1,000 resampled values of species-specific thermal niches and abundance data, as described above. Then, for each value of the distribution, we computed a kernel density estimate and extracted the probability density at the corresponding T clim value using linear interpolation. The change in probability density between the baseline survey and resurvey indicates whether the local climate has become more or less aligned with the realized thermal niche of the community, and a decline in density suggests an increasing mismatch between climate and community thermal preferences. The results showed significant and marginal declines in the probability density at T clim from the baseline survey to resurvey in forests and grasslands, respectively (Extended Data Fig. 5a), which indicated an increasing climatic mismatch. However, no change in the probability density at T clim was detected on alpine summits.

Furthermore, to estimate the proportion of thermally mismatched species, we derived 2.5th and 97.5th percentiles of the thermal niche distribution for each species. For each plot and survey time, we identified all recorded species and classified each as ‘thermally mismatched’ if the background T clim fell outside its thermal niche bounds (that is, between the 2.5th and 97.5th percentiles of thermal niche distribution). We then calculated the proportion of mismatched species per plot as the number of mismatched species divided by the total number of species present, and estimated the change between the resurvey and baseline survey. This metric therefore captures the prevalence of species experiencing thermal stress relative to their historical niche distribution. We found that the proportion of mismatched species showed a significant increase in forests, but not in grasslands or on alpine summits (Extended Data Fig. 5b). These results suggest that thermal mismatch is progressing unevenly across ecosystems. In forests, the observed increase in mismatched species, together with declining probability density, indicates that climatic conditions are moving away from the realized thermal preferences of the plant community. By contrast, the lack of detectable change in grasslands and summits may reflect either reduced climatic displacement relative to forests or species assemblages with broader thermal affinities. Overall, these results suggest that comparing the difference in ΔT clim and ΔT plant may have overestimated climatic debts in grasslands and alpine summits. The data also refine our understanding of climate debt by highlighting both directional responses and ecosystem-specific characteristics in tracking climate change.

Alpine summits are characterized by rugged topographic gradients that generate strong microclimatic heterogeneity, which may have strong impacts on vegetation thermophilization 36 . To assess the potential impacts of alpine summit topography on thermophilization, we derived TRI and slope data from a 30-m resolution digital elevation model using the terrain function from the terra R package (v.1.8.29) and extracted snow cover data from CHELSA (v.2.1) with a resolution of 30 arcsec (around 1 km) 75 . We found that neither TRI nor slope showed significant relationships with thermophilization rates (Extended Data Fig. 3a), which may be because in our database, only summit-level coordinates were available, and metrics such as TRI fail to capture the fine-scale microclimate that more strongly influence community responses. However, the thermophilization rate was significantly lower in north-facing slopes compared with south-facing slopes, and the snow cover change rate was lowest on northern aspects (Extended Data Fig. 3b). These findings indicate that microclimatic habitats, generated by rugged topographic gradients and slope aspect, could be a key mechanism in shaping thermophilization patterns on alpine summits.

To further evaluate whether the use of 20 × 20 km temperature data would bias the estimation of thermophilization rates on alpine summits, given that data at such a resolution may not adequately capture the fine-scale thermal conditions experienced by alpine plants, we extracted the climatic data of monthly mean temperature of the growing season from WorldClim (v.2) as a raster file with 30 arcsec resolution (about 1 × 1 km) 64 and re-estimated thermophilization rates on alpine summits using the same method and compared the results with those estimated using 20 × 20 km resolution temperature data. The results showed that the two sets of estimates were strongly correlated (Extended Data Fig. 8a), and only a small difference was found between the thermophilization rates of the 95th percentile as estimated using Bayesian mixed-effects models (Extended Data Fig. 8b). Nevertheless, this difference will not affect the interpretation of the overall patterns of thermophilization on alpine summits, which demonstrated a high degree of consistency between the estimation of thermophilization derived from fine-resolution and coarse-resolution climate data, which indicates that the use of the 20 × 20 km resolution data remains appropriate and our core conclusions are robust to the choice of climate data resolution.

To evaluate whether spatial autocorrelation might bias the ecosystem effects on thermophilization rates, we incorporated a spatial Gaussian process (GP) term into our Bayesian mixed-effects models (model 3) 76 :

Thermophilization rate ~ ecosystem type + plot size + T clim baseline + GP(longitude, latitude) + (1|site).

We then compared the posterior distributions of the fixed-effect estimates for ecosystem type between the spatial and non-spatial model specifications. The results showed a high degree of consistency, with substantial overlap in the 95% CIs (Extended Data Fig. 8c). This strong agreement between the models indicates that spatial autocorrelation did not influence the estimated ecosystem effects, which is probably because the random intercept for ‘site’ captured most of the large-scale spatial structure. Given the robustness of the estimates and for the sake of model parsimony, we present the results from the simpler, non-spatial models.

Sources of uncertainty

Although we document general patterns of vegetation thermophilization and climatic debts here, uncertainties remain, which reflect inherent limitations in the data. First, declines and increases in species abundance may reflect sampling variation (for example, shifts in plant locations or inexact placement of quasi-permanent plots) rather than actual community changes between baseline surveys and resurveys, and the inclusion of quasi-permanent plots may therefore influence the overall patterns of vegetation response to warming. In our study, all the plots on alpine summits are permanent, whereas most of the plots are quasi-permanent in forests and grasslands. Regardless of their formal designation, all plots were located again using high-precision GPS coordinates or clearly documented map references, which ensured accurate resurveys across time. Although metadata did not always distinguish between permanent and quasi-permanent plots, the accuracy of finding the same location was closely comparable across plot types. Importantly, we restricted our inferences to averages across plots and at the community level rather than focusing on results from individual plots. Hence, small relocation errors only added some uncertainty to results rather than introducing any systematic bias in the magnitude or direction of effects at the community level (thermophilization). Moreover, many plants move over time, which means that even permanent plots that were exactly found at the correct location show shifts in composition and abundance. Our key conclusions are robust to these sources of noise as supported by a previous study 77 using a subset of the forest plots in our dataset. That study showed that the impact of relocation errors is limited when focusing on overall temporal changes in vegetation communities 77 . As such, the inclusion of both types of plots does not introduce systematic bias or influence the observed ecosystem-level patterns 77 , 78 .

Second, although our study provides a unique opportunity to compare vegetation responses to climate warming across ecosystems using harmonized plot-level resurveys and standardized climate metrics, we acknowledge that full methodological standardization is not possible. The underlying datasets were originally collected by different monitoring networks using different survey protocols, plot sizes and time intervals. These differences may introduce some uncertainty into cross-ecosystem comparisons. Nevertheless, by applying consistent analytical methods and filtering to include only plots with comparable abundance data, we minimized potential bias and ensured the robustness of our main conclusions.

Third, another source of uncertainty may stem from the variation in time intervals between the baseline survey and resurvey across plots (Extended Data Fig. 1a). Some plots span multiple decades and are therefore more likely to capture clear climatic changes, whereas others have shorter intervals with potentially weaker climate signals. This temporal heterogeneity, which also differs among ecosystems, may influence the observed rates of thermophilization. Although such variation is inherent in resurvey-based studies and reflects logistical and historical constraints, it represents an important consideration when interpreting the strength and consistency of observed patterns.

Last, the distribution of our plots was skewed towards Northern and Central Europe, with comparatively fewer sites in Southern Europe (Fig. 1). Although the dataset still covers a broad range of European climatic and ecological gradients and is therefore suitable for continental-scale analysis, this spatial imbalance would still represent a source of uncertainty. Future research should aim to increase plot coverage in underrepresented southern regions of Europe to ensure more robust and spatially comprehensive assessments of thermophilization patterns.