Benoit Richard

Jean-Luc Dupouey

Emmanuel Corcket

Didier Alard

Frédéric Archaux

Michaël Aubert

Vincent Boulanger

François Gillet

Estelle Langlois

Sébastien Macé

| Didier Alard

| Frédéric Archaux

Pierre Montpied

Thérèse Beaufils

Carole Begeot

| Patrick Behr

Jean-Michel Boissier

Sylvaine Camaret

| Richard

Guillaume Decocq

Yann Dumas

Richard Eynard-Machet

Jean- Claude Gégout

Sandrine Huet

| Valéry Malécot

| Pierre Margerie

Arnaud Mouly

Thierry Paul

Benoît Renaux

Pascale Ruffaldi

Fabien Spicher

| Erwin Thirion

Erwin Ulrich

Manuel Nicolas

Jonathan Lenoir
email: jonathan.lenoir@u-picardie.fr

The climatic debt is growing in the understorey of temperate forests: Stand characteristics matter

Keywords: 43 Biodiversity, climate change, climatic debt, anthropogenic disturbances, long-term 44 monitoring, plant community, thermophilization

come   L'archive ouverte pluridisciplinaire

 , 2018). Some previous studies even relied on non-permanent or quasi-91 permanent plots to assess changes in CTI values over time [START_REF] Bertrand | Changes in plant community composition lag behind climate 608 warming in lowland forests[END_REF]De 92 Frenne et al., 2013;[START_REF] Becker-Scarpitta | Four decades of plant 601 community change along a continental gradient of warming[END_REF]Zellweger et al., 2020).

93 Unfortunately, such studies relying on non-permanent plots tend to infer temporal 94 changes blurred by space-for-time substitutions, which are known to overestimate 95 the magnitude of thermophilization of communities in response to contemporary 96 climate warming [START_REF] Elmendorf | Experiment, monitoring, and gradient methods 668 used to infer climate change effects on plant communities yield consistent patterns[END_REF]. Even when studies rely on strictly 97 permanent plots, it requires more than one or two resurveys to build a detailed time 98 series of community composition changes from long-term monitoring programs.

99 Strictly permanent monitoring programs relying on frequent resurveys (e.g. from 100 one-year to five-year intervals) during a sufficiently long time period (e.g. at least 101 two decades) are needed to capture a temporal dynamic in the velocity of community 102 thermophilization and confirm linearity over time. Although such long-term time 103 series of biodiversity monitoring already exist (Dornelas et al., 2018), we still lack a 104 good knowledge and understanding of the temporal dynamic of the climatic debt.

105

Whether the climatic debt of plant communities is steadily growing over time

106 and what are the main underlying determinants of this dynamic remain open 107 questions in global change ecology [START_REF] Essl | Delayed biodiversity change: no time to waste[END_REF][START_REF] Bertrand | Ecological constraints increase the climatic debt in forests[END_REF]Zellweger 108 et al., 2020). For instance, a previous study focusing on temperate forests has 109 suggested an increasing climatic debt in the response of understory plant 110 communities from the lowlands [START_REF] Bertrand | Changes in plant community composition lag behind climate 608 warming in lowland forests[END_REF]. A follow-up study has 111 revealed that this growing climatic debt involved a complex combination of abiotic 114 occurring within initially warmer conditions) [START_REF] Bertrand | Ecological constraints increase the climatic debt in forests[END_REF]. Noteworthy, 115 the same authors also found that anthropogenic disturbances contributed very little 116 and unexpectedly to the climatic debt. Yet, previous and more recent findings 117 suggest the opposite with local disturbances amplifying forest community responses 118 to climate change [START_REF] De Frenne | Microclimate moderates plant responses to 644 macroclimate warming[END_REF][START_REF] Stevens | Forest disturbance 758 accelerates thermophilization of understory plant communities[END_REF][START_REF] Perring | Global environmental change effects on plant 744 community composition trajectories depend upon management legacies[END_REF]119 Brice et al., 2019). A more recent study has even demonstrated that the climatic lag 120 in forest plant communities is primarily controlled by forest microclimate dynamics 121 suggesting forest stand characteristics to be a strong determinant of the climatic 122 debt (Zellweger et al., 2020). These discrepancies call for more in-depth 123 investigations on the temporal dynamic of the climatic debt in forest understory 124 plant communities and the relative contribution of anthropogenic and natural 125 disturbances.

126

In addition to baseline environmental conditions and stand characteristics 127 modulated by anthropogenic and natural disturbances, ungulates directly impact 128 vegetation dynamics [START_REF] Post | Opposing plant community responses to warming 749 with and without herbivores[END_REF][START_REF] Olofsson | Herbivores inhibit climate-driven shrub expansion on the tundra[END_REF]Boulanger et al., 129 2018). On the one hand, [START_REF] Olofsson | Herbivores inhibit climate-driven shrub expansion on the tundra[END_REF] 150

Materials and Methods

A long-term national monitoring network of forest ecosystems

152 To assess changes in understory plant communities, we used vegetation data 153 collected every five years, and even every year for some sites, throughout a long- [START_REF] Ulrich | Le réseau RENECOFOR: objectifs et réalisation[END_REF][START_REF] Nicolas | How monitoring networks contribute to 732 the understanding and to the management of soil and forest ecosystems? Revue 733 forestière française[END_REF]. All plots 162 are located in mature and even-aged pure stands, and cover an area of 2 ha with a 163 central fenced zone of 0.5 ha to exclude large herbivores (Fig. 1) (for more details, 164 see Boulanger et al., 2018).

165

Since 1995, vegetation surveys have been conducted every five years (baseline 166 survey in 1995 and resurveys in 2000, 2005, 2010 and 2015). From the total pool of 167 102 permanent plots, we discarded 58 plots which suffered from major data gaps 168 throughout the studied time period (e.g. the resurvey of 2010 was systematically 169 missing for 53 plots due to a shortage in funding). This yielded a total of 44 plots for 170 which a full time series, at a 5-yr interval, was available between 1995 and 2015, 171 hereafter referred to as the five-year dataset (Fig. 1). Among these 44 plots, 8 plots 172 were monitored annually without discontinuity between 1995 and 2017, hereafter 173 referred to as the yearly dataset (Fig. 1). The whole workflow of analyses was 174 carried out separately for both the five-year and yearly datasets, as these two 175 datasets are very complementary (Supporting Information App. Fig. 1). 1.

182

During each monitoring campaign, the vegetation was surveyed by expert 183 botanists according to standardized protocols under a quality assurance and quality 184 control procedure. The quality assurance and quality control procedure involved 185 intercalibration exercises organized the year before each monitoring campaign with 186 all expert botanists surveying the same plots in the same forest (Camaret & 187 Brêthes, 2004;[START_REF] Archaux | Can we reliably estimate species richness with large plots? an 595 assessment through calibration training[END_REF]. During each monitoring campaign, all plots 188 were surveyed both at spring and in summer/autumn, to cover the entire growing 189 season. All vascular plants were recorded in eight 2 m 50 m subplots: four inside × 190 and four outside the fenced area of 0.5 ha (Fig. 1). The cover-abundance of all taxa 191 occurring in four vegetation layers (tree: > 7 m; high shrub: 2-7 m; low shrub: 0.3-2 192 m; herb: < 0.3 m) was recorded. Species names were standardized following the 193 nomenclature of Flora Europaea [START_REF] Tutin | Flora Europaea 5 Volume Set and CD-ROM Pack[END_REF]. A total of 877 taxa were 228 from -0.56 to 1.54 over the period 1992-2017), except for three years (1993, 1996 and 229 2010), for which anomalies were negative.

230

To infer temperature conditions from local species assemblages, we used the 231 community temperature index (CTI) [START_REF] Devictor | Birds are tracking climate 650 warming, but not fast enough[END_REF][START_REF] Devictor | Differences in the climatic debts of birds and butterflies at a 654 continental scale[END_REF] (see Equation 1). 𝐼𝑉 𝑇𝑒𝑚𝑝𝑒𝑟𝑎𝑡𝑢𝑟𝑒 𝑆𝑝 𝑖 𝑆 240 Species' IVs for temperature were extracted from both the CATMINAT (Julve, 1998) 241 and EcoPlant (Gégout et al., 2005) databases, and we used them in combination to 242 take advantage of their respective strengths (Supporting Information App. Fig. 4).

F

Statistical analyses

309 We used linear mixed-effects models (LMMs) to assess the velocity of MAT changes 310 and thermophilization of understory plant communities. For CTI values, we used 311 the continuous variable "year" of the monitoring as the main fixed effect while 312 accounting for potentially diverging temporal dynamics between the outside and 313 inside of the fenced area, by adding the interaction effect: year exclosure. The × 314 RENECOFOR plot ID (44 or 8 levels depending on the dataset used) and the 315 vegetation subplot ID (8 levels), nested within the plot ID, were used as a random 316 term to account for the hierarchical structure of the sampling design and the 317 repeated surveys over time. For MAT values, we only used the "year" of the 318 monitoring as a fixed effect and the RENECOFOR plot ID as a random term.

319 Indeed, for MAT, the subplot level was not used in the random component of LMMs 320 since all 8 subplots within a given plot have similar MAT values. To allow between-321 site variation in the temporal dynamics of MAT and CTI values, we set the random 322 term (e.g. plot ID) of our LMMs so that it modifies the slope coefficient (i.e. random 323 slope term) of the variable "year" for each level of the random term variable. The 324 coefficient estimate of the "year" variable in our LMMs was used as a proxy for 325 estimating the velocities of temperature changes, both for MAT and CTI values.

326 Assumptions of normality and homoscedasticity of the residuals of our models were 327 visually checked and no departure from normality and homoscedasticity were 328 detected. Significance of the predictors in LMMs were assessed using Wald tests, 𝜒 2 329 after fitting models with the restricted maximum likelihood estimation. To analyze the temporal dynamic of the difference between MAT and CTI 331 values over time (i.e. the climatic debt), we used LMMs with the value of the 332 difference (MAT CTI) as the response variable and the continuous variable "year" -333 of monitoring as the main fixed effect variable. In addition to "year", we added 334 several abiotic and biotic covariates as fixed effect variables to assess their relative 335 contribution in modulating the climatic debt. These covariates were: baseline 336 temperature conditions; sum of natural disturbances; sum of anthropogenic 337 disturbances; disturbance anteriority; stand basal area; stand age; type of dominant 338 tree species in forest stand (deciduous vs. coniferous); and the exclosure effect. From 339 the full model including all covariates together with time as fixed effect terms, we 340 ran all possible combinations of models and ranked all candidate models using the 341 corrected Akaike information criteria (AICc). Based on the best candidate model 342 with the lowest AICc value, we then computed the difference in AICc values ( AICc) 𝛥 343 for all models relative to the best candidate model. Finally, from the set of selected 344 candidate models with AICc values < 2, we used a model averaging approach to 𝛥 345 calculate the conditional mean value of each coefficient and its associated 95%

346 confidence interval (Burnham & Anderson, 2002). All continuous predictor variables 347 were standardized (z-score) before analyses to improve their interpretability on a 348 comparable scale [START_REF] Schielzeth | Simple means to improve the interpretability of regression 756 coefficients[END_REF]. None of the fixed-effect variables included in 349 our models showed multicollinearity issues (VIF < 3;Zuur et al., 2007). Note that 350 correlation values among all our explanatory variables did not exceed a Spearman's 351 r coefficient of 0.6 (Supporting Information App. Fig. 8). Again, we used the 352 RENECOFOR plot ID (44 or 8 levels depending on the dataset used) and the Irrespective of the dataset considered in our analyses, we found very 395 consistent and similar trends in the factors contributing to the climatic debt (Fig. 3).

396 For instance, the year of sampling had a positive effect on the magnitude of the 397 difference, confirming an increasing climatic debt of understory plant communities

411

In contrast, forest disturbance had a negative effect on the climatic debt (Fig. 420 albeit its tendency towards a negative effect did not differ from zero (-0.04 0.04 ± 421 SE, Fig. 3a).

Discussion

423 Velocities of thermophilization 424 Irrespective of the temporal resolution in the data (i.e. the yearly vs. five-year 425 dataset), we show a very similar increase in temperature conditions inferred from 426 understory plant communities over the period 1995-2017, consistent with a 427 thermophilization process. This finding adds to a growing body of evidence 428 describing thermophilization of plant communities across many forest ecosystems 429 around the globe (e.g. in temperate, boreal and tropical ecosystems; Lenoir et al., 430 2010;[START_REF] Bertrand | Changes in plant community composition lag behind climate 608 warming in lowland forests[END_REF][START_REF] De Frenne | Microclimate moderates plant responses to 644 macroclimate warming[END_REF][START_REF] Savage | Elevational shifts, biotic homogenization and time 753 lags in vegetation change during 40 years of climate warming[END_REF]431 Fadrique et al., 2018;Zellweger et al., 2020), as well as across other ecosystems and 432 taxonomic groups (e.g. in mountain plant communities, [START_REF] Gottfried | Continent-wide response of mountain vegetation to 688 climate change[END_REF] 433 marine fishes and invertebrates, [START_REF] Cheung | Signature of ocean warming in 633 global fisheries catch[END_REF]birds and butterflies, Devictor 434 et al., 2012;bumblebee communities, Fourcade et al., 2018) 442 Nevertheless, these rates are still much lower than those found in marine 443 ecosystems, where environmental constraints may weakly limit the dispersal of 444 species in response to rising temperatures [START_REF] Lenoir | Species better track climate warming 720 in the oceans than on land[END_REF].

445

Our findings consist in the most recent updates supporting and consolidating 446 previous reports of plant community thermophilization which happened across the 447 same study area (France) [START_REF] Bertrand | Changes in plant community composition lag behind climate 608 warming in lowland forests[END_REF][START_REF] Martin | Short-term 724 climate-induced change in French plant communities[END_REF] (Bertrand et al., 482 2016). Considering the velocity of the predicted future climate change (Loarie et al., 483 2009), this is an important result because it raises the question whether the lag will 484 continue to steadily accumulate and approach a potential critical breakpoint for 485 understory biodiversity and forest ecosystem functioning. Interestingly, the linear 486 increase of the climatic debt over time is clearer and more pronounced (higher effect 487 size) in the yearly dataset than in the five-year dataset. This is probably because the 488 long-term increasing trend in the climatic debt was better separated from inter-489 annual variations in the yearly dataset than in the five-year dataset. Inter-annual 490 variations tend to play a bigger role in blurring the long-term response of understory 491 plant communities when the frequency between resurveys is lower.

492

Importantly, we demonstrate that forest stand characteristics and 493 anthropogenic disturbances significantly matter in modulating the climatic debt, by 494 either increasing (positive effect) or mitigating (negative effect) the lag. For 495 instance, both the anteriority of the last disturbance event and the cumulative 496 number of anthropogenic disturbances have a negative effect on the climatic debt.

497 Such findings are in accordance with recent works indicating that local disturbances 498 speed up the process of thermophilization in plant communities (Stevens et al., 499 2015;[START_REF] Brice | Disturbances amplify tree 619 community responses to climate change in the temperate-boreal ecotone[END_REF]Zellweger et al., 2020). However, it remains to be 500 investigated in more details whether the thermophilization of understory plant [START_REF] Olofsson | Herbivores inhibit climate-driven shrub expansion on the tundra[END_REF]Boulanger et al., 2018). The indirect effects of herbivores 550 on vegetation responses to climate warming can be highly dependent on the 551 ecological context, as pointed out by [START_REF] Vuorinen | Cool as a moose: How can browsing counteract 769 climate warming effects across boreal forest ecosystems[END_REF].

501

552

Together with previous studies, our results show that the climatic debt in 553 understory plant communities seem to be much more pronounced in situations 554 where temperature conditions were warmer during the period preceding the 555 baseline surveys [START_REF] Bertrand | Ecological constraints increase the climatic debt in forests[END_REF], such as in the lowlands (Bertrand et al., 556 2011). However, the fact that warmer baseline conditions increase the climatic debt 557 of understory plant communities could also stem from a bias towards an (Julve, 1998) and EcoPlant (Gégout et 83 al., 2005). These two databases were chosen, because CATMINAT contains almost 84 all of species recorded in our dataset and EcoPlant, poorer in species, has the 85 advantage to allow the direct computation of CTI values in degrees Celsius.

558
86 CATMINAT (Julve, 1998) is based on expert knowledge similarly to Ellenberg IVs 87 (Ellenberg et al., 1992), with temperature preferences of plant species coded along 88 an ordinal scale from cold (1) to warm (9) environments. While the CATMINAT 89 database contains almost all of the species sampled in the RENECOFOR monitoring 90 program (858 out of the 877 species recorded within all strata of vegetation; App.

91 Fig. 5a), CTI values estimated with these species IVs are unitless and direct 92 comparisons with macroclimatic temperatures are not possible. On the other hand, 93 EcoPlant allows the computation of CTI values in degrees Celsius. EcoPlant is a 94 phyto-ecological database which provides information about ecological optimum for 95 key variables of forest plant species based on their distribution over the French 96 territory. In EcoPlant, species IV for temperature are expressed in degrees Celsius 97 as they are estimated using logistic regression models based on climatic data from 98 AURELHY dataset for the historical period . Nevertheless, the EcoPlant 99 database covers much less species from our dataset (396 out of 877 species).

100 1965-1986and 1987[START_REF] Brice | hence 307 we hypothesized a[END_REF]periods 1965-1986and 1987-2008 Comparison between Comparison between

  factors with the greatest effect size captured by environmental conditions 113 during the baseline surveys (i.e. a greater climatic debt for plant communities

  have shown that herbivores inhibit 130 the climate-driven shrub expansion towards high latitudes, indicating that 131 herbivores could contribute to the climatic debt. On the other hand, Post & Pedersen 132 (2008) have shown that plant community composition (e.g. dwarf birch, willow, 133 graminoids, forbs and mosses) on warmed plots (open-top chambers) did not differ on ambient plots where herbivores grazed, suggesting no effect of 135 herbivory on the climatic debt. Altogether, these findings show contrasting impacts 136 of herbivory on the climatic debt. 137 Here, we rely on a long-term (1995-2017) national (France) monitoring 138 network of permanent forest plots organized in a paired design of exclosure (fenced 139 area excluding ungulates) vs. control plots -the RENECOFOR monitoring program 140 -with multiple resurveys (every year to every five years) to generate a high-141 resolution and unprecedented time series of CTI values for understory plant 142 communities inside and outside the fenced area. Using this information-rich time 143 series, we aimed to: (i) compare the rates of CTI changes to the velocity of MAT 144 changes during the period 1995-2017; (ii) test whether the climatic debt of 145 understory plant communities increased steadily over time; and (iii) assess the 146 relative contribution of several abiotic and biotic factors known to be influential 147 from the scientific literature, including baseline environmental conditions, stand 148 characteristics that are influenced by anthropogenic and natural disturbances, and 149 biotic disturbances from large herbivores.

194

  recorded in all vegetation layers, including 751 occurring in the herb layer 195 (hereafter understory plants). Here, we restricted our analyses to understory plants, 196 including tree seedlings, because they represent a relevant component of forest 197 biodiversity which is not directly (but indirectly via the upper canopy layers) 198 impacted by forest management practices, and are expected to be the most each of the 44 studied plots (by averaging maximum and minimum 222 temperature values for each month), and then computed mean annual temperature 223 (hereafter MAT) for each year between 1958 and 2017. Between 1958 and 2017, 224 MAT increased across all 44 studied plots, leading to a shift from negative to only 225 positive anomalies towards the end of the reference period (1958-2017) (Supporting 226 Information App. Fig. 3). From the very beginning (1992) of the RENECOFOR 227 monitoring program until 2017, anomalies were mostly positive (values ranging

2 382

 2 ID (8 levels), nested within the plot ID, as a random slope term 354 modifying the effect of "year" of monitoring. The explained variance of LMMs was 355 determined by calculating the marginal ( m, associated with fixed effects only) and 𝑅 2 356 conditional ( c, associated with both fixed and random effects) values 𝑅 2 𝑅 2 357 (Nakagawa & Schielzeth, 2013) using the r.squared.GLMM function from the MuMIn 358 package.359All statistical analyses were performed with the R software version 3.6.0 (R 360 Core Team, 2019). Models were performed with the nlme package v.3.1-140 (Pinheiro 361 et al., 2019), and the model selection procedure was assessed with the MuMIn package 362 v.1.43.6 (Bartoń, 2019), using the dredge function. Other packages were used for 363 data handling extraction and graphical visualization: dplyr v.0.8.3 (Wickham et al., 364 2019); raster v.2.9-23 (Hijmans, 2019); ggplot2 v.3.2.0 (Wickham, 2016); cowplot 365 v.1.0.0 (Wilke, 2019); and sf v.0.7-7 (Pebesma, 2018). 366 Results 367 Temporal trends in macroclimate and community temperatures 368 Using the 5-yr interval dataset, mean annual temperature (MAT) during 1995-2015 369 increased by 0.23 0.02 °C/decade (95 % Confidence Intervals (CI): 0.18-0.28, t = ± 370 9.05, P < 0.001; Fig.2a). In the yearly dataset, MAT increased by 0.31 0.05 ± 371 °C/decade (95 % CI: 0.20-0.42, t = 5.51, P < 0.001; Fig.2b) during 1995-2017. Over 372 the same time periods, the velocities of CTI changes were weaker, but similar trends 373 and magnitudes were observed in both the five-year (0.10 0.03 °C/decade, 95 % ± 0.04-0.15, t = 3.62, P < 0.001; Fig.2a) and yearly (0.08 0.03 °C/decade, 95 % ± 375 CI: 0.01-0.14, t = 2.15, P = 0.031; Fig.2b) datasets. No significant differences were 376 observed in CTI trends between inside and outside of the fenced area. Further, the 377 difference between MAT and CTI values showed a tendency to increase steadily over 378 time (Fig.2): the lag in the response of herbaceous plant communities to climate 379 warming increased linearly over time. Although we found a significant effect of time, 380 it contributed weakly to the overall variation in MAT and CTI values: m 0.01 𝑅 2 ≤ 381 while c ranged between 0.92 to 0.96. 𝑅 Factors contributing to the climatic debt 383 For both the yearly and five-year datasets, the model averaging procedure to explain 384 the climatic debt (MAT CTI) selected 4 and 3 candidate models ( AICc < 2), -𝛥 385 respectively. The best candidate models explained 55% and 70% of the total 386 variation in the difference between MAT and CTI values for the five-year ( m = 𝑅 2 387 0.33; c = 0.55) and yearly ( m = 0.46; c = 0.70) datasets, respectively. The best 𝑅 2 𝑅 2 𝑅 2 388 candidate model for the five-year dataset contained 6 variables including the year of 389 sampling, stand age, the anteriority of the last disturbance, the cumulative number 390 of anthropogenic disturbances, stand basal area, and baseline temperature 391 conditions. For the yearly dataset, the best candidate model included the same 6392 variables together with the type of dominant tree species (Supporting Information 393 App.Tables 3 & 4 and App. Fig. 9 & 10).

  Fig.3band Fig.4a).400Baseline temperature conditions had the strongest positive influence on the401climatic debt (five-year dataset = 0.77 0.04 SE, Fig. 3a; yearly dataset = 0.72 ± ± 402 0.10 SE, Fig. 3b), indicating greater lags in plots with warmer baseline temperature 403 conditions during the period 1961-1990. We also found a positive effect of stand 404 characteristics, through the variables stand age in 1995 and basal area during the 405 year of monitoring, on the magnitude of the lag. More specifically, the lag was larger 406 in denser (five-year dataset = 0.22 0.04 SE, Fig. 3a; yearly dataset = 0.35 0.05 ± ± 407 SE, Fig. 3b and Fig. 4b) and older (five-year dataset = 0.15 0.04 SE, Fig. 3a; ± 408 yearly dataset = 0.46 0.07 SE, Fig. 3b and Fig. 4c) forests. We found no clear ± 409 differences in the climatic debt between coniferous and deciduous forests as well as 410 between inside and outside of the fenced area.

412 3 )

 3 . Notably, the magnitude of the lag was negatively correlated with the anteriority 413 of the last disturbance event (five-year dataset = -0.12 0.02 SE, Fig.3a; yearly ± 414 dataset = -0.10 0.03 SE, Fig.3b), revealing that lags decrease in plots recently ± 415 disturbed. Moreover, the cumulative number of anthropogenic disturbances were to the climatic debt (five-year dataset = -0.17 0.04 SE, Fig. ± 417 3a; yearly dataset = -0.19 0.07 SE, Fig. 3b), indicating that lags decrease in the ± 418 most frequently disturbed plots, following logging events. The cumulative number of 419 natural disturbances was only retained in the models based on the five-year dataset,

807 Figure 2 .808815 Figure 3 .816Figure 1 .

 231 Figure 1. Workflow implemented in the study. The 102 plots of the RENECOFOR network are mapped, and the two datasets used in our analyses are represented in color (in red and blue for the five-year and yearly datasets, respectively). Analyses of temporal trends in mean annual temperature (MAT) and the community temperature index (CTI) and factors contributing to the lag between MAT and CTI (MAT-CTI) are conducted on the two datasets separately. For more details, see the Materials and Methods section.

Figure 2 .Figure 3 .Figure 4 . 2 .7 1 . 2 .

 234212 Figure 2. Temporal trends for mean annual temperature (MAT) and the community temperature index (CTI) in (a) the five-year dataset from 1995 to 2015 and (b) the yearly dataset from 1995 to 2017. Linear trends (dashed lines) with standard errors (ribbons) of estimated MAT or CTI values are shown. For CTI, the values for exclosure or enclosure are plotted separately. Observed values for temperatures (MAT or CTI) are represented with their associated standard errors and symbol sizes are proportional to sample size (number of plots). 926x1296mm (72 x 72 DPI)

Figure 4 .

 4 Workflow conducted to calculate the Community Temperature 66 Indices (CTIs) through the transfer function. 67 App. Figure 5. Supporting information for the construction of Community 68 Temperature Indices (CTIs) using the transfer function. 69 a) Venn diagram depicting the number of taxa recorded in the CATMINAT and 70 EcoPlant databases for species indicator value of temperature. Our dataset 71 includes a total of 877 taxa (all strata of vegetation combined). 72 b) Relationship between the IVs of temperature of CATMINAT and EcoPlant at 73 the species level. 74 c) Map of the sites where the NFI floristic surveys have been performed, and 75 which have been used to model the transfer function. 76 d) Hexagonal binning plot between the CTIs calculated with CATMINAT and 77 EcoPlant at the community level, and the transfer function fitted to obtain CTI 78 values from CATMINAT in degrees Celsius. 79 Text details for App. Figure 5 : To infer temperature from plant species 80 assemblages, we used an approach based on Community Temperature Indices 81 (CTIs) with species Indicator Values (IVs), and a transfer function to model the

  //www1.onf.fr/renecofor). The RENECOFOR network belongs to the European 156 ICP Forests Level II monitoring program (ICP Forests, http://www.icp-forests.org) 157 and aims at monitoring the health and functioning of forest ecosystems. This 158 monitoring network was established in 1992 by the French National Forest Service 159 (ONF, Office National des Forêts). It includes 102 permanents plots distributed in 160 state and communal forests throughout the whole metropolitan French territory 161 (www1.onf.fr/renecofor/sommaire/sites)

	155 http:
	F o r
	P e
	e r
	R e
	v
	i e
	w
	9

154 term French monitoring network of forest ecosystems (RENECOFOR,

  France, its temporal resolution is rather coarse. Per contra, the yearly 178 dataset covers a restricted spatial extent but provides a very detailed time series at 179 high temporal resolution. More information on the spatial location and dominant 180 tree species for each of 44 study plots are provided in Supporting Information App.

	177 forests in 181 Table
	F o r
	P e
	e r
	R e
	v
	i e
	w
	10

Although

176 the five-year dataset covers a large spatial extent representative of all temperate

260 The main determinants of the climatic debt

	243 On the one hand, the CATMINAT database is based on expert knowledge similarly 308
	244 to Ellenberg's IVs (Ellenberg et al., 1992), with temperature preferences of plant
	245 species coded along an ordinal scale from cold (1) to warm (9) conditions. It contains
	246 almost all of species recorded in our dataset (858 out of the 877 species recorded in
	247 all vegetation layers), but CTI values are unitless and direct comparisons with
	248 macroclimatic temperatures are not possible. On the other hand, EcoPlant IVs for
	249 temperature are directly expressed in degrees Celsius (Gégout et al., 2005) allowing
	o r 250 the direct computation of CTI values in degrees Celsius, but much less species from 251 our dataset are recorded (396 out of 877 species). Based on these two databases, we F o r F o r F o r
	P e e r 252 calibrated a transfer function to produce CTI values directly expressed in degrees P e P e P e 253 Celsius and transform the unitless CTI values from the CATMINAT database into a 254 CTI value expressed in degrees Celsius (Supporting Information App. Fig. 4 & 5; e r e r e r
	R e 255 and see Lenoir et al., 2013 for a similar approach). R e R e R e
	v i e w v v v Based on the time series of MAT values and the time series of CTI values, i e i e i e 257 both expressed in degrees Celsius, we computed the lag as the difference between 256 258 the two (MAT CTI), thus getting a new time series of the difference: the climatic -w w w
	259 debt (Bertrand et al., 2016).
	261 Although we were primarily interested in assessing the temporal dynamic of the
	262 difference between MAT and CTI values, thereby testing the effect of time itself on
	263 the direction and magnitude of this difference -as we assumed the climatic debt to
	13 16

  . The velocities of 435 thermophilization we report here (ranging from 0.08 to 0.10 °C/decade) are broadly 436 in the range of those reported in the scientific literature (from 0.03 to 0.38 Supporting Information App. Table. 5), but also are among the highest 438 rates recorded for plant communities (maximum value of 0.27 °C/decade;Duque et 439 al., 2015). Specifically, we show velocities that are twice as fast as the average 440 thermophilization rate estimated for understory plant communities across 441 temperate deciduous forests in Europe (i.e. 0.041 °C/decade; De[START_REF] De Frenne | Microclimate moderates plant responses to 644 macroclimate warming[END_REF].

	437 °C/decade;
	F o r
	P e
	e r
	R e
	v
	i e
	w
	22

climatic debt and its potential determinants 472

  are thus likely to be more accurate than those previously reported 460[START_REF] Elmendorf | Experiment, monitoring, and gradient methods 668 used to infer climate change effects on plant communities yield consistent patterns[END_REF], which matters for informing our society on future scenarios 461 of biodiversity changes(IPBES, 2019). Furthermore, our yearly resurveys over a 462 period of 25 years confirm that the thermophilization rate of understory plant 463 communities in temperate forests is following a linear trajectory over time. This is 464 an important and novel finding that helps us to understand the temporal dynamic of 465 the biotic responses to anthropogenic climate change. Yet, time alone had a very 466 minor contribution to explain the variation in temperature conditions inferred from 467 understory plant communities. This suggests that other abiotic and biotic 468 determinants, such as the macroclimatic context as well as local stand 469 characteristics influenced by anthropogenic disturbances may potentially explain 470 local variations in temperature conditions. One of our most striking results is that the velocity at which atmospheric air 473 temperatures are rising is twice faster, on average, than the velocity at which 474 understory plant communities are responding, generating a climatic debt that is 475 steadily growing over time, even after accounting for the relative contribution of 476 several abiotic and biotic determinants modulating it. Although the idea of an 477 increasing climatic debt over time has already been suggested(Bertrand et al., 478 2011), we provide the first evidence that the climatic debt of understory plant 479 communities in response to the ongoing climate change is steadily increasing over 480 time. This contrasts with a previous report showing no clear temporal signal in the

	F o r r o 459 substitutions and F 481 climatic debt of herbaceous plant communities in French forests
	P e P e
	e r e r
	R e R e
	v v
	i e i e
	w w
	23 24

. Focusing on a 448 completely independent dataset from the period

1965[START_REF] Brice | hence 307 we hypothesized a[END_REF][START_REF] Bertrand | Changes in plant community composition lag behind climate 608 warming in lowland forests[END_REF]

449 were the first to highlight a thermophilization of understory plant communities in 450 the French forests. More recently,

[START_REF] Martin | Short-term 724 climate-induced change in French plant communities[END_REF]

have described a similar 451 thermophilization trend extending to other plant communities at national and 452 regional scales in France over the period 2009-2017. Unlike these previous studies 453 relying either on massive collection of opportunistic surveys over time

(Bertrand et 454 al., 2011)

or citizen science data

[START_REF] Martin | Short-term 724 climate-induced change in French plant communities[END_REF]

, our study is the first that used 455 a long-term monitoring program based on strictly permanent plots intensively 456 monitored over time (up to every year), providing a very detailed time series of 457 changes in understory plant community composition. Hence, we argue that the 458 velocities we report here do not suffer from biases due to space-for-time 471 The

  [START_REF] Lenoir | Climatic microrefugia under anthropogenic 716 climate change: implications for species redistribution[END_REF][START_REF] De Frenne | Global buffering of temperatures under forest 648 canopies[END_REF] Zellweger et al., 2020). In contrast, 528 repeated disturbances are likely to create canopy openings in forest stands, 529 modifying the microhabitat and associated microclimate that strongly exacerbate 530 the impact of macroclimate temperature warming for understory vegetation 531[START_REF] Vanderwel | How do disturbances and environmental 765 heterogeneity affect the pace of forest distribution shifts under climate change? 766[END_REF] Zellweger et al., 2020). This can in turn lead to drastic 532 changes in the overall composition of the plant community (i.e. by increasing the 533 proportion of warm-adapted species through either species turnover or an increase 534 in local species richness) that will reduce the climatic debt. Paying off the climatic Contrary to our initial expectation that large herbivores should contribute to 543 the climatic debt by further increasing the lag in the response of understory plant 544 communities to macroclimate warming, we found no difference in the observed 545 climatic debt between outside and inside the fenced area. This suggests that the 546 climatic debt generated by the buffering effect of microclimate in the understory 547 may mostly stem from the canopy layer and independently from the shrub layer, as

	526 plant communities thus providing thermal microrefugia for species persistence 548 large herbivores are mostly impacting woody species from the understory shrub
	F o r F o r r o 527 (F 549 layer
	P e P e P e
	e r e r e r
	R e R e R e
	v v v
	i e i e i e
	w w w
	26 27

communities following anthropogenic disturbances is due to warm-adapted species 502 that are: (i) replacing more mesophilous species (pure turnover); (ii) increasing the 535 debt of understory plant communities in temperate forests implies to break off this 536 microclimate stability, which can consequently reduce the microrefugia capacity of 537 forest habitats for the survival of forest-dwelling species under anthropogenic 538 climate change

[START_REF] Lenoir | Climatic microrefugia under anthropogenic 716 climate change: implications for species redistribution[END_REF]

. Accepting the idea that the climatic debt of 539 understory plant communities might be a positive sign of local adaptation has 540 significant implications on the concept of climatic debt itself, which might be 541 interpreted as a macroclimatic debt but not a microclimatic debt.

542

Table 1 |

 1 Informations on sites characteristics.

	33 14. App. Table 2 | List of 751 forest species found in the herbaceous layer in our
	34	study.

Table 2

 2 Therefore, we calibrated a transfer function to transform the unitless CTI values 101 from the CATMINAT database into a CTI value expressed in degrees Celsius. By 102 regressing CTI values from EcoPlant against CTI values from CATMINAT, we can 103 obtain predicted CTI values from CATMINAT expressed in degrees Celsius (see these two pools are used to calculate temporal -diversity between two 𝛽169 surveys of a given plot[START_REF] Legendre | A temporal beta-diversity index to identify sites that have 248 changed in exceptional ways in space-time surveys[END_REF]. Then, we computed the average of IVs for 170 Temperature within the gained and lost species. List of 751 forest species inventoried in the herbaceous layer. Species 184 retained for analyses are marked with a "1" in the éponymous column.

	F o r
	P e
	e r
	R e
	v
	i e
	w