# Digging deeper: deep joint species distribution modeling reveals environmental drivers of Earthworm Communities

Si-Moussi Sara, Thuiller Wilfried, Galbrun Esther, Decaëns Thibaud, Gérard Sylvain, Marchán Daniel F., Marsden Claire, Capowiez Yvan, Hedde Mickaël

► To cite this version:

Si-Moussi Sara, Thuiller Wilfried, Galbrun Esther, Decaëns Thibaud, Gérard Sylvain, et al. Digging deeper: deep joint species distribution modeling reveals environmental drivers of Earthworm Communities. Soil Biology and Biochemistry, 2026, 212, pp.110021. 10.1016/j.soilbio.2025.110021. hal-05336331

HAL Id: hal-05336331

https://hal.inrae.fr/hal-05336331v1

Submitted on 29 Oct 2025

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
Digging deeper: deep joint species distribution modeling reveals
environmental drivers of Earthworm Communities

Si-moussi Sara<sup>a, ID</sup>, Thuiller Wilfried <sup>a, ID</sup>, Galbrun Esther <sup>b, ID</sup>, Decaëns Thibaud <sup>c, ID</sup>, Gérard Sylvain <sup>d, ID</sup>, Marchán Daniel F. <sup>e, ID</sup>, Marsden Claire <sup>d, ID</sup>, Capowiez Yvan <sup>d, f</sup>, Hedde Mickaël <sup>d,*, ID</sup>

ARTICLE INFO

Keywords:

Earthworms

Soil biodiversity

Joint species distribution models

Deep learning

Explainable AI

Land cover

Climate variability

A B S T R A C T

Earthworms are key drivers of soil function, influencing organic matter turnover, nutrient cycling, and soil structure. Understanding the environmental controls on their distribution is essential for predicting the impacts of land use and climate change on soil ecosystems. While local studies have identified abiotic drivers of earth- worm communities, broad-scale spatial patterns remain underexplored. We developed a multi-species, multi-task deep learning model to jointly predict the distribution of 77 earthworm species across metropolitan France, using historical (1960–1970) and contemporary (1990–2020) records. The model integrates climate, soil, and land cover variables to estimate habitat suitability. We applied SHapley Additive exPlanations (SHAP) to identify key environmental drivers and used species clustering to reveal ecological response groups. The joint model achieved high predictive performance (TSS >0.7) and improved predictions for rare species compared to traditional species distribution models. Shared feature extraction across species allowed for more robust identification of common and contrasting environmental responses. Precipitation variability, temperature seasonality, and land cover emerged as dominant predictors of earthworm distribution but differed in ranking across species and functional groups. Species clustering into response groups to climatic, land use and soil revealed distinct ecological strategies including a gradient of sensitivity to precipitation seasonality, differential habitat prefer- ences in terms of vegetation cover and wetness and trade-offs between soil acidity and organic matter quality. Our study advances both the methodological and ecological understanding of soil biodiversity. We demonstrate the utility of interpretable deep learning approaches for large-scale soil fauna modeling and provide new insights into earthworm habitat specialization. These findings highlight land cover and seasonal climate variability as efficient proxies for soil biodiversity, providing actionable indicators for global monitoring initiatives and helping to identify habitat requirements of earthworm species to guide emerging earthworm conservation strategies in the face of global environmental change.

# 1. Introduction

Earthworms are among the most ecologically significant in-
vertebrates in terrestrial ecosystems, contributing substantially to
overall animal biomass and playing an essential role in soil functioning

(Blouin et al., 2013). These organisms are pivotal in driving key soil
processes such as organic matter decomposition, nutrient cycling, and
the formation of soil structure (Vidal et al., 2023). Through their bur-
rowing activities, earthworms increase soil porosity, enhance water
infiltration, and promote plant growth, making them crucial for

<footnote><sup>a</sup> Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LECA, Grenoble, France</footnote>

<footnote><sup>b</sup> School of Computing, University of Eastern Finland, Kuopio, Finland</footnote>

<footnote><sup>c</sup> CEFE, Univ. Montpellier, CNRS, EPHE, IRD, Montpellier, France</footnote>

<footnote><sup>d</sup> Eco&Sols, INRAE, IRD, CIRAD, Institut Agro, Montpellier, France</footnote>

<footnote><sup>e</sup> Department of Biodiversity, Ecology and Evolution, Faculty of Biology, Universidad Complutense de Madrid, Spain</footnote>

<footnote><sup>f</sup> Emmah, INRAE, Univ Avignon, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: sara.si-moussi@univ-grenoble-alpes.fr (S.-m. Sara), Wilfried.thuiller@univ-grenoble-alpes.fr (T. Wilfried), esther.galbrun@inria.fr (G. Esther), Thibaud.Decaens@cefe.cnrs.fr (D. Thibaud), sylvain.gerard@ens-paris-saclay.fr (G. Sylvain), danifermch@gmail.com (M. Daniel F.), claire.marsden@supagro.fr (M. Claire), yvan.capowiez@inrae.fr (C. Yvan), mickael.hedde@inrae.fr (H. Mickaël).</footnote>
ecosystem productivity and soil health (Arrázola-Vásquez et al., 2022; Jouquet et al., 2006). Additionally, as key components of soil food webs, earthworms connect various soil organisms, influencing nutrient flows and energy transfer across trophic levels below and above ground (Goury et al., 2025; Wan et al., 2025). Understanding their spatial distributions is therefore crucial for predicting ecosystem responses to environmental change such as climate shifts or land use transformations.

While local-scale studies have identified several abiotic drivers shaping earthworm distributions (Decaëns et al., 2008; Eggleton et al., 2009; Eisenhauer et al., 2014; Pelosi et al., 2014), large spatial patterns remain less understood (Mathieu and Jonathan Davies, 2014; Phillips et al., 2019; Zeiss et al., 2022). This gap in knowledge is particularly critical in the context of the proposed EU Directive on Soil Monitoring and Resilience, which emphasizes the urgent need for robust indicators of soil health and biodiversity at large scale (European Commission, 2021). Improving our understanding of soil organism responses to environmental changes is therefore essential. While traditional earthworm species classifications based on morphological traits (Bouché, 1977) have long guided ecological interpretations of species' occurrences, they can be subjective and inconsistently applied (Bottinelli and Capowiez, 2021). A more objective approach could involve the modeling of species responses to environmental gradients, enabling data-driven ecological classifications and improved predictive accuracy.

Single-species species distribution models (SDMs) are commonly used to predict species distributions based on environmental predictors (Guisan et al., 2017). However, many earthworm species are rather rare, often with sparse observations, making single-species SDMs prone to overfitting or underfitting. Also, in many single-species modeling workflows, predictions for each species are later “stacked” to infer community composition (Grenié et al., 2020). This process can lead to overestimation of species richness when summed across all single-species predictions (Calabrese et al., 2014; Deschamps et al., 2023). Finally, extracting common species responses to environmental gradients from single SDMs is not straightforward since different models can have selected different variables. Joint species distribution models (jSDM) offer interesting avenues to address these issues (Pollock et al., 2014; Warton et al., 2015). By jointly modeling multiple species, jSDMs can borrow strength from more common or data-rich species that share similar environmental responses. This approach provides more robust parameter estimates for each species, including rare ones (Norberg et al., 2019). If many species show parallel responses to certain environmental gradients (e.g., precipitation), these gradient-response relationships can be learned more accurately by pooling data, and reveal ecological clusters or guilds—groups of species responding similarly to temperature, precipitation, or land cover. This reduces uncertainty in model estimates and can improve predictions for each individual species. Additionally, since jSDMs explicitly model the joint occurrence of species, they are supposed to yield more accurate estimations of community composition. Traditional jSDMs can be computationally cumbersome, especially those using purely Bayesian or likelihood-based approaches and latent variable implementations (LV-jSDMs) (Ovaskainen and Abrego, 2020). In addition, most jSDM are constrained by linear assumptions, requiring manual feature engineering to capture nonlinear relationships.

To overcome these limitations, we have designed the MTEC (Multi-Task modeling of Ecological Communities), a deep latent variable jSDM integrating multi-task neural networks to capture complex species-environment relationships (Hu et al., 2025) and variational autoencoders for scalable inference of latent factors (Tang et al., 2018). Using large-scale earthworm data from mainland France including Corsica, we compare the predictive performance of MTEC to that of traditional SDMs. We then identify groups of species with shared habitat preferences, by using the explainable AI framework of Shapley additive explanations SHAP (Lundberg and Lee, 2017) to attribute earthworm distribution patterns to specific environmental drivers. Combining joint species modeling with explainable AI enhances predictive accuracy

while providing insights into the ecological processes shaping earthworm assemblages.

## 2. Materials and methods

### 2.1. Overview

We modeled earthworm distributions using historical calibration datasets covering diverse environmental conditions and evaluated the models on contemporary datasets (Fig. 1). The calibration dataset originated from Marcel Bouché's extensive surveys across mainland France including Corsica (Bouché, 1972), recently digitized and taxonomically updated by (Gérard et al., 2025). A separate dataset of recent georeferenced observations served as independent model evaluation data. Climate, soil, hydrology, and land cover variables were integrated as environmental predictors to capture major abiotic gradients. We implemented a deep joint species distribution model (MTEC), enabling simultaneous modeling of multiple earthworm species. Using the Shapley Additive Explanations (SHAP) framework, we identified and assessed key environmental drivers shaping species distributions. Lastly, unsupervised clustering was applied to categorize species into ecological response groups based on their environmental affinities.

### 2.2. Earthworm observation data

#### 2.2.1. Calibration dataset

During the 1960s, Marcel Bouché collected extensive data on earthworm communities in mainland France (including Corsica), documented in detail along with habitat characteristics (geographical coordinates, vegetation, land use, soil measurements) in his seminal work (Bouché, 1972). The earthworm occurrence data (presence/absence) used here were digitized by (Gérard et al., 2025), resulting in 34,737 individual earthworm records from 1399 sampling sites, covering 126 species and subspecies.

For modeling, we selected taxa occurring in at least five sites with detailed habitat descriptions. The refined dataset included 1346 sites and 77 taxa (species and subspecies sensu Bouché, 1972).

#### 2.2.2. Evaluation dataset

We aggregated a dataset of geotagged earthworm observations from multiple sources sampled in the [1998–2019] period. This dataset contains 1126 observations of 37 out of the 77 taxa occurring in the calibration dataset. Due to the multiplicity of sources and the heterogeneous sampling effort, we treated this collection as presence-only data.

### 2.3. Environmental covariates

We assembled environmental predictors related to climate, soil chemical and structural properties, hydrology, and land cover type. Two distinct datasets were prepared for spatial and temporal alignment with earthworm occurrences for model training and projection:

The first dataset (1960–1970, calibration) integrated in-situ soil and land cover data from (Bouché, 1972; Gérard et al., 2025), supplemented by GIS-derived environmental layers. Climate data were sourced from WorldClim (1960–1990; 1 km resolution), while soil and hydrological data were extracted from the European Soil Database at 1 km resolution (see Table 1 for more details about the specific variables we selected).

The second dataset (1990–2019, projection and evaluation) covering metropolitan France used climate variables from WorldClim, soil data from the European Soil Database, and land cover information from Corine Land Cover 2012–2018, maintaining the remaining predictors consistent.
<figure><img src="image_2.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(c)</figcaption></figure>

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

(caption on next page)
**Fig. 1.** Overview of the earthworm species distribution modeling study. (a) Training: a deep-learning joint species distribution model (MTEC) is trained on earthworm communities sampled across France. The model integrates in-situ field descriptions (land cover, soil physico-chemical properties) with environmental GIS databases (climate, hydrography, pedology). These inputs are embedded by a neural network into a shared environmental representation, along with latent factors capturing missing predictors, unmeasured features and interspecific associations to optimize multi-species predictions. (b) Projection and evaluation: The trained model is applied to continuous environmental layers to predict habitat suitability for earthworm species across France. The resulting habitat suitability maps are evaluated using independent observations to assess the accuracy of the predictions. (c) Explainability: Model predictions are analyzed using SHAP (SHapley Additive exPlanations). For a given set of observations (calibration dataset), SHAP values quantify the positive or negative impact of each environmental feature on the habitat suitability of each species within each site. Statistical summaries of these SHAP values help identify local and global drivers (SHAP importance), habitat preferences (SHAP summaries), and response groups (SHAP clustering).

<table><caption>**Table 1**

List of environmental predictors, definition, units and sources (calibration and projection).</caption><thead><tr><th>Group</th><th>Variable</th><th>Description</th><th>Unit, Type</th><th>Source (calibration)</th><th>Source (projection)</th></tr></thead><tbody><tr><td rowspan="10">Temperature</td><td>TMeanY</td><td>Annual mean temperature</td><td rowspan="10">Celsius</td><td rowspan="10">WorldClim (1960–1990) (Fick and Hijmans, 2017)</td><td rowspan="10">CHELSA (1979–2013) (Karger et al., 2017)</td></tr><tr><td>TMeanRngD</td><td>Mean diurnal range</td></tr><tr><td>Tiso</td><td>Isothermality</td></tr><tr><td>Tseason</td><td>Temperature seasonality</td></tr><tr><td>TMaxWarmM</td><td>Maximum of Warmest Month</td></tr><tr><td>TMinColdM</td><td>Minimum of Coldest Month</td></tr><tr><td>TRngY</td><td>Annual range</td></tr><tr><td>TMeanWetQ</td><td>Mean temperature of Wettest/ Driest/Warmest/Coldest Quarters</td></tr><tr><td>TMeanDryQ</td><td></td></tr><tr><td>TMeanWarmQ</td><td></td></tr><tr><td>TMeanColdQ</td><td></td></tr><tr><td>TMinColdM</td><td>Min/Max temperature of Coldest/ Warmest Months</td></tr><tr><td>TMaxWarmM</td><td></td></tr><tr><td rowspan="7">Precipitation</td><td>PTotY</td><td>Annual precipitation</td><td rowspan="7">mm</td><td rowspan="7"></td><td rowspan="7"></td></tr><tr><td>PSeason</td><td>Precipitation seasonality</td></tr><tr><td>PWetM</td><td>Precipitation of Wettest/Driest</td></tr><tr><td>PDryM</td><td>Month</td></tr><tr><td>PWetQ</td><td>Precipitation of Wettest/Driest/ Warmest/Coldest Quarters</td></tr><tr><td>PDryQ</td><td></td></tr><tr><td>PWarmQ</td><td></td></tr><tr><td>PColdQ</td><td></td></tr><tr><td>Land cover</td><td>LandCover</td><td>Land cover class</td><td>Categorical (16 levels)</td><td>(Bouché, 1972)</td><td>Corine land cover (Copernicus Land Monitoring Service, 2020)</td></tr><tr><td rowspan="3">Soil chemical properties</td><td>pH</td><td>Water acidity level (pH)</td><td>None</td><td rowspan="3">(Bouché, 1972)</td><td rowspan="3">ISRIC soilgrids (Poggio et al., 2021)</td></tr><tr><td>C</td><td>Organic carbon (C)</td><td>mg/kg soil</td></tr><tr><td>C/N</td><td>Carbon to Nitrogen ratio (C/N)</td><td>No unit</td></tr><tr><td rowspan="4">Hydrological properties</td><td>ProxiWater</td><td>Proximity to watercourses of at least 50 m width</td><td>Binary</td><td>(Botella et al., 2019)</td><td>(Botella et al., 2019)</td></tr><tr><td>AWC</td><td>Available water capacity of topsoil</td><td>mm</td><td rowspan="4">European soil database v2 raster library (ESDB) (Panagos et al., 2022)</td><td rowspan="4">European soil database v2 raster library (ESDB) (Panagos et al., 2022)</td></tr><tr><td>BaseSat</td><td>Base saturation of topsoil</td><td>%</td></tr><tr><td>CEC</td><td>Cation exchange of topsoil</td><td>cmol/kg soil</td></tr><tr><td>WaterRegime</td><td>Water regime</td><td>Categorical - 5 levels</td></tr><tr><td rowspan="5">Soil structural properties</td><td>Clay</td><td>Soil clay content</td><td>%</td><td rowspan="5">(Bouché, 1972)</td><td rowspan="5">ISRIC Soilgrids rasters (Poggio et al., 2021)</td></tr><tr><td>Silt</td><td>Soil silt content</td><td>%</td></tr><tr><td>Sand</td><td>Soil sand content</td><td>%</td></tr><tr><td>Crusting</td><td>Soil crusting class</td><td>Categorical (6 levels)</td></tr><tr><td>Erodi</td><td>Erodibility class</td><td>Categorical (6 levels)</td></tr><tr><td rowspan="3">Soil structural properties</td><td>DepthGleyH</td><td>Depth to gley horizon</td><td>cm</td><td rowspan="3">ESDB (Panagos et al., 2022)</td><td rowspan="3">ESDB (Panagos et al., 2022)</td></tr><tr><td>DepthImpL</td><td>Depth to impermeable layer</td><td>cm</td></tr><tr><td>ParentMat</td><td>Dominant parent material</td><td>Categorical (9 levels)</td></tr><tr><td>SBD</td><td></td><td>Soil bulk density</td><td>g/cm3</td><td></td><td></td></tr></tbody></table>

## 2.4. Earthworm distribution models

### 2.4.1. Joint model of earthworm species

To model earthworm species distributions, we designed MTEC to simultaneously predict probabilities of species presence based on environmental covariates and latent factors representing unmeasured ecological processes (e.g., missing environmental variables, interspecies interactions). MTEC leverages neural network architectures to efficiently capture nonlinear species-environment relationships and

improve multi-species prediction accuracy and computational efficiency simultaneously (Fig. 2). See Supplementary Information for details on model formulation (SI, section A).

### 2.4.2. Training

We implemented MTEC in Python using Keras (Chollet, 2015). We trained MTEC (Fig. 1a) using the Adam (Kingma and Ba, 2014) gradient-based optimization algorithm, minimizing a total loss composed of the species prediction loss (weighted binary cross-entropy,
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** **Architecture of MTEC.** MTEC consists of three main modules: a probabilistic encoder, an environmental feature encoder, and a multi-task decoder. The probabilistic encoder is a fully connected neural network that learns a variational posterior distribution *q* over latent factors *h*, capturing unobserved ecological influences. The environmental feature encoder processes raw environmental covariates *e* through a separate fully connected network, extracting a shared embedding *x* relevant to all species. This introduces non-linearity in the response and reduces the need for feature engineering. Finally, the multi-task decoder predicts species occurrence probabilities *θ* as a function of the shared environmental embedding (response coefficients *β*) and the latent factors (latent factor loadings *α*). During training (a), latent factors are inferred from the approximate posterior distribution *q* given known community data *y*. During prediction (b) on new sites, the latent factors are sampled from the prior distribution *p*.</figcaption></figure>

Goodfellow et al., 2016), the latent factor inference loss (Kullback-Leibler divergence, Kullback and Leibler, 1951) and a regularization term to control model complexity and prevent overfitting (Hastie et al., 2009; Hoerl and Kennard, 1970).

To improve the prediction accuracy of rare taxa, we used an

imbalance-aware partitioning strategy (SI, section B2) to ensure their
adequate representation in the training dataset. Model selection was
based on 5 × 2 cross-validation (Dietterich, 1998), with Area Under the
Receiver Operating Characteristic Curve (ROC-AUC) and the True Skill
Statistics (TSS) as the primary evaluation metrics (SI, section B1).
### 2.4.3. Stacked species distribution models

For comparison with MTEC, single species distribution models based on Generalized Linear Models (Nelder and Wedderburn, 1972), Gradient Boosting Machines (Elith et al., 2008), Random Forests (Breiman, 2001) were fitted independently using biomod2 in R (Thuiller et al., 2009, 2016) on the same calibration dataset as MTEC. Community composition predictions from these single-species models were combined through stacking (S-SDMs).

### 2.4.4. Projection and comparative evaluation

Using environmental data from the projection dataset, we predicted community composition across France at a 1 km resolution using MTEC (Fig. 1b) and S-SDMs, producing habitat suitability maps for all species. Model extrapolation was evaluated by comparing predicted suitability to occurrence points in the evaluation dataset using recall scores.

Finally, we compared the predictive performances of the three S-SDMs to MTEC in terms of ROC-AUC and True Skill Statistics (TSS) on the calibration dataset and in terms of recall on the evaluation dataset using the Wilcoxon rank-sum tests (Wilcoxon, 1945).

### 2.5. MTEC interpretability

To identify the key environmental factors shaping earthworm distribution, we applied the Shapley Additive exPlanations (SHAP) framework (Lundberg and Lee, 2017) using the R package iml (Molnar, 2018). SHAP is model-agnostic and provides a detailed breakdown of how each environmental variable contributes to the predicted species probabilities in every single pixel of interest (Fig. 1c). A positive SHAP value indicates that an environmental variable increases habitat suitability, while a negative value suggests a reduction compared to the average prediction. The effect of a variable can vary across its range: positive in some conditions, negative or neutral in others. This analysis was restricted to observed environmental covariates. We analyze ecological structure captured by the latent factors and resulting earthworm association networks in (SI, section C2).

#### 2.5.1. Attributing the environmental drivers of earthworm distribution

We computed the SHAP values for all environmental covariates for each species across samples of the calibration dataset (Fig. 1c). From these, we generated various statistical summaries to derive ecological insights into earthworm distribution through the lenses of the model. The absolute value of the SHAP statistic reflects the local importance of an environmental factor at a specific site. Mapping these values geographically revealed spatial variability in influential factors. Averaging SHAP values yielded global measures of variable importance for each species (global importance).

#### 2.5.2. Comparing earthworm species responses to the environment

Pairwise distances among species based on SHAP values were calculated, followed by hierarchical clustering (Ward, 1963). The optimal number of clusters was identified using the average of the GAP and WSS statistics (Tibshirani et al., 2001). Principal Component Analysis (PCA) was conducted on mean SHAP values per species, facilitating visualization and interpretation of species response groups. We particularly focused on precipitation and land cover, given their known significance in structuring earthworm communities and relevance under global environmental changes (Phillips et al., 2019).

## 3. Results

### 3.1. Predictive performances

Comparison between MTEC and single-species distribution models revealed significantly improved predictive performance of MTEC on the calibration dataset (Table 2). Predictive performance was generally lower for widespread taxa and variable for less prevalent taxa (SI,

<caption>Table 2</caption>

Median and standard deviation of the predictive performance metrics on the cross-validated calibration dataset (TSS, ROC-AUC) and on the independent evaluation dataset (Recall) across the 77 modeled taxa, comparing the proposed multi-task architecture (MTEC) with single species distribution models (GBM, GLM, and RF).

<table><thead><tr><td>Model</td><td>TSS</td><td>ROC AUC</td><td>Recall (Evaluation)</td></tr></thead><tbody><tr><td>MTEC</td><td>0.746 ± 0.243</td><td>0.912 ± 0.102</td><td>0.914 ± 0.214</td></tr><tr><td>GLM</td><td>0.411 ± 0.234</td><td>0.681 ± 0.186</td><td>0.519 ± 0.389</td></tr><tr><td>GBM</td><td>0.552 ± 0.281</td><td>0.753 ± 0.238</td><td>0.302 ± 0.267</td></tr><tr><td>RF</td><td>0.646 ± 0.247</td><td>0.816 ± 0.145</td><td>0.250 ± 0.314</td></tr></tbody></table>

section C). Specifically, the lowest ROC-AUC values (~0.69) of MTEC were obtained for widely distributed taxa such as *Lumbricus castaneus disjonctus* (Tetry, 1936) and *Octolasion cyaneum* (Savigny, 1826).

For highly prevalent species such as *Aporrectodea caliginosa* (Savigny, 1826) and *Lumbricus terrestris* (Linnaeus, 1758), performance differences between MTEC and more advanced SDMs (RF and GBM) were marginal or slightly unfavorable to MTEC. This indicates that common species benefit less from shared environmental information than rare or less common taxa. However, MTEC demonstrated a higher ability (+0.4) to predict contemporary occurrences from the evaluation dataset (Table 2), suggesting that single-species SDMs may be more prone to overfitting.

### 3.2. Habitat suitability

#### 3.2.1. Abiotic drivers of earthworm habitat suitability

Land cover, precipitation extremes, and temperature seasonality emerged as the most influential abiotic factors across taxa (Fig. 3). Interestingly, total annual precipitation had limited predictive value, whereas specific precipitation metrics such as rainfall during the wettest month (typically autumn in France), and precipitation during the coldest (winter) and warmest (summer) quarters were critical, along with precipitation seasonality. Among soil properties, soil pH was ranked among the most important variables, while organic matter quality (C/N ratio) influenced suitability more than the quantity of soil organic carbon. The soil water regime strongly affected habitat suitability for many taxa, whereas proximity to water bodies was largely insignificant, likely due to the coarse scale of hydrological rasters. Additionally, soil bulk density had a greater importance than soil texture.

The contribution of environmental variables can be assessed spatially. As an illustration, Fig. 4 shows the geographic distribution of the SHAP values, which can be interpreted as the local partial effect, of the most important climate variables for *Allolobophora icterica*. (See SI, section C4 for other modeled taxa).

#### 3.2.2. Response groups of earthworms to environmental gradients

Focusing on precipitation and land cover, two major ecological drivers highlighted in previous studies, we visualized species-specific responses through multi-species SHAP summary plots.

The PCA on SHAP values related to precipitation explained 60 % of the total variation (Fig. 5). Axis 1 (43 % of variance) differentiated species based on their affinity to climates with varying precipitation seasonality, with high seasonality associated with negative PCA scores. Using hierarchical clustering informed by the GAP statistic, we identified four response groups:

- The first cluster predominantly included Mediterranean species, including *Scherotheca porotheca* (Bouché, 1972), *Hormogaster samnitica* (Cognetti, 1914), *Kritodrilus calarensis* (Tetry, 1944), or species within the genus *Diporodrilus*.

- The second cluster comprised numerous endemics of the South of France and the Massif Central, from genera *Zophoscolex*, *Prosellodrilus*, *Scherotheca*, *Gatesona* and *Hemigastrodrilus*, along with
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Global environmental feature contributions. For each taxon, the importance of each environmental factor was quantified as the average (SHAP) contribution of that variable to habitat suitability probability across all sites. (a) The distribution of environmental features importance among taxa is presented in descending order, highlighting the relative influence of each variable on the model's predictions. Meaning of the codes for environmental features are given in Table 1. (b) Individual feature contributions were aggregated at the feature group level and normalized to obtain a variation partitioning for each earthworm species.</figcaption></figure>

eurytopic species such as *Allolobophora chlorotica* (Savigny, 1826) or *Aporrectodea rosea* (Savigny, 1826).

- The third cluster included taxa like *Eiseniella tetraedra* (Savigny, 1826), *Satchellius mammalis* (Savigny, 1826), *Dendrobaena octaedra* (Savigny, 1826), *Aporrectodea giardi* (Ribaucourt, 1901), and *Lumbricus centralis* Bouché (1972).

- The fourth cluster consisted primarily of widespread species (e.g., *Ap. caliginosa*, *O. cyaneum*) and a species known for riparian preferences (*Haplotaxis gordioides* (Hartman, 1821)).

Regarding land cover, the PCA on SHAP values accounted for approximately 80 % of the overall variation (Fig. 6). Axis 1 (62 % of the variance) distinguished species preferring closed canopy habitats (negative scores) from those favoring open habitats (positive values) whereas axis 2 (17.5 % of the variance) relates to species affinity to soil wetness. Clustering identified six distinct ecological response groups, corresponding to specialists of forests (cluster 6), shrublands (cluster 3), wetlands (cluster 5) and three grassland clusters: meadows (cluster 1), wet grasslands (cluster 2) and pastures (cluster 4). Taxa were unevenly distributed among these groups. Shrubland and grassland specialists each comprised 15–20 taxa, while forest and pasture specialists were less numerous (~10 taxa). Only four taxa were classified as wetland specialists. However, some clusters presented affinities with different habitats (e.g. cluster 2).

A similar analysis was conducted for earthworm response groups to soil properties (SI, section C5.3).

## 4. Discussion

In this study, we advance earthworm species distribution modeling by jointly predicting the distributions of 77 taxa across France, revealing more precise environmental drivers of their distribution. Using MTEC, an innovative joint Species Distribution Model (jSDM) powered by neural networks, we confirmed that multi-task neural networks notably outperform traditional single-species SDMs, especially in modeling rare taxa (SI, section C1). Additionally, our analysis identified primary environmental drivers shaping earthworm taxa distributions and used clustering to illustrate distinct ecological response groups.

### 4.1. *A deep joint species distribution model*

MTEC represents an evolution in jSDM techniques, combining neural networks with latent variables in a probabilistic framework akin to generalized linear latent variable models (GLLVM) (see Harris, 2015; Knape, 2016). Unlike conventional jSDMs employing generalized linear models, MTEC leverages neural networks that combine the functional representation capabilities of neural networks with the inferential capacity of hierarchical models to accommodate complex features and
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Geolocalized contribution of environmental factors to explaining the occurrence of *Aporrectodea icterica* (Savigny, 1826). This species is mostly influenced by precipitation and temperature factors, then structural, physico-chemical factors and land cover. For instance, in the mediterranean region, *A. icterica* presence is limited by the mean diurnal temperature range (TMeanRngD), the minimum temperature of the coldest month (TMinColdM), the precipitation of the warmest Quarter (PWarmQ) and the precipitation of the driest month (PDryM) and positively influenced by the precipitation during the coldest quarter (PColdQ). On the Channel facade, the factors with most positive importance were the mean temperature of the warmest quarter (TMeanWarmQ), the precipitation of the coldest Quarter (PColdQ) while those with negative influence were the isothermality (TIso), the mean daily temperature range (TMeanRngD) and the temperature seasonality (TSeason).</figcaption></figure>

model non-linear species-environment relationships. The use of variational autoencoders (VAEs) within MTEC facilitates efficient latent variable modeling with fixed complexity relative to dataset size. Although VAEs (Kingma and Welling, 2013) inherently assume independent latent factors, recent developments propose spatially and temporally structured latent priors (e.g., Gaussian Process Priors, Casale et al., 2018), potentially enhancing ecological modeling applications. Beyond undirected association networks, VAEs have also been applied in domains outside ecological research to learn the structures of Bayesian Networks (Zhang et al., 2019) and Gaussian Copulas (Mazoure, 2019).

## 4.2. Predictive performances

Contrary to previous reviews indicating limited advantages of jSDMs over stacked SDMs (Poggiato et al., 2021; Wilkinson et al., 2019), our findings demonstrate that MTEC frequently surpasses single-species models. MTEC's superior predictive performance, particularly for rare taxa, is likely due to shared environmental feature extraction, increased statistical power, and targeted loss weighting strategies addressing data imbalance. However, a few widespread taxa showed relatively low predictive performance across all fitted methods, suggesting that forecasting habitat suitability is more challenging for species with high adaptive capacities. This could be due to missing relevant environmental drivers (e.g. soil disturbance). For such widespread species which also occur beyond our study extent (France), it is likely that a larger scale modeling is required to cover their full niche distribution. Modeling relative abundances, rather than occurrences, could provide a deeper understanding of these species' ecological niches. In some instances, poor model performance might also indicate the presence of cryptic species creating taxonomic confusion in available data (e.g. among *L. terrestris* and *Lumbricus herculeus* (Savigny, 1826) (James et al., 2010), which further complicates accurate predictions (Marchán et al., 2020). Integrative taxonomy, which combines morphological, molecular, ecological, and geographical data, offers a promising approach for

resolving complex species boundaries, especially in soil biodiversity. Coupled with AI, it can accelerate species identification and pattern discovery across large, multidimensional datasets (Karbstein et al., 2024).

Beyond improving taxonomic discrimination, future research should integrate finer-scale environmental data that describe microhabitats to improve soil biodiversity modeling. Examples include proximal sensing data like soil electrical conductivity (EC) and near-infrared absorbance (NIR) (Schirrmann et al., 2016). Proxies for vegetation structure and temporal variability, such as NDVI-based metrics, along with soil depth and other topographic features like slope, topographic wetness index, and solar radiation, could be obtained via airborne or satellite remote sensing (Skidmore et al., 2021).

## 4.3. Abiotic drivers of habitat suitability for earthworm species

Consistently with prior research, land cover, precipitation, and temperature emerged as key determinants of earthworm habitat suitability (Phillips et al., 2019; Rutgers et al., 2016). Although our study operates at a mesoscopic scale (France) compared to Phillips et al. (2019), the geographic breadth of our work encompasses diverse clients and land uses, allowing for a clear hierarchy to emerge: at fine scales (landscape, catchment), soil properties dominate, whereas at broader scales, where climate and land use shift significantly, climate and land use become the primary drivers of species distributions.

Land cover reflects vegetation structure and management which shape litter inputs and microclimate, which in turn regulate resource availability and soil habitat conditions for earthworms (Curry, 2004; Edwards et al., 1996; Lavelle and Spain, 2001).

Precipitation regimes influence soil moisture, which in turn affects the activity and development of earthworms. Eisenhauer et al. (2014) demonstrated through experimental warming that earthworm development is constrained by temperature only when it leads to drought conditions. Intra-annual variability, seasonal metrics and extreme
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Response of earthworms to precipitation. (a) Response groups to precipitation: SHAP contributions to precipitation covariates across all observations (sites) for all species were organized into a matrix (species × site-covariate). Hierarchical clustering (Ward's algorithm) was applied to group species with similar responses across sites and precipitation covariates, with the GAP statistic identifying four optimal clusters. To aid visualization, Principal Component Analysis (PCA) was applied to reduce the dimensionality of the input matrix. The projection of species onto the first two principal components (explaining 56.9 % of the variation) is shown, along with their assigned clusters. (b) Multi-species SHAP summary plot of Precipitation Seasonality (PSeason): This plot visualizes how Precipitation Seasonality influences the habitat suitability of different earthworm taxa. Species label colors correspond to cluster colors (panel a). The y-axis shows the SHAP contribution, indicating whether PSeason has a positive or negative effect on predictions. A color gradient represents the range of PSeason values, helping to identify species preferred/limiting PSeason levels. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

events had stronger effects than mean temperature or total precipitation. Wet-season rainfall promotes surface foraging and activity after rain pulses, whereas hot-dry periods cause sharp declines (Curry, 2004; Holmstrup and Hornum, 2012). Temperature seasonality further constrains communities because earthworms rely on cutaneous respiration and osmoregulation, making them highly sensitive to thermal and moisture stress. Accordingly, activity, abundance, and biomass increase with temperature only under moist conditions, while drought or flooding strongly reduce performance (Curry, 2004; Singh et al., 2020). While some species developed strategies for coping with high climate variability both during the day and across months, others clearly prefer more uniform and optimal climatic conditions, such as observed in

regions like northwestern France.

At the edaphic level, soil pH was ranked among the most important variables, aligning with previous knowledge (Bouché, 1972). Soil pH influences not only plant growth and litter production but also the availability of metals (Al, Zn, Cu, Cd) that can reach toxic levels for earthworms (Spurgeon and Hopkin, 1996). The soil water regime was highly significant for most taxa, while proximity to water courses had a negligible impact, either due to the coarse scale of the raster layers, or its redundancy with soil water regime and wetland land cover classes. Additionally, soil bulk density had a greater effect at the assemblage level compared to texture, a result that is not surprising given that clay, silt and sand proportions were converted from discrete classes
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Response of earthworms to land cover. (a) Earthworm response Groups to Land Cover – SHAP contributions across sites were structured into a species × site matrix, with hierarchical clustering identifying six response groups (GAP statistic). To aid visualization, PCA was applied, and the projection of species onto the first two principal components (explaining 79.5 % of the variation) is shown along with their assigned clusters. (b) **Multi-Species SHAP Summary Plot** – This plot visualizes the effect of land cover categories on different earthworm species. The y-axis represents SHAP contributions, indicating whether a land cover category (scatter point colors) has a positive or negative effect on species predictions helping to identify species preferential habitats. Species label colors correspond to the cluster colors (panel a). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)</figcaption></figure>

(Hiederer, 2013). This transformation may have resulted in some loss of information and could have reduced the discriminative power of the texture features. Nevertheless, texture remained the most distinguishing factor for species such as *Dendrodrilus subrubicundus* (Eisen, 1874) and *E. tetraedra*, which are dependent on sandy soils near water bodies.

4.4. Differences across functional groups

Differences among functional groups (SI, section C3) highlight how environmental gradients filter species according to both their habitat requirements and their ecological roles. Litter dwellers proved most sensitive to wet-season precipitation and vegetation cover, a pattern consistent with their surface activity where litter moisture regulates food availability and refuge. By contrast, deep-burrowing and anecic
species were constrained mainly by temperature seasonality and dry-
quarter heat, reflecting the energetic costs and structural challenges of
maintaining galleries under fluctuating moisture regimes. Shallow bio-
turbators exhibited intermediate responses, combining sensitivity to
vegetation cover with temperature variation, in line with their role in
topsoil mixing (Bottinelli and Capowiez, 2021). Importantly, these
contrasting responses illustrate how habitat niches map onto functional
niches (Devictor et al., 2010; Violle et al., 2007): shifts in precipitation
or soil structure will not only redistribute species, but also alter
decomposition, mixing, and infiltration processes that underpin soil
functioning (Blouin et al., 2013; Lavelle and Spain, 2001).

4.5. *Building groups of response to environmental gradients*

In the context of European efforts to monitor and safeguard soil life (e.g., EU soil strategy for 2030), there is growing recognition of the need for robust, interpretable metrics that can inform management and policy decisions (Eisenhauer et al., 2022). One promising approach is to classify taxa into ecological response groups, such as habitat specialists vs. generalists or disturbance-sensitive vs. tolerant species. This strategy, already highlighted by Gérard et al., (2025a) in the development of new biodiversity indicators, can help detect shifts in community composition that are not visible through richness-based metrics alone. Similarly, Jupke et al. (2024) have emphasized the relevance of identifying typical earthworm assemblages across land-use gradients, offering a more integrative view of how soil fauna respond to environmental change. By building such response groups, we aimed to move beyond taxonomic inventories and towards ecological interpretations that are both actionable and scalable. We illustrated this approach with response groups to precipitation and land use.

Regarding precipitation regimes, amongst the four identified species
groups, the first cluster primarily consisted of Mediterranean species
that exhibit a strong ability to thrive in warm ecosystems with highly
seasonal rainfall, such as *H. samnitica* and *K. calarensis*. Species within
this cluster demonstrated a positive response to precipitation season-
ality, likely avoiding adverse climate conditions by their ability to enter
diapause in the deeper layers of the soil (Bouché, 1972). The second
cluster included many species from the genera *Zophoscolex*, *Prosello-
drilus*, *Scherotheca*, *Gatesona* and *Hemigastrodrilus*, which are endemic to
Southwestern Europe. It also comprised some eurytopic, widespread
species like *All. chlorotica* and *Ap. rosea*. We interpreted this cluster as a
group of mesophilous, warmth-tolerant species. The third cluster
included species such as *E. tetraedra*, *S. mammalis*, *D. octaedra*, *Ap. giardi*,
and *L. centralis*. These species require moist conditions for most of the
year to thrive (Bouché, 1972). The final cluster consisted of a few spe-
cies, some of which are known to be riparian. These species share a high
affinity for a precipitation regime that is uniformly distributed during
the year. Climate models project an increase in the frequency of extreme
precipitation events, while the total amount of rainfall is expected to
remain comparable to its current state, with greater variability. Given
earthworm sensitivity to changes in precipitation seasonality, these al-
terations are likely to impact their distributions. Rising temperatures,
drought, and increased winter rainfall will influence soil moisture and
temperature regimes, each affecting earthworms differently. The com-
bination of elevated temperatures, which heighten metabolic demands,
and increased extreme precipitation events may alter earthworm life
cycles. Additionally, extreme rainfall and changes in land use can make
soils more vulnerable to erosion, destroying earthworm habitats (Singh
et al., 2020).

The classification procedure on land cover response revealed sig-
nificant differences in earthworm habitat preferences. We identified six
response groups, five of which corresponded to open ecosystems and
three of which were associated with very distinct land cover classes:
forests, pastures, and inland wetlands. Taxa with affinities for shrub and
herbaceous vegetation were divided into two groups along the forest-to-
pasture gradient: scrub woodlands and meadows. An additional

intermediate group was identified for species from wet grasslands. Most
taxa exhibited a neutral or negative response to forest soils, except for a
limited number of litter-dwelling species, such as *D. octaedra*. A few ri-
parian species, such as *E. tetraedra* and *H. gordioides*, and several endo-
geics, e.g. *O. cyaneum* or *Helodrilus oculatus* (Hoffmeister, 1845),
demonstrated a strong preference for inland wetlands. Other litter-
dwelling species were positively associated with pastures, which may
support a higher species richness (Decaëns et al., 2008). Shrubs and
meadows were co-dominated by endogeic and anecic worms. Forests,
pastures, and shrub/herbaceous land cover classes accounted for
approximately 80 % of the calibration dataset. In contrast, arable and
permanent crop lands were underrepresented in the training dataset.
Their composition included widespread species. The lack of epigeics is
currently reported to be linked to farming practices such as tillage that
impact surface soil fauna (Briones and Schmidt, 2017; Decaëns et al.,
2008). There were no species with a strong affinity for agricultural sites,
although a few demonstrated a strong aversion toward these systems
(negative SHAP values). These conclusions may vary when considering
relative abundances rather than presence/absence data.

4.6. Implications for biodiversity monitoring and conservation planning

Our results have several implications for monitoring and conservation of earthworm, and more generally of soil fauna. The identification of land cover and seasonal climate variability as major predictors shows that simple, easily available products of remote sensing can capture both habitat preferences and functional roles of earthworm species, providing an efficient basis for initiatives such as SoilBON (Ganault et al., 2024; Guerra et al., 2021). At the same time, the classification into response groups makes it possible to identify which parts of the community are most at risk: for instance, litter dwellers in forests sensitive to canopy and moisture regimes, or deep-burrowers constrained by temperature seasonality and soil structure. This perspective goes beyond taxonomic inventories by linking habitat niches to ecosystem effects, thereby highlighting that shifts in environmental drivers will reorganize not only species distributions but also the processes they sustain. Such an approach echoes recent calls to move from generic indicators such as richness or biomass towards ecological attributes that better reflect vulnerability, specialization, or functional importance (Gérard et al., 2025). In this sense, conservation planning should explicitly integrate projected changes in precipitation and land cover in order to secure the functions provided by soil fauna as much as their diversity.

4.7. *Model transferability*

Our findings on how earthworm species respond to climate, land use,
and soil should be seen as a template for comparative analysis rather
than a universal rule. We recommend that extending these patterns to
other species, ecosystems, or biomes be treated as a testable hypothesis,
as differences in community composition and environmental ranges can
alter these relationships.

Transferring these findings to other locations (spatial transfer) or
time periods (temporal transfer) constitutes extrapolation. Such trans-
fers may be affected by factors not included in our models, such as
dispersal, biotic interactions (e.g., predation), or anthropogenic
pressures.

The MTEC framework is based on a multi-task neural network ar-
chitecture which allows transfer learning. This means the trained model
can be fine-tuned on new, related tasks, such as predicting new species
or community-level properties, by reusing the existing learned repre-
sentations. This process is most effective when the new data falls within
similar environmental ranges, biomes, and ecosystem types as the
original training data. Applying the model outside of this range would
require additional training data for reliable results.
## 4.8. Conclusions and future directions

Understanding the drivers of earthworm distribution is essential for predicting their response to ongoing global changes and estimating the spatial emergence of their functional benefits as ecosystem engineers and trophic resources for wildlife. So far, most studies have utilized single species distribution models that focus on a limited number of widespread taxa or, conversely, aggregated indicators such as biomass or richness, without considering species identities. Studying each species separately using dedicated species distribution models (SDMs) can be time-consuming and may be impractical for species with few available observations due to their intrinsic rarity or an insufficient sampling effort. The multi-task model presents new opportunities to study earthworm responses to global change. Beyond predicting distributions, we demonstrated the potential of classifying species into response groups along major environmental gradients, offering a framework to interpret biodiversity patterns more effectively. Such classifications could be extended to encompass other critical constraints influencing earthworm populations, providing further ecological insight.

Future improvements in earthworm distribution modeling should incorporate biotic interactions and microhabitat characteristics, alongside advanced remote sensing methods, to better capture complex ecological patterns (Hu et al., 2025; Poggiato et al., 2025). Integrating species dispersal capabilities through spatially explicit modeling will further enhance the accuracy of predicted distributions under global environmental changes. Ultimately, multi-species distribution models like MTEC represent valuable tools for understanding earthworm biodiversity dynamics and guiding effective conservation and management policies in soil ecosystems.

### CRediT authorship contribution statement

**Si-moussi Sara:** Writing – review & editing, Writing – original draft, Methodology, Data curation, Conceptualization. **Thuiller Wilfried:** Writing – review & editing, Supervision. **Galbrun Esther:** Supervision, Methodology. **Decaëns Thibaud:** Writing – review & editing, Validation. **Gérard Sylvain:** Writing – review & editing, Validation, Data curation. **Marchán Daniel F.:** Writing – review & editing, Validation. **Marsden Claire:** Writing – review & editing. **Capowiez Yvan:** Writing – review & editing. **Hedde Mickaël:** Writing – original draft, Validation, Supervision, Funding acquisition, Data curation, Conceptualization.

### Data and code availability statement

All data supporting the findings of this study are available within the article and its supplementary materials. Code used to run the analysis is available through a GitHub repository: https://github.com/bettasimouss/DeepLombrics.git.

### Ethics approval statement

This study did not involve human participants or experimental animals and therefore did not require ethics approval.

### Funding statement

This research is a by-product of the IMPACTS and LandWorm groups funded by the Centre for the Synthesis and Analysis of Biodiversity (CESAB) of the Foundation for Research on Biodiversity (FRB) and the Ministry of Ecological Transition. WT and SS-M. also acknowledge support from the HorizonEurope OBSGESSION (N°101,134,954) and NaturaConnect (N°101060429) projects and the MIAI@Grenoble Alpes (ANR-19-P3IA-0003). SG PhD thesis was supported by the ENS and by the INRAE Agroecosystem division through the GloWorm project.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.soilbio.2025.110021.

### References

Arrázola-Vásquez, E., Larsbo, M., Capowiez, Y., Taylor, A., Sandin, M., Iseskod, G., Keller, T., 2022. Earthworm burrowing modes and rates depend on earthworm species and soil mechanical resistance. Applied Soil Ecology 178, 104568. https://doi.org/10.1016/j.apsoil.2022.104568.

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. European Journal of Soil Science 64, 161-182. https://doi.org/10.1111/ejss.12025.

Botella, C., Servajean, M., Bonnet, P., Joly, A., 2019, September. Overview of geolifeclef 2019: plant species prediction using environment and animal occurrences. In CLEF 2019-Conference and Labs of the Evaluation. Forum 2380 (257).

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., 1977. Strategies lombriciennes. Ecological Bulletins 122-132.

Bouché, M.B., 1972. Lombriciens de France. Ecologie Et Systématique. Inra Editions.

Breiman, L., 2001. Random forests. Machine Learning 45, 5-32.

Briones, M.J.I., Schmidt, O., 2017. Conventional tillage decreases the abundance and biomass of earthworms and alters their community structure in a global meta-analysis. Global Change Biology 23, 4396-4419. https://doi.org/10.1111/gcb.13744.

Calabrese, J.M., Certain, G., Kraan, C., Dormann, C.F., 2014. Stacking species distribution models and adjusting bias by linking them to macroecological models. Global Ecology and Biogeography 23, 99-112. https://doi.org/10.1111/geb.12102.

Casale, F.P., Dalca, A., Saglietti, L., Listgarten, J., Fusi, N., 2018. Gaussian process prior variational autoencoders. In: Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., Garnett, R. (Eds.), Advances in Neural Information Processing Systems, vol. 31. Curran Associates, Inc. Available at: https://proceedings.neurips.cc/paper_files/paper/2018/file/1c336b8080f82bcc2cd2499b4c57261d-Paper.pdf.

Chollet, F., 2015. Keras.

Copernicus Land Monitoring Service, 2020. CORINE Land Cover 2018 (raster 100 m), Europe, 6-yearly. In: version 2020.20u1, May 2020. European Environment Agency. https://doi.org/10.2909/960998c1-1870-4e82-8051-6485205ebbc.

Curry, J.P., 2004. Factors affecting the abundance of earthworms in soils. In: Edwards, C. A. (Ed.), Earthworm Ecology. CRC Press, Boca Raton, FL, pp. 91-113.

Decaëns, T., Margerie, P., Aubert, M., Hedde, M., Bureau, F., 2008. Assembly rules within earthworm communities in North-Western France—A regional analysis. Applied Soil Ecology 39, 321-335. https://doi.org/10.1016/j.apsoil.2008.01.007.

Deschamps, G., Poggiato, G., Brun, P., Galiez, C., Thuiller, W., 2023. Predict first-assemble later versus assemble first-predict later: revisiting the dilemma for functional biogeography. Methods in Ecology and Evolution 14, 2680-2696. https://doi.org/10.1111/2041-210X.14203.

Devictor, V., Clavel, J., Julliard, R., Lavergne, S., Mouillot, D., Thuiller, W., Venail, P., Villéger, S., Mouquet, N., 2010. Defining and measuring ecological specialization. Journal of Applied Ecology 47, 15-25. https://doi-org.inee.bib.cnr.sfr/10.1111/j.1365-2664.2009.01744.x.

Dietterich, T.G., 1998. Approximate statistical tests for comparing supervised classification learning algorithms. Neural Computation 10, 1895-1923. https://doi.org/10.1162/08997669830001797.

Edwards, C.A., Bohlen, P.J., Hendrix, P., Arancon, N., 1996. Biology and ecology of earthworms, 3. Chapman & Hall, London.

Eggleton, P., Inward, K., Smith, J., Jones, D.T., Sherlock, E., 2009. A six year study of earthworm (Lumbricidae) populations in pasture woodland in southern England shows their responses to soil temperature and soil moisture. Soil Biology and Biochemistry 41, 1857-1865. https://doi.org/10.1016/j.soilbio.2009.06.007.

Eisenhauer, N., Bender, S.F., CalderónSanou, I., Vries, F.T., Lembrechts, J.J., Thuiller, W., et al., 2022. Frontiers in soil ecology—Insights from the World Biodiversity Forum 2022. Journal of Sustainable Agriculture and Environment 1, 245-261. https://doi.org/10.1002/sae2.12031.

Eisenhauer, N., Stefanski, A., Fischelli, N.A., Rice, K., Rich, R., Reich, P.B., 2014. Warming shifts 'warming': effects of experimental warming on invasive earthworms in northern North America. Scientific Reports 4, 6890. https://doi.org/10.1038/srep06890.

Elith, J., Leathwick, J.R., Hastie, T., 2008. A working guide to boosted regression trees. Journal of Animal Ecology 77, 802-813.

European Commission, 2021. EU Soil Strategy for 2030: Reaping the benefits of healthy soils for people, food, nature and climate. Communication from the Commission to the European Parliament, the Council, the European Economic and Social Committee
and the Committee of the Regions. Brussels, November 2021. Available at: https://
eur-lex.europa.eu/legal-content/EN/TXT/?uri=CELEX%3A520212DC069.

Fick, S.E., Hijmans, R.J., 2017. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. International Journal of Climatology 37, 4302–4315. https://doi.org/10.1010/joc.s0586.

Ganault, P., Ristok, C., Phillips, H.R.P., Hedde, M., Capowiez, Y., Bottinelli, N., Decaëns, T., Marchan, D., Gérard, S.O.E., Mathieu, J., Potapov, A., Cameron, E.K., Brown, G., Bartz, M., Zeiss, R., Zi, Y., Tsiafouli, M., Russell, D.J., Guerra, C., Eisenhauer, N., 2024. Soil BON Earthworm - A global initiative on earthworm distribution, traits and spatiotemporal diversity patterns. Soil Organisms 96 (1), 47–60. https://doi.org/10.25674/362.

Gérard, S., Decaëns, T., Butt, K.R., Briones, M.J.I., Capowiez, Y., Cluzeau, D., Hoeffner, K., Le Bayon, R.-C., Marchán, D.F., Marsden, C., Muys, B., Pelosi, C., Pérès, G., Phillips, H.R.P., Santini, L., Thuiller, W., Hedde, M., 2025. Putting earthworm conservation on the map: Shortfalls and solutions for developing earthworm conservation. Biological Conservation 302, 110911. https://doi.org/10.1016/j.biocon.2024.110911.

Gérard, S., Decaëns, T., Marchán, D.F., Beauchesne, M., Berlioz, L., Capowiez, Y., Clause, J., Decaëns, L., Dellavedova, R., Duhaut, C.-B., others, 2025. # Vers2022: 50-year resurvey data of French earthworm assemblages obtained after resampling Bouché's historical sites. Ecology 106, e70109.

Goodfellow, I., Bengio, Y., Courville, A., 2016. Deep Learning. MIT Press.

Goury, R., Münkemüller, T., Bonfanti, N., Calderón-Sanou, I., Poulenard, J., Consortium, ORCHAMP, Thuiller, W., 2025. The role of earthworms in mediating environmental impacts on mountain ecosystem functioning. Global Ecology and Biogeography 34 (8), e70093. https://doi.org/10.1111/geb.70093.

Grenié, M., Violle, C., Munoz, F., 2020. Is prediction of species richness from stacked species distribution models biased by habitat saturation? Ecological Indicators 111, 105970. https://doi.org/10.1016/j.ecolind.2019.105970.

Guerra, C.A., Bardgett, R.D., Caon, L., Crowther, T.W., Delgado-Baquerizo, M., Montanarella, L., Navarro, L.M., Orgiazzi, A., Singh, B.K., Tedersoo, L., Vargas-Rojas, J., Briones, M.J.I., Buscot, P., Cameron, E.K., Zarass, S., Chatzinotas, A., Cowan, D.A., Djukic, I., van den Hoogen, J., Lehmann, A., Maestre, F.T., Marin, C., Reitz, T., Rillig, M.C., Smith, L.C., de Vries, F.T., Weigelt, A., Eisenhauer, N., Wall, D. H., 2021. Tracking, targeting, and conserving soil biodiversity. Science 371 (6526), 239–241. https://doi.org/10.1126/science.abd7926.

Guisan, A., Thuiller, W., Zimmermann, N.E., 2017. Habitat Suitability and Distribution Models: with Applications in R, first ed. Cambridge University Press. https://doi.org/10.1017/9781139028271.

Harris, D.J., 2015. Generating realistic assemblages with a joint species distribution model. Methods in Ecology and Evolution 6, 465–473. https://doi-org.inee.bib.cnrs.fr/10.1111/2041-210X.12332.

Hastie, T., Tibshirani, R., Friedman, J., 2009. The Elements of Statistical Learning: Data Mining, Inference, and Prediction, second ed. Springer.

Hiederer, R., et al., 2013. EN, Scientific and Technical Research Series, 26082. Publications Office of the European Union, Luxembourg.

Hoerl, A.E., Kennard, R.W., 1970. Ridge regression: biased estimation for nonorthogonal problems. Technometrics 12, 55–67. https://doi.org/10.1080/00401706.1970.10488634.

Holmstrup, M., Hormum, H.D., 2012. Earthworm colonisation of abandoned arable soil polluted by copper. Pedobiologia 55, 63–65. https://doi.org/10.1016/j.pedobi.2011.08.005.

Hu, Y., Si-Moussi, S., Thuiller, W., 2025. Introduction to deep learning methods for multi-species predictions. Methods in Ecology and Evolution 16, 228–246. https://doi.org/10.1111/2041-210X.14466.

James, S.W., Porco, D., Decae ns, T., Richard, B., Rougerie, R., et al., 2010. DNA
Barcoding Reveals Cryptic Diversity in Lumbricus terrestris L., 1758 (Clitellata):
Resurrection of L. herculeus (Savigny, 1826). PLoS ONE 5 (12), e15629. https://doi.org/
org/10.1371/journal.pone.0015629.

Jouquet, P., Dauber, J., Lagerlöf, J., Lavelle, P., Lepage, M., 2006. Soil invertebrates as ecosystem engineers: intended and accidental effects on soil and feedback loops. Applied Soil Ecology 32, 153–164. https://doi.org/10.1016/j.apsoec.2005.07.004.

Jupke, J.F., Scheu, S., Cameron, E.K., Eisenhauer, N., Phillips, H.R.P., Römble, J., Rutgers, M., Schäfer, R.B., Entling, M.H., 2024. Typical earthworm assemblages of European ecosystem types. European Journal of Soil Science 75 (5), e13584. https://doi.org/10.1111/ejs.13584.

Karbstein, K., Kösters, L., Hodač, L., Hofmann, M., Hörandl, E., Tomasello, S., Wagner, N. D., Emerson, B.C., Albach, D.C., Scheu, S., Bradler, S., De Vries, J., Irisarri, T., Li, H., Soltis, P., Mäder, P., Wäldchen, J., 2024. Species delimitation 4.0: integrative taxonomy meets artificial intelligence. Trends in Ecology & Evolution 39, 771–784. https://doi.org/10.1016/j.tree.2023.11.002.

Karger, D., Conrad, O., Böhner, J., et al., 2017. Climatologies at high resolution for the earth's land surface areas. Scientific Data 4, 170122. https://doi.org/10.1038/sdata.2017.122.

Kingma, D.P., Ba, J., 2014. Adam: a method for stochastic optimization. https://doi.org/10.48550/ARXIV.1412.6980.

Kingma, D.P., Welling, M., 2013. Auto-encoding variational bayes. arXiv preprint arXiv: 1312.6114.

Knape, J., 2016. Decomposing trends in Swedish bird populations using generalized additive mixed models. Journal of Applied Ecology 53 (6), 1852–1861. https://doi.org/10.1111/1365-2664.12720.

Kullback, S., Leibler, R.A., 1951. On information and sufficiency. The Annals of Mathematical Statistics 22, 79–86. https://doi.org/10.1214/aoms/117772964.

Lavelle, P., Spain, A.V., 2001. Soil ecology. Springer Netherlands, Dordrecht.

Lundberg, S.M., Lee, S.-I., 2017. A unified approach to interpreting model predictions. Advances in Neural Information Processing Systems 30.

Marchán, D.F., Fernández, R., Domínguez, J., Díaz Cosín, Novo, M., 2020. Genome-informed integrative taxonomic description of three cryptic species in the earthworm genus *Carpetania* (Oligochea, Hym Rogadistria). Systematics and Biodiversity 18 (3), 203–215. https://doi.org/10.1080/14772000.2020.1730474.

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

Mazoure, B., 2019. Estimating the structure of probabilistic graphical models through a Gaussian copula with discrete marginals. McGill University (Canada).

Molnar, C., 2018. Im: an R package for Interpretable Machine Learning. Journal of Open Source Software 3, 786. https://doi.org/10.21105/joss.000786.

Nelder, J.A., Wedderburn, R.W., 1972. Generalized linear models. J. R. Stat. Soc. Ser. A Stat. Soc. 135, 370–384.

Norberg, A., Abrego, N., Blanchet, F.G., Adler, F.R., Anderson, B.J., Anttila, J., Araújo, M. B., Dallas, T., Dunson, D., Elith, J., Foster, S.D., Fox, R., Franklin, J., Godsoe, W., Guisan, A., O'Hara, B., Hill, N.A., Holt, R.D., Hui, F.K.C., Husby, M., Kålås, J.A., Lehikoinen, A., Luoto, M., Mod, H.K., Newell, G., Renner, I., Roslin, T., Soininen, J., Thuiller, W., Vanhatalo, J., Warton, D., White, M., Zimmermann, N.E., Gravel, D., Ovaskainen, O., 2019. A comprehensive evaluation of predictive performance of 33 species distribution models at species and community levels. Ecological Monographs 89, e01370. https://doi.org/10.1002/ecm.1370.

Ovaskainen, O., Abrego, N., 2020. Joint Species Distribution Modelling: with Applications in R. Cambridge University Press.

Panagos, P., Van Liedekerke, M., Borrelli, P., Köninger, J., Ballabio, C., Orgiazzi, A., Lugato, E., Liakos, L., Hervas, J., Jones, A., Montanarella, L., 2022. European Soil Data Centre 2.0: Soil data and knowledge in support of the EU policies. European Journal of Soil Science 73 (6), e13315. https://doi.org/10.1111/ejss.13315.

Pelosi, C., Pey, B., Hedde, M., Caro, G., Capowiez, Y., Guerniction, 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.apsoec.2013.10.005.

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., Mathieu, J., Mulder, C., Van Der Putten, W.H., Ramirez, K.S., Rillig, M. C., Russell, D., Rutgers, M., Thakur, M.P., De Vries, F.T., Wall, D.H., Wardle, D.A., Arai, M., Ayuke, F.O., Baker, G.H., Beauséjour, R., Bedano, J.C., Birkhofer, K., Blanchart, E., Fragsey, B., Bolger, T., Bradley, R.L., Callaham, M.A., Capowiez, Y., Caulfield, M.E., Choi, A., Crotty, F.V., Crumsey, J.M., Davalos, A., Diaz Cosin, D.J., Dominguez, A., Duhour, A.E., Van Eekeren, N., Emmerling, C., Falco, L.B., Fernández, R., Fonte, S.J., Fragoso, C., Franco, A.L.C., Fugère, M., Fusilero, A.T., Gholami, S., Gundale, M.J., López, M.G., Hackenberger, D.K., Hernández, L.M., Hishi, T., Holdsworth, A.R., Holmstrup, M., Hopfensperger, K.N., Lwanga, E.H., Huhta, V., Hurisso, T.T., Iannone, B.V., Iordache, R., Matela, R., Moos, J.H., Moreno, G., Morón-Ríos, A., Muys, B., Neiryntc, J., Norgrove, L., Novo, M., Nuutinen, V., Nuzzo, V., Mujeeb, Rahman P., Pansu, J., Paudel, S., Pérès, G., Pérez-Camacho, L., Piñeiro, R., Ponge, J.-F., Rashid, M.I., Rebollo, S., Rodeiro-Iglesias, J., Rodríguez, M.A., Roth, A.M., Rousseau, G.X., Rozen, A., Sayad, E., Van Schaik, L., Scharenbroch, B.C., Schirrmann, M., Schmidt, O., Schröder, B., Seeger, J., Shashkov, M.P., Singh, J., Smith, S.M., Steinwandter, M., Talavera, J.A., Trigo, D., Tsukamoto, J., De Valença, A.W., Vanek, S.J., Virto, I., Wackett, A.A., Warren, M.W., Wehr, N.H., Whalen, J.K., Wironen, M.B., Wolters, V., Zenkova, I.V., Zhang, W., Cameron, E.K., Eisenhauer, N., 2019. Global distribution of earthworm diversity. Science 366, 480–485. https://doi.org/10.1126/science.aax4851.

Poggiatto, G., Andréoletti, J., Pollock, L.J., Thuiller, W., 2025. Integrating food webs in species distribution models can improve ecological niche estimation and predictions. Ecography 2025, e07546. https://doi.org/10.1111/ecog.07546.

Poggiatto, G., Münkemüller, T., Bystrova, D., Arbel, J., Clark, J.S., Thuiller, W., 2021. On the interpretations of joint modeling in community ecology. Trends in Ecology & Evolution 36 (5), 391–401. https://doi.org/10.1016/j.tree.2021.01.002.

Poggio, L., de Sousa, Batjes, N.H., Heuvelink, G.B.M., Kemen, P., Ribeiro, E., Rossiter, D., 2021. SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty. Soil 7, 217–240. https://doi.org/10.5194/soil-7-217-2021.

Pollock, L.J., Tingley, P., Morris, W.K., Golding, N., O'Hara, R.B., Parris, K.M., Vesk, P. A., McCarthy, M.A., 2014. Understanding co-occurrence by modelling species simultaneously with a Joint Species Distribution Model (JSDM). Methods in Ecology and Evolution 5, 397–406. https://doi.org/10.1111/2041-210X.tj tree.12180.

Rutgers, M., Orgiazzi, A., Gardi, C., Römble, J., Jänsch, S., Keith, A.M., Neilson, R., Boag, B., Schmidt, O., Murchie, A.K., Blackshaw, R.P., 2016. Mapping earthworm communities in Europe. Applied Soil Ecology 97, 98–111.

Schirrmann, M., Giebel, A., Gleiniger, F., Pflanz, M., Lentschke, J., Dammer, K.-H., 2016. Monitoring Agronomic Parameters of Winter Wheat Crops with Low-Cost UAV Imagery. Remote Sensing 8 (9), 706. https://doi.org/10.3390/rs8090706.

Singh, J., Schädler, M., Demetrio, W., Brown, G.G., Eisenhauer, N., 2020. Climate change effects on earthworms - a review. Soil Organisms 91 (3), 113–137. https://doi.org/10.25674/so91iss3pp114.

Skidmore, A.K., Coops, N.C., Neinavaz, E., et al., 2021. Priority list of biodiversity metrics to observe from space. Nature Ecology & Evolution 5, 896–906. https://doi.org/10.1038/s41559-021-01451-x.

Spurgeon, D.J., Hopkin, S.P., 1996. The effects of metal contamination on earthworm populations around a smelting works: quantifying species effects. Applied Soil Ecology 4 (2), 147–160. https://doi.org/10.1016/0929-1393(96)00109-6.
Tang, L., Xue, Y., Chen, D., Gomes, C., 2018. Multi-Entity dependence learning with rich context via conditional variational auto-encoder. In: Proc. AAAI Conf. Artif. Intell., vol. 32. https://doi.org/10.1609/aaai.v32i1.11335.

Thuiller, W., Georges, D., Engler, R., Breiner, F., Georges, M.D., Thuiller, C.W., 2016. Package 'biomod2.' Species Distrib. Model. Ensemble Forecast. Framew. 10, 1600-0587.

Thuiller, W., Lafourcade, B., Engler, R., Araújo, M.B., 2009. BIOMOD – a platform for ensemble forecasting of species distributions. Ecography 32, 369–373. https://doi.org/10.1111/j.1600-0587.2008.05742.x.

Tibshirani, R., Walther, G., Hastie, T., 2001. Estimating the number of clusters in a data set via the gap statistic. J. R. Stat. Soc. Ser. B Stat. Methodol. 63, 411–423. https://doi.org/10.1111/1467-9868.00293.

Vidal, A., Blouin, M., Lubbers, I., Capowiez, Y., Sanchez-Hernandez, J.C., Calogiuri, T., Van Groenigen, J.W., 2023. The role of earthworms in agronomy: consensus, novel insights and remaining challenges. In: Advances in Agronomy. Elsevier, pp. 1–78. https://doi.org/10.1016/bs.agron.2023.05.001.

Violle, C., Navas, M.-L., Vile, D., Kazakou, E., Fortunel, C., Hummel, I., Garnier, E., 2007. Let the concept of trait be functional! Oikos 116, 882–892. https://doi.org/10.1111/j.0030-1299.2007.15559.x.

Wan, B., Barnes, A.D., Li, M., Song, Y., Yang, Q., Chen, X., Hu, F., Liu, M., 2025. Earthworms buffer the impacts of nitrogen enrichment on energy dynamics of soil

micro-food webs. Soil Biology and Biochemistry 202, 109705. https://doi.org/10.1016/j.soilbio.2024.109705.

Ward, J.H., 1963. Hierarchical grouping to optimize an objective function. Journal of the American Statistical Association 58, 236–244. https://doi.org/10.1080/01621459.1963.10500845.

Warton, D.I., Blanchet, F.G., O'Hara, R.B., Ovaskainen, O., Taskinen, S., Walker, S.C., Hui, F.K.C., 2015. So many variables: joint modeling in community ecology. Trends in Ecology & Evolution 30, 766–779. https://doi.org/10.1016/j.tree.2015.09.007.

Wilcoxon, F., 1945. Individual comparisons by ranking methods. Biomedical Bulletin 1, 80–83. https://doi.org/10.2307/3001968.

Wilkinson, D.P., Golding, N., Guillera-Arroita, G., Tingley, R., McCarthy, M.A., 2019. A comparison of joint species distribution models for presence-absence data. Methods in Ecology and Evolution 10, 198–211. https://doi-org.inee.bib.cnrs.fr/10.1111/2041-210X.13106.

Zeiss, R., Eisenhauer, N., Orgiazzi, A., Rillig, M., Buscot, F., Jones, A., Lehmann, A., Reitz, T., Smith, L., Guerra, C.A., 2022. Challenges of and opportunities for protecting European soil biodiversity. Conservation Biology 36, e13930. https://doi.org/10.1111/cobi.13930.

Zhang, Y., Wang, L., Duan, Z., Sun, M., 2019. Structure learning of Bayesian network based on adaptive thresholding. Entropy, 21 (7), 665. https://doi.org/10.3390/e21070665.