Gabrielle Martin¹
email: gabrielle.martin@edu.mnhn.fr

Vincent Devictor

Eric Motard

Nathalie Machon¹

Emmanuelle Porcher¹

Short-term climate-induced change in French plant communities

Keywords: Plants, monitoring, mean community thermal preference, climate change, lifespan, citizen science

Latitudinal and altitudinal range shifts in response to climate change have been reported for numerous animal species, especially those with high dispersal capacities. In plants, the impact of climate change on species distribution or community composition has been documented mainly over long periods (decades) and in specific habitats, often forests. Here, we broaden the results of such long-term, focused studies by examining climate-driven changes in plant community composition over a large area (France) encompassing multiple habitat types and over a short period (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017). To this end, we measured mean community thermal preference, calculated as the community-weighted mean of the Ellenberg temperature indicator value, using data from a standardized participatory monitoring scheme. We report a rapid increase in the mean thermal preference of plant communities at national and regional scales, which we relate to climate change. This reshuffling of plant community composition corresponds to a relative increase in the abundance of warm-vs. cold-adapted species.

However, support for this trend was weaker when considering only the common species, including common annuals. Our results thus suggest for the first time that the response of plant communities to climate change involves subtle changes affecting all species rare and common, which can nonetheless be detected over short time periods. Whether such changes are sufficient to cope with the current climate warming remains to be ascertained.

INTRODUCTION

Climate change is considered a significant global threat to biodiversity [1][2][3], with impacts from individuals to ecosystems. In ecological communities, climate change drives compositional change, in particular via latitudinal and altitudinal shifts, especially in species with high dispersal capacities, including birds, insects or marine invertebrates [3][4][5]. In plants, which are less mobile, elevation has shifted upward during the 20 th century, e.g. in highland French temperate forests [6], although these shifts lag temperature change. In contrast, lowland forest plant communities have not shifted northward, highlighting their greater inertia [7] or higher shading mitigating climate warming [8]. Whether and how plant communities are able to track climate change in all habitats is still unclear.

Owing to the scarcity of available data, most studies demonstrating an effect of climate change on plant communities used heterogeneous data collected over decades and focused on a particular habitat. To understand better how climate change and other drivers influence biodiversity, real-time monitoring data using standardized schemes are required across large spatio-temporal scales [9]. Citizen science, the involvement of volunteers in research [10], represents a powerful tool that can contribute to investigating the impacts of climate change on biodiversity [4,11]. Such monitoring makes it possible to report subtle responses to climate change.

In this study, we examine nine years of change in plant community composition across France, using data collected through a participatory monitoring scheme, to assess whether such data can be used to detect short-term effects of climate warming beyond alpine or forest areas. The systematic sampling of monitoring sites particularly encompasses numerous lowland open habitats, which have received limited attention in studies of the response of plant communities to climate change. We address the following questions: (1) has the mean plant community thermal preference changed over the past nine years? (2) Can this change be related to climatic variables? (3) Are the temporal changes in mean community thermal preference and in abundances of individual species related to species lifespan?

MATERIALS AND METHODS

Monitoring data

We used data from Vigie-flore, a French citizen science program monitoring wild flora. The data were collected yearly by 321 skilled amateur botanists between 2009 and 2017 in 586 1kmx1km squares sampled from a systematic grid (one square every 10 km), which ensures representative sampling of habitats. Each square contains eight systematically distributed 10 m² plots (Figure S1) divided into ten 1 m²-quadrats. In each quadrat, the presence of all vascular plants was recorded, as well as habitat type following the CORINE biotope nomenclature, a European hierarchical classification of vegetation types [12]. Most squares were in open habitats (27% artificial land cover, 29% farmland, 16% meadows vs. 22% in forests, Electronic Supplementary Material I). For each plot, the number of quadrats in which a species was observed provides a proxy for species abundance. Individual plots were recorded at different intensities over time, depending on the behaviour of recorders, such that on average there were 2.9 years of observation per plot, and 5.3 plots sampled per 1 km² square, for a total of 3,118 plots.

Species characteristics

We collected information for two species-specific attributes most likely to be involved in plant species response to climate change: (1) the Ellenberg temperature indicator value (hereafter ETIV), characterizing the optimum temperature class for growth and survival of a species (nine classes), and ( 2) lifespan (Electronic Supplementary Material II and Table S1).

ETIV and lifespan were available for 1,709 and 1,780 out of the 2,428 species sampled, representing 85.6% and 89.3% of total observations, respectively. We calculated the community-weighted Mean Thermal Preference (hereafter MTP), i.e. the abundance-weighted sum of ETIV of all species [13]. We also calculated MTP with presence/absence data instead of abundance, to test whether temporal trends could be detected with lower resolution data [14]. Finally, we calculated MTP within two extreme lifespan classes: annuals vs. perennials.

The survey sites are not distributed homogeneously in space and time, which may cause spurious temporal trends in MTP due to a spatial displacement of sampled sites. To control for this, we performed the same analysis at a regional level, in Île-de-France, the region with the highest site density (188 squares sampled throughout the nine years, Figure S1).

Climate data

For each plot, we collected daily precipitation sum, mean, minimum and maximum temperature between 2009 and 2017, from which we calculated 10 climatic variables per plot.

Two variables used annual data: mean annual temperature and annual temperature anomaly;

while the remaining eight used information from January to May (the growing period of most plant species): mean temperature, temperature anomaly, mean maximum temperature, number of heatwave days (T°C>27°C), freezing days (T°C<0°C), rainy days, drought days, and pluviometry. We used a random forest analysis to extract the most important climatic variable explaining MTP, which was mean annual temperature (Electronic Supplementary Material III and Table S2).

Statistical analyses

We first checked for a potential bias attributable to uneven sampling effort by verifying that there was no temporal trend in the average latitude of sampled sites. We then used Bayesian hierarchical models to estimate the temporal trend in MTP (1) at a national scale and within Île-de-France, (2) using either abundance data or presence/absence data and (3) in communities where only annuals vs. only perennials species were retained (Electronic Supplementary Material IV). Note that this approach captures changes in the relative abundance of species with fixed thermal preference rather than changes in species-specific adaptation to new climatic conditions. We then examined the relationship between (1) the temporal change in MTP in a plot between 2009 and 2017, and (2) the temporal change in mean annual temperature over the same period, using Spearman's correlation coefficient.

To further investigate the role of species lifespan in the community-level trends, we used Bayesian models to estimate the temporal trend in abundance for the 550 most common species (75 annuals and 349 perennials) observed at least four out of nine years and in at least ten squares (Electronic Supplementary Material IV). We tested whether species trends were related to lifespan (annuals vs. perennials), to ETIV, or to the interaction of both variables by a linear regression in which each species was weighted by the inverse of the standard error associated with the estimated temporal trend. All analyses were performed using R [15] (Electronic Supplementary Material V).

RESULTS

At the community level, MTP increased through time at the national scale, both with abundance and presence/absence data (Figure 1a,c), particularly in northern France (Electronic Supplementary Material IV). This increase was also observed within Île-de-France (Figures 1c andS2). The increase in MTP of plant communities occurred simultaneously with an increase in mean annual temperature (Electronic Supplementary Material III, Table S3, Figure S3), such that there was a significant positive but weak correlation between MTP change and the mean annual temperature change over time (Figure 1d): plots in which the temperature has increased tended to be also plots for which MTP has increased. Moreover, we found a significant increase in MTP of annual, but not of perennial species, at national and regional scales (Figure 1b,c).

Changes at the community level were weakly explained by differences in temporal trends for common species considered independently. We found a significant relationship between species temporal trends and lifespan (Estimate=0.194, n=424, df=1.402, F=3.978, P=0.047, Figure 2). The interaction between ETIV and lifespan was marginally significant (Estimate=-0.028, n=424, df=1.402, F=3.013, P=0.083), but no relation was found between temporal trends and ETIV.

DISCUSSION

Our work complements longer-term, single habitat studies by documenting a small but steady increase in MTP of plant communities over all habitat types at national and regional scales, detectable over nine years only, and most likely driven by temperature change. This increase is measurable even with presence/absence data, suggesting a true reshuffling of plant communities, as opposed to mere abundance fluctuations of common species. In contrast to Bertrand and collaborators [7], who observed no response of lowland forest plant communities to climate change over a longer period, we show with a dataset encompassing mostly low-altitude sites that lowland plant communities can be modified by climate change, despite their limited dispersal capabilities. The difference may arise because our dataset contained a high proportion of open habitats, in which plant communities may experience higher turnover than in forests [8]. MTP changes resulted from local increases in the relative abundance of warm-vs. cold-adapted species, which, over nine years, are detectable in annuals but not perennials. Yet, examination of the temporal trends of the 424 most common annual or perennial species did not reveal this pattern. Non-congruence of the community vs. species-level analysis, the latter covering only 55 % of the observations, suggests that all species, including rarer species, are affected and calls for a more thorough analysis of species traits likely to drive the observed changes.

The general trend of increased relative abundance of warm-vs. cold-adapted species is consistent with longer-term studies documenting responses to climate change via species latitudinal and altitudinal range shifts, or local changes in abundance (including local decline and extinction of cold-adapted species) [6][7][8]16,17]. Yet, the causal role of climate change in the observed increase in the MTP of plant communities still has to be confirmed, for example by checking that this increase is sustained in the near future, or by examining possible confounding drivers, such as eutrophication or urbanization. Urbanization is unlikely to be responsible for increased MTP, because only a quarter of plots are located in urban areas.

Although the Ellenberg nitrogen indicator value is weakly but positively correlated with ETIV, the mean community nitrogen preference or the abundance of nitrophilous species did not increase over the same period (Electronic Supplementary Material VI and Figure S4). Therefore, eutrophication should have a limited role in the observed increase in MTP.

The detection of this subtle but rapid change in plant communities was made possible by citizen science monitoring and a large sampling effort in space and time. Such data can be of great value [18] but come with several caveats deserving further attention. For example, possible biases could be related to uneven sampling effort, identification errors or uncomplete trait data. The effect of spatial heterogeneity in the sampling effort was tested by examining trends within a small region. Identification errors and missing trait data are unlikely to be systematically biased in favour of high or low species thermal preference, such that these problems should not influence our results. Finally, the observed trend in MTP could be caused partly by temporal changes in species phenology. Because the Vigie-flore protocol recommends a single annual survey, there may be a phenology-influenced bias in species detection and identifiability against undeveloped or non-flowered species. Climate change has led to phenological shifts in flowering plants [19], which could increase the detectability of warm-adapted species as the climate warms.

Whether such rapid community changes are sufficient a response to current climate warming, or instead reflect mostly a maladaptive loss of cold-adapted species, remains to be ascertained. Future research efforts could also be focused on analysing trends in MTP across habitats or protection regimes and on identifying possible consequences of changes for biotic homogenization or plant-pollinator interactions.

I. Habitat of sampled plots

A total of 3,118 plots were sampled in 586 1kmx1km squares between 2009 and 2017 as part of the Vigie-flore program. Habitat type following the CORINE biotope nomenclature [1] was reported by observers for 2,985 plots (ca. 96% of plots). These plots were distributed among the first levels of the CORINE biotope classification as follows: 21 plots in coastal habitats, 13 plots in wetland, 488 plots in grassland and shrub, 682 plots in woodland, 24 plots in marsh and bog, 36 plots in rocky, habitats, 895 plots in agricultural habitats and 826 plots in artificial land use.

II. Species attributes

We collected information for two species-specific attributes: (1) the Ellenberg temperature indicator values (hereafter ETIV), characterizing the optimum temperature class for growth and survival of a species (nine classes), from the database Baseflor [2], and (2) lifespan, described using seven broad classes obtained by combining classes from the Biolflor [3] and LEDA [4] databases (Table S1).

Table S1. Procedure used to combine lifespan information from Biolflor and LEDA traitbases.

III. Climate data and Random Forest analysis

For each plot, we collected daily mean temperature, daily minimum temperature, daily maximum temperature and daily precipitation sum between 2009 and 2017 from the ENSEMBLES dataset [5], using data from the nearest E-OBS station in France on a 0.25 degree grid. Climate data were extracted using the R package climateExtract [6], and devtools [7]. We then calculated 10 climatic variables per plot, two using annual data (mean annual temperature and annual temperature anomaly) and eight using information from January to May (growing period of many plants): mean temperature, temperature anomaly, mean maximum temperature, number of heatwave days (T°C>27°C), freezing days (T°C<0°C), rainy days, drought days, and pluviometry. The temperature anomaly is the sum of the differences between the observed daily temperature and the mean temperature over a period on a site.

We then performed a Random Forest regression analysis to identify the climatic variable(s) explaining the greatest variability in MTP of plant communities. The Random Forest analysis is a nonparametric technique derived from classification and regression trees. This method, which has been applied in various ecological studies [8,9], shows high accuracy and the ability to model complex interactions between variables. To perform the model, we defined two essential parameters: mtry, i.e. the number of variables to try at each split, and ntree, i.e.

the number of trees to run. To maximize the accuracy of the Random Forest analysis, i.e.

minimize the measure of the prediction error (Out Of Bag (OOB) error), mtry and ntree can be varied. We determined the parameter mtry thanks to the internal RF function TuneRF [R package randomForest, 10]; this function computes the optimal number of variables starting from the default (total number of variables/3 for regression) and tests half and twice this threshold for the value with the minimum OOB error rate. It is already known that a mtry of 1 can produce good accuracy [10,11], but it is recommended to include at least two variables to avoid using also the weaker regressors as splitters [12]. Here, even if the best accuracy was obtained for a mtry of 1, we performed the Random Forest regression analysis with a mtry of 2 following Grömping [12]. The ntree parameter was set to 1000 to obtain stable results.

Variable importance of the 10 climatic variables (i.e. the mean annual temperature, the annual temperature anomaly, and eight variables estimated from January to May: mean temperature, temperature anomaly, mean maximum temperature, number of heatwave days, freezing days, rainy days, drought days, and pluviometry) was estimated as the total decrease in node impurities, measured by residual sum of squares from splitting on the variable, averaged over all trees.

Finally, as a result of the Random Forest regression analysis, we selected the annual mean temperature, the most relevant climatic variable by ranking the variables according to their importance measure (Table S2).

We performed a linear regression of the annual mean temperature against years to describe the temporal trend of climate change on the sampled plots of the study. We included geographical coordinates in the model to account for residual spatial autocorrelation. We recovered one slope per square, since climate data were only available at that scale; all plots within a scale were thus attributed the same temporal trend in climate data, describing the temporal trend of annual mean temperature.

Table S2. Importance measures of the 10 climatic variables resulting from a Random Forest regression analysis and explaining the share of variability of the MTP of plant communities. where MTPijk is the mean temperature preference of the plant community in plot k of square j and year i, m is the intercept, lat and long are the latitude and longitude of the plot, Square is the random effect of square with variance ²s and ²ijk is the residual variance. We chose to model non-linear spatial patterns in MTP squared latitude and longitude terms over more complex models, because on visual examination of the relationship between MTP and latitude or longitude, we expected mostly a linear effect, potentially with a plateau at high or low latitudes that a second order polynomial model correctly represents.

Climatic variables Importance

Note that we performed a linear regression to test for an interaction between latitude and temporal trend of the MTP in order to detect spatial differences at the national scale (not detailed here). At the national scale, the temporal trend increased significantly with latitude, i.e. the trend was significantly larger in northern vs. southern France (significant latitude:year interaction, not shown).

We also performed Bayesian hierarchical models to generate estimates of the temporal trend in common species abundance, as estimated at the 10m² plot level by the number of 1m² quadrats in which these species were observed. We selected the 550 most common species that were observed at least four out of nine years and in at least ten squares. We modelled species abundance assuming a binomial distribution with n = 10 trials 1 and a logittransformed probability p depending linearly on year, survey plot, square and geographical

FIGURESFigure 1 .

 1 FIGURES

Figure 2 .

 2 Figure 2. Relationship between the temporal trends of the most common species and their

  ) perennials perennial p (perennial-pollakanthic) perennials perennial 396 p (perennial-pollakanthic) perennials + poly-annuals < 5 years (shortlived perennials) perennial 28 p (perennial-pollakanthic) perennials + poly-annuals < 5 years (shortlived perennials) + poly-annuals > 50 years (long-lived perennials) perennial p (perennial-pollakanthic) perennials + poly-annuals < 5 years (shortlived perennials) + poly-annuals > 50 years (long-lived perennials) + polyannuals 5-annuals > 50 years (long-lived perennials) + poly-annuals 5-50 years (medium-lived perennials) perennial a (annual) perennials + poly-annuals 5-50 years (medium-lived perennials) annual/perennial a (annual) -p (perennialpollakanthic) perennials + poly-annuals 5-50 years (medium-lived perennials) annual/perennial p (perennial-pollakanthic) perennials + poly-annuals 5-50 years perennial (medium-lived perennials) monocarpic bi-annuals and poly-annuals + poly-annuals < 5 years (short-lived perennials) biennial/perennial b (biennial) -p (perennialpollakanthic) perennials + strict monocarpic bi-annuals and poly-annuals + poly-annuals < 5 years (short-lived perennials) biennial/perennial a (annual) -b (biennial) -h (perennial-hapaxanthic)p (perennial-pollakanthic) perennials + strict monocarpic bi-annuals and poly-annuals + poly-annuals < 5 years (short-lived perennials) + poly-annuals 5-50 years (medium-lived perennials) annual/biennial/perennial p (perennial-pollakanthic) perennials + strict monocarpic bi-annuals and poly-annuals + poly-annuals < 5 years (short-lived perennials) + poly-annuals 5-50 years (medium-lived perennials) biennial/perennial b (biennial) -p (perennialpollakanthic) perennials + strict monocarpic bi-annuals and poly-annuals + poly-annuals 5-50 years (medium-lived perennials) biennial/perennial p (perennial-pollakanthic) perennials + strict monocarpic bi-annuals and poly-annuals + poly-annuals 5-50 years (medium-lived perennials) ) -b (biennial) -p (perennial-pollakanthic) perennials + summer annuals + polyannuals < 5 years (short-lived perennials) annual/biennial/perennial a (annual) -p (perennialpollakanthic) perennials + summer annuals + polyannuals < 5 years (short-lived perennials) + poly-annuals 5-50 years (medium-lived perennials) annuals + strict monocarpic bi-annuals and poly-annuals + poly-annuals < 5 years (short-lived perennials) annual/biennial/perennial p (perennial-pollakanthic) perennials + summer annuals + strict monocarpic bi-annuals and poly-annuals + poly-annuals < 5 years (short-lived perennials) annual/biennial/perennial h (perennial-hapaxanthic) -p (perennialpollakanthic) perennials + summer annuals + strict monocarpic bi-annuals and poly-annuals + poly-annuals 5-50 years (medium-lived annual/biennial/perennial annuals + winter annuals + strict monocarpic bi-annuals and poly-annuals annual/biennial/perennial p (perennial-pollakanthic) perennials + winter annuals annual/perennial p (perennial-pollakanthic) perennials + winter annuals + poly-annuals > 50 years (long-lived perennials) annual/perennial a (annual) -b (biennial) perennials + winter annuals + strict monocarpic bi-annuals and poly-annuals annual/biennial/perennial b (biennial) perennials + winter annuals + strict monocarpic bi-annuals and poly-annuals annual/biennial/perennial b (biennial) -h (perennialhapaxanthic) perennials + winter annuals + strict monocarpic bi-annuals and poly-annuals annual/biennial/perennial a (annual) -b (biennial) poly-annuals < 5 years (short-lived perennials) annual/biennial/perennial b (biennial) poly-annuals < 5 years (short-lived perennials) biennial/perennial b (biennial) -p (perennialpollakanthic) poly-annuals < 5 years (short-lived perennials) biennial/perennial p (perennial-pollakanthic) poly-annuals < 5 years (short-lived perennials) perennial p (perennial-pollakanthic) poly-annuals < 5 years (short-lived perennials) + poly-annuals 5-50 years (medium-lived perennials) perennial p (perennial-pollakanthic) poly-annuals > 50 years (long-lived perennials) perennial p (perennial-pollakanthic) poly-annuals > 50 years (long-lived perennials) + poly-annuals 5-50 years (medium-lived perennials) strict monocarpic bi-annuals and polyannuals annual/biennial/perennial a (annual) -b (biennial) strict monocarpic bi-annuals and polyannuals annual/biennial/perennial b (biennial) strict monocarpic bi-annuals and poly-) -b (biennial) -h (perennial-hapaxanthic) strict monocarpic bi-annuals and polyannuals + poly-annuals < 5 years (shortlived perennials) annual/biennial/perennial b (biennial) -h (perennialhapaxanthic) strict monocarpic bi-annuals and polyannuals + poly-annuals < 5 years (shortlived perennials) biennial/perennial p (perennial-pollakanthic) strict monocarpic bi-annuals and polyannuals + poly-annuals < 5 years (shortlived perennials) summer annuals + strict monocarpic biannuals and poly-annuals annual/biennial/perennial a (annual) -b (biennial) summer annuals + strict monocarpic biannuals and poly-annuals annual/biennial/perennial a (annual) -b (biennial) -h (perennial-hapaxanthic) summer annuals + strict monocarpic biannuals and poly-annuals annual/biennial/perennial

  Annual mean temperature 725 Number of freezing day from January to May 365 Mean temperature from January to May 342 Mean maximum temperature from January to May 243 Number of drought day from January to May 204 Pluviometry from January to May 160 Annual temperature anomaly 135 Temperature anomaly from January to May 133 Number of rainy day from January to May 99 Number of heatwave day from January to May 55 IV. Estimation of temporal trends in mean temperature preference and species occurrence We performed Bayesian hierarchical models to estimate the temporal trend in MTP of plant communities (1) at a national scale and within Île-de-France, (2) using either abundance data or presence/absence data and (3) of annuals vs. perennials species. Note that MTP is comparable to the Community Temperature Index [13], but is not a direct measure of mean temperature, because ETIV consists in classes of temperature preference. The underlying model describes the MTP value as a function of year for all plant communities. We assumed that MTP of plant communities was normally distributed, with a mean that depended linearly on year, geographical coordinates and survey square, and constant variance. The linear model was the following: MTPijk = m + a*yeari + b*latijk + c*longijk + d*latijk² + e*longijk² + Squarej + ²ijk

ACKNOWLEDGEMENTS

We thank the Vigie-flore volunteers (List S1) who made this work possible, the E-OBS dataset from the EU-FP6 project ENSEMBLES (http://ensembles-eu.metoffice.com), the data providers in the ECA&D project (http://www.ecad.eu), François Duchenne, Nicolas Dubos, Karine Princé, Jeanne Vallet, Sandrine Pavoine, Guillaume Fried, François Chiron, Bruno Colas for their constructive comments and advice. Datasets and code are available at https://sandbox.zenodo.org/record/245806#.W7OIVqKaHvY .

FUNDING

This project was supported by grants from the Région Île-de-France (DIM ASTREA) and the GDR Ecostat.

where pijk is the probability to observe the species in each 1m² square of plot k of square j and year i, m is the intercept, lat and long are the latitude and longitude of the plot, Square is the random effect of square with variance ²s, Plot is the random effect of plot with variance ²p.

For all Bayesian hierarchical models, we chose only non-informative priors, sampled from Gaussian distributions with a mean of 0 and a variance of 10 10 for most parameters, and gamma distributions with shape parameter k = 0.001 and rate  = 0.001 for precision (inverse of variance) parameters. Parameters were estimated by fitting the model in a Bayesian mode of inference using JAGS [14], through R [15] using the package R2jags [16]. Four chains were run using priors to draw 100,000 MCMC iterations, discarding the first 50,000 iterations. Samples were thinned by 1 in 10 to reduce autocorrelation and obtain mean trends from the remaining posterior distributions. Convergence was checked using the Rubin statistics (Rhat), comparing within-chain variance to between-chain variance [17] (all Rhat < 1.2) and visual inspection of the traceplots (superposition of the 4 chains) produced through the R package mcmcplots [18]. Posterior mean community and species trends were interpreted as different from zero when 95 % of the highest posterior density interval did not include zero.

V. R packages used

All analyses were performed using R version 3.2.4 [15] and R packages Taxonstand version 2.1 [19], R packages TR8 [START_REF] Bocci | TR8: an R package for easily retrieving plant species traits[END_REF], R packages climateExtract version 0.2.0 [6], R packages randomForest [10], R packages caret version 6.0-78 [START_REF] Kuhn | caret: Classification and Regression Training[END_REF], R packages plotmo version 3.3.4 [START_REF] Milborrow | plotmo: Plot a Model's Response and Residuals[END_REF], R packages R2jags version 0.5-7 [16], R packages mcmcplots version 0.4.2 [18], R packages car [START_REF] Fox | An {R} Companion to Applied Regression[END_REF], R packages lme4 [START_REF] Bates | Fitting Linear Mixed-Effects Models Using lme4[END_REF], R packages nlme version 3.1-131 [START_REF] Pinheiro | nlme: Linear and Nonlinear Mixed Effects Models[END_REF], R packages dplyr version 0.7.3 [START_REF] Wickham | dplyr: A Grammar of Data Manipulation[END_REF], R packages reshape2 [START_REF] Wickham | Reshaping Data with the reshape Package[END_REF], R packages devtools version 1.13.4 [7], R packages effects version 2.1.0 [START_REF] Fox | Effect Displays in R for Generalised Linear Models[END_REF], R packages ggplot2 version 2.1.0 [START_REF] Wickham | Ggplot2: Elegant Graphics for Data Analysis[END_REF], R packages ggthemes version 3.4.0 [START_REF] Arnold | ggthemes: Extra Themes, Scales and Geoms for 'ggplot2[END_REF], R packages gridExtra version 2.3 [START_REF] Auguie | gridExtra: Miscellaneous Functions for 'Grid' Graphics. List S1. The 321 Vigie-flore volunteers involved in the wild flora monitoring[END_REF]. Plant functional traits were collected using the R package TR8 [START_REF] Bocci | TR8: an R package for easily retrieving plant species traits[END_REF].

VI. Relationship between Ellenberg Temperature and Ellenberg Nitrogen Indicators

Values

Across the 1,473 species for which we could obtain Ellenberg indicator values (out of the 2,428 species of the Vigie-flore dataset), we observed that the Ellenberg temperature indicator value was positively correlated with the Ellenberg nitrogen indicator value (i.e. nitrogen preference, Pearson's correlation = 0.27), although the correlation was not very strong. There is thus a possibility that some of the changes in plant community composition we observed are also driven by ecosystem eutrophication, in addition to climate warming. To rule out this possibility, we tested whether mean nitrogen preference changed over time in the plant communities we monitored. To this end, we calculated the Community Weighted-Mean Ellenberg nitrogen indicator value (hereafter ENIV for the species characteristic, and MNP for the community-level Mean Nitrogen Preference), in the same way as the Mean Temperature Preference. We used a Bayesian hierarchical model to estimate the temporal trend in MNP in plant communities. The linear model was similar to the one we used for Mean Temperature Preference:

where MNPijk is the mean nitrogen preference of the plant community in plot k of square j and year i, m is the intercept, lat and long are the latitude and longitude of the plot, Square is the random effect of square with variance  ² s and  ² ijk is the residual variance.

The temporal trend in MNP was not different from zero (Figure S4, -0.0047, CI95%=[-0.013, 0.0036]), i.e. the mean Ellenberg value for nitrogen preference of plant communities has not changed significantly over the past nine years, which suggests that the increase in Mean Temperature Preference we observed is unlikely to be driven by the (weak) correlation between the thermal preference and the nitrophily of plant species. In other words, while we do observe an increase in the relative abundance of thermophilic species (increase in MTP), the relative abundance of nitrophilous species remains stable. Table S3. Type II Anova results for the linear mixed models testing for a temporal change in the mean annual temperature over time. Plant species are sampled in 1 km X 1 km squares. In each square, plant species were recorded in several plots with a standardized protocol. On the map, square colours indicate the average (over all plots) temporal trend in Mean Temperature Preference (MTP) of plant communities on the study period. Six squares had average temporal trends in MTP of plant communities greater than 0.5 (0.51, 0.55, 0.62 and 1) or below -0.5 (-0.67 and -0.79), they were capped at 0.5 and -0.5 respectively to improve readability. Plots and squares that were surveyed only once, for which no temporal trend could be calculated, are plotted as circles.