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 earthworm communities, broad-scale spatial patterns remain underexplored. We developed a multi-species, multi-task deep latent variable model to jointly predict the distribution of 77 earthworm species across metropolitan France, using historical (1960)(1961)(1962)(1963)(1964)(1965)(1966)(1967)(1968)(1969)(1970) and contemporary 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. Species clustering revealed distinct ecological strategies tied to climatic and land use gradients. Our study advances both the methodological and ecological understanding of soil 3 biodiversity. We demonstrate the utility of interpretable deep learning approaches for largescale soil fauna modeling and provide new insights into earthworm habitat specialization. These findings support improved soil biodiversity monitoring and conservation planning in the face of global environmental change. Earthworms are among the most ecologically significant invertebrates 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 burrowing activities, earthworms increase soil porosity, enhance water infiltration, and promote plant growth, making them crucial for ecosystem productivity and soil health (Arrazola-Vasquez 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 (Decaens 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, 2023). Improving our understanding of soil organism responses to environmental changes is therefore essential. While traditional earthworm species classifications based on morphological traits (Bouche, 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 (Grenie 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. We modelled earthworm distributions using historical calibration datasets covering diverse environmental conditions and evaluated the models on contemporary datasets . The calibration dataset originated from Marcel Bouche's extensive surveys across mainland France including Corsica (Bouche, 1972), recently digitized and taxonomically updated by (Gerard 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. During the 1960s, Marcel Bouche 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 (Bouche, 1972). The earthworm occurrence data (presence/absence) used here were digitized by (Gerard et al., 2025), resulting in 34,737 individual earthworm records from 1,399 sampling sites, covering 126 species and subspecies. For modelling, we selected taxa occurring in at least five sites with detailed habitat descriptions. The refined dataset included 1,346 sites and 77 taxa (species and subspecies sensu Bouche, 1972). We aggregated a dataset of geotagged earthworm observations from multiple sources sampled in the period. This dataset contains 1,126 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. 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 Bouche et al. (1972) and Gerard et al. (2025b), supplemented by GIS-derived environmental layers. Climate data were sourced from WorldClim (1960WorldClim ( -1990; 1 km resolution), while soil and hydrological data were extracted from the European Soil Database at 1km resolution . 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. 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 . See Supplementary Information for details on model formulation (SI, section A). We implemented MTEC in Python using Keras (Chollet, 2015). We trained MTEC 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, 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 x 2 cross-validation (Dietterich, 1998), with Area Under the (ROC-AUC) and the True Skill Statistics (TSS) as the primary evaluation metrics (SI, section B1). 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., 2016(Thuiller et al., , 2009) ) on the same calibration dataset as MTEC. Community composition predictions from these single-species models were combined through stacking (S-SDMs). Using environmental data from the projection dataset, we predicted community composition across France at a 1 km resolution using MTEC 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). 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 . 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. We computed the SHAP values for all environmental covariates for each species across samples of the calibration dataset . 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). 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). Comparison between MTEC and single-species distribution models revealed significantly improved predictive performance of MTEC on the calibration dataset . Predictive performance was generally lower for widespread taxa and variable for less prevalent taxa (SI, 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 andOctolasion 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 , suggesting that single-species SDMs may be more prone to overfitting. Land cover, precipitation extremes, and temperature seasonality emerged as the most influential abiotic factors across taxa . 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, Figure 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 C3 for other modeled taxa). 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 . 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 (Bouche 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 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 Bouche, 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 . 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 temperature and soil properties (SI, section C4). 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 multitask neural networks notably outperform traditional single-species SDMs, especially in modelling 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. 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 Knape (2016) and Harris ( 2015 Contrary to previous reviews indicating limited advantages of jSDMs over stacked SDMs (Pichler and Hartig, 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 (Marchan 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 finerscale 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 (Ou et al., 2021). The relative importance of environmental factors and the nature of species' responses differed among taxa. 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 climates 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 becomes the primary driver of species distributions. 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. Rather than total precipitation, specific seasonal metrics (e.g., wettest month precipitation) proved crucial. While extreme precipitation events are known to affect earthworms (Singh et al., 2019), the significance of seasonality has led some species to develop strategies for coping with high climate variability both during the day and across months. Conversely, several species, primarily epigeic ones, 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 (Bouche, 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 variables. Additionally, soil bulk density had a greater effect at the assemblage level compared to texture, a result that is not surprising given that bulk density was estimated from soil texture and structure classes using pedo-transfer rules (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. 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 Gerard 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 seasonality, likely avoiding adverse climate conditions by their ability to enter diapause in the deeper layers of the soil (Bouche, 1972). The second cluster included many species from the genera Zophoscolex, Prosellodrilus, 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 (Bouche, 1972). The final cluster consisted of a few species, 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 alterations 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 combination 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., 2019). The classification procedure on land cover response revealed significant 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, with the exception of a limited number of litter-dwelling species, such as D. octaedra. A few riparian species, such as E. tetraedra and H. gordioides, and several endogeics, 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 (Decaens 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;Decaens 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. Both climate conditions and land cover emerged as key predictors in our model, underscoring their combined influence on earthworm distributions. Anticipating future shifts will therefore require the development of scenarios that integrate the interplay between changing precipitation regimes and land cover dynamics. Beyond this, we demonstrated the potential of classifying species into response groups along two 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. 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. 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. 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: Tables In section 1.1, we describe the mathematical details of MTEC. In section 1.2, we describe the neural network architecture used to fit the model to environmental and community data. This represents variables that are observed in the dataset. A deterministic neural network (feature extractor) f ext transforms raw input covariates into a more informative vector representation (eq. 1), shared across species and referred to as environmental embedding, which captures relevant features driving the community composition while reducing noise and irrelevant variability in the raw input. x i = f ext (e i ) (1) At a given site i, the presence/absence of a given species j follows a Bernoulli distribution parameterized with a mean probability of occurrence  ij . The occurrence probability of a given species j varies across sites. At a given site i, it is modeled as the combined response (eq. 3) of the species j to the observed environment x i (environment effect) and unobserved or hidden factors h i (residual effect) which represent unmeasured confounders (e.g. microbial activity) or interspecies association drivers (e.g. resource for competition, microhabitat for facilitation). The parameters  and  represent the coefficients of species response to the environment and their latent factor loadings respectively. The link function g maps the linear response to the mean. We assume that species are independent (eq. 4), conditional on the environmental and hidden factors. y ij  y ij  |x i , h i (4) We further assume that the hidden factors follow an independent normal distribution (prior) with mean vector  p and variance vector  p 2 (eq. 5). h i ~ N( p ,  p 2  I) (5) The goal of inference is to maximize the marginal loglikelihood (eq. 6) of observed community composition conditional to the observed environment P(Y|E) w.r.t the generative parameters  = (W e , , ), denoted as p  (Y). log p  (Y ) = log  p  (Y, H) dH ( 6) Given that hidden factors are de facto unobserved, this requires computing an intractable integral. To circumvent this issue, we introduce a variational approximation q(H|Y) of the posterior distribution of the latent variables p(H|Y) within a class Q of distributions. Here, we choose the family of factorized Gaussian distributions (eq. 7). From the previous equation, the number of posterior parameters (means and variances) increases linearly with N the number of observations O(2LN). To ensure the scalability of the inference, we follow the AEVB algorithm and use a sufficiently complex parametric function f ENC  y s  ( s ,  s 2 ). In practice, f ENC is a neural network parameterized with weights , referred to as recognition or probabilistic encoder network. The combination of the encoder and the generative model (namely, decoder), leads to a variational autoencoder (VAE) [Kingma and Welling, 2013]. This transforms the inference problem to an optimization problem where we maximize the likelihood J w.r.t the global parameters  minus the information loss D induced by the posterior approximation w.r.t the encoder parameters . The KL divergence between the prior and approximate posterior of the latent factors, both normally distributed, can be computed in closed form. We further introduce a regularization term on the model parameters in the final optimization objective (eq. 9) to prevent overfitting.
The optimization objective is to minimize the total cost J reg , which includes the reconstruction loss (measuring the difference between the predicted and actual species presence), the KL divergence loss (penalizing the divergence between the approximate and true posterior distributions of latent factors), and a regularization term (preventing overfitting by controlling model complexity). In the context of this paper, the reconstruction loss is the balanced binary crossentropy loss as the reconstruction loss, whereas we use Elasticnet regularization loss for the parameter regularization term. We define the following posterior statistics over the latent factors H, assuming the dimension of the latent space is L and the number of data samples is N From there, the species-by-species residual covariance matrix can be formulated as: We use the Graphical lasso 2 algorithm to estimate a sparse inverse for the posterior residual covariance matrix  r =  r 1 with a sample size of N. This corresponds to the precision matrix which encodes conditional dependencies between species in the latent space. Gaussian Random Field that we refer to as the species association network. It is an undirected graphical model in which the absence of a link represents pairwise independence conditional to other species. An environmental feature encoder: a fully connected neural network with parameters (weights and biases) which takes as input the raw environmental covariates and transforms them to extract a shared embedding with relevant features for all species. We use category embedding layers for categorical features and dense layers for numerical features. A multi-task decoder: a fully connected neural network that predicts multispecies probabilities of presences as a function of the shared environmental embedding and the latent factors. The latter are sampled from the posterior distribution during training on known communities and from the prior distribution on new sites. We fitted MTEC to our dataset of earthworm communities and environmental covariates across training sampling sites. We used a probit link g =  to accommodate the binary presence/absences. Given the limited size of the calibration dataset, we used a single layer in the probabilistic encoder and fixed the number of latent factors to 3 to facilitate the visualization. Furthermore, we applied a combined ridge and lasso (elastic net) regularization to the generative  and variational  parameters to control model complexity. We experimented with different settings to process environmental features both at the data preprocessing stage (dimensionality reduction) and within the feature encoder (architecture complexity). (3) preprocessing via principal component analysis. To find the optimal configuration of network complexity and transferability amongst prediction tasks, we varied the shared feature encoder network depth, and the width (number of neurons) of each layer. Combining learning architectures with the preprocessing modes yielded 6 configurations. We trained all 6 models on the full calibration dataset while varying the regularization strength in the range {1,2,5,10}  104. We used a 5x2 cross-validation on the calibration dataset coupled with a paired t-test (Dietterich, 1998) to compare models in pairs. We report average performances on the calibration dataset (ROC-AUC) and the evaluation dataset (recall) as well as optimal regularization settings for each configuration in We fit the model using the Adam optimization algorithm for at most 400 iterations (epochs), using a batch size of 32. To prevent overfitting, we used early stopping to stop the training whenever the objective function of the validation set stopped improving after more than 10 epochs. Like most species' occurrence data, the earthworm dataset used in this study exhibits a long-tailed species prevalence distribution where a vast majority of species occur in less than 10% of the community samples while. In this section, we report species-wise predictive performances both on the calibration dataset (cross-validation) and on the evaluation dataset. Here, we examine how the cross-validation predictive performances in the calibration dataset vary along a gradient of species prevalences from rarest to most prevalent. We find that both MTEC and machine learning based SDMs exhibit generally decreasing ROC-AUC with prevalence , except GLMs. Furthermore, MTEC obtained better scores for rare species, likely due to imbalance correction strategies during training. Additionally, single-species SDM models struggled to converge for some species with very low prevalence (fewer than 10 records). For more prevalent species, performances were variable, and this can potentially be attributed to the difference in niche complexity between species. In section A1.4, we described the process to infer the earthworm association network from the latent factor estimation. Figure 6 shows the resulting association network. This network illustrates associations between earthworm taxa that are not accounted for by their shared responses to environmental variables. Instead, these associations are structured by latent factors, which may reflect unmeasured environmental influences (e.g., microhabitat conditions) or biotic interactions (e.g., microbial activity). (epigeic, endogeic, anecic). Edges between two taxa nodes represent significant associations. Association strength (partial correlation) normalized between -1 and 1, is shown on each node in red (resp. blue) for positive (resp.negative) associations. The applied method identified 31 significant pairwise associations, resulting in a network density of 0.015. These associations involved only 31 of the 77 modeled taxa, and the network visualization was accordingly restricted to these taxa. Notably, the partial correlations-representing the strength of associations and visualized as edge attributes-were relatively weak and uniform, ranging from 0.14 to 0.22 in absolute value. This suggests that, overall, earthworm distributions are primarily governed by abiotic factors. The resulting network comprised three connected components. The largest component encompassed 80% of the nodes and 87% of the edges, and featured a hub formed by four of the most prevalent taxa: Lumbricus castaneus, Satchellius mammalis, and two subspecies of Aporrectodea caliginosa-A. caliginosa caliginosa, which is widespread, and A. caliginosa meridionalis, which is restricted to southern regions. These two subspecies were rarely co-occurring, as evidenced by their negative association. Additionally, a positive ternary association was observed among Lumbricus rubellus, Aporrectodea longa ripicola, and Dendrobaena rubidus rubidus-three epigeic to epianecic species predominantly found in northern areas. Most species represented in the network were classified as epigeic, epi-endogeic, or epi-anecic.
This pattern may reflect the influence of an unmodeled above-ground factor, such as vegetation type, which affects litter quality (e.g., leaf carbon-tonitrogen ratio) and its seasonality (e.g., deciduous vs. evergreen trees). Interestingly, species positioned toward the "epigeic" end of Bouche's ecological categories triangle (including epigeics, epi-endogeics, and epi-anecics) were predominantly involved in positive associations with endogeic and anecic species. This suggests a potential facilitative interaction, where surface-dwelling species may benefit from the use of burrows and galleries constructed by soil-dwelling taxa. These subterranean structures could serve as refuges from above-ground stressors such as cold temperatures, wind, drought, or predation by birds, mammals, or reptiles. Notably, endemic or narrow-range species were almost entirely absent from the association network, which was predominantly composed of widespread, peregrine taxa. This pattern suggests that the success of eurytopic species in dispersing and colonizing previously inaccessible habitats may be attributed not only to their inherent competitive and dispersal capacities [Mathieu and Jonathan Davies, 2014] 3 but also to positive interspecific interactions, as evidenced by the observed association patterns. We computed the SHAP values for all environmental covariates for each species across samples of the calibration dataset. 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 the effect of the environmental factors. We provide in the supplementary file LimitingFactorMaps.zip all limiting factor maps for modeled taxa. Pairwise distances among species based on SHAP values were calculated, followed by hierarchical clustering (Ward, 1963) to unravel response groups for each feature group (precipitation, land cover and soil physico-chemical properties). The optimal number of clusters was identified using the GAP statistic (Tibshirani et al., 2001). MH designed the study. SS-M developed the analytical framework with guidance from EG and WT. MH, TD, YC, and SG provided the earthworm datasets. SL and colleagues contributed the environmental data layers. TD, MH, DM, YC, SG and CM participated in the interpretation and ecological contextualization of the results. MH, and SS-M drafted the initial manuscript. All authors contributed to manuscript revisions and approved the final version. 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 deg 101134954) and NaturaConnect (N deg 101060429) projects and the MIAI@Grenoble Alpes (ANR--19--P3IA--0003). SG PhD thesis was The authors declare no conflict of interest. This study did not involve human participants or experimental animals and therefore did not require ethics approval. We ensured all taxa appeared in each training set and weighted likelihood presences (positive loss) by their odds relative to absences (negative loss). To prevent bias fitting, we initialized the mean prediction intercepts to each taxon's logitprevalence, while other parameters used the Glorot-uniform initializer (Glorot & Bengio, 2010). -off between tolerance to acidity and ability to degrade the carbon. The first two components of the PCA of the responses to pH and C:N ratio explained  70% of the variation and ordinated species according to tolerance to soil acidity and organic matter quality respectively. At one end, we find the typical forest assemblages (cluster 3). They include epigeics and epi-endogeics that survive in spite of acidic soils by taking advantage of the rich organic matter in the litter. In this group, we also found endogeics of the Zophoscolex genus, described by Qiu and Bouche [1998] as acid-tolerant. At the other end are taxa that prefer neutral, slightly alkaline environments (cluster 1) despite poorer organic matter. An intermediate group (cluster 2) is made of species which showed little sensitivity to pH values (SHAP values closed to zero), or with preference for the average pH value. Riparian taxa stand out in a distinct cluster (cluster 4) associated to neutral wetland soils. MH designed the study. SS-M developed the analytical framework with guidance from EG and WT. MH, TD, YC, and SG provided the earthworm datasets. SL and colleagues contributed the environmental data layers. TD, MH, DM, YC, SG and CM participated in the interpretation and ecological contextualization of the results. MH, and SS-M drafted the initial manuscript. All authors contributed to manuscript revisions and approved the final version. 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 deg 101134954) and NaturaConnect (N deg 101060429) projects and the MIAI@Grenoble Alpes (ANR--19--P3IA--0003). SG PhD thesis was The authors declare no conflict of interest. This study did not involve human participants or experimental animals and therefore did not require ethics approval. We ensured all taxa appeared in each training set and weighted likelihood presences (positive loss) by their odds relative to absences (negative loss). To prevent bias fitting, we initialized the mean prediction intercepts to each taxon's logitprevalence, while other parameters used the Glorot-uniform initializer (Glorot & Bengio, 2010). -off between tolerance to acidity and ability to degrade the carbon. The first two components of the PCA of the responses to pH and C:N ratio explained  70% of the variation and ordinated species according to tolerance to soil acidity and organic matter quality respectively. At one end, we find the typical forest assemblages (cluster 3). They include epigeics and epi-endogeics that survive in spite of acidic soils by taking advantage of the rich organic matter in the litter. In this group, we also found endogeics of the Zophoscolex genus, described by Qiu and Bouche [1998] as acid-tolerant. At the other end are taxa that prefer neutral, slightly alkaline environments (cluster 1) despite poorer organic matter. An intermediate group (cluster 2) is made of species which showed little sensitivity to pH values (SHAP values closed to zero), or with preference for the average pH value.
Riparian taxa stand out in a distinct cluster (cluster 4) associated to neutral wetland soils.