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

Yoan Fourcade, Mathias Vercauteren

► To cite this version:

Yoan Fourcade, Mathias Vercauteren. Predicted changes in the functional structure of earthworm assemblages in France driven by climate change. Diversity and Distributions, 2022, 28 (5), pp.1050-1066. 10.1111/ddi.13505. hal-03607878

HAL Id: hal-03607878

https://hal.u-pec.fr/hal-03607878

Submitted on 14 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Predicted changes in the functional structure of earthworm assemblages in France driven by climate change

Yoan Fourcade ID |

Mathias Vercauteren

Univ. Paris Est Creteil, Sorbonne Université, CNRS, INRAE, IRD, Université de Paris, Institute of Ecology and Environmental Sciences Paris iEES, Créteil, France

**Correspondence**

Yoan Fourcade, Univ. Paris Est Creteil, Sorbonne Université, CNRS, INRAE, IRD, Université de Paris, Institute of Ecology and Environmental Sciences Paris iEES, 94010 Créteil, France.

Email: yoan.fourcade@u-pec.fr

Editor: Andrew Barnes

**Abstract**

**Aim:** 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.

**Location:** Metropolitan France.

**Methods:** We fitted boosted regression trees to earthworm data collected using a standardized protocol across France in the 1960s. We used model projections constrained by a macroecological model of species richness to predict the composition of earthworm assemblages in the present and in two scenarios of climate change and two future time periods. We coupled these results with a large set of species traits to calculate predicted changes in functional diversity, which we summarized by ecoregion.

**Results:** Models predicted a clear decline in functional richness between the period of sampling and nowadays which are expected to continue in the future, with substantial differences depending on ecoregions and on whether species will be able to disperse or not. However, predicted changes in functional evenness and divergence are much weaker, suggesting that climate change will not affect all facets of functional diversity in the same way.

**Main Conclusions:** Our results mostly pointed to a potential reduction of the functional richness of earthworm communities in the future, but this predicted loss of diversity could be weaker if species are able to colonize new suitable sites or to persist in microclimate refugia. There are concerns, though, that these changes lead to an alteration of soil processes and of the ecosystem services they provide.

**KEYWORDS**

climate change, earthworms, ecological niche modelling, functional divergence, functional evenness, functional richness, species distribution modelling, species traits
# 1 | INTRODUCTION

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 (Cardinale et al., 2012; 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 (Capowiez et al., 2014), 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 (Huang et al., 2020). 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.

## 2 | METHODS

### 2.1 | 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² 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

<figure><img src="image_2.png" /><figcaption>**FIGURE 1** Location of sampling sites of earthworms conducted by Bouché (1972) across France, with the 14 French ecoregions represented as coloured polygons</figcaption></figure>
in a way that is likely relevant for ecosystem functioning. Site locations, presence/absence and abundance of earthworms, as well as species traits, were all compiled by Mathieu and Davies (2014), and are available from the Dryad repository (https://doi.org/10.5061/dryad.g7046).

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 order to describe climatic conditions during the sampling period, and 2000-2016 to describe current climate. In addition, we downloaded future projections of the same variables for the years 2041-2060 and 2061-2080 (respectively referred to as 2050 and 2070 later on) according to two representative concentration pathways (RCP 4.5 and RCP 8.5), modelled under the MIROC5 global circulation model. RCP 8.5 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.copernicus.eu/pan-european/corine-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.

## 2.2 | 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).

## 2.3 | 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 location, using a macroecological model of species richness (MEM). For this, we fitted a BRT model to the number of species observed at each sampling point, using the same variables as for SDMs but with a Poisson distribution. This model was then projected in space and time in the same way as the SDMs, providing spatially explicit predictions of species richness in the 1960s, in the present and in future scenarios of climate change. Predicted species richness was rounded to the upper integer to ensure that all grid cells have at least one species. We estimated the predictive ability of the MEM by calculating Spearman's correlation between observed and predicted species richness.

We combined the outputs of MEM and SDMs to predict the identity of co-occurring species. Because the calculations of functional diversity were computationally intensive, all species distribution and species richness maps were aggregated at a resolution of 2.5 arc-min, and later projected to ETRS89-LAEA coordinate system (EPSG:3035). In each time period and scenario, we assumed that SDMs define the potential species pool constrained by abiotic drivers, and combined SDM predictions by applying a probability ranking rule such that the number of species was equal to MEM predictions, and species were selected based on decreasing suitability (D'Amen et al., 2015, 2018).

## 2.4 | 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 (FEve), which measures if species are regularly distributed within functional space or if there exist clusters of species sharing similar traits. Third, functional divergence (FDiv) is an estimate of whether species are mostly located at the centre or at the edge of functional space. Altogether, these three indices describe the functional diversity of earthworm assemblages in all its aspects.

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 interpreted 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 “picate” 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 SDM predictions constrained by the MEM, we calculated for each time period and scenario the average maximum weight of all the species predicted to be present in each grid cell. This is analogous to a community-weighted mean where all species are given the same weight, since we have no information about their respective abundance. This way, we described changes in the mean size of earthworm communities in different scenarios of climate change compared to the size structure of communities in the 1960s.

We overlapped species richness and functional diversity maps with the map of French ecoregions, and extracted the mean values within each ecoregion, in addition to calculating mean values across the whole country. For the raw functional indices, we calculated the average change (in %) between values predicted in the 1960s (the period of sampling) and each temporal projection of diversity. This way, we estimated how much our models predicted changes in earthworm's functional diversity in each ecoregion of Metropolitan France, between the 1960s and the present, and between the 1960s and different scenarios of future climate change. For standardized functional indices, we analysed directly the values of ses.FRic, ses.
FEve and ses.FDiv estimated at each time period in order to test whether some regions are predicted to experience a shift from functional clustering to overdispersion, or vice versa.

## 3 | RESULTS

Species distribution models varied in predictive performance, with a mean AUC of ca. 0.796 (minimum = 0.610 and maximum = 0.995), and a mean TSS of ca. 0.589 (minimum = 0.223 and maximum = 0.985) (Figure S1a-c). In total, 14 species were excluded from the analyses of functional diversity because they did not reach the thresholds of AUC = 0.7 and TSS = 0.4. There was a clear negative link between predictive performance and species' prevalence (Figure S1d,e). Overall, the most important variable during SDM fitting was land cover, which had an average permutation importance of 19% (SEM = 1.67) and was ranked as the most important variable for 16 species out of 44. The second most important variable appeared to be BIO 3 (isothermality), which had a mean importance of 18% (SEM = 3.17) and ranked first for eight species. It was followed by potential evapotranspiration during the coldest quarter (mean importance = 14.5, SEM = 2.91, ranked first for six species; see Figure S2).

Predictions from the macroecological model of species richness were correlated with observed richness (Spearman's correlation: $\rho$ = 0.47, $p$ < .001) and errors were centred on zero (mean = 0.003), which is similar to what has been found in comparable studies (D'Amen et al., 2015, 2018). Predicted species richness in the 1960s ranged from 2 to 6 species per grid cell, with a mean richness of 4.27 species. Species richness was generally lower along the Atlantic coast and in the Mediterranean region, with maximum values along the English Channel (Figure 2). Compared to the 1960s, the MEM predicted in the present and in all future scenarios a decrease in species richness overall, especially in the regions that already exhibited lower richness (western half of France and Mediterranean region) (Figure 2).

As expected, FRic was positively correlated with species richness, especially in scenarios of no dispersal, although this relationship remained weak ($R^2$ < .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 north-eastern 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 and mountainous ecoregions: an increase of FEve in both Corsican ecoregions, while a decrease of FEve is predicted in all other southern ecoregions, especially in the RCP 8.5 climate change scenario (Figure 4). Standardized functional evenness was distributed very similarly as ses.Fric in the 1960s (Figure 5): mostly clustered (ses. 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 co-nifer 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

<figure><img src="image_3.png" /><figcaption>FIGURE 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</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 3 Earthworm functional diversity predicted by species distribution models 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. Three complementary descriptors of functional diversity were calculated within each 2.5 arc-min grid cell: (a) functional richness, (b) functional evenness and (c) functional divergence. For each metrics, projections either assume that species fill their predicted climatic niche (top row) or that they cannot disperse (bottom row), while the maximum number of co-occurring species is limited following the macroecological model of species richness shown in Figure 2</figcaption></figure>

predicted that FDiv will differ only little compared to the 1960s,
with changes < ±10% (Figures 3 and 4), although it will likely de-
crease compared to the 1960s, especially if species are unable to
disperse (Table 1). Interestingly, models predicted an increase in
FDiv in the mountainous ecoregions ("Alps conifer and mixed for-
ests" and "Pyrenees conifer and mixed forests") in almost all disper-
sal and climate change scenarios. Mediterranean ecoregions show
a clear tendency for a decrease in FDiv in future climate. Changes

in the largest ecoregions, though, are more uncertain and will likely
remain of limited magnitude (Figure 4). In the 1960s, standardized
FDiv showed evidence of clustering (ses.FDiv <0) in the western
and Mediterranean areas ("Western European broadleaf forests,"
"Northeastern Spain and Southern France Mediterranean forests"
and "Italian sclerophyllous and semi-deciduous forests" ecoregions),
while eastern ecoregions were generally characterized by ses.FDiv
>0 ("Alps conifer and mixed forests" and "Atlantic mixed Forests;"
<figure><img src="image_5.png" /><figcaption>FIGURE 4 Predicted changes in functional diversity, separated into FRic (a), functional evenness (b) and functional divergence (c), compared to the 1960s. Values are presented as mean percent change within each French ecoregion, coloured as in Figure 1</figcaption></figure>

ecoregions) (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).

4 | DISCUSSION

Climate-driven range shifts are responsible for changes in the local composition of assemblages (Devictor et al., 2008; Parmesan, 2006),
<table><caption>TABLE 1 Summary of the mean values (and standard deviation) of all diversity metrics computed across France</caption><thead><tr><th rowspan="2">Period/Model</th><th colspan="2">Species richness</th><th colspan="2">Functional richness</th><th colspan="2">Functional evenness</th><th colspan="2">Functional divergence</th></tr><tr><th>Mean</th><th>SD</th><th>Raw Mean</th><th>SD</th><th>Raw Mean</th><th>SD</th><th>Raw Mean</th><th>SD</th></tr></thead><tbody><tr><td colspan="9"><strong>Full dispersal</strong></td></tr><tr><td>Historical (1960s)</td><td>4.190</td><td>0.671</td><td>0.135</td><td>0.077</td><td>-0.014</td><td>0.086</td><td>0.814</td><td>0.064</td><td>-0.001</td><td>-0.001</td><td>0.102</td></tr><tr><td>Present</td><td>3.760</td><td>0.654</td><td>0.127</td><td>0.087</td><td>-0.016</td><td>0.098</td><td>0.810</td><td>0.076</td><td>0.005</td><td>0.000</td><td>0.101</td></tr><tr><td>RCP 4.5–2050</td><td>3.460</td><td>0.550</td><td>0.119</td><td>0.087</td><td>-0.019</td><td>0.090</td><td>0.808</td><td>0.083</td><td>-0.008</td><td>-0.014</td><td>0.103</td></tr><tr><td>RCP 4.5–2070</td><td>3.360</td><td>0.551</td><td>0.109</td><td>0.082</td><td>-0.019</td><td>0.085</td><td>0.814</td><td>0.091</td><td>-0.005</td><td>-0.011</td><td>0.108</td></tr><tr><td>RCP 8.5–2050</td><td>3.550</td><td>0.644</td><td>0.109</td><td>0.082</td><td>-0.029</td><td>0.085</td><td>0.813</td><td>0.078</td><td>-0.010</td><td>-0.016</td><td>0.102</td></tr><tr><td>RCP 8.–2070</td><td>3.330</td><td>0.568</td><td>0.105</td><td>0.083</td><td>-0.025</td><td>0.088</td><td>0.817</td><td>0.092</td><td>0.000</td><td>-0.016</td><td>0.106</td></tr><tr><td colspan="9"><strong>No dispersal</strong></td></tr><tr><td>Present</td><td>0.135</td><td>0.150</td><td>-0.014</td><td>0.097</td><td>0.810</td><td>0.079</td><td>0.004</td><td>0.103</td><td>-0.007</td><td>-0.012</td><td>0.102</td></tr><tr><td>RCP 4.5–2050</td><td>0.114</td><td>0.138</td><td>-0.026</td><td>0.086</td><td>0.808</td><td>0.089</td><td>-0.007</td><td>0.104</td><td>-0.009</td><td>-0.001</td><td>0.103</td></tr><tr><td>RCP 4.5–2070</td><td>0.105</td><td>0.137</td><td>-0.028</td><td>0.079</td><td>0.812</td><td>0.096</td><td>-0.009</td><td>0.106</td><td>-0.013</td><td>-0.001</td><td>0.111</td></tr><tr><td>RCP 8.5–2050</td><td>0.112</td><td>0.146</td><td>-0.036</td><td>0.081</td><td>0.812</td><td>0.082</td><td>-0.013</td><td>0.102</td><td>-0.006</td><td>-0.006</td><td>0.106</td></tr><tr><td>RCP 8.5–2070</td><td>0.096</td><td>0.135</td><td>-0.035</td><td>0.081</td><td>0.813</td><td>0.097</td><td>-0.004</td><td>0.105</td><td>-0.009</td><td>-0.009</td><td>0.108</td></tr></tbody></table>

Note: Functional diversity indices were computed assuming unlimited dispersal of species (top) or keeping species within their historical range (bottom), and we report both their raw and standardized values. Within each column, highest values are presented in green and lowest values in red. Functional diversity indices were computed assuming unlimited dispersal of species (top) or keeping species within their historical range (bottom), and we report both their raw and standardized values. Within each column, values are presented from green (highest) to red (lowest), with intermediate values appearing, in decreasing order, from light green to yellow to orange.
<figure><img src="image_6.png" /><figcaption>FIGURE 5 Earthworm-standardized functional diversity indices [(a) functional richness, (b) functional evenness, (c) functional divergence] predicted by species distribution models in the 1960s, and projected in the present and in future climate following two scenarios of climate change and two time periods. Compared with Figure 3, these measures are independent from species richness and can be interpreted as functional clustering (values <0) or overdispersion (values >0). For each metrics, projections either assume that species fill their predicted climatic niche (top row) or that they cannot disperse (bottom row), while the maximum number of co-occurring species is limited following the macroecological model of species richness shown in Figure 2</figcaption></figure>

potentially affecting the distribution of functional traits within bio-
logical communities (Wieczynski et al., 2019). In this study, we used
a comprehensive sampling of earthworms carried out in metropoli-
tan France to predict potential changes in the functional diversity of
earthworms' assemblages in future climate. We observed that diver-
sity 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 perfor-
mance metrics and species' prevalence, a pattern that has already
been described (Lobo et al., 2008) and that suggests that both AUC
<figure><img src="image_7.png" /><figcaption>FIGURE 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</figcaption></figure>

and TSS provide a poor measure of the actual predictive perfor-
mance of SDMs, even trained with presence-absence data obtained
from high-quality sampling. In any case, we emphasize that the ap-
proach 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
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 com-
munities 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' distri-
butions 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 perfor-
mance, 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 dis-
tributions. Here, we observed instead that most species were asso-
ciated 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 colo-
nization 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' dis-
tribution 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 earth-
worms' 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 im-
pact 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 projec-
tions 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 cold-
est quarter) revealed the importance of temperature for determin-
ing 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 communi-
ties (Uvarov, 2009). However, again, this effect must be limited in
space to the fine scale relevant to processes of interspecific com-
petition. 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 vari-
ables 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 cli-
mate 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 compo-
sition within a local sampling plot.

The predicted changes in diversity we obtained were partly de-
pendent 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 ex-
cluded from our analyses all species that were observed in less than
10 sites has possibly influenced the difference (or lack thereof) be-
tween full and no dispersal scenarios. By removing narrowly distrib-
uted 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 defi-
nitely 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 commu-
nity 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.

## 5 | CONCLUSION

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, in-
cluding in the present. Field-validation of our models may be pos-
sible in the future thanks to a resampling of Bouché's (1972) sites
that is currently underway.

ACKNOWLEDGEMENT

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.

CONFLICT OF INTEREST

The authors declare that they have no conflict of interest.

DATA AVAILABILTY STATEMENT

All data used in this study are freely available online from the
CHELSA website (climate variables, https://chelsa-climate.org/),
from the Zenodo repository (soil pH, https://zenodo.org/recor
d/2525664; soil organic carbon content, https://zenodo.org/recor
d/2525553), from the Copernicus website (CORINE land cover,
https://land.copernicus.eu/pan-european/corine-land-cover)
and from the Dryad repository (earthworm data, https://doi.
org/10.5061/dryad.g7046).

ORCID

Yoon Fourcade https://orcid.org/0000-0003-3820-946X

REFERENCES

Bardgett, R. D., & van der Putten, W. H. (2014). Belowground biodiver-
sity and ecosystem functioning. *Nature*, *515*, 505–511. https://doi.
org/10.1038/nature13855

Biber, M. F., Voskamp, A., Niamir, A., Hickler, T., & Hof, C. (2019). A comparison of macroecological and stacked species distribution models to predict future global terrestrial vertebrate richness. *Journal of Biogeography*, *47*(1), 114–129. https://doi.org/10.1111/jbi.13696

Blouin, M., Hodson, M. E., Delgado, E. A., Baker, G., Brussaard, L., Butt, K. R., Dai, J., Dendooven, L., Peres, G., Tondoh, J. E., Cluzeau, D., & Brun, J.-J. (2013). A review of earthworm impact on soil function and ecosystem services: Earthworm impact on ecosystem services. *European Journal of Soil Science*, *64*, 161–182. https://doi.org/10.1111/ejss.12025

Booth, T. H., Nix, H. A., Hutchinson, M. F., & Busby, J. R. (1987). Grid matching: A new method for homoclime analysis. *Agricultural and Forest Meteorology*, *39*, 241–255. https://doi.org/10.1016/0168-1923(87)90041-4

Bottinelli, N., & Capowiez, Y. (2021). Earthworm ecological categories are not functional groups. *Biology and Fertility of Soils*, *57*, 329–331. https://doi.org/10.1007/s00374-020-01517-1

Bouché, M. B. (1972). *Lombriciens de France: écologie et systématique*.
INRA.

Bouché, M. B. (1977). Stratégies lombriciennes. *Ecological Bulletins*, *25*, 122-132

Briones, M. J. I., & Álvarez-Otero, R. (2018). Body wall thickness as a po-
tential functional trait for assigning earthworm species to ecolog-
ical categories. *Pedobiologia*, *67*, 26–34. https://doi.org/10.1016/j.
pedobi.2018.02.001

Burgess, M. G., Ritchie, J., Shapland, J., & Pielke, R. (2020). IPCC baseline scenarios have over-projected CO₂ emissions and economic growth. *Environmental Research Letters*, *16*, 014016. https://doi.org/10.1088/1748-9326/abcdd2

Cadotte, M. W., Carscadden, K., & Mirotnick, N. (2011). Beyond species: Functional diversity and the maintenance of ecological processes and services. *Journal of Applied Ecology*, *48*, 1079–1087. https://doi.org/10.1111/j.1365-2664.2011.02048.x

Cameron, E. K., Martins, I. S., Lavelle, P., Mathieu, J., Tedersoo, L., Gottschall, F., Guerra, C. A., Hines, J., Patoine, G., Siebert, J., Winter, M., Cesarz, S., Delgado-Baquerizo, M., Ferlian, O., Fierer, N., Kreft, H., Lovejoy, T. E., Montanarella, L., Orgiazzi, A., Eisenhauer, N. (2018). Global gaps in soil biodiversity data. *Nature Ecology & Evolution*, *2*, 1042–1043. https://doi.org/10.1038/s41559-018-0573-8

Capowiez, Y., Sammartino, S., & Michel, E. (2014). Burrow systems of endogeic earthworms: Effects of earthworm abundance and con-sequences for soil water infiltration. *Pedobiologia*, *57*, 303–309. https://doi.org/10.1016/j.pedobi.2014.04.001

Cardinale, B. J., Duffy, J. E., Gonzalez, A., Hooper, D. U., Perrings, C., Venail, P., Narwani, A., Mace, G. M., Tilman, D., Wardle, D. A., Kinzig, A. P., Daily, G. C., Loreau, M., Grace, J. B., Larigauderie, A., Srivastava, D. S., & Naeem, S. (2012). Biodiversity loss and its impact on humanity. *Nature*, *486*, 59–67. https://doi.org/10.1038/nature11148

Carmona, C. P., Tamme, R., Pärtel, M., de Bello, F., Brosse, S., Capdevila, P., González-M, R., González-Suárez, M., Salguero-Gómez, R., Vásquez-Valderrama, M., & Toussaint, A. (2021). Erosion of global functional diversity across the tree of life. *Science Advances*, *7*, eabf2675. https://doi.org/10.1126/sciadv.abf2675

Caro, G., Decaëns, T., Lecarpentier, C., & Mathieu, J. (2013). Are dispersal behaviours of earthworms related to their functional group? *Soil Biology and Biochemistry*, *58*, 181–187. https://doi.org/10.1016/j.soilbio.2012.11.019

Ceballos, G., Ehrlich, P. R., & Dirzo, R. (2017). Biological annihilation via the ongoing sixth mass extinction signaled by vertebrate population losses and declines. *Proceedings of the National Academy of Sciences of the United States of America*, *114*(30), E6089–E6096. https://doi.org/10.1073/pnas.1704949114

Chan, K.-Y., & Mead, J. A. (2003). Soil acidity limits colonisation by *Aporrectodea trapezoides*, an exotic earthworm. *Pedobiologia*, *47*, 225–229. https://doi.org/10.1078/0031-4056-00186

Chen, I.-C., Hill, J. K., Ohlemüller, R., Roy, D. B., & Thomas, C. D. (2011). Rapid range shifts of species associated with high levels of climate warming. *Science*, *333*, 1024–1026. https://doi.org/10.1126/scien ce.1206432

Chen, M., Vernon, C. R., Graham, N. T., Hejazi, M., Huang, M., Cheng, Y., & Calvin, K. (2020). Global land use for 2015–2100 at 0.05° resolution under diverse socioeconomic and climate scenarios. *Scientific Data*, *7*, 320. https://doi.org/10.1038/s41597-020-00669-x

Chen, S.-Y., Hsu, C.-H., & Soong, K. (2021). How to cross the sea: Testing the dispersal mechanisms of the cosmopolitan earthworm *Pontodrilus litoralis*. *Royal Society Open Science*, *8*, 202297. https://doi.org/10.1098/rsos.202297

Costa, G. C., Nogueira, C., Machado, R. B., & Colli, G. R. (2009). Sampling bias and the use of ecological niche modeling in conservation planning: A field evaluation in a biodiversity hotspot. *Biodiversity and Conservation*, *19*, 883–899. https://doi.org/10.1007/s1053 1-009-9746-8

D'Amen, M., Dubuis, A., Fernandes, R. F., Pottier, J., Pellissier, L., & Guisan, A. (2015). Using species richness and functional traits predictions to constrain assemblage predictions from stacked species distribution models. *Journal of Biogeography*, *42*, 1255–1266. https://doi.org/10.1111/jbi.12485
D'Amen, M., Mateo, R. G., Pottier, J., Thuiller, W., Maiorano, L., Pellissier, L., Ndiribe, C., Salamin, N., & Guisan, A. (2018). Improving spatial predictions of taxonomic, functional and phylogenetic diversity. *Journal of Ecology*, **106**, 76-86. https://doi.org/10.1111/1365-2745.12801

Daniel, O., Kohli, L., & Bieri, M. (1996). Weight gain and weight loss of the earthworm *Lumbricus terrestris* L. at different temperatures and body weights. *Soil Biology and Biochemistry*, **28**, 1235-1240. https://doi.org/10.1016/0038-0717(96)00121-6

de Bello, F., Lavorel, S., Díaz, S., Harrington, R., Cornelissen, J. H. C., Bardgett, R. D., Berg, M. P., Cipriotti, P., Feld, C. K., Hering, D., Martins da Silva, P., Potts, S. G., Sandin, L., Sousa, J. P., Storkey, J., Wardle, D. A., & Harrison, P. A. (2010). Towards an assessment of multiple ecosystem processes and services via functional traits. *Biodiversity and Conservation*, **19**, 2873-2893. https://doi.org/10.1007/s10531-010-9850-9

De Wandeler, H., Sousa-Silva, R., Ampoorter, E., Bruelheide, H., Carnol, M., Dawud, S. M., Dänilä, G., Finer, L., Hättenschwiler, S., Hermy, M., Jaroszewicz, B., Joly, F.-X., Müller, S., Pollastrini, M., Ratcliffe, S., Raulund-Rasmussen, K., Selvi, F., Valladares, F., Van Meerbeek, K., ... Muys, B. (2016). Drivers of earthworm incidence and abundance across European forests. *Soil Biology and Biochemistry*, **99**, 167-178. https://doi.org/10.1016/j.soilbio.2016.05.003

Devictor, V., Julliard, R., Couvet, D., & Jiguet, F. (2008). Birds are tracking climate warming, but not fast enough. *Proceedings of the Royal Society B: Biological Sciences*, **275**(1652), 2743-2748. https://doi.org/10.1098/rspb.2008.0878

Devictor, V., Mouillot, D., Meynard, C., Jiguet, F., Thuiller, W., & Mouquet, N. (2010). Spatial mismatch and congruence between taxonomic, phylogenetic and functional diversity: The need for integrative conservation strategies in a changing world. *Ecology Letters*, **13**, 1030-1040. https://doi.org/10.1111/j.1461-0248.2010.01493.x

Dormann, C. F., Elith, J., Bacher, S., Buchmann, C., Carl, G., Carré, G., Marquéz, J. R. G., Gruber, B., Lafourcade, B., Leitão, P. J., Münkemüller, T., McClean, C., Osborne, P. E., Reineking, B., Schröder, B., Skidmore, A. K., Zurell, D., & Lautenbach, S. (2013). Collinearity: A review of methods to deal with it and a simulation study evaluating their performance. *Ecography*, **36**, 27-46. https://doi.org/10.1111/j.1600-0587.2012.07348.x

Elith, J., & Leathwick, J. R. (2009). Species distribution models: Ecological explanation and prediction across space and time. *Annual Review of Ecology Evolution and Systematics*, **40**, 677-697. https://doi.org/10.1146/annurev.ecolsys.110308.120159

Elith, J., Leathwick, J. R., & Hastie, T. (2008). A working guide to boosted regression trees. *Journal of Animal Ecology*, **77**, 802-813. https://doi.org/10.1111/j.1365-2656.2008.01390.x

Eriksen-Hamel, N. S., & Whalen, J. K. (2006). Growth rates of *Aporrectodea caliginosa* (Oligochaetae: Lumbricidae) as influenced by soil temperature and moisture in disturbed and undisturbed soil columns. *Pedobiologia*, **50**, 207-215. https://doi.org/10.1016/j.pedobi.2005.10.008

Lafourcade, Y., Åström, S., & Öckinger, E. (2021). Decline of parasitic and habitat-specialist species drives taxonomic, phylogenetic and functional homogenization of sub-alpine bumblebee communities. *Oecologia*, **196**, 905-917. https://doi.org/10.1007/s00442-021-04970-3

Galego de Oliveira, A., Bailly, D., Cassemiro, F. A. S., Couto, E. V. D., Bond, N., Gilligan, D., Rangel, T. F., Agostinho, A. A., & Kennard, M. J. (2019). Coupling environment and physiology to predict effects of climate change on the taxonomic and functional diversity of fish assemblages in the Murray-Darling Basin, Australia. *PLoS One*, **14**, e0225128. https://doi.org/10.1371/journal.pone.0225128

Gauzere, P., Jiguet, F., & Devictor, V. (2015). Rapid adjustment of bird community compositions to local climatic variations and its functional consequences. *Global Change Biology*, **21**, 3367-3378. https://doi.org/10.1111/gcb.12917

Gotelli, N. J. (2000). Null model analysis of species co-occurrence patterns. *Ecology*, **81**, 2606-2621. https://doi.org/10.1890/0012-9658(2000)081[2606:NMAOSC]2.0.CO;2

Heemsbergen, D. A., Berg, M. P., Loreau, M., van Hal, J. R., Faber, J. H., & Verhoef, H. A. (2004). Biodiversity effects on soil processes explained by interspecific functional dissimilarity. *Science*, **306**, 1019-1020. https://doi.org/10.1126/science.1101865

Hengl, T. (2018) Soil pH in H₂O at 6 standard depths (0, 10, 30, 60, 100 and 200 cm) at 250 m resolution. https://zenodo.org/record/2525664

Hengl, T., & Wheeler, I. (2018) Soil organic carbon content in ×5 g/kg at 6 standard depths (0, 10, 30, 60, 100 and 200 cm) at 250 m resolution. https://zenodo.org/record/2525553

Hijmans, R. J., Phillips, S., Leathwick, J., & Elith, J. (2020). *dismo: Species distribution modeling*. Version R package version 1.3-3. https://CRAN.R-project.org/package=dismo

Hoffmann, A. A., & Sgro, C. M. (2011). Climate change and evolutionary adaptation. *Nature*, **470**, 479-485. https://doi.org/10.1038/natur e09670

Hooper, D. U., Adair, E. C., Cardinale, B. J., Byrnes, J. E. K., Hungate, B. A., Matulich, K. L., Gonzalez, A., Duffy, J. E., Gamfeldt, L., & O'Connor, M. I. (2012). A global synthesis reveals biodiversity loss as a major driver of ecosystem change. *Nature*, **486**, 105-108. https://doi.org/10.1038/nature11118

Huang, W., González, G., & Zou, X. (2020). Earthworm abundance and functional group diversity regulate plant litter decay and soil organic carbon level: A global meta-analysis. *Applied Soil Ecology*, **150**, 103473. https://doi.org/10.1016/j.apoilso.2019.103473

Hurtt, G. C., Chini, L., Sahajpal, R., Frolking, S., Bodirsky, B. L., Calvin, K., Doelman, J. C., Fisk, J., Fujimori, S., Klein Goldewijk, K., Hasegawa, T., Havlik, P., Heinimann, A., Humpenöder, F., Jungclaus, J., Kaplan, J. O., Kennedy, J., Krisztin, T., Lawrence, D., ... Zhang, X. (2020). Harmonization of global land use change and management for the period 850-2100 (LUH2) for CMIP6. *Geoscientific Model Development*, **13**, 5425-5464. https://doi.org/10.5194/gmd-13-5425-2020

Jiménez, J.-J., Decaëns, T., Amézquita, E., Rao, I., Thomas, R. J., & Lavelle, P. (2011). Short-range spatial variability of soil physico-chemical variables related to earthworm clustering in a neotropical gallery forest. *Soil Biology and Biochemistry*, **43**, 1071-1080. https://doi.org/10.1016/j.soilbio.2011.01.028

Jordan, D., Miles, R. J., Hubbard, V. C., & Lorenz, T. (2004). Effect of management practices and cropping systems on earthworm abundance and microbial activity in Sanborn Field: A 115-year-old agricultural field. *Pedobiologia*, **48**, 99-110. https://doi.org/10.1016/j.pedobi.2003.06.001

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., Zimmermann, N. E., Linder, H. P., & Kessler, M. (2017). Climatologies at high resolution for the earth's land surface areas. *Scientific Data*, **4**, 170122. https://doi.org/10.1038/sdata.2017.122

Karger, D. N., Nobis, M. P., Normand, S., Graham, C. H., & Zimmermann, N. E. (2021). CHELSA-TraCE21k v1.0. Downscaled transient temperature and precipitation data since the last glacial maximum. *Climate of the Past Discussions*, **1-27**. https://doi.org/10.5194/cp-2021-30

Karger, D. N., & Zimmermann, N. E. (2018). CHELSAcruts - High resolution temperature and precipitation timeseries for the 20th century and beyond. *EnviDat*. https://doi.org/10.16904/envidat.159

Kembel, S. W., Cowan, P. D., Helmus, M. R., Cornwell, W. K., Morlon, H., Ackerly, D. D., Blomberg, S. P., & Webb, C. O. (2010). Picante: R tools for integrating phylogenies and ecology. *Bioinformatics*, **26**, 1463-1464. https://doi.org/10.1093/bioinformatics/btq166

Kramer-Schadt, S., Niedballa, J., Pilgrim, J. D., Schröder, B., Lindenborn, J., Reinfelder, V., Stillfried, M., Heckmann, I., Scharf, A. K., Augeri, D. M., Cheyne, S. M., Hearn, A. J., Ross, J., Macdonald, D. W., Mathai, J., Eaton, J., Marshall, A. J., Semiadi, G., Rustam, R., ... Wilting, A. (2013). The importance of correcting for sampling bias in MaxEnt
species distribution models. *Diversity and Distributions*, 19, 1366-1379. https://doi.org/10.1111/ddi.12096

Laliberte, E., & Legendre, P. (2010). A distance-based framework for measuring functional diversity from multiple traits. *Ecology*, 91, 299-305. https://doi.org/10.1890/08-2244.1

Le Bagousse-Pinguet, Y., Soliveres, S., Gross, N., Torices, R., Berdugo, M., & Maestre, F. T. (2019). Phylogenetic, functional, and taxonomic richness have both positive and negative effects on ecosystem multifunctionality. *Proceedings of the National Academy of Sciences of the United States of America*, 116(17), 8419–8424. https://doi.org/10.1073/pnas.1815727116

Lembrechts, J. J., Aalto, J., Ashcroft, M. B., De Frenne, P., Kopecký, M., Lenoir, J., Luoto, M., Maclean, I. M. D., Roupsard, O., Fuentes-Lillo, E., García, R. A., Pellissier, L., Mitteloud, C., Alatalo, J. M., Smith, S. W., Björk, R. G., Muffler, L., Ratier Backes, A., Cesarz, S., ... Nijs, I. (2020). SoilTemp: A global database of near-surface temperature. *Global Change Biology*, 26, 6616–6629. https://doi.org/10.1111/gcb.15123

Lenoir, J., & Svenning, J. C. (2014). Climate-related range shifts - a global multidimensional synthesis and new research directions. *Ecography*, 38, 15-28. https://doi.org/10.1111/ecog.00967

Liu, C., White, M., & Newell, G. (2013). Selecting thresholds for the prediction of species occurrence with presence-only data. *Journal of Biogeography*, 40, 778–789. https://doi.org/10.1111/jbi.12058

Liu, T., Chen, X., Gong, X., Lubbers, I. M., Jiang, Y., Feng, W., Li, X., Whalen, J. K., Bonkowski, M., Griffiths, B. S., Hu, F., & Liu, M. (2019). Earthworms coordinate soil biota to improve multiple ecosystem functions. *Current Biology*, 29(20), 3420–3429.e5. https://doi.org/10.1016/j.cub.2019.08.045

Lobo, J. M., Jimenez-Valverde, A., & Real, R. (2008). AUC: A misleading measure of the performance of predictive distribution models. *Global Ecology and Biogeography*, 17, 145-151. https://doi.org/10.1111/j.1466-8238.2007.00358.x

Marchán, D. F., Refoyo, P., Novo, M., Fernández, R., Trigo, D., & Díaz Cosín, D. J. (2015). Predicting soil micro-variables and the distribution of an endogeic earthworm species through a model based on large-scale variables. *Soil Biology and Biochemistry*, 81, 124-127. https://doi.org/10.1016/j.soilbio.2014.10.023

Marichal, R., Praxedes, C., Decaëns, T., Grimaldi, M., Oszwald, J., Brown, G. G., Desjardins, T., da Silva, M. L., Feijoo Martinez, A., Oliveira, M. N. D., Velasquez, E., & Lavelle, P. (2017). Earthworm functional traits, landscape degradation and ecosystem services in the Brazilian Amazon deforestation arc. *European Journal of Soil Biology*, 83, 43-51. https://doi.org/10.1016/j.ejsobi.2017.09.003

Mason, N. W. H., de Bello, F., Mouillot, D., Pavoine, S., & Dray, S. (2013). A guide for using functional diversity indices to reveal changes in assembly processes along ecological gradients. *Journal of Vegetation Science*, 24, 794-806. https://doi.org/10.1111/jvs.12013

Mason, N. W. H., Mouillot, D., Lee, W. G., & Wilson, J. B. (2005). Functional richness, functional evenness and functional divergence: The primary components of functional diversity. *Oikos*, 111, 112-118. https://doi.org/10.1111/j.0030-1299.2005.13886.x

Mathieu, J., Barot, S., Blouin, M., Caro, G., Decaëns, T., Dubs, F., Dupont, L., Jouquet, P., & Nai, P. (2010). Habitat quality, conspecific density, and habitat pre-use affect the dispersal behaviour of two earthworm species, *Aporrectodea icterica* and *Dendrobaena veneta*, in a mesocosm experiment. *Soil Biology and Biochemistry*, 42, 203-209. https://doi.org/10.1016/j.soilbio.2009.10.018

Mathieu, J., & Davies, J. T. (2014). Glaciation as an historical filter of below-ground biodiversity. *Journal of Biogeography*, 41, 1204-1214. https://doi.org/10.1111/jbi.12284

Meehl, G. A., & Tebaldi, C. (2004). More intense, more frequent, and longer lasting heat waves in the 21st century. *Science*, 305, 994-997. https://doi.org/10.1126/science.1098704

Morelli, T. L., Barrows, C. W., Ramirez, A. R., Cartwright, J. M., Ackerly, D. D., Eaves, T. D., Ebersole, J. L., Krawchuk, M. A., Letcher, B.

H., Mahalovich, M. F., Meigs, G. W., Michalak, J. L., Millar, C. I., Quiñones, R. M., Stralberg, D., & Thorne, J. H. (2020). Climate-change refugia: Biodiversity in the slow lane. *Frontiers in Ecology and the Environment*, 18, 228-234. https://doi.org/10.1002/fee.2189

Muscarella, R., Galante, P. J., Soley-Guardia, M. et al (2014). ENMeval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. *Methods in Ecology and Evolution*, 5, 1198-1205. https://doi.org/10.1111/2041-210X.12261

O'Gorman, P. A. (2015). Precipitation extremes under climate change. *Current Climate Change Reports*, 1, 49-59. https://doi.org/10.1007/s40641-015-0009-3

Palm, J., van Schaik, N. L. M. B., & Schröder, B. (2013). Modelling distribution patterns of anecic, epigeic and endogeic earthworms at catchment-scale in agro-ecosystems. *Pedobiologia*, 56, 23-31. https://doi.org/10.1016/j.pedobi.2012.08.007

Parmesan, C. (2006). Ecological and evolutionary responses to recent climate change. *Annual Review of Ecology Evolution and Systematics*, 37, 637-669. https://doi.org/10.1146/annurev.ecols.ys.37.091305.110100

Pelosi, C., Pey, B., Hedde, M., Caro, G., Capowiez, Y., Guernion, M., Peigné, J., Piron, D., Bertrand, M., & Cluzeau, D. (2014). Reducing tillage in cultivated fields increases earthworm functional diversity. *Applied Soil Ecology*, 83, 79-87. https://doi.org/10.1016/j.apsoil.2013.10.005

Perkins-Kirkpatrick, S. E., & Lewis, S. C. (2020). Increasing trends in regional heatwaves. *Nature Communications*, 11, 3357. https://doi.org/10.1038/s41467-020-16970-7

Phillips, H. R. P., Guerra, C. A., Bartz, M. L. C., Briones, M. J. I., Brown, G., Crowther, T. W., Ferlian, O., Gongalsky, K. B., van den Hoogen, J., Krebs, J., Orgiazzi, A., Routh, D., Schwarz, B., Bach, E. M., Bennett, J. M., Brose, U., Decaëns, T., König-Ries, B., Loreau, M., ... Eisenhauer, N. (2019). Global distribution of earthworm diversity. *Science*, 366, 480-485. https://doi.org/10.1126/science.aax4851

Piearce, T. G. (1983). Functional morphology of lumbricid earthworms, with special reference to locomotion. *Journal of Natural History*, 17, 95-111. https://doi.org/10.1080/00222938300770071

Pimm, S. L., Jenkins, C. N., Abell, R., Brooks, T. M., Gittleman, J. L., Joppa, L. N., Raven, P. H., Roberts, C. M., & Sexton, J. O. (2014). The bio-diversity of species and their rates of extinction, distribution, and protection. *Science*, 344, 1246752. https://doi.org/10.1126/Science.1246752

Pradervand, J. N., Pellissier, L., Randin, C. F., & Guisan, A. (2014). Functional homogenization of bumblebee communities in alpine landscapes under projected climate change. *Climate Change Responses*, 1, 1. https://doi.org/10.1186/s40665-014-0001-5

Quillin, K. J. (2000). Ontogenetic scaling of burrowing forces in the earthworm *Lumbricus terrestris*. *Journal of Experimental Biology*, 203, 2757-2770. https://doi.org/10.1242/jeb.203.18.2757

Radosavljevic, A., & Anderson, R. P. (2014). Making better Maxent models of species distributions: Complexity, overfitting and evaluation. *Journal of Biogeography*, 41, 629-643. https://doi.org/10.1111/jbi.12227

Rutgers, M., Orgiazzi, A., Gardi, C., Römbke, J., Jänsch, S., Keith, A. M., Neilson, R., Boag, B., Schmidt, O., Murchie, A. K., Blackshaw, R. P., Pérès, G., Cluzeau, D., Guernion, M., Briones, M. J. I., Rodeiro, J., Piñeiro, R., Cosín, D. J. D., Sousa, J. P., ... Zwart, D. D. (2016). Mapping earthworm communities in Europe. *Applied Soil Ecology*, 97, 98-111. https://doi.org/10.1016/j.apsoil.2015.08.015

Santorufo, L., Cortet, J., Arena, C., Goudon, R., Rakoto, A., Morel, J.-L., & Maisto, G. (2014). An assessment of the influence of the urban environment on collembolan communities in soils using taxonomy-and trait-based approaches. *Applied Soil Ecology*, 78, 48-56. https://doi.org/10.1016/j.apsoil.2014.02.008
Schröder, B. (2008). Challenges of species distribution modeling below-ground. *Journal of Plant Nutrition and Soil Science*, **171**, 325-337. https://doi.org/10.1002/jpln.200700027

Schwalm, C. R., Glendon, S., & Duffy, P. B. (2020). RCP8.5 tracks cumulative CO₂ emissions. *Proceedings of the National Academy of Sciences of the United States of America*, **117**, 19656-19657. https://doi.org/10.1073/pnas.2007117117

Setälä, H., Laakso, J., Mikola, J., & Huhta, V. (1998). Functional diversity of decomposer organisms in relation to primary production. *Applied Soil Ecology*, **9**, 25-31. https://doi.org/10.1016/S0929-1393(98)00049-3

Sims, R. W., & Gerard, B. M. (1985). *Earthworms: Keys and notes for the identification and study of the species*. Brill Archive.

Singh, J., Schädler, M., Demetrio, W. et al (2019). Climate change effects on earthworms - a review. *Soil Organisms*, **91**, 113-137. https://doi.org/10.2567/so91iss3pp114

Smith, R., Mcswiney, C., Grandy, A., Suwanwaree, P., Snider, R., & Robertson, G. (2008). Diversity and abundance of earthworms across an agricultural land-use intensity gradient. *Soil and Tillage Research*, **100**, 83-88. https://doi.org/10.1016/j.still.2008.04.009

Soberon, J., & Nakamura, M. (2009). Niches and distributional areas: Concepts, methods, and assumptions. *Proceedings of the National Academy of Sciences of the United States of America*, **106**(Supplement_2), 19644-19650. https://doi.org/10.1073/pnas.0901637106

Spurgeon, D. J., Keith, A. M., Schmidt, O., Lammertsma, D. R., & Faber, J. H. (2013). Land-use and land-management change: Relationships with earthworm and fungi communities and soil structural properties. *BMC Ecology*, **13**, 46. https://doi.org/10.1186/1472-6785-13-46

Thomson, A. M., Calvin, K. V., Smith, S. J., Kyle, G. P., Volke, A., Patel, P., Delgado-Arias, S., Bond-Lamberty, B., Wise, M. A., Clarke, L. E., & Edmonds, J. A. (2011). RCP4.5: A pathway for stabilization of radiative forcing by 2100. *Climatic Change*, **109**(77), https://doi.org/10.1007/s10584-011-0151-4

Title, P. O., & Bemmels, J. B. (2018). envirem: An expanded set of bioclimatic and topographic variables increases flexibility and improves performance of ecological niche modeling. *Ecography*, **41**, 291-307. https://doi.org/10.1111/ecog.02880

Toro, I. D., Silva, R. R., & Ellison, A. M. (2015). Predicted impacts of climatic change on ant functional diversity and distributions in eastern North American forests. *Diversity and Distributions*, **21**, 781-791. https://doi.org/10.1111/ddi.12331

Uvarov, A. V. (2009). Inter- and intraspecific interactions in lumbricid earthworms: Their role for earthworm performance and ecosystem functioning. *Pedobiologia*, **53**, 1-27. https://doi.org/10.1016/j.pedobi.2009.05.001

Vignali, S., Barras, A. G., Arlettaz, R., & Braunisch, V. (2020). SDMtune: An R package to tune and evaluate species distribution models. *Ecology and Evolution*, **10**, 11488-11506. https://doi.org/10.1002/ece3.6786

Villéger, S., Mason, N. W. H., & Mouillot, D. (2008). New multidimensional functional diversity indices for a multifaceted framework in functional ecology. *Ecology*, **89**, 2290-2301. https://doi.org/10.1890/07-1206.1

Vincent, Q., Leyval, C., Beguiristain, T., & Auclerc, A. (2018). Functional structure and composition of Collembola and soil macrofauna communities depend on abiotic parameters in derelict soils. *Applied Soil Ecology*, **130**, 259-270. https://doi.org/10.1016/j.apsoil.2018.07.002

Wang, S., Olatunji, O. A., Guo, C., Zhang, L., Sun, X., Tariq, A., Wu, X., Pan, K., Li, Z., Sun, F., & Song, D. (2020). Response of the soil macrofauna abundance and community structure to drought stress under agroforestry system in southeastern Qinghai-Tibet Plateau. *Archives of Agronomy and Soil Science*, **66**, 792-804. https://doi.org/10.1080/03650340.2019.1639154

Wieczynski, D. J., Boyle, B., Buzzard, V., Duran, S. M., Henderson, A. N., Hulshof, C. M., Kerkhoff, A. J., McCarthy, M. C., Michaletz, S. T., Swenson, N. G., Asner, G. P., Bentley, L. P., Enquist, B. J., & Savage, V. M. (2019). Climate shapes and shifts functional biodiversity in forests worldwide. *Proceedings of the National Academy of Sciences of the United States of America*, **116**(2), 587-592. https://doi.org/10.1073/pnas.1813723116

Wikramanayake, E., Dinerstein, E., Loucks, C., Olson, D., Morrison, J., Lamoreux, J., McKnight, M., & Hedao, P. (2002). Ecoregions in ascendance: Reply to Jepson and Whittaker. *Conservation Biology*, **16**, 238-243. https://doi.org/10.1046/j.1523-1739.2002.01403.x

Zurell, D., Franklin, J., König, C., Bouchet, P. J., Dormann, C. F., Elith, J., Fandos, G., Feng, X., Guillera-Arroita, G., Guisan, A., Lahoz-Monfort, J. J., Leitão, P. J., Park, D. S., Peterson, A. T., Rapacciuolo, G., Schmatz, D. R., Schröder, B., Serra-Diaz, J. M., Thuiller, W., ... Merow, C. (2020). A standard protocol for reporting species distribution models. *Ecography*, **43**, 1261-1277. https://doi.org/10.1111/ecog.04960

## BIOSKETCHES

Yoan Fourcade is an associate professor at the Univ. Paris-Est Creteil and the Institute of Ecology and Environmental Sciences of Paris; his research focuses on macroecological biodiversity trends in a climate change context. **Mathias Vercauteren** is currently a master candidate at Univ. Paris-Est Creteil.

Author contributions: Yoan Fourcade conceived the idea, analysed the data and wrote the manuscript. Mathias Vercauteren prepared and analysed the data, wrote a first draft in French language. Both authors contributed to the final version of the manuscript.

## SUPPORTING INFORMATION

Additional supporting information may be found in the online version of the article at the publisher's website.

**How to cite this article:** Fourcade, Y., & Vercauteren, M. (2022). Predicted changes in the functional structure of earthworm assemblages in France driven by climate change. *Diversity and Distributions*, **00**, 1-17. https://doi.org/10.1111/ddi.13505
SUPPORTING INFORMATION

<table><caption>**Table S1:** Summary of species traits used to compute functional diversity indices. For each trait, we show the mean, minimum and maximum value across all species, or the proportion of species in each category.</caption><thead><tr><td>Trait</td><td>Summary of trait values</td><td>Type of trait</td></tr></thead><tbody><tr><td>Maximum length (mm)</td><td>174.69 [18 - 1050]</td><td>General morphology</td></tr><tr><td>Maximum width (mm)</td><td>5.867 [1 - 25]</td><td>General morphology</td></tr><tr><td>Maximum weight (mg)</td><td>6566 [20 - 106000]</td><td>General morphology</td></tr><tr><td>Body shape</td><td>cylindrical (90%), cylindrical to subtrapezoidal (1%), cylindrical with puberculus enlargement (1%), quadrangular (2%), subtrapezoidal (6%)</td><td>General morphology</td></tr><tr><td>Prostomium shape</td><td>absent (1%), epilobic (76%), prolobic (13%), tanylobic (9%)</td><td>General morphology</td></tr><tr><td>Flat tail</td><td>no (54%), yes (46%)</td><td>General morphology</td></tr><tr><td>Heart position (segment no.)</td><td>4 (1%), 5 (2%), 6 (62%), 7 (31%), 9 (1%), 10 (3%)</td><td>Anatomy</td></tr><tr><td>Heart size</td><td>2 (1%), 3 (3%), 4 (3%), 5 (90%), 6 (1%), 7 (2%)</td><td>Anatomy</td></tr><tr><td>Body color</td><td>none (34%), melanic (35%), green (1%), red (30%)</td><td>General morphology</td></tr><tr><td>Longitudinal color gradient</td><td>yes (52%), no (48%)</td><td>General morphology</td></tr><tr><td>Transversal color gradient</td><td>yes (51%), no (49%)</td><td>General morphology</td></tr><tr><td>Epidermis thickness (ranked 1-3)</td><td>2.17 [1 - 3]</td><td>Anatomy</td></tr><tr><td>Epidermis smoothness (ranked 1-3)</td><td>1.73 [1 - 3]</td><td>Anatomy</td></tr><tr><td>Mucus quantity</td><td>normal (85%), abondant (15%)</td><td>Anatomy</td></tr><tr><td>Muscle type</td><td>elementary (22%), intermediate (12%), squamous (7%), pinnate (47%), radial (11%)</td><td>Anatomy</td></tr><tr><td>Septa thickness (ranked 0-6)</td><td>3.44 [0 - 6]</td><td>Anatomy</td></tr><tr><td>Longitudinal furrows</td><td>yes (41%), no (59%)</td><td>General morphology</td></tr><tr><td>Transverse furrows</td><td>yes (91%), no (9%)</td><td>General morphology</td></tr><tr><td>Male pores position</td><td>6.63 [0 - 13.5]</td><td>Reproduction</td></tr><tr><td>Setae proximity</td><td>0.16 [0 - 1]</td><td>General morphology</td></tr><tr><td>Crop presence</td><td>yes (6%), no (94%)</td><td>Digestion</td></tr><tr><td>Morren gland presence</td><td>yes (91%), small or absent (2%), no (7%)</td><td>Digestion</td></tr><tr><td>Typhlosole location (segment no.)</td><td>21.23 [0 - 28]</td><td>Digestion</td></tr><tr><td>Typhlosole complexity (ranked 1-3)</td><td>1.89 [1 - 3]</td><td>Digestion</td></tr><tr><td>Clitellum shape</td><td>annular (5%), annular to saddle (1%), saddle-shaped (94%)</td><td>Reproduction</td></tr><tr><td>Clitellum size</td><td>9.51 [2 - 27]</td><td>Reproduction</td></tr><tr><td>Flat clitellum</td><td>yes (55%), little flat (1%), slightly flat (1%), no (43%)</td><td>Reproduction</td></tr><tr><td>Tubercula shape</td><td>gutter (42%), band (29%), other (29%)</td><td>Reproduction</td></tr><tr><td>Tubercula presence</td><td>yes (90%), no (10%)</td><td>Reproduction</td></tr><tr><td>Seminal vesicles complexity (ranked 1-3)</td><td>1.15 [1 - 3]</td><td>Reproduction</td></tr><tr><td>Seminal vesicles location</td><td>9.58 [0 - 14]</td><td>Reproduction</td></tr><tr><td>Spermathecal pore size (ranked 0-3)</td><td>2.20 [0 - 3]</td><td>Reproduction</td></tr><tr><td>Female pore location</td><td>14.54 [14.17 - 17.55]</td><td>Reproduction</td></tr><tr><td>Spermathecal pore location</td><td>15.60 [13.5 - 22.3]</td><td>Reproduction</td></tr><tr><td>Spermathecal pore area</td><td>A [2%], B [98%]</td><td>Reproduction</td></tr><tr><td>Nephridiopore distribution</td><td>aligned (65%), not aligned (35%)</td><td>Reproduction</td></tr></tbody></table>
**Methods S1:** ODMAP protocol for reporting details about the SDM procedure

## OVERVIEW

**Title of study:** Potential decrease in earthworm functional diversity caused by climate change

**Authorship:** Yoan Fourcade, Mathias Vercauteren

**Contact:** Univ. Paris Est Creteil, Sorbonne Université, CNRS, INRAE, IRD, Université de Paris, Institute of Ecology and Environmental Sciences Paris iEES, 94010 Créteil, France
yoan.fourcade@u-pec.fr

## MODEL OBJECTIVE

**Model objective:** Forecast and transfer

**Target output:** Functional diversity of earthworm assemblages

**Focal Taxon:** Earthworms (44 species)

## LOCATION

**Location:** Metropolitan France

## SCALE OF ANALYSIS

**Spatial extent:** -4.77, 9.57, 41.38, 51.09 (xmin, xmax, ymin, ymax)

**Spatial resolution:** 30 arc-sec

**Temporal extent:** 1960s, present (2000-2016), future (2040-2060 and 2060-2080)

**Boundary:** Political

## BIODIVERSITY DATA

**Observation type:** Standardized monitoring data

**Response data type:** Presence/absence

## PREDICTORS

**Predictor types:** Climatic, edaphic, land cover

## ASSUMPTIONS

**Model assumptions:** We assumed that at the scale of the study, earthworms' species distributions are mostly controlled by climate, although we still include soil and land cover predictors. We adopted two alternative hypotheses regarding dispersal: no dispersal or full dispersal, and combined individual models according to the species richness estimated by a macroecological model.

## ALGORITHMS

**Modelling techniques:** Boosted regression trees (BRT)

**Model complexity:** All possible combinations of the following parameters were tested: number of trees = 50, 100, 150, 200, 250, 300, 350, 400, 450, 500, interaction depth = 1, 2, 3, 4, shrinkage = 0.005, 0.01, 0.05, 0.1, 0.15, 0.2. The final models were the ones that provided the highest AUC.

**Model averaging:** No

## WORKFLOW

**Model workflow:** (1) Variables were cropped to the extent of France, (2) presence/absence data were obtained from a standardized protocol of soil sampling conducted in the
1960s, (3) BRTs were fitted to each earthworm species with hyperparameters chosen based on AUC, (4) models were projected in the present and in two future scenarios of climate change and two time periods, (5) functional diversity indices were calculated from stacked predictions at 2.5 arc-min resolution, with a probability ranking rule to limit the number of co-occurring species to the predictions of a macroecological model of species richness.

## SOFTWARE

**Software:** SDMtune package in R

**Code availability:** Can be provided upon request

**Data availability:** Climate variables: https://chelsa-climate.org/
Soil pH: https://zenodo.org/record/2525664
Soil organic carbon content: https://zenodo.org/record/2525553
Land cover: https://land.copernicus.eu/pan-european/corine-land-cover
Earthworm data: https://doi.org/10.5061/dryad.g7046

## DATA

### BIODIVERSITY DATA

**Taxon names:** 44 earthworm species

**Ecological level:** Species

**Data sources:** Data from Bouché MB (1972) Lombriciens de France : écologie et systématique. INRA, Paris. Available at: https://doi.org/10.5061/dryad.g7046

**Sampling design:** 1366 sites regularly spaced across France

**Sample size:** Prevalence from 0.008 to 0.702 (mean = 0.086)

**Clipping:** Metropolitan France

**Cleaning:** Only species with > 10 presence locations were used

**Absence data:** Input data included absence data

**Background data:** NA

### Data partitioning

**Training data:** Spatial blocks in a double 'checkerboard' grid of 1.5 arc-min and 2.5° resolution

**Validation data:** Spatial blocks in a double 'checkerboard' grid of 1.5 arc-min and 2.5° resolution

### Predictor variables

**Predictor variables:** 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, annual precipitation (BIO12) during the last glacial maximum, land cover.

The number of predictors was reduced by removing variables whose importance was < 10% if the AUC of the model was improved.

**Data sources:** CHELSA (climate): https://chelsa-climate.org/, OpenLandMap (soil): https://zenodo.org/record/2525664, https://zenodo.org/record/2525553, Copernicus (land cover): https://land.copernicus.eu/pan-european/corine-land-cover
Spatial extent: -4.77, 9.57, 41.38, 51.09 (xmin, xmax, ymin, ymax)

Spatial resolution: 30 arc-sec

Coordinate reference system: Initially EPSG 4326, model predictions were projected to EPSG 3035

Temporal extent:
Climate: 1960s
Past climate: last glacial maximum
Soil: 1950-2017
Land cover: 1990

Data processing: Bioclimatic variables have been created from rasters of monthly minimum temperature, maximum temperature and precipitation. Variables were cropped to the extent of Metropolitan France.

## Transfer data

Data sources: CHELSA (climate): https://chelsa-climate.org/

Spatial extent: -4.77, 9.57, 41.38, 51.09 (xmin, xmax, ymin, ymax)

Spatial resolution: 30 arc-sec

Temporal extent: Climate: 2000-2016 (present), 2041-2060, 2061-2080
Land-cover: 2012 (present), 2018 (future)

Models and scenarios:
Global circulation model: MIROC 5
Climate change scenarios: RCP 4.5 and RCP 8.5

Data processing: Rasters were cropped to the extent of Metropolitan France, including Corsica

Quantification of Novelty: none

## MODEL

### Multicollinearity

Multicollinearity: We removed all variables that had a variance inflation factor > 4

### Model settings

Distribution: Bernoulli

nTrees: 20-500 depending on species

interactionDepth: 2-4 depending on species

Shrinkage: 0.005-0.2 depending on species

bagFraction: 0.2

### Threshold selection

Threshold selection: Based on the threshold that maximizes the sum of sensitivity and specificity

## ASSESSMENT

### Performance statistics

Performance on validation data: Area under the ROC curve (AUC) calculated from spatial block partitioning following a double 'checkerboard' grid of 2.5 arc-min and 2.5° resolution

### Plausibility check
Response shapes: Not extracted

Expert judgement: Visualization of output binary maps

## PREDICTION

### Prediction output

Prediction unit: Continuous suitability maps

Binary presence/absence maps used for limiting dispersal

Post-processing: We stacked models in 2.5 arc-min grid cells following probability ranking rule to limit the number of co-occurring species to the predictions of a macroecological model of species richness, then calculated functional diversity (raw and standardized values) from species traits within each grid cell

### Uncertainty quantification

Scenario uncertainty: Projection into two scenarios and two time periods; two dispersal hypotheses

Novel environments: Not quantified
**Figure S1:** Distribution of AUC (a) and TSS (b) values obtained from the species distribution models trained on 44 species, and relationship between AUC and TSS (c). Panels (d-e) show the relationship between AUC (d) or TSS (e) and species prevalence, with the regression line and its standard-error obtained from a linear model between fitted between AUC or TSS and the logarithm of prevalence.

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>
**Figure S2:** Summary of variable importance obtained from permutation for the 44 modelled species. In (a) is shown the mean and standard-error of variable importance across all species; in (b) is shown the number of species for which each variable was ranked as the most important in model training.

<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /><figcaption>**Figure S3:** Scatterplots of the relationship between species richness and the three indices of functional diversity, for 5000 random grid cells. Plots are shown with regression lines, along with the corresponding R².</figcaption></figure>
**Figure S4:** (a) Mean community values of maximum species weight as predicted by stacked BRT species distribution models 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. Models are projected in space and time by allowing species to fill their predicted climatic niche (top row), or assuming that they cannot disperse (bottom row), while maximum species richness is constrained by a macroecological model (see Figure 2). (b) Predicted changes in mean community values of maximum species weight compared to the 1960s. Values are presented as mean percent change within each French ecoregion, coloured as in Figure 1.

a

<figure><img src="image_15.png" /></figure>

b

<figure><img src="image_16.png" /></figure>