Predicted changes in the functional structure of earthworm assemblages in France driven by climate change

Species shift their ranges as a consequence of climate change, hence modifying the structure of local assemblages. This may have important consequences for ecosystem functioning in the case of ecosystem engineers such as earthworms, especially when community restructuring leads to an alteration of their functional diversity. Here, we aimed to model the potential modification of the functional diversity of French earthworm assemblages in a context of climate change.


| INTRODUC TI ON
Studies of large-scale biodiversity trends often focus on the taxonomic identity of species to describe modifications in populations' abundance, range areas or community composition (Ceballos et al., 2017;Pimm et al., 2014). One reason is that it is believed that these changes will translate into an alteration of the functionality of ecosystems and the services they provide to human well-being Hooper et al., 2012). However, not all species have equal actions on ecosystem processes. In fact, the role of organisms in ecosystems is closely related to their functional traits (Cadotte et al., 2011). Therefore, the impact of anthropogenic disturbance on biodiversity would often be more efficiently characterized by the changes in the richness, redundancy or complementarity of species traits (Le Bagousse-Pinguet et al., 2019), that is the various facets of functional diversity (Mason et al., 2005). There is evidence that taxonomic diversity-the number and relative abundance of species as defined by taxonomic classification -is not necessarily a good proxy for functional diversity (Devictor et al., 2010). Because of this, many assessments of biodiversity changes in an era of rapid human-induced global transformations lack the functional perspective that would make them more useful to anticipate the impact they may have on ecosystems (Carmona et al., 2021).
Community composition and diversity in the soil compartment have proved to be a key component of ecosystem processes (Bardgett & van der Putten, 2014). In this regard, a relationship between soil invertebrates' traits and several ecosystem services such as climate and water regulation, soil stability and fertility, and primary production, has been established (de Bello et al., 2010;Setälä et al., 1998). At the same time, the functional diversity of soil communities is often affected by environmental gradients and human disturbance (Pelosi et al., 2014;Santorufo et al., 2014;Vincent et al., 2018). Among the many organisms that compose the soil biota, earthworms play an essential role in soil function and ecosystem services (Blouin et al., 2013;Liu et al., 2019). Some of these roles are directly dependent on species traits (Marichal et al., 2017). For instance, burrowing ability, which influences soil water retention and physical properties , can be linked to earthworm body size (Piearce, 1983;Quillin, 2000). Earthworms' main ecological categories-anecics, epigeics and endogeics, which presumably describe how earthworms influence soil functioning [Bouché (1977), but see Bottinelli and Capowiez (2021)]-are also closely related to a number of morphological and life-history traits, for example variations in generation time or tegument thickness (Briones & Álvarez-Otero, 2018). In addition, some evidence shows that the diversity of functional traits within earthworm assemblages regulate some aspects of soil functioning such as plant litter decay and soil organic carbon level . For these reasons, future changes in the functional structure of earthworm communities may be among the most striking examples of the consequences of global changes on ecosystem functioning.
Although the impact of land use changes, such as agricultural practices or pollution, on earthworm communities has been extensively studied (e.g. Jordan et al., 2004;Pelosi et al., 2014), climate change also has the potential to contribute to the restructuring of macrofauna functional diversity. Climate change affects the environmental conditions of organisms in such a way that, unless they adapt or are preadapted to their new climatic conditions (Hoffmann & Sgro, 2011), species must shift their distribution towards newly suitable habitats to remain in the climatic niche they are adapted to (Chen et al., 2011). Locally, it means that some species will go locally extinct if temperature conditions become unsuitable, while new species may colonize sites that have become warm enough for them to survive (Parmesan, 2006). This process will thus result in a turnover of species within local assemblages, which reshuffles species according to their climatic niche (Devictor et al., 2008). However, species that go extinct or that colonize also harbour other non-climatic traits, leading to changes in the local composition of species traits and, potentially, in the functionality of the reshuffled communities (Fourcade et al., 2021;Wieczynski et al., 2019). Soil species mostly respond to the microclimatic conditions they experience locally rather than to macroclimate (Lembrechts et al., 2020). In this regard, the metabolism of earthworms is known to be largely dependent on soil temperature variation (e.g. Daniel et al., 1996;Eriksen-Hamel & Whalen, 2006). Nevertheless, large changes in macroclimatic regimes may affect soil properties through, for example, droughts that can affect soil macrofauna (Wang et al., 2020). There is also evidence that distribution pattern of earthworm communities has a large-scale spatial structure that is partially mediated by climatic factors (Phillips et al., 2019). Therefore, it is reasonable to assume that community assembly processes of soil organisms such as earthworms will be impacted by climate change (Singh et al., 2019).
Species range shifts in response to climate change can be investigated using methodological tools known as species distribution models (SDMs), or ecological niche models (ENMs). The general principle of these approaches is to establish the statistical relationships between species' presence and environmental predictors to model species' niche and project it in space, producing a spatially-explicit assessment of environmental suitability (Elith & Leathwick, 2009). If predictors of the niche can be projected in time under various scenarios of climate change, models can be similarly projected in the future to predict range shifts in a changing climate (Elith & Leathwick, 2009). Fitting SDMs to soil biota can be challenging, because of the heterogeneity of soil habitats at various spatial scales and because of the existing feedbacks between the metabolism of soil organisms and their habitat (Schröder, 2008). Such approach, though, has been successfully carried out for modelling the distributions of a few earthworm species and to identify their environmental drivers (De Wandeler et al., 2016;Marchán et al., 2015;Palm et al., 2013).
Aggregating multiple single-species SDMs (a technique known as stacked SDMs) allows to model the current and future richness of assemblages (Biber et al., 2019), providing that their species richness is effectively calibrated (D'Amen et al., 2015). By associating the predicted composition of assemblages to the functional traits of the modelled species, stacked SDMs allow to predict current and future functional diversity of assemblages under various scenarios of global change (e.g. Oliveira et al., 2019;Pradervand et al., 2014;Toro et al., 2015). Such functional approach, though, is relatively rare compared to the profusion of SDM studies in the literature; to our knowledge, it has never been employed in the context of soil macrofauna.
To provide reliable predictions, SDMs must be based on solid field data reporting the presence, and optionally the absence, of the modelled species in a way that is as unbiased as possible (Costa et al., 2009;Kramer-Schadt et al., 2013). Generally, the amount of data available for soil organisms is scarce compared to above-ground biodiversity (Cameron et al., 2018). However, in France, a remarkable inventory of earthworms has been compiled in the 1960s by Marcel Bouché, who sampled over 1300 sites evenly located across the whole country using a standardized protocol (Bouché, 1972).
Because Bouché also reports a large number of morphological and anatomical traits, this constitutes an invaluable source of data to investigate the functional diversity of earthworm assemblages at a country-scale. Here, we fitted SDMs to Bouché's (1972) data, using soil variables and bioclimatic predictors describing climatic conditions in the 1960s. Models were projected in time to map the potential composition of earthworm assemblages-defined as 2.5 arc-min grid cells-in the present and under different scenarios of climate change. Using the database of species traits, we then calculated three complementary descriptors of functional diversity within these assemblages. Evidence shows that earthworms are able to disperse rapidly (moving decision being triggered within days) into new suitable habitats in response to changes in population density and habitat quality (Mathieu et al., 2010), and can occasionally perform long-distance movements through passive dispersal (e.g. crossing water bodies on floating wood; see Chen et al., 2021). However, it is also observed that French earthworm communities are not perfectly at equilibrium with current climate conditions, suggesting that colonization can be a long process at biogeographic scales (Mathieu & Davies, 2014). Therefore, we proposed predictions under two hypotheses: complete dispersal, where species are allowed to fill their entire potential niche, and no dispersal, where changes in climatic conditions can only result in local species extirpations. These results were used to describe potential changes in earthworm functional diversity in each French ecoregion depending on climate change scenarios.

| Data collection
We used data on the presence and absence of earthworm species in France that were collected in the 1960s by Bouché (1972).
More precisely, 1366 sites were sampled between 1966 and 1969, evenly spaced across the entire French European territory ( Figure 1).
Sampling sites were selected in order to represent the most common habitat types within each region, avoiding anthropogenic soils such as crops. Sampling was performed mostly in autumn and spring, which correspond to the period of maximum activity of adult earthworms. In each site, three or four soil samples were collected, each corresponding to a 1 × 1 m 2 of ca. 30-cm depth. Soil fauna was extracted, earthworms isolated from other organisms, and all individuals were identified at the species level and counted twice, included juveniles that were reared until maturity. For our analyses, we used only the presence or absence of species, and discarded species that were found at less than 10 locations because their modelling would be unreliable. In total, 105 species were identified, which were reduced to 44 after filtering rare species.
In addition to distributional data, Bouché (1972) reported several information regarding the morphology and internal anatomy of the species he collected. Complemented with data from Sims and Gerard (1985), this constitutes a large database of 36 species traits related to size, mobility, ecology and reproduction (Table   S1), suitable to describe the diversity of earthworm communities F I G U R E 1 Location of sampling sites of earthworms conducted by Bouché (1972)  We selected climate predictors from the 19 standard bioclimatic variables (Booth et al., 1987) and the 18 ENVIREM variables that were developed recently to complement the former (Title & Bemmels, 2018). Because many of them are highly correlated, and because we had no a priori expectation regarding which could be the most influential for describing earthworms' distributions, we retained only the eight variables that allowed for variance inflation factor to remain <4: isothermality (BIO3), mean temperature of the wettest quarter (BIO 8), precipitation seasonality (BIO 15), precipitation of coldest quarter (BIO19), number of months with mean temperature >10°C, potential evapotranspiration of the driest quarter, potential evapotranspiration of the coldest quarter and potential evapotranspiration of the warmest quarter. Using only a selection of those variables ensures that we avoid overfitting in SDMs that can be caused by multicollinearity (Dormann et al., 2013). Moreover, because there is evidence for a legacy of past climate in the structure of French earthworms' communities (Mathieu & Davies, 2014), we also considered mean annual temperature (BIO1) and annual precipitation (BIO12) during the last glacial maximum (LGM). BIO1 during LGM was discarded because it was highly correlated with the contemporary climate variables.
All climatic variables were computed with functions in the "dismo" (Hijmans et al., 2020) and "envirem" (Title & Bemmels, 2018) R packages, using temperature and precipitation data obtained for different time periods from the CHELSA project (Karger et al., 2017). The CHELSA-TraCE21k dataset (Karger et al., 2021) was used to produce LGM variables, while the CHELSAcruts dataset (Karger & Zimmermann, 2018) served as a source of contemporary climate data : 1950-1960, in  is a very pessimistic scenario of climate change, which assumes that no mitigation measures will limit greenhouse gases emissions during the 21st century (Burgess et al., 2020;Schwalm et al., 2020). RCP 4.5, on the other hand, is a more optimistic scenario in which temperature rise starts to level-off before the end of the century (Thomson et al., 2011).
In addition to climate, we included soil variables as predictors in our SDMs. Specifically, we used soil pH in water (Hengl, 2018) and soil organic carbon content (Hengl & Wheeler, 2018), obtained from the OpenLandMap project (openlandmap.org) where they were estimated through machine-learning. Both variables are available at different soil depths; to account for the uncertainty and variation of earthworms' habitat, we averaged values of pH and soil organic carbon between the surface (0 cm) and 30-cm depth. Finally, we obtained land cover data from the CORINE project (https://land. coper nicus.eu/pan-europ ean/corin e-land-cover), which provides land cover maps of Europe derived from photo-interpretation of satellite images. Here, we used data from the year 1990 for training SDMs because it is the earliest dataset in the CORINE database. We projected models using land cover data for the year 2012 to represent present conditions, and using land cover in the year 2018 to represent future conditions because this is the last dataset in the CORINE database. All environmental predictors were cropped to the extent of Metropolitan France (mainland France and Corsica) and processed at a resolution of 30 arc-second, which corresponds to ca. 800 m in France.
In order to describe regional changes in diversity under the various modelling scenarios, we obtained a map of ecoregions from the World Wildlife Fund (WWF), which are biogeographical units describing distinct biota and habitats (Wikramanayake et al., 2002).
Fourteen ecoregions are present in France, mostly separating western from eastern France, as well as highlighting mountainous and Mediterranean areas ( Figure 1). Most of the country is included in the "Atlantic mixed forests" and "Western European broadleaf forests ecoregions." Results for these ecoregions will thus have a broader impact on earthworm communities at the country scale, compared to smaller ecoregions.

| Species distribution modelling
We fitted SDMs using boosted regression trees (BRTs), a statistical method suitable for training SDMs with presence/absence data, which can fit nonlinear relationships and handle interactions between predictors (Elith et al., 2008). There are several parameters that must be chosen in BRTs that control model complexity. To select the most appropriate settings, we fitted for each species several BRTs with all possible combinations of the following parameters: number of trees = from 50 to 500 in increments of 50, interaction depth = 1, 2, 3, 4, shrinkage = 0.005, 0.01, 0.05, 0.1, 0.15, 0.2. We retained among the set of models the one that yielded the highest predictive performance. The performance of models was evaluated through the area under the receiver operating curve (AUC) and the True Skill Statistics (TSS). AUC and TSS were estimated from partitioning data in spatial blocks following a double "checkerboard" grid of 2.5 arc-min and 2.5° resolution, producing evaluation bins that were spatially independent from the training data (Muscarella et al., 2014;Radosavljevic & Anderson, 2014). For later analyses, we kept only species with AUC > 0.7 and TSS > 0.4 to make predictions from reliably modelled species only.
We extracted variable importance using a permutation procedure. It was further used to reduce the number of predictors when possible: from the 12 variables initially included in the SDMs (eight bioclimatic variables + precipitation during LGM + two soil variables + land cover), those whose importance <10% were removed iteratively, and the new model was retained if its AUC was improved.
Models were projected in space over the entire area of Metropolitan France, and in time for the three periods for which we have climate data (present, 2050, 2070, in addition to the 1960s), assuming that soil did not change, and using land cover data for the year 2012 (present) and 2018 (future). Future projections were computed for the RCP 4.5 and RCP 8.5 scenarios, resulting in a total of six outputs in the form of suitability maps with values ranging between 0 and 1. Here, we considered two hypotheses regarding the dispersal of species. First, assuming full dispersal of species following their climatic niche, we did not transform model predictions. Second, assuming that species are unable to disperse within the time frame of analysis, we restricted future projections of species distributions to the locations already predicted as suitable under the conditions of sampling (1960s). For this, suitability maps for the 1960s were converted to binary outputs of predicted presence or absence using the threshold that maximizes the sum of sensitivity and specificity (Liu et al., 2013).
BRTs and parameter tuning were computed using the "SDMtune" R package (Vignali et al., 2020) and data partitioning used functions from the "ENMeval" R package (Muscarella et al., 2014). The entire procedure was reported as an ODMAP protocol (Methods S1) as recommended by Zurell et al. (2020).

| Macroecological constraints of species richness
Generally, stacked SDMs overestimate local species richness, because they ignore many local factors such as biotic interactions or dispersal limitations (D'Amen et al., 2018). Therefore, we imposed a limit to the number of species that can co-occur at a given loca-

| Functional diversity and post-hoc analyses
From the predicted species composition, we calculated within each grid cell and for each time period and scenario three complementary measures of functional diversity (Mason et al., 2005;Villéger et al., 2008), estimating functional distance with Gower's distance based on the matrix of species' traits. First, we computed functional richness (FRic), which corresponds to the volume occupied by a convex hull polygon drawn in the species' functional space. It is thus a measure of the total amount of trait variation within an assemblage. Second, we estimated the functional evenness of assemblages Functional diversity indices are known to be correlated with species richness (Mason et al., 2013). Therefore, in order to obtain measures of functional richness, evenness and divergence that are independent of species richness, we also calculated standardized effect sizes using a randomization procedure (ses.FRic, ses.FEve, ses.FDiv), where we randomized assemblages 500 times with an independent swap algorithm (Gotelli, 2000). Values of ses.FRic, ses.
FEve and ses.FDiv <0 can be interpretated as functional clustering while indices >0 correspond to assemblages that are functionally overdispersed (i.e. less-or more, respectively-variation in species traits than expected given the number of species). Functional diversity indices were calculated using the "FD" R package (Laliberte & Legendre, 2010), and randomizations were carried out with the "picante" R package (Kembel et al., 2010).
In order to illustrate changes in individual traits, we chose to focus on species' maximum weight, because it was strongly correlated to width and length, hence synthetizing species' size. Using As expected, FRic was positively correlated with species richness, especially in scenarios of no dispersal, although this relationship remained weak (R² < .15, Figure S3). Model projections predicted a strong increase (up to +90%) in FRic compared to the 1960s in some areas, especially in the "Italian sclerophyllous and semi-deciduous forests" ecoregion (Figures 3 and 4). However, at the country-scale and in the largest ecoregions ("Atlantic mixed forests" and "Western European broadleaf forests"), FRic is predicted to decline in the future (Table 1; Figure 4). Interestingly, in the no dispersal hypothesis, we predicted that FRic will concentrate in a small region of northeastern France (Figure 3). FRic was clustered in the 1960s (ses.Fric

<0) in all regions of France except in both Corsican ecoregions and
in a narrow strip connecting the Atlantic and Mediterranean coasts where ses.FRic was >0 ( Figure 5). Models predict in future climate a decrease in ses.Fric in most ecoregions, including the largest ones, leading to a functional richness that will be even more clustered compared to the present and the 1960s (Table 1; Figure 6).
Generally, we observed that ses.Fric is predicted to become more homogeneous across France in the future, especially in the "no dispersal" scenario ( Figure 6).
Functional evenness was not related to species richness in the 1960s, but we predicted that they will become slightly negatively correlated in future scenarios ( Figure S3). Changes in functional evenness always remained <±15%, and models predicted that FEve will be mostly stable, or even slightly increasing, in the largest ecoregions ("Atlantic mixed forests" and "Western European broadleaf forests") ( Table 1; Figures 3 and 4). We noted, however, that the strongest changes were predicted to occur in the southern FEve <0) across the country, except for a strip characterized by overdispersion (ses.FEve >0) in southwestern France, located in the "Atlantic mixed forests" ecoregion. There, as well as in the "Alps conifer and mixed forests" ecoregion where ses.FEve was also overdispersed in the 1960s, models predicted a decrease of ses.FEve that in some scenarios will lead to a switch towards clustered functional evenness (ses.FEve <0). In the south-eastern ecoregions, ses.FEve is predicted to increase (up to ses.FEve >0). Overall, ses.FEve is predicted to decrease substantially in the "no dispersal" hypothesis (Table 1), although in the largest ecoregion ("Western European broadleaf forests"), it is predicted to remain approximately stable in all scenarios ( Figure 6).
We observed a weak negative correlation between FDiv and species richness in the 1960s and in the present, which is predicted to become slightly positive in the future ( Figure S3). Generally, we F I G U R E 2 Earthworm species richness predicted by a macroecological model trained in the 1960s (the time of data collection), and projected in the present and in future climate following two scenarios of climate change and two time periods predicted that FDiv will differ only little compared to the 1960s, with changes < ±10% (Figures 3 and 4), although it will likely decrease compared to the 1960s, especially if species are unable to disperse (Table 1) (Figures 5 and 6). In the largest ecoregions ("Atlantic mixed Forests" and "Western European broadleaf forests"), models predicted either a stable or a slightly decreasing ses.FDiv in future climate. Overall, ses.FDiv is predicted to decrease more strongly in the "full dispersal" hypothesis ( Table 1). The most notable predicted changes were a strong increase in ses.FDiv in the "Alps conifer and mixed forests" ecoregion and a strong decline in the "Northeastern Spain and Southern France Mediterranean forests" ecoregion ( Figure 6).
Predictions of average species' weight indicated that in the 1960s, the heaviest earthworm communities were clustered in a small area in southern France ( Figure S4a). Future predictions revealed a shift towards heavier communities on average in all ecoregions in the strongest scenarios of climate change ( Figure S4b). In this case, though, the heaviest communities are predicted to remain in roughly the same, albeit larger, region of southern France ( Figure S4a).

| DISCUSS ION
Climate-driven range shifts are responsible for changes in the local composition of assemblages (Devictor et al., 2008;Parmesan, 2006),    potentially affecting the distribution of functional traits within biological communities (Wieczynski et al., 2019). In this study, we used a comprehensive sampling of earthworms carried out in metropolitan France to predict potential changes in the functional diversity of earthworms' assemblages in future climate. We observed that diversity is expected to have changed, both in terms of raw numbers and regarding its spatial distribution in France, between the period of sampling (1960s) and nowadays. Projections in scenarios of climate change predicted that this modification of the functional diversity of assemblages is expected to continue in the future, with stronger changes in the most extreme scenarios of climate change. There were substantial differences depending on the region considered and on whether species will be able to disperse and track climate change.
The first noticeable results we obtained concern the estimated accuracy of SDMs. We noticed a clear relationship between performance metrics and species' prevalence, a pattern that has already been described (Lobo et al., 2008) and that suggests that both AUC  and TSS provide a poor measure of the actual predictive performance of SDMs, even trained with presence-absence data obtained from high-quality sampling. In any case, we emphasize that the approach used in this study is designed to provide general trajectories under some specific hypotheses and scenarios, and should not be interpreted as a prediction of the true diversity of earthworm communities at a precise location. Given the performance and data quantity thresholds that we imposed to incorporate species into F I G U R E 6 Predicted changes in standardized functional diversity obtained from randomization (mean within each French ecoregion, coloured as in Figure 1), separated into standardized functional richness, standardized functional evenness and standardized functional divergence. Models were fitted with data from the 1960s, and projected in the present and in the future climate (2050 and 2070) following two scenarios of climate change (RCP 4.5: plain lines; RCP 8.5: dotted lines). For each metrics, projections either assume that species can or cannot disperse, while the maximum number of co-occurring species is limited following a macroecological model of species richness our predictions of functional diversity, only a handful of earthworm species (30/105) were actually included in the analyses. Therefore, it must be kept in mind that the true diversity of earthworm communities that will be present in the field may differ strongly from our predictions. It is also important to acknowledge that we most likely missed several important predictors of earthworms' distributions because they are relevant at a finer scale than our study.
Incorporating microclimate and chemical conditions within the soil layer, for example, would have certainly improved model performance, although such variables are difficult to map at fine scale and to project in the future.
We fitted here SDMs using climate, soil and land cover variables, with the drawback that we did not have future projections of soil pH, soil organic carbon and land cover changes. Since earthworms are soil-dwelling organisms, that moreover are ecosystem engineers that contribute to many soil processes, we know that feedbacks exist between earthworms' activity and soil structure (Blouin et al., 2013;Liu et al., 2019). This could have resulted in an apparent strong relationship between soil variables and species' presence/absence, making it difficult to infer the effect of climate change on their distributions. Here, we observed instead that most species were associated with variations in land cover or climate variables. It is probable that a relationship between soil properties and the distribution of earthworm species could have been detected at a smaller scale. For example, it is known that soil pH can be a limiting factor for the colonization of some earthworm species (e.g. Chan & Mead, 2003). There is also evidence that soil environmental variability-including pH and soil organic carbon concentration-contributes to earthworms' distribution at a scale of a few meters (Jiménez et al., 2011), far finer that the 30 arc-second resolution of our country-wide variables.
The fact that land cover appeared as the primary driver of earthworms' distribution in our models is consistent with the results of Rutgers et al. (2016), who observed that earthworm abundance and diversity were related to the presence or absence of some land cover types (grasslands, croplands, forests, heathlands, vineyards) across Europe. This reflects the habitat preference of earthworm species, as well as the fact that land use intensity is known to have a large impact on earthworm communities (Smith et al., 2008;Spurgeon et al., 2013). Currently, this variable accounted for country-scale variation in habitat types and anthropogenic activity. Our predictions may have been refined by incorporating models of land use change under various socioeconomic and climate scenarios; however, such projections are currently available at a coarser scale than the predictors we employed (Chen et al., 2020;Hurtt et al., 2020).
The most influential climate variables (BIO3, PET of the coldest quarter) revealed the importance of temperature for determining earthworms' distribution at the scale of France. An increase in the frequency and intensity of extreme temperature events is one of the most often observed (Perkins-Kirkpatrick & Lewis, 2020) and predicted (Meehl & Tebaldi, 2004) consequence of climate change. Similarly, precipitation extremes are predicted to intensify in response to a warming climate (O'Gorman, 2015). Therefore, it seems likely that climate change will play (and probably have already played) a role in the distribution of earthworm species and hence in the structure of communities. The projections we produced confirm that, providing that the association between earthworm species and climate we detected is legit, climate change has the potential to be responsible for strong shifts in the functional diversity of earthworm assemblages.
Besides abiotic factors, it is evident that species interactions are partly responsible for the structure of earthworm communities (Uvarov, 2009). However, again, this effect must be limited in space to the fine scale relevant to processes of interspecific competition. Looking at the scale of whole species ranges, large-scale factors such as climate are usually recognized as the primary drivers shaping range limits (Soberon & Nakamura, 2009). In this regard, a number of studies have reported an effect of temperature variation on earthworms' activity and survival (Singh et al., 2019), and a global mapping of earthworm diversity also concluded that climate variables were more important in shaping earthworm communities than soil properties (Phillips et al., 2019). Therefore, we are confident that our models provide useful projections of the potential impact of climate change on the diversity of French earthworm's assemblages.
However, these predictions must be understood at the relevant scale of analysis: they represent the expected diversity within 2.5 arc-min grid cells and cannot be used to predict community composition within a local sampling plot.
The predicted changes in diversity we obtained were partly dependent on the imposition of a dispersal limit or not. There were, on average, larger changes (especially declines) in functional diversity when we constrained species to remain in their historical ranges, most notably for functional richness. However, differences were rather limited, as evidenced by prediction maps that looked largely similar in both hypotheses (see Figures 3 and 5). The fact that we excluded from our analyses all species that were observed in less than 10 sites has possibly influenced the difference (or lack thereof) between full and no dispersal scenarios. By removing narrowly distributed and rare species, we could not identify a potential expansion of these species that could occur if conditions become more suitable for them in the future. Here, we simply assumed that species either moved freely to fill their entire niche, or on the contrary were so slow dispersers that they could not colonize new grid cells under the time frame of the study. Although active or passive dispersal definitely occurs in earthworms (see e.g. Caro et al., 2013), we know for sure that large-scale colonizations take time, as evidenced by the imprint glaciation has left of earthworm community structuring in France (Mathieu & Davies, 2014). Thus, the real changes in community composition at a scale of a few decades must be closer to the "no dispersal" scenario, even though knowing the actual dispersal ability of each species would certainly improve predictions.
We predicted a rather strong decrease in both species and FRic in the future, in most regions of France. In the most extreme scenario and in the absence of dispersal, climate change is expected to be responsible for a reduction of FRic of ca. 50% in the eastern half of the country ("Western European broadleaf forests" and "Alps conifer and mixed forests" ecoregions). Given the many effects earthworms have on ecosystem services (Blouin et al., 2013;Liu et al., 2019), this may profoundly change the functioning of soil processes in the areas affected by this decline, potentially with deleterious consequences on natural and agroecosystems. Experimental removal of earthworm species in mesocosms could provide empirical evidence of this effect in controlled conditions (e.g. Heemsbergen et al., 2004). Generally, climate change is expected to lead to a loss of rare species traits and to a functional homogenization. Loss of specialist species and functional homogenization of communities driven by climate change has already been empirically observed in other taxa (e.g. Fourcade et al., 2021;Gauzere et al., 2015). Functional differences among soil species is generally a factor that promotes biodiversity effects on many aspects of ecosystem functionality (Heemsbergen et al., 2004;Huang et al., 2020). Therefore, there is a risk that this loss of diversity result in the loss-or decrease-of some key aspects of ecosystem processes, for example because a key but rare trait will be missing from climate-altered earthworm communities. Since species and FRic are known to be correlated (Villéger et al., 2008), we produced standardized estimates that are independent of species richness. This revealed that FRic was mostly clustered in the 1960s, indicating that functionally similar species were generally present in the same grid cells (whether they actually co-occurred at the plot scale is something we could not estimate here). In most scenarios of climate change, FRic is predicted to become even more clustered, that is there will be even less variation of functional traits than expected given species richness.
Results for functional evenness and divergence give a different picture than richness. First of all, it appears that these indices of diversity will be less affected by climate change, since the predicted changes in the present and future conditions, compared to the 1960s, remain limited. We observe that, in the largest French ecoregions, functional evenness may slightly increase in the scenarios of strongest climate change (RCP 8.5 in 2070), while it is predicted to decrease in mountainous regions. On the contrary, FDiv should decrease in most of the country, except in the mountains where an increase in future climate is predicted. In the latter regions, ses.FDiv is also predicted to increase, showing that, even in a declining functional richness, climate change will lead in mountainous areas to a greater dispersion of functional traits around the mean trait values.
Since higher FDiv is a sign of greater niche differentiation, it may contribute to increase ecosystem functioning thanks to a more efficient use of resources (Mason et al., 2005). Standardized functional evenness should remain largely stable in all ecoregions, showing that despite changes in species and functional richness, functional space is predicted to remain evenly utilized.
The predicted changes in diversity we describe in this study are not spatially homogeneous across France but exhibit clear regional differences. Some of this heterogeneity corresponds to known patterns in climate change ecology. Here, we predicted a clear reduction in species richness in southern France, which likely reflects the retraction of species' ranges at their southern edge (Lenoir & Svenning, 2014). Training models at a larger scale could have allowed to identify the colonization of climate-tracking species coming from southern Europe (Chen et al., 2011;Devictor et al., 2008;Parmesan, 2006), but only under the hypothesis that earthworms are able to perform long-distance range shifts in this time frame. We also observed that some mountainous regions show a slower decrease in species richness than in lowlands (see e.g. the Alps and Massif Central in Figure 2). This is typical of climate refugia, which are areas where climatic conditions will remain relatively colder than elsewhere, even in a global trajectory of warming. These regions have the potential to buffer the effects of climate change because they allow the persistence or the colonization of species that would have been extirpated otherwise (Morelli et al., 2020). Since declines in species richness drive the large changes we predicted for functional richness, it is important to consider these climate refugia for predicting future functional diversity. We note here in this regard that the "Alps conifer and mixed forests" ecoregion is the only one where standardized FDiv is predicted to increase in the future, despite a general trend of decreasing species and functional richness. There are other possible climate refugia that are not accounted for here, though. Colder microclimates may be created by the interaction between landscape structure, biological processes, topographical features and weather (Lembrechts et al., 2020), in such a way that the temperature conditions experienced at fine scale may be different from those observed from macroclimate variables. Therefore, there probably exist microclimatic refugia such as in shady landscapes or north-oriented slopes where species may persist longer than predicted by our models. Similarly, soil conditions may provide a buffering effect from macroclimatic warming that increases with depth; in this case, endogeic species may be less sensitive to climate change that species that live closer to the surface.

| CON CLUS ION
In conclusion, we showed here how the functional diversity of a key group of ecosystem engineers may change in the future as a consequence of climate change. Our results mostly pointed to a strong potential reduction of species and functional richness.
Although we were able to model only ca. one third of the species present in the country, it is likely that earthworm communities will be significantly altered by climate change in terms of their richness, size structure and distribution of traits. Since distinct functional groups of earthworms differ in their activity in the soil, there are concerns that changes in community structure may lead to an alteration of soil processes and of the ecosystem services they provide (Blouin et al., 2013;Heemsbergen et al., 2004). However, we also predicted that, relatively to species richness, standardized functional diversity indices may increase in some specific regions of France. In addition to the observation that functional evenness and divergence will remain relatively stable overall, this indicates that declines in species and FRic may be partially compensated.
The approach of using stacked SDMs coupled with species' traits to predict functional diversity in a changing climate is not new (see e.g. D'Amen et al., 2018; Pradervand et al., 2014), but this is the first time it is applied to belowground taxa. Here, we used data from the 1960s (Bouché, 1972) and projected models in time, including in the present. Field-validation of our models may be possible in the future thanks to a resampling of Bouché's (1972) sites that is currently underway.

ACK N OWLED G EM ENT
The authors thank Jérôme Mathieu for having compiled Bouché's data and having made it available for future research. The study did not benefit from any financial support.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.

DATA AVA I L A B I LT Y S TAT E M E N T
All data used in this study are freely available online from the