Florian Zellweger
email: florian.zellweger@wsl.ch

Pieter De Frenne

Jonathan Roger Michel

Henri Lenoir

Pieter Vangansbeke

Kris K Verheyen

Markus Bernhardt-Römermann

Lander Baeten

Radim Hédl

Imre Berki

Jörg Brunet

Jonathan Lenoir

Hans Van Calster

Markéta Chudomelová

Guillaume Decocq

Thomas Dirnböck

Tomasz Durak

Thilo Heinken

Bogdan Jaroszewicz

Martin Kopecký

František Máliš

Martin Macek

Malicki Marek

Tobias Naaf

Thomas A Nagel

Adrienne Ortmann- Ajkai

Petr Petřík

Remigiusz Pielech

Kamila Reczyńska

Wolfgang Schmidt

Tibor Standovár

Krzysztof Świerkosz

Balázs Teleki

Ondřej Vild

Monika Wulf

David Coomes

Adrienne Ortmann-Ajkai

Forest microclimate dynamics drive plant responses to warming

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Short title: Forest plant responses to microclimate warming

Abstract: Climate warming is causing a shift in biological communities in favor of warmaffinity species (i.e., thermophilisation). However, species responses often lag behind climate warming, but the reasons for such lags remain largely unknown. Here we analyze multidecadal understorey microclimate dynamics in European forests and show that thermophilisation and the climatic lag in forest plant communities are primarily controlled by microclimate. Increasing tree canopy cover reduces warming rates inside forests, but loss of canopy cover leads to increased local heat that exacerbates the disequilibrium between community responses and climate change.

Reciprocal effects between plants and microclimates are key to understanding the response of forest biodiversity and functioning to climate and land-use changes.

One Sentence Summary: The cooling effect of tree canopies can reduce forest biodiversity loss in times of climate warming.

Main Text

Climate warming is having profound effects on ecological processes and biodiversity thus seriously affecting ecosystem functioning and human well-being (1)(2)(3)(4). Our knowledge and predictions about biotic responses to anthropogenic climate warming are largely based on air temperature data measured at official meteorological stations, which record free-air (macroclimate) temperature in open areas at 1.2 to 2 m above short grass (5,6). However, most organisms on Earth experience temperature conditions that differ from the macroclimate, mainly because the topography and vegetation create heterogeneous microclimates near the ground via interception of solar radiation, air mixing and evapotranspiration (7,8). Local microclimates may explain why responses of biological communities and ecosystem processes are often partially uncoupled from macroclimate warming (6,(9)(10)(11)(12)(13)(14).

Range shifts towards higher latitudes and elevations are now commonly observed for many species and systems, as organisms shift their geographical distributions in order to track their thermal requirements (15). With rising temperatures at a location, the presence or abundance of species adapted to higher temperatures is therefore expected to increase, whereas species adapted to lower temperatures may decline and eventually become excluded. Such directional shifts in community composition in favour of warm-affinity species are referred to as "thermophilisation", a phenomenon that is increasingly documented in terrestrial and marine plants and animals (12-14, 16, 17). Yet, the thermophilisation rate of many biological communities is not keeping pace with the velocity of contemporary macroclimate change (18,19), leading to a climatic lag or debt in community responses to macroclimate warming (10)(11)(12)(13). Climatic debt effects may be the inevitable consequence of habitat fragmentation, slow dispersal and long lifespans (20), but the magnitude of the climatic debt may also be affected by different warming rates of localised microclimates. However, we know very little about how microclimates have changed over time and it is unclear how any such change has modulated the temporal thermophilisation rate and climatic debt observed in plant and animal communities (12)(13)(14)17). Effects of changes in vegetation cover on microclimates near the ground could have either accelerated or counteracted the effects of macroclimate warming on biological communities, but a long-term, large-scale and multi-taxa assessment of these effects is currently missing.

Microclimates are perhaps nowhere more evident than in forests, owing to their threedimensional canopy structure that drives shading, air mixing and evapotranspirative cooling (7,21). The tree canopy buffers forest-floor temperatures against extreme heat (9), and this buffering capacity constantly changes with tree species, growth and mortality, leading to highly dynamic microclimates across space and over time (22). Accounting for changes in canopy cover and the associated microclimate dynamics is therefore important to better understand the response of forest biodiversity to climate change. Here we provide multi-decadal evidence of forest sub-canopy temperature changes, enabling the comparison between anthropogenic climate change, as measured by weather stations (macroclimate), and forest microclimate dynamics triggered by canopy cover changes over time. To this end, we combined sub-canopy temperature measurements in 100 forest stands in temperate forest in Europe with 2955 permanent vegetation plots from 56 regions, where each plot has been resurveyed over a period of 12 to 66 years (23) (Fig. 1A, Fig. S1). Using a continental-scale analysis of forest microclimates based on in situ empirical temperature and canopy cover data we then predict changes in understorey temperature during the growing season, building upon the relationship between canopy cover and the buffering of macroclimate temperatures (21) (Fig. 1B &C). (purple circles, scaled to the number of plots as indicated in the top right) across the temperate forest biome (green area) in Europe. We representatively sampled microclimate temperature in 100 forest stands, i.e., in 10 stands in each of 10 regions (black circles; effective n = 96), to estimate the maximum (macroclimate) temperature buffering during the growing season as a function of canopy cover (23) (Fig. S2). (B), Schematic overview of the method used to approximate microclimate change in the forest understorey. In this example, canopy cover at the time of the baseline survey was higher than during the resurvey, resulting in a decrease in macroclimate temperature buffering from 2 to 1 °C, which in turn led to a relatively larger increase in microclimate warming (20 to 23 °C) compared to macroclimate warming (22 to 24 °C). The relationship between canopy cover and the buffering of maximum macroclimate temperature was empirically assessed across the study area (Fig. S2) (21). (C), The rate of macroclimate change is plotted against the rate of microclimate change, with the black bisecting line representing the 1:1 relationship. Micro-and macroclimate have both significantly warmed (see text for statistical results). The distributions of values in the rates of micro-and macroclimate change are indicated by grey rugs on each axis. Microclimate change rates are 45 % more variable than macroclimate change rates and macroclimate change rates only accounted for about half of the variation in microclimate change rates, as indicated by the marginal (conditional) R-squared value of 0.52 (0.69). All statistical results are based on mixed-effects models with region as a random-effect (intercept) term.

We found that temporal changes in canopy cover varied greatly across the 56 European regions studied, ranging from -110 % (significant canopy opening) to +113 % (strong densification of the canopy) (1st and 99th percentile of distribution, respectively) with a mean canopy cover change not significantly different from zero (+2.6 %; mixed-effects models P = 0.426, Fig. S3). To predict how the microclimate in the understorey of each plot had changed between the baseline survey and resurvey, we applied a previously published statistical model to estimate temperature buffering as a function of canopy cover (21, 23) (Fig. S2). The predicted maximum temperatures in the forest understories have significantly warmed over the past decades, with mean (+/-standard error of the mean, SEM) rates of 0.40 +/-0.04 °C and 0.38 +/-0.03 °C per decade for micro-and macroclimate warming, respectively (both estimates of warming rates are based on mixed-effects models: P < 0.001). Yet, the rate of microclimate change was 45 % more variable (1 st and 99 th percentiles: -0.32 -1.36 °C per decade) than the rate of macroclimate change (1 st and 99 th percentiles: -0.08 -1.08 °C per decade) (Fig. S4). The rate of macroclimate change was significantly (P < 0.001) related to the rate of microclimate change, but left 48 % of the total variation in microclimate change unexplained (slope: 1.05, Rsquared: 0.52, P < 0.001) (Fig. 1C).

To quantify the thermophilisation, we inferred the thermal affinity for each vascular plant species present in our dataset from its current distribution ranges. Using these species-specific temperature affinity values, we calculated the rate of change in the community-based maximum temperature affinity values between the resurvey and baseline survey (14, 23) (Fig. S6). We expected changes in maximum temperature affinity values to be most closely related to changes in micro-and macroclimate maximum temperatures during the growing season (23). This biotic reconstruction of temperature changes based on the observed changes in the composition of species assemblages has been widely used to assess community-level climate-change impacts in a variety of terrestrial and marine taxa (12)(13)(14). The resulting thermophilisation rates across the 2955 permanent plots ranged from -0.84 to 1.05 °C per decade, with a mean of 0.01 +/-0.01 (SEM) °C that was not significantly different from zero (P = 0.09) (23). Remarkably, the thermophilisation rate of forest understorey vegetation was positively linked to the rate of microclimate warming (scaled slope estimate: 0.02, 95 th confidence interval (CI): 0.01 -0.03, P < 0.001), but not to macroclimate warming (scaled slope estimate: -0.002, CI: -0.01 -0.01, P = 0.70) (Fig. 2). To quantify how forest microclimate affected the observed climatic debt accumulated by a plant community in a given plot, we subtracted the thermophilisation rate (∆Tplant) per unit of time (∆t) from the rate of microclimate change (∆Tmicro) per unit of time (i.e., microclimate debt: (∆Tmicro/∆t) -(∆Tplant/∆t)) and from macroclimate change (∆Tmacro) per unit of time (i.e., macroclimate debt: (∆Tmacro/∆t) -(∆Tplant/∆t)), in each focal plot. Despite surprisingly similar means for the microclimatic debt (0.38 +/-0.04 (SEM) °C per decade) and macroclimatic debt (0.37 +/-0.04 (SEM) °C per decade), the climatic debts calculated using macroclimate data underrepresent the variability in microclimatic debt (Fig. S7). Importantly, we found the greatest microclimate warming in areas where canopy cover and thus the temperature buffering declined, and these were also areas where the microclimatic debt was greatest (Fig. 3). Despite higher thermophilisation rates with increasing microclimate warming (Fig. 2), locally increased heat due to a reduction of canopy cover impedes the ability of understorey plant communities to respond to such high rates of warming. On the contrary, we found lower microclimatic debts in sites with increased canopy cover; there, temperature buffering led to a cooling effect during the growing season. These patterns remain hidden when analysing climatic debts based on macroclimate data (Fig. 3). Realistic assessments of the current pressures on communities due to climate warming thus requires long-term data on microclimate change. Canopy cover dynamics have triggered microclimate changes over time in forest interiors that can differ considerably from macroclimate changes outside forests. This has important implications for predicting biodiversity responses to climate and land use (e.g., forest management) change, which interactively drive the emergence of novel thermal environments. With the predicted increase of heat-waves (4), many species and communities may greatly suffer from loss of canopy cover, e.g. following tree harvesting or dieback (24). The resulting impacts are serious, because forests harbour most of the terrestrial biodiversity, and because many ecosystem services and livelihoods critically depend on forest biodiversity (1,25). Forest managers and policymakers should therefore consider the effects of different forest management practices on local microclimates in their endeavours to safeguard forest biodiversity in a warming world.

Our results support the hypothesis that the thermophilisation rate in forest understorey plant communities is primarily driven by the rate of sub-canopy microclimate change (10,12) and not by the rate of macroclimate change. This finding provides empirical evidence that microclimate change ultimately drive organismal responses to climate change, a frequently ignored fact when using macroclimate data to study biotic responses to climate change (8,26,27).

Increasing climatic debts in community responses to climate change mean that a growing number of species are occurring in sub-optimal climatic conditions, potentially accelerating the loss of biodiversity. Our results suggest that microclimates can amplify as well as decrease the disequilibrium between community responses and macroclimate change, suggesting that climatic debts based on macroclimate data (13,20) should be revisited and interpreted with caution. Microclimate data, therefore, considerably improves the local relevance of the climatic debt concept for climate-change impact assessments on biodiversity, a field that will benefit from emerging datasets and methods to quantify microclimatic variability in space and over time (28,29). In fact, high rates of microclimate warming can greatly exceed the capacity of understorey plant species to spatially track their thermal niche, suggesting that other factors limiting species establishment, such as plant-water relations (30), habitat fragmentation as well as dispersal limitation may impede or severely delay community responses (11). Such effects may outweigh remedial effects of microclimate variability to reduce the pressures of climate change on biological communities, e.g. by providing thermal refuges and facilitating short-distance thermal niche tracking (27,31,32).

In sum, we provide evidence that forest community responses to climate change are most closely related to microclimate and not to macroclimate change. Despite widespread evidence for thermophilisation trends in plant communities (14,17), many community responses are strongly lagging behind warming, thereby accumulating a climatic debt (10). Growing pressures from woody biomass extraction and the increasing vulnerability of forests to climate change will lead to frequent canopy cover disturbance and tree dieback (33), which will severely intensify the emergence of adverse thermal habitat conditions for many species, impeding the ability of communities to keep track with anthropogenic environmental changes. Our findings also show that climate-change impacts on forest plant communities have been reduced by higher standing stocks and associated cooling following increases in thermal buffering (34). Accounting for microclimate in global-change impact assessments on forest biodiversity and functioning is crucial if we are to better understand and counteract the increasing pressures imposed on forests.

Materials and Methods Methods

Vegetation surveys We used data from a resurveyed vegetation plot network across deciduous temperate forests in Europe (www.forestreplot.ugent.be) (35) (Fig. 1A). Each of the 2955 permanent or quasi-permanent plots belongs to one of 56 studied regions and has been surveyed at two time points, i.e., during the baseline and resurvey years between 1934 and 2017 (Fig. 1A; Supplementary Fig. S1). The median time interval between the baseline survey and the resurvey was 38, ranging from 12 to 66 years, which is sufficient to detect directional change in the forest understorey layer (14). All plots are located in ancient (never cleared for any other land use than forest since the oldest available sources such as maps, typically at least 150-300 years ( 36)) semi-natural forests (37) (e.g., excluding fertilized and heavily managed plantations). Between the two surveys, none of the plots have been subject to clearcutting and/or replanting. All plots were either permanently marked (i.e. physical mark permanently visible in the field, such as labelled trees) or quasi-permanent (i.e. with known and precise coordinates for subsequent relocation) and the plot sizes ranged between 25 and 1300 m 2 (for more information about the plots see Refs (14,35,38,39)). Plot size was neither related to the rate of microclimate change, nor to the thermophilisation rate, nor to the climatic debt (mixed effects models; all P > 0.05), suggesting that plot size did not affect our findings and conclusions. All vegetation surveys were performed during the peak of the vegetation season. We used complete understorey (below 1m) vascular plant species lists and determined canopy cover (sum of the visually estimated tree (> 7m) and shrub (1 -7m) species' canopy cover in percentage per species) changes over time for each plot (Supplementary Fig. S3) (14). Changes in canopy cover occurred due to both natural as well as anthropogenic causes, and showed no geographic pattern, which is in line with the fact that the forests considered here are mostly subject to small-scale disturbances. Note that our estimate of canopy cover is commonly larger than 100 % because tree and shrub canopies tend to significantly overlap (14). Taking into account this overlap is important to calculate the microclimate buffering (21). The total number of vascular plant (tree, shrub, graminoid, forb and fern) species recorded in the 2955 plots and over both the baseline survey and the resurvey was 1292.

Macroclimate temperature data

To quantify the rate of macroclimatic temperature change between the baseline survey and the resurvey, we obtained monthly mean maximum temperature data during summer (June, July, August). Two sources of interpolated data were used: TerraClimate (40) and data from the Climatic Research Unit (CRU, TS v.4.03) (41). The TerraClimate dataset provides monthly climate data for global terrestrial surfaces over the period 1958-2017 at a spatial resolution of 1/24 ° (ca. 4 km at the equator), and has been shown to improve coarser resolution gridded dataset (40). We therefore here use TerraClimate data rather than the CRU data. However, in 251 (7.3%) out of the 2955 permanent or quasi-permanent plots, the baseline survey was performed prior to 1958, when no TerraClimate data was available. In these cases only, we had to use the CRU data, which provides updated gridded climate datasets at a spatial resolution of 0.5 ° (ca. 50 km at the equator) since 1901, based on monthly observations at meteorological stations across the world's land areas. For each plot, we calculated the rate of maximum temperature change between the baseline survey and the resurvey (expressed in degrees Celsius per decade) as the difference between the average summer maximum temperature during the 5 years preceding the baseline survey and the resurvey, divided by the time interval between the baseline survey and the resurvey (14). The resulting distribution of rates in macroclimate change is presented in Supplementary Fig. S4. The use of the coarser-gridded CRU data (i.e. analysing plot-level changes with complete TerraClimate data coverage only) did not affect our results (see section Sensitivity analyses).

Microclimate temperature measurements and buffering effect

To representatively sample microclimate in our resurveyed vegetation plot network, we selected 10 regions in each of which we selected 10 forest stands along a gradient of canopy cover, from dense and closed canopy to open. The resulting sample size was n = 96 (four out of 100 loggers got damaged or stolen). In each stand, we recorded subcanopy air temperature every hour for a full year (22 February 2017 to 21 February 2018), using Lascar EasyLog EL-USB-1 temperature logger devices with an accuracy of +/-0.5° Celsius. The loggers were shielded in white, ventilated PVC shields and attached to the north side of a tree trunk (DBH > 25 cm) to prevent exposure to direct sunlight, at 1 m above the ground. The hourly temperature data were then aggregated to daily maximum temperatures, and retained for the summer months June, July and August. More details on the methods can be found in Zellweger et al. (21) Maximum temperatures during summer are often considerably lower inside forests than outside, owing to the temperature buffering (or offset) effect of the trees and shrubs (9). To quantify the degree to which the canopy cover buffers free-air (macroclimate) maximum temperatures in our 96 stands, we compared our microclimate data to the daily maximum temperatures measured at an official weather station closest to each plot. The temperature buffering is thus expressed as the microclimate maximum temperature minus the macroclimate maximum temperature. This analysis revealed a positive, non-linear relationship between temperature buffering and canopy cover (Supplementary Fig. S2) (21). Based on cross-validation procedure, we found that combining canopy cover with distance to the coast to predict the maximum temperature buffering resulted in an improved predictive accuracy with R-squared and Root Mean Square Error (RMSE) values of 0.33 and 0.92, respectively (the R-squared value of the full model, without splitting the data for independent cross-validation, was 0.38). Distance to the nearest coastline was extracted for each plot, using free vector map data 1 . Adding topographic variables such as elevation, topographic position or slope did not improve the predictive accuracy of the model, which is attributable to the fact that local canopy cover and associated shading may mainly control the local maximum temperature buffering in the primarily lowland forests considered here (median, 1 st and 99 th percentile elevation: 325, 16 and 1110 meters above sea level) (21). The predictive accuracy with R-squared and Root Mean Square Error (RMSE) of the model solely based on canopy cover was 0.16 and 1.02, respectively. We applied a conservative cross-validation procedure with blocked data splitting, accounting for our hierarchical sampling design ('region' as a random effect) (42). To this end, we calibrated ten different models, each of which was calibrated using the data from nine regions, and validated based on the predictions made to the 10th, left-out region.

Following this cross-validation procedure, we then used our model to predict the plot-level maximum temperature buffering in each of our plots across the 56 regions. Because both the canopy cover data and the macroclimate temperatures were available for the baseline and the resurvey years, we were able to predict subcanopy microclimate at both time steps, and thus calculate the rates of both micro-and macroclimate change (i.e., °C per decade) between the two surveys. This implies that the model used for predicting microclimate change rates was calibrated and validated with data collected after the vegetation plots were surveyed. We found that the microclimate warmed most in areas with decreasing canopy cover (Supplementary Fig. S5). It should be noted that including distance to the coast did not affect the temporal dynamics of microclimate in a plot, because distance to the coast is constant between the baseline survey and the resurvey for a given plot. Therefore, all the changes in microclimate reported here are due to changes in canopy cover and the associated changes in macroclimate temperature buffering. The data and methods used to develop the predictive model for maximum temperature are explained in detail in Zellweger et al. (21). Our focus lies on maximum temperatures dynamics during summer (i.e., the growing season), because they are most affected by canopy cover in deciduous forests and are among the most important climatic determinants of plant growth and survival (9,30,43). Dynamics in other microclimate temperature parameters, such as mean and minimum temperatures, have not been addressed here because previous results have shown that minimum temperatures are less affected by forest attributes such as canopy cover and tree densities, and more by landscape properties (e.g. distance to the coast and topography) (21). These landscape properties do not change over the time span considered here and thus we expect that the small-scale temporal change of maximum temperatures to be much larger than minimum temperatures. However, in order to arrive at a more complete understanding of the role of microclimate in driving the response of biodiversity and ecosystems to climate change, further studies should address the effects of a full range of ecologically relevant microclimate temperature parameters, including soil, air and surface temperatures (26,44,45). Promising opportunities in this respect arise from latest developments in process-based microclimate modelling (46), which is assumed to be more generalizable over time, an issue we were not able to test within our statistical modelling framework.

Our approach of calculating the rates of microclimate change in forests assumes a linearly occurring change in canopy cover and temperature buffering over time. This will indeed apply to the large majority of forests considered here, as most of our plots have not been affected by heavy disturbances (heavily managed forests, e.g. plantations, have been excluded in our data). Only 45 or 1.5 % of our 2955 plots have been subject to a canopy cover decrease by 100 % or more. However, we acknowledge that sudden losses of canopy cover due to e.g. windthrow will have triggered instantaneous losses of temperature buffering, potentially leading to a non-linear increase of warming rates but such effects are nearly impossible to capture at high spatiotemporal resolution. Whether canopy cover change was due to anthropogenic or natural causes did not affect our calculation of the microclimate change rates, because the relationship between canopy cover and maximum temperature buffering is assumed to remain the same between the two causes of canopy cover change.

Thermophilisation rates

To calculate the thermophilisation rate in the understorey layer of each plot, we adopted the data and methods used by De Frenne et al. (14). First, we extracted gridded (20 km by 20 km) data on monthly averages of maximum temperatures during April to September, based on long-term climatic averages from WorldClim (version 2) (47). Next, we estimated species' thermal response curves by random sampling with replacement of 1000 grid cells within the current distribution range of the species (see Supplementary Fig. S6). We used 20 x 20 km² resolution climate data because the continental-scale species range maps are not available at a higher level of spatial detail. The distribution maps were digitized from two renowned atlases (48,49) and were available for 993 out of 1292 species in our dataset. A total of 133 species, out of the 1292 studied species, were not identified to the species level, but recorded as grouped or aggregated species (e.g., Rosa sp.). After removing the grouped species from the 1292 species, the 993 species for which we had distribution maps thus accounted for 85.68 % of the species, which corresponded to a share of 98.39 % of the total cover of the vegetation surveys. These thermal response curves are interpreted as species-specific thermal affinities and were used to test for directional shifts in community composition over time as a response to increasing temperatures. We calculated the floristic temperature (in degrees Celsius) for each plot by sampling the thermal response curves of all species that were present at the time of each survey. To account for variability and uncertainty in species' thermal preferences and niche widths, the distribution of plot-level floristic temperatures at each survey was constructed by resampling 1,000 times from the species' thermal response curves, as explained in Rodríguez-Sánchez et al. (50). In line with our microclimate data representing changes in maximum temperatures during summer, we were specifically interested in analysing shifts in warm-affinity species due to species-specific affinities to heat. Shifts in relative dominance in warm-affinity species within the community are expected to be best represented in shifts at the upper end (right tail) of the floristic temperature distribution (14). Therefore, we quantified the thermophilisation rate as the difference between the 95 th percentiles of the floristic temperature distribution between the baseline survey and resurvey (∆Tplant), divided by the time interval (∆t, in decades) between the two surveys (see Sensitivity analyses section for more details).

We expected the time interval between the baseline surveys and resurveys to be long enough to allow for species replacements due to range shifts. Therefore, we calculated all thermophilisation rates based on changes in species presence or absence, as opposed to changes in species cover (abundance). In 12 out of all baseline surveys and resurveys, no vascular plant species occurred in the understorey-layer, which prevented the calculation of a floristic temperature and resulted in the exclusion of 12 plots from the 2955 plots, leading to a sample size of n = 2943 plots for the thermophilisation analyses.

Climatic debt

To calculate the climatic debt accumulated by a understorey plant community, i.e. the degree to which the rates of compositional changes in vegetation are lagging behind the rates of micro-and macroclimate warming, we subtracted the thermophilisation rates (i.e. ∆Tplant/∆t) from the rates of micro-and macroclimate temperature change (12,20) (i.e., macroclimatic debt: (∆Tmacro/∆t) -(∆Tplant/∆t); vs. microclimatic debt: (∆Tmicro/∆t) -(∆Tplant/∆t). Large differences thus indicate a large disequilibrium between vegetation dynamics and climate change (51), indicating that a relatively high number of species are occurring under sub-optimal thermal conditions. In line with our previously applied statistical framework, we tested whether our findings presented in the main text change when we consider (1) thermophilisation rates calculated using the shifts in the 50 th and 5 th percentiles of the floristic temperature distributions rather than the 95 th percentiles, and (2) when we exclude rare species. We found that none of these two approaches changed our main findings and conclusions (Supplementary Figs. S8 andS9 & Sensitivity analyses section).

Statistical analyses

To report summary statistics of the variables we analysed, i.e. the rates of micro-and macroclimate change, the thermophilisation rate, as well as the micro-and macroclimatic debts, we applied a multilevel-modelling framework with group-level effects for 'region' to account for the spatial non-independence of plots within the same region, using linear mixed-effects models (LMMs). For each of the above variables we calculated a separate LMM, without fixed predictor variables but using 'region' as a random intercept term. The intercept and standard error of these maximum temperatures are more closely related to shifts at the warm end of the floristic temperature distribution, i.e., that community shifts have been caused by differences in species-specific heat affinities and tolerances.

To test whether the presence of rare species affected our findings and conclusions, we re-run the analyses using a dataset where all rare plant species with a cover abundance of less than 1 % at the plot level were excluded. Doing so led to a reduction of the sample size for the thermophilisation analysis, because in 109 surveys out of all baseline surveys or resurveys, no vascular plant species in the herb-layer occurred with a cover of 1 % or more. This prevented the calculation of a floristic temperature and resulted in the exclusion of 109 plots (3.7%) from the 2955 plots, leading to a sample size of n = 2846 plot for the thermophilisation analysis. The resulting thermophilisation rates as calculated based on shifts in the 95 th percentile of floristic temperatures ranged from -1.03 to 1.54 (mean: 0.03 +/-0.01 (SEM) °C per decade, P = 0.015). Relating the thermophilisation rate to the rates of micro-and macroclimate change revealed that thermophilisation was primarily linked to the rate of microclimate warming, and not to the rate of macroclimate warming (scaled slope estimate for microclimate change: 0.036, CI: 0.025 -0.048, P < 0.001; scaled slope estimate for macroclimate change: 0.017, CI: -0.002 -0.037, P = 0.076) (Supplementary Fig. S9). Thus, excluding rare species did not change our main findings and reinforced our conclusions.

It is worth noting that we found a large variation in the thermophilisation rate, which can also be caused by factors other than sub-canopy temperature. For example, plant species distribution may be constrained by biotic interactions, soil moisture or soil acidity or nutrient concentrations, all of which can affect community responses to climate change, and thus also the thermophilisation rate (11).

Variation in community composition over time resulting in thermophilisation could be confounded by changes in light availability brought about by changes in canopy cover. To test for potentially confounding effects of light availability we calculated the plot-level changes in canopy cover between the baseline survey and the resurvey as response ratios, i.e., log(baseline cover/resurvey cover), following the approach in De Frenne et al. (14). We then calculated a univariate LMM with canopy cover response ratios as a fixed effect variable and 'region' as a random intercept term. This model showed that thermophilisation increases with decreasing canopy cover response ratios, i.e. increasing light availability (slope: -0.01, P = 0.051) (Supplementary Fig. S10, see also De Frenne et al. ( 14)), a finding that we attribute to the correlation between the rate of microclimate change and changes in light availability, as indicated by an increase in the rate of microclimate change (warming) with decreasing canopy cover response ratios (slope: -0.24, R-squared: 0.21, P < 0.001). Based on niche theory we expected that microclimate change be more directly related to thermophilisation than light availability, because microclimate temperature is directly related to species thermal niche requirements. In support of this assumption we found that our microclimate model explains thermophilisation rates better than the light availability model, as indicated by AIC values of -2478 and -2465 (P < 0.001), respectively. Moreover, the rate of microclimate change was significantly and positively related to the residuals of the light availability model, suggesting that our microclimate temperature change data holds important additional information for explaining the thermophilisation rate (slope: 0.02, P = 0.049). The combination of the rate of microclimate change and light availability did not significantly improve the model including the rate of microclimate change alone (difference in AIC values < 2, P > 0.05). Adding the interaction term between the rate of microclimate change and light availability improved the model, as indicated by AIC values of -2483 and -2478 (P = 0.016), respectively, suggesting that the effect of microclimate warming on thermophilisation increased with increasing light availability, a finding that corroborates previous experimental findings (59). Taken together, these results suggest that our main findings and conclusions are not confounded by the effects of changes in light availability per se.

We further tested if the exclusion of the coarser-gridded CRU data (i.e. analysing plot-level changes with complete TerraClimate data coverage only) affected our results, i.e. the relationship between the rates of micro-and macroclimate change and the thermophilisation rate, and found that it did not: the thermophilisation rate was only statistically related to the rate of microclimate change (scaled slope estimate: 0.02, 95 th confidence interval (CI): 0.01 -0.03, P < 0.001), and not to the rate of macroclimate change (scaled slope estimate: -0.01, CI: -0.03 -0.00, P = 0.12).

Fig. 1 .

 1 Fig. 1. Forest microclimate change following canopy cover changes over time is considerably more variable than macroclimate change. (A), The distribution of the 2955 resurveyed forest plots (black dots) in 56 regions(purple circles, scaled to the number of plots as indicated in the top right) across the temperate forest biome (green area) in Europe. We representatively sampled microclimate temperature in 100 forest stands, i.e., in 10 stands in each of 10 regions (black circles; effective n = 96), to estimate the maximum (macroclimate) temperature buffering during the growing season as a function of canopy cover (23) (Fig.S2). (B), Schematic overview of the method used to approximate microclimate change in the forest understorey. In this example, canopy cover at the time of the baseline survey was higher than during the resurvey, resulting in a decrease in macroclimate temperature buffering from 2 to 1 °C, which in turn led to a relatively larger increase in microclimate warming (20 to 23 °C) compared to macroclimate warming (22 to 24 °C). The relationship between canopy cover and the buffering of maximum macroclimate temperature was empirically assessed across the study area (Fig.S2)(21). (C), The rate of macroclimate change is plotted against the rate of microclimate change, with the black bisecting line representing the 1:1 relationship. Micro-and macroclimate have both significantly warmed (see text for statistical results). The distributions of values in the rates of micro-and macroclimate change are indicated by grey rugs on each axis. Microclimate change rates are 45 % more variable than macroclimate change rates and macroclimate change rates only accounted for about half of the variation in microclimate change rates, as indicated by the marginal (conditional) R-squared value of 0.52 (0.69). All statistical results are based on mixed-effects models with region as a random-effect (intercept) term.

Fig. 2 .

 2 Fig. 2. Thermophilistaion in forest understorey plant communities is related to microclimate change, not to macroclimate change. Thermophilisation rates increase with increasing microclimate warming of maximum temperatures during the growing season (Tmax), as shown by the regression slope and 95 % confidence intervals for microclimate. The thermophilisation rate was not statistically related to the rate of macroclimate warming (see text for statistical results).

Fig. 3 .

 3 Fig. 3. Temperature buffering by canopy cover explains the climatic debt in forest plant communities. (A), The climatic debt calculated based on microclimate temperature change (i.e., microclimatic debt) increases with decreasing maximum temperature (Tmax) buffering following a reduction of canopy cover (slope: -1.88, marginal (conditional) R-squared: 0.51 (0.75), P < 0.001). Negative values on the x-axis represent a warming effect (reduced canopy cover and thus less Tmax buffering), positive values represent a cooling effect (increased canopy cover and thus more Tmax buffering). (B), Climatic debts calculated using macroclimate temperature change (i.e., macroclimate debt) are only weakly related to differences in temperature buffering (slope: 0.13, marginal (conditional) R-squared: 0 (0.25), P: 0.06). The linear regression lines are plotted including the 95 % confidence intervals (grey bands). The green dots indicate exemplified micro-and macroclimate debts following the illustrated change (reduction) in Tmax buffering in Fig. 1B.

Fig. S1 .

 S1 Fig. S1. Frequency distribution of forest plot survey years. Vegetation surveys in each plot were carried out twice between the years 1934 and 2017, with a median (range) time difference between the baseline surveys and resurveys of 38 (12 -66) years. 5

Fig. S2 .

 S2 Fig.S2. Macroclimate maximum temperatures (Tmax) are buffered by canopy cover. To quantify the degree to which the canopy cover buffers free-air (macroclimate) maximum temperatures during summer (i.e., the growing season; June to August), we compared ground-truth microclimate sensor measurements in 96 forest stands to daily maximum temperatures measured at the weather station that was the closest to each stand. The temperature 5buffering values represent the daily microclimate maximum temperature minus the daily macroclimate maximum temperature, averaged during the growing season. This analysis revealed a positive, non-linear relationship between temperature buffering and canopy cover, as indicated by the smoothed curve with 95 % confidence intervals (light red polygons), with an R-squared of 0.26. We then used this relationship in combination with the distance to coast to predict the temperature buffering for each plot and survey year. Adding topographic variables such as elevation, 10 topographic position or slope did not improve the predictive accuracy of the model (see Methods).

Fig. S3 .

 S3 Fig. S3. Frequency distribution of temporal canopy cover changes in temperate European forests. Canopy cover changes were highly variable with a mean not significantly different from zero (2.6 %; P = 0.426). The plotlevel differences in canopy cover over time was calculated by subtracting the canopy cover observed in the baseline 5 survey from the canopy cover observed in the resurvey; mean and s.e.m. are based on mixed-effects models with region as a random-effect term.Difference in canopy cover (%) over time

Fig. S5 .

 S5 Fig. S5. The microclimate warmed most in areas with decreasing canopy cover. Microclimate change was negatively correlated with canopy cover differences (slope: -0.005, R-squared: 0.33, P < 0.001). The canopy cover difference was calculated by subtracting the canopy cover observed in the baseline survey from the canopy cover observed in the resurvey. Positive values on the x-axis thus indicate an increase in canopy cover over time, and vice-5 versa.

Fig. S7 .

 S7 Fig. S7. The mean micro-and macroclimatic temperature debts are similar, but climatic debts calculated using macroclimate data underrepresent the variability in microclimatic debt. (A), Histograms of microclimatic and macroclimatic debts, calculated by subtracting the thermophilisation rate from the respective micro-and macroclimate temperature change rate. The resulting microclimatic debt ranged from -0.86 to 2.03 5 (mean: 0.38 +/-0.04 (SEM)) °C per decade; the range for macroclimatic debt was -0.67 to 1.72 (mean: 0.37 +/-0.04 (SEM)) °C per decade. The two dashed vertical lines represent the means of the distributions. (B), Scatterplot between microclimatic and macroclimatic debt; macroclimatic debt explained only 56 % of the variation in microclimate debt (slope = 0.92, marginal (conditional) R-squared = 0.56 (0.72), P < 0.001). 10

Fig. S8 .

 S8 Fig. S8. Results derived from shifts of the median of the community-level floristic temperature distributions.(A), The results are presented analogous to Fig.2in the main text, but representing thermophilisation rates calculated using the shifts in the 50 th percentile of the community-level floristic temperature distributions rather than the shifts in the 95 th percentile as used in the main text (see Methods). Thermophilisation rates increase with increasing microclimate warming, as shown by the regression slopes and 95 % confidence intervals (CI) (scaled slope estimate: 0.006, CI: 0 -0.011, P = 0.036). The macroclimate model, however, was not statistically significant, with a scaled slope estimate of -0.001 (CI: -0.008 -0.007, P = 0.850). (B), Climatic debt effects based on thermophilisation rates representing shifts in the 50 th percentiles of the community-level floristic temperature distributions (see Methods), analogous to Fig.3in the main text. The climatic debt based on microclimate data (i.e, microclimatic debt) is negatively related to the difference in maximum temperature (Tmax) buffering brought about by canopy cover dynamics (slope: -2.393, marginal (conditional) R-squared: 0.581 (0.869), P < 0.001). Negative values on the x-axis represent a warming effect (reduced canopy cover and thus less Tmax buffering), positive values represent a cooling effect (increased canopy cover and thus more Tmax buffering). (C), Climatic debt calculated using macroclimate data (i.e., macroclimatic debt) are not significantly related to differences in temperature buffering (slope: -0.065, marginal (conditional) R-squared: 0 (0.242), P = 0.489). Thus, using the thermophilisation rates based on shifts in the 50 th percentiles of the floristic temperature distributions reveal the same conclusion as reported in the main text: The climatic debt in forest plant communities depend on microclimate dynamics mediated by canopy dynamics.

Fig. S9 .

 S9 Fig.S9. Rare species did not affect our main results and conclusions. (A), Relationships between thermophilisation rates and micro-, as well as macroclimate change, based on the dataset excluding rare species with a cover abundance of less than 1 % at the plot level, as shown by the regression slopes and 95 % confidence intervals (CI). The thermophilisation rate was primarily linked to the rate of microclimate warming, and not to the rate of macroclimate warming: the effect of microclimate change on thermophilisation was stronger (scaled slope estimate: 0.036, CI: 0.025 -0.048, P < 0.001) than the effect of macroclimate change (scaled slope estimate: 0.017, CI: -0.002 -0.037, P = 0.076). (B), The climatic debt based on microclimate data (i.e., microclimatic debt) is negatively related to the difference in maximum temperature (Tmax) buffering brought about by canopy cover dynamics (slope: -1.358, marginal (conditional) R-squared = 0.379 (0.586), P < 0.001). Negative values on the xaxis represent a warming effect (reduced canopy cover and thus less Tmax buffering), positive values represent a cooling effect (increased canopy cover and thus more Tmax buffering). (C), The climatic debt calculated using macroclimate data (i.e., macroclimatic debt) was only weakly related to differences in temperature buffering (slope = 0.197, marginal (conditional) R-squared = 0.007 (0.256), P < 0.001). The linear regression lines are plotted including the 95 % confidence intervals (grey bands). Rare species did thus not affect our conclusions reported in the main text: The climatic debt in forest plant communities depend on microclimate dynamics mediated by canopy dynamics. Note that the sample size for the analysis and results presented here was n = 2846. The reason for the reduced sample size is that in 109 surveys out of all baseline surveys or resurveys, no vascular plant species in the herb-layer occurred with a cover of 1 % or more. This prevented the calculation of a floristic temperature and resulted in the exclusion of 109 plots (3.7%) from the 2955 plots.

Fig. S10 .

 S10 Fig. S10. Thermophilisation tends to be higher in forests where light availability has increased following canopy opening (slope: -0.01, P = 0.051). Please refer to the Sensitivity analyses section above for more details.

www.naturalearthdata.com, last accessed 18 th February

Acknowledgments Funding: F.Z. received funding from the Swiss National Science Foundation (project 172198) and the Isaac Newton Trust; P.D.F. and P.V. from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (ERC Starting Grant FORMICA 757833); K.V. from ERC Consolidator Grant PASTFORWARD (614839). The Czech Science Foundation and Czech Academy of Sciences funded M.C., O.V. and R.H. (17-09283S, RVO 67985939), and M.M. and M.K. (17-13998S, RVO 67985939), and P.P. (RVO 67985939). T.N. received funding from the Slovenian Research Agency (J4-1765); T.B. from the Ministry for Innovation and Technology in Hungary; B.I. was funded by Agrárklíma II. (VKSZ_12-1-2013-0034); R.P. from the Polish Ministry of Science and Higher Education (DS-3421); F.M. from the Slovak government (APVV-15-0270).

Author contributions: F.Z., P.D.F. and D.C. conceived and designed research. F.Z. analyzed the data with help from P.D.F. P.V. and P.D.F. calculated the thermal response curves. J.L. contributed data and valuable methodological comments. All other authors contributed data. F.Z. wrote the manuscript together with P.D.F. and D.C. All other authors commented, edited and approved the manuscript.

Competing interests: Authors declare no competing interests.

Data and materials availability: Data and code used in this study are available in the online DRYAD repository with the identifier doi:10.5061/dryad.r7sqv9s83.

Supplementary Materials:

Materials and Methods

Figures S1-S10

References 35 -59 are only called out in the SM Supplementary Materials for models represent the mean and standard error of the mean (SEM) of each variable. When fitting our intercept-only LMMs, we used the restricted maximum likelihood method in the lmer function from the lme4 package (52) as recommended by Zuur et al. (53). We used χ 2 tests to determine whether the means were significantly different from zero. The mean and SEM of each variable were then complemented by the range and quantiles drawn from the raw data distribution.

Next, we assessed how the rates of macroclimate change predict variation in the rates of microclimate change. We therefore fitted a LMM with the rate of macroclimate change as a fixed effect variable and 'region' as a random intercept term variable using restricted maximum likelihood in the lmer function from the lme4 package (52). Moreover, χ 2 tests were performed by comparing the univariate LMM including the rate of macroclimate change with the intercept-only model (53). Goodness-of-fit was determined by calculating the marginal coefficient of determination (R-squared) as proposed by Nakagawa & Schielzeth (54) using the r.squaredGLMM function in the MuMIn package (55). The marginal R-squared describes the variation explained by the fixed factor only, as opposed to the conditional R-squared, which describes the variation explained by the fixed and random factors together (54).

The same approach was chosen to analyse the relationship between climate change (i.e., the rates of microand macroclimate change) and the thermophilisation rate. We fitted an individual LMM with either the rate of micro-or macroclimate change as a fixed effect variable and 'region' as a random intercept term, and thermophilisation rate as the response variable. Prior to model fitting, we centered and standardised (scaled) the rates of micro-and macroclimate change (mean = 0, standard deviation = 1), so that the scaled slope estimates can be interpreted as relative effect sizes (56). The rates of microclimate change were significantly and positively related to the residuals of the macroclimate change model (slope: 0.025; P < 0.01), suggesting that microclimate data holds complementary information compared to macroclimate data for explaining biotic responses to climate change. We also tested the linearity of the relationship between the rate of micro-or macroclimate change and the thermophilisation rate, by running general additive mixed-effects models using the gamm function in the mgcv package (57) and 'region' as a random term. However, we did not find any evidence of a non-linear relationship (P = 0.60). Lastly, we fitted an individual LMM with the change in Tmax buffering as a fixed effect variable and 'region' as a random intercept term, and either the micro-or macroclimatic debts as the response variables. All analyses were performed in R version 3.5.0 (58).

Sensitivity analyses

We performed a thorough sensititvity analysis to test: (1) for effects from different approaches to calculate the thermophilisation rate; (2) if our main results and conclusions are affected by the presence of rare species; (3) whether our results are confounded by changes in light availability brought about by changes in canopy cover, and (4) whether the use of coarser-gridded CRU climate data affected our results.

As stated above, we were particularly interested in testing for shifts in the relative dominance in warmaffinity species and therefore quantified the thermophilisation rate as the difference (shifts) between the 95 th percentiles of the floristic temperature distribution between the baseline survey and the resurvey (14). Alternatively, thermophilisation rates are often derived from shifts in the 50 th percentiles of floristic temperatures, which indicate a shift in median floristic temperatures, which appears most adequate when related to shifts in mean temperatures. To test if using the median floristic temperatures to calculate thermophilisation would have resulted in different findings and conclusions, we re-run all major analysis using the 50 th percentile (i.e., median) of floristic temperatures. The resulting thermophilisation rates ranged from -0.74 to 0.92 (mean: 0.02 +/-0.005 (SEM) °C per decade, P < 0.001), indicating that thermophilisation has occurred in our plots. The relationships between the thermophilisation rate and the rate of micro-as well as macroclimate change were similar but statistically weaker than the relationships reported in the main text (Supplementary Fig. S8). Moreover, thermophilisation may be caused by shifts in relative abundances of cold-affinity species, which in turn should be best represented in shifts in the 5 th percentile of floristic temperatures (14). Re-running the analyses using shifts in the 5 th percentiles revealed thermophilisation rates ranging from -1.42 to 1.48 (mean: 0.04 +/-0.01 (SEM) °C per decade, P < 0.001). This finding indicates that thermophilisation may result from community filtering effects at cold range edges of species distributions. Testing how temperature changes may have triggered shifts in cold-affinity species provides an interesting avenue requiring data on floristic temperature distributions calculated from minimum temperatures, as well as actual changes in minimum temperatures. Here, we focus on maximum temperatures and indeed, we found that the thermophilisation rate based on the shifts in 5 th percentile was neither related to the rate of microclimate change (scaled slope estimate: 0.00, CI: -0.01 -0.01, P = 0.35), nor to the rate of macroclimate change (scaled slope estimate: 0.00, CI: -0.01 -0.02, P = 0.67). In summary, these additional analyses support our assumption that changes in microclimate