# Effects of climate on the distribution and conservation of commonly observed European earthworms

Romy Zeiss, Maria J I Briones, Jérome Mathieu, Angela Lomba, Jessica Dahlke, Laura-fiona Heptner, Gabriel Salako, Nico Eisenhauer, Carlos A Guerra

► To cite this version:

Romy Zeiss, Maria J I Briones, Jérome Mathieu, Angela Lomba, Jessica Dahlke, et al. Effects of climate on the distribution and conservation of commonly observed European earthworms. Conservation Biology, 2024, 38, 10.1111/cobi.14187. hal-04946106

HAL Id: hal-04946106

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

Submitted on 13 Feb 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
<figure><img src="image_2.png" /></figure>

# Effects of climate on the distribution and conservation of commonly observed European earthworms

Romy Zeiss<fnref n="1" />,<fnref n="2" /> ID

Maria J. I. Briones<fnref n="3" /> ID

Jérôme Mathieu<fnref n="4" /> ID

Angela Lomba<fnref n="5" />,<fnref n="6" />,<fnref n="7" /> ID

Jessica Dahlke<fnref n="1" />,<fnref n="8" /> ID

Laura-Fiona Heptner<fnref n="1" />,<fnref n="2" /> ID

Gabriel Salako<fnref n="9" />,<fnref n="10" /> ID

Nico Eisenhauer<fnref n="1" />,<fnref n="2" /> ID

Carlos A. Guerra<fnref n="1" />,<fnref n="2" />,<fnref n="8" /> ID

**Correspondence**

Romy Zeiss, German Centre for Integrative Biodiversity Research (iDiv), Puschstraße 4, 04103 Leipzig. Email: romyzeiss@gmail.com

**Article impact statement:** Predicted distribution of commonly observed earthworm species indicates climate vulnerability and insufficient conservation

**Funding information**

Fundação para a Ciência e a Tecnologia, Grant/Award Number: in the context of the Transitory Norm DL57/2016/CP1440/CT0001; Deutsche Bundesstiftung Umwelt, Grant/Award Number: 20021/752; Deutsches Zentrum für integrative Biodiversitätsforschung Halle-Jena-Leipzig, Grant/Award Numbers: DFG FZT 118, 202548816; sDiv, Grant/Award Number: SFW9.02; France-Stanford Center for Interdisciplinary Studies

**Abstract**

Belowground biodiversity distribution does not necessarily reflect aboveground biodiver-
sity patterns, but maps of soil biodiversity remain scarce because of limited data availability.
Earthworms belong to the most thoroughly studied soil organisms and—in their role as
ecosystem engineers—have a significant impact on ecosystem functioning. We used species
distribution modeling (SDMs) and available data sets to map the spatial distribution of
commonly observed (i.e., frequently recorded) earthworm species (Annelida, Oligochaeta)
across Europe under current and future climate conditions. First, we predicted potential
species distributions with commonly used models (i.e., MaxEnt and Biomod) and esti-
mated total species richness (i.e., number of species in a 5 × 5 km grid cell). Second, we
determined how much the different types of protected areas covered predicted earthworm
richness and species ranges (i.e., distributions) by estimating the respective proportion of
the range area. Earthworm species richness was high in central western Europe and low
in northeastern Europe. This pattern was mainly associated with annual mean temperature
and precipitation seasonality, but the importance of predictor variables to species occur-
rences varied among species. The geographical ranges of the majority of the earthworm
species were predicted to shift to eastern Europe and partly decrease under future climate
scenarios. Predicted current and future ranges were only poorly covered by protected areas,
such as national parks. More than 80% of future earthworm ranges were on average not
protected at all (mean [SD] = 82.6% [0.04]). Overall, our results emphasize the urgency of
considering especially vulnerable earthworm species, as well as other soil organisms, in the
design of nature conservation measures.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is
properly cited.

© 2024 The Authors. *Conservation Biology* published by Wiley Periodicals LLC on behalf of Society for *Conservation Biology*.

<footnote n="1">¹German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig, Leipzig, Germany</footnote>

<footnote n="2">²Institute of Biology, Leipzig University, Leipzig, Germany</footnote>

<footnote n="3">³Departamento de Ecologia y Biología Animal, Universidad de Vigo, Vigo, Spain</footnote>

<footnote n="4">⁴Sorbonne Université, CNRS, IRD, INRAE, Université Paris Est Créteil, Université de Paris Cité, Institute of Ecology and Environmental Sciences of Paris (IEES-Paris), Paris, France</footnote>

<footnote n="5">⁵CIBIO, Centro de Investigação em Biodiversidade e Recursos Genéticos, InBIO Laboratório Associado, Campus de Vairão, Universidade do Porto, Vairão, Portugal</footnote>

<footnote n="6">⁶Departamento de Biologia, Faculdade de Ciências, Universidade do Porto, Porto, Portugal</footnote>

<footnote n="7">⁷BIOPOLIS Program in Genomics, Biodiversity and Land Planning, CIBIO, Campus de Vairão, Vairão, Portugal</footnote>

<footnote n="8">⁸Martin Luther University Halle-Wittenberg (MLU), Naturwissenschaftliche Fakultät 1, Halle (Saale), Germany</footnote>

<footnote n="9">⁹Soil Zoology Division, Senckenberg Museum of Natural History, Görlitz, Germany</footnote>

<footnote n="10">¹⁰Department of Environmental Management and Toxicology, Kwara State University, Malete, Nigeria</footnote>
<figure><img src="image_3.png" /></figure>

KEYWORDS

soil biodiversity, SDM, protected area, Lumbricidae, climate change, vulnerability, modeling, future

Efectos del clima sobre la distribución y conservación de la lombriz de tierra europea

**Resumen:** La distribución de la biodiversidad del subsuelo no refleja necesariamente los patrones de biodiversidad, pero los mapas de la biodiversidad del suelo aún son escasos debido a la disponibilidad limitada de datos. Las lombrices son uno de los organismos del suelo más estudiados a detalle—en su papel de ingenieros del ecosistema—y tienen un impacto significativo sobre el funcionamiento de ecosistema. Usamos modelos de distribución de especies (MDE) y conjuntos de datos disponibles para mapear la distribución espacial de las especies (Annelida, Oligochaeta) de lombrices más observadas (es decir, registradas con frecuencia) en toda Europa bajo el clima actual y el futuro. Primero pronosticamos la distribución potencial de las especies con modelos de uso común (MaxEnt y Biomod) y estimamos la riqueza total de especies (número de especies en una cuadrícula de 5 × 5 km). Después determinamos cuánto pronosticaban los diferentes tipos de áreas protegidas contempladas la riqueza de lombrices y la distribución de las especies mediante la estimación de la proporción respectiva del rango del área. La riqueza de especies fue alta en el occidente central y baja en el noreste de Europa. Este patrón estuvo asociado principalmente con la temperatura media anual y la estacionalidad de la precipitación, aunque la importancia de las variables de pronóstico para la presencia de la especie varió entre especies. Se pronosticó que la distribución geográfica de la mayoría de las especies cambiaría al este de Europa y disminuiría parcialmente bajo los escenarios climáticos futuros. El pronóstico de la distribución actual y futura contaba con una cobertura deficiente de las áreas protegidas, como los parques nacionales. En promedio, más del 80% de la distribución futura de las lombrices no estaba protegido (promedio [SD] = 82.6% [0.04]). En general, nuestros resultados destacan la urgencia por considerar a las especies vulnerables de lombrices, así como a otros organismos del suelo, en el diseño de las medidas de conservación.

PALABRAS CLAVE

área protegida, biodiversidad del suelo, cambio climático, futuro, Lumbricidae, MDE, modelo, vulnerabilidad

INTRODUCTION

Soil harbors an important portion of the world's biodiversity, including organisms of various sizes (Bardgett & van der Putten, 2014; Decaëns et al., 2006; FAO et al., 2020; Guerra et al., 2020). This multidimensionality drives many ecosystem functions (e.g., bioturbation, soil respiration, aggregate stability) that are key for the sustainability of terrestrial ecosystems (Bardgett & van der Putten, 2014; Eisenhauer et al., 2020; FAO et al., 2020). Nevertheless, the distribution and main drivers of most soil taxa are relatively unknown when compared with other biodiversity groups, such as vertebrates, birds, or plants (Phillips et al., 2017), and information about the effects and coverage of nature conservation areas in relation to soil biodiversity is scarce (Cameron et al., 2018; Ciobanu et al., 2019; Guerra et al., 2020; Zeiss et al., 2022). Attempts to map soil biodiversity have shown that its complex spatial patterns do not resemble those of aboveground organisms (Cameron et al., 2019; Guerra et al., 2022; Mathieu & Davies, 2014). As a result, protecting aboveground species does not necessarily provide protection to belowground organisms (Zeiss et al., 2022). Recent large-scale studies describe the main drivers of large taxonomic groups distributions (e.g., earth-

worms, nematodes, fungi, bacteria) (Delgado-Baquerizo et al., 2018; Köninger et al., 2023; Phillips et al., 2019; Tedersoo et al., 2014; van den Hoogen et al., 2019) and highlight the effects of climate (in particular temperature and precipitation). Regional and local-scale studies confirm the generality of these main drivers of distribution (Orgiazzi et al., 2016; Desie et al., 2020).

Soil biodiversity mapping is mostly done for single taxonomic
groups (e.g., fungi and nematodes) (Tedersoo et al., 2014; van
den Hoogen et al., 2019) and functional clusters (e.g., bacteria)
(Delgado-Baquerizo et al., 2018) but not at the species level or
systematically within each taxonomic group. Beyond the bio-
geographical interest, species-level information allows for more
effective conservation actions and the evaluation of the conser-
vation status and threats (Phillips et al., 2017; Rodrigues et al.,
2006). From all soil animal groups, earthworms are one of the
best surveyed invertebrates in Europe (Phillips et al., 2017),
with abundant species-level data available across several regions
(for example, Edaphobase.org, 2021; GBIF.org, 2021). Because
earthworms create biogenic structures and environmental con-
ditions that are suitable for other soil organisms, they are
considered ecosystem engineers (Lavelle, 2011; Le Bayon et al.,
2017). In agroecosystems, for example, earthworms strongly
affect crop yield, plant growth (i.e., aboveground biomass), and other plant traits (Bardgett & van der Putten, 2014; Cunha et al., 2016; Eisenhauer et al., 2019; Kooch & Jalilvand, 2008; van Groenigen et al., 2014). As for most soil taxa, data availability of most earthworm species is often limited (e.g. German earthworms; Lehmitz et al., 2016), something that can be overcome, to some extent, through systematic standardization of all available information (e.g., GBIF.org, 2021; Jetz et al., 2019) and ecological modeling, such as with species distribution models (SDMs) (Navarro et al., 2017).

We tackled the aforementioned knowledge gaps by investigating the distribution (as potential species' geographical range) under climate change of the commonly observed European earthworm species. We focused on commonly observed species, that is species with a sufficient number of occurrences, to avoid model overfitting. Europe is especially interesting because of comparable previous studies in which richness values were extrapolated instead of single-species models being built (Rutgers et al., 2016); its broad range of environmental conditions (e.g., mean annual temperatures from -14 to 21°C and elevation between from -102 to 4454 m asl [Supporting Information Appendix S1]); and the availability of environmental and species occurrence datasets. To properly assess potential changes in the conservation status of each earthworm species, we also assessed their distribution across multiple future scenarios of climate change (Pörtner & Roberts, 2022) and examined their distribution overlap with current protected areas. Based on the protection coverage and future predictions, we also identified species whose ranges only sparsely overlap with current protected areas and that are strongly affected by climatic factors.

## METHODS

Our methodology included 4 steps (Figure 1): building the species distribution database; importing and harmonizing environmental data; calibrating distribution models to select the main predictors; and predicting the potential distribution of a set of well-documented European earthworm species with an SDM. We used a uniform grid system for Europe that covered the 27 member states of the European Union and Great Britain. Models were built at a resolution of ~2 km and predicted at a ~5-km resolution, a trade-off that allowed optimization of computational power and time (one-fourth of time compared with ~1-km resolution) but maintained sufficient fine spatial resolution to identify regional environmental drivers and sufficient numbers of species records. All analyses were performed in R 4.2.0 (R Core Team, 2020) and RStudio 2022.02.3 (RStudio Team, 2020), unless otherwise indicated.

### Species occurrence data

We collected occurrence data of earthworm species from the Global Biodiversity Information Facility (GBIF) (GBIF.org, 2022), Edaphobase (2021), and other available large datasets (Phillips et al., 2021) (details in Supporting Information Appendix S2 and S3).

Data sets are available at the iDiv data repository (https://doi.org/10.25829/idiv.3524-gqvs4z). The GBIF records provided by iNaturalist were not considered in the final dataset because species' identification has proven unreliable (Di Cecco et al., 2021; McMullin & Allen, 2022) (Supporting Information Appendix S3). In addition to the 3 data sets, we used 2 unpublished datasets from European researchers, Jérôme Mathieu and Carlos A. Guerra (detailed description in Supporting Information Appendix S3, data available on request). All species names were harmonized by taxonomic expert knowledge (Maria J.I. Briones) and based on previous standardization (Phillips et al., 2019).

After data collection, only records collected from 1970 to March 2022 were retained to estimate recent rather than past species occurrence patterns. The R package CoordinateCleaner was used to exclude data with common spatial and temporal errors (Zizka et al., 2019). We accounted for country capitals, country centroids, equal numbers for latitude and longitude, GBIF headquarters, plain zeros, and default sea areas. We manually checked whether the exclusion of a record was appropriate. Most of the flagged records were duplicates (i.e., Edaphobase data in GBIF) or were missing coordinates. We excluded records from GBIF with coordinate uncertainty of >1 km, that were not observational (i.e., but including records of living specimens, human observations, or preserved specimens); that were based on <1 observation; and that described taxonomic levels higher than species. The combined data set contained 98,732 occurrence records of 142 unique earthworm species collected across 45 European countries (Supporting Information Appendices S2 & S4, Figure 1).

The GBIF provided most of the records (54,883 records of 69 species) followed by Mathieu et al. 2022 (25,592 records, 77 species) and Edaphobase (13,054 records, 53 species). The occurrence records were spatially thinned to adjust for sampling bias caused by varying sampling density and resolution (e.g., Kramer-Schadt et al., 2013; Steen et al., 2021). Spatial thinning was done by transferring the records into a ~2 × 2 km grid system in WGS84 (Hijmans, Bivand et al., 2022) and back transforming the grid into point data. After spatial thinning, our final dataset contained 26,389 records of 127 earthworm species.

To avoid overfitting and increase the robustness of the models, we focused on the 41 species that were observed in at least 10 grid cells. We generated background data for the focal species (i.e., number of records ≥10) with the R package biomod2 (Thuiller et al., 2016) by randomly sampling 10,000 grid cells within the spatial extent of the environment for each species (Barbet-Massin et al., 2012). Such a procedure has been used widely with presence-only data sets because it allows implementation of all modeling algorithms available in the biomod2 package.

### Environmental data

Earthworm distributions are driven by multiple factors, including soil moisture, climate, chemical and physical soil properties, and topographic factors (Desie et al., 2020; Mathieu et al.,
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**FIGURE 1** Workflow of study of earthworm distribution and protection in Europe under current and future climatic conditions (small squares, predictor variables [sources in Table 1]; *n*, number of occurrences at 4-km² resolution per species). MaxEnt models identify the 10 most important drivers of species distribution, and species distribution models (SDMs) in biomod2 predict current and future distributions based on the most important 10 drivers only. Step 4 does not show the evaluation of the effect of future precipitation and temperature, which resulted in 3 additional datasets for each species and climate scenario (i.e., 45 future models per species).</figcaption></figure>

2022b; Mueller et al., 2016; Phillips et al., 2019; Singh et al., 2019). Soil moisture is related to the availability of water through precipitation and can be expressed as aridity (proportion of precipitation to evapotranspiration). Chemical and physical soil properties encompass variables such as pH, base cations, heavy metal concentrations, soil texture, and amount of organic matter, a food resource for earthworms. Topographic factors drive the large-scale distribution of many organisms, including earthworms (Phillips et al., 2019). We included distance to the coast in addition to commonly used predictors of earthworm distributions because only a few earthworm species are tolerant of high salinity and can thus exist near the sea (e.g., *Eisenia fetida*, Owojori et al., 2008).

We considered 4 groups of environmental descriptors: climatic, land-use, topographic, and edaphic factors (n = 29) (Table 1). Given that most of the data on earthworm occurrences were collected around 2007 (mean = 2004, median = 2007), we used the CORINE Land Cover dataset from 2012 as a land-cover baseline. While other land-cover datasets exist (Buchhorn et al., 2020; Ellis et al., 2013; ESDAC, 2020),

we used CORINE given its spatial and coverage accuracy (Aune-Lundberg & Strand, 2021; Caetano et al., 2006; European Environmental Agency, 2019; Torma & Harma. 2004). Environmental variables were harmonized when necessary (Supporting Information Appendix S1). All environmental data were reprojected or resampled onto a 2- and 5-km grid systems and standardized by dividing mean-centered values by their respective standard deviations. We used ArcGIS 10.7.1 (Bajjali, 2018) to reproject variables into WGS84, if necessary, but merged them into one table in R to avoid conflicts with missing data. Areas with ≥1 missing value in environmental variables were excluded.

## Model calibration

We used the R package usdm (Naimi, 2017) to check for correlation and collinearity between variables. Variables with r>0.8 (correlation coefficient of Pearson) and variable inflation factor (VIF) ≥10 were excluded from further analyses (n = 4)
<table><caption>TABLE 1 Environmental predictors used for species distribution modeling (SDM) of European earthworm species.ª</caption><thead><tr><th>Category</th><th>Abbreviation</th><th>Description</th><th>Unit</th><th>Time scale</th><th>Reference</th></tr></thead><tbody><tr><td>Climate</td><td>Aridity<sup>b</sup></td><td>aridity index</td><td>Value</td><td>1981–2010</td><td>Karger et al. (2017)</td></tr><tr><td>Climate</td><td>MAP</td><td>annual precipitation (bioclim 12)</td><td>mm year<sup>-1</sup></td><td>1981–2010</td><td>Karger et al. (2017)</td></tr><tr><td>Climate</td><td>MAP_Seas<sup>c</sup></td><td>precipitation seasonality (coefficient of variation) (bioclim 15)</td><td>kg m<sup>-2</sup></td><td>1981–2010</td><td>Karger et al. (2017)</td></tr><tr><td>Climate</td><td>MAT<sup>c</sup></td><td>annual mean temperature (bioclim 1)</td><td>°C *10</td><td>1981–2010</td><td>Karger et al. (2017)</td></tr><tr><td>Climate</td><td>MAT_Seas<sup>b</sup></td><td>temperature seasonality (standard deviation) (bioclim 4)</td><td>°C *10</td><td>1981–2010</td><td>Karger et al. (2017)</td></tr><tr><td>Climate</td><td>Snow<sup>b</sup></td><td>number of days with snow cover, based on gap-filled fractional snow cover (GFSC)</td><td>% + 100</td><td>2018–2021</td><td>Copernicus<sup>d</sup></td></tr><tr><td>LandCover</td><td>Agriculture<sup>c</sup></td><td>proportion of area covered by agricultural land based on CORINE land cover (CLC) 2012: class 2 except 231 (pasture)</td><td>Proportion</td><td>2012</td><td>Copernicus<sup>d</sup></td></tr><tr><td>LandCover</td><td>Dist_Urban</td><td>distance to urban areas based on CLC 2012: class 1 except 122 (roads) and 131 (mineral extraction)</td><td>m</td><td>2012</td><td>Copernicus<sup>d</sup></td></tr><tr><td>LandCover</td><td>Forest_Coni</td><td>proportion of area covered by coniferous forests based on CLC 2012: class 312 and 334 (burnt)</td><td>Proportion</td><td>2012</td><td>Copernicus<sup>d</sup></td></tr><tr><td>LandCover</td><td>Forest_Deci</td><td>proportion of area covered by deciduous forests based on CLC 2012: class 311 and 313 (mixed)</td><td>Proportion</td><td>2012</td><td>Copernicus<sup>d</sup></td></tr><tr><td>LandCover</td><td>NDVI</td><td>normalized difference vegetation index</td><td>Value</td><td>1999–2019</td><td>Copernicus Service Information 2022</td></tr><tr><td>LandCover</td><td>Pastures</td><td>proportion of area covered by pastures based on CLC 2012: class 231 (pasture) and 321(natural grassland) and 333 (sparsely open vegetation)</td><td>Proportion</td><td>2012</td><td>Copernicus<sup>d</sup></td></tr><tr><td>LandCover</td><td>Pop_Dens<sup>c</sup></td><td>UN WPP-adjusted population density</td><td>number of people/km<sup>2</sup></td><td>2000–2020</td><td>Doxsey-Whitfield et al. (2015)</td></tr><tr><td>LandCover</td><td>Shrubland</td><td>proportion of area covered by shrubland based on CLC 2012: class 32 except 321 (natural grassland)</td><td>Proportion</td><td>2012</td><td>Copernicus<sup>d</sup></td></tr><tr><td>Location</td><td>Aspect</td><td>GISCO aspect (i.e., direction of slope)</td><td>Value</td><td>2000</td><td>Copernicus<sup>d</sup></td></tr><tr><td>Location</td><td>Dist_Coast<sup>c</sup></td><td>distance from the coast based on EU-Hydro - River Network Database: class 52 (marine water)</td><td>m</td><td>2006–2012</td><td>Copernicus<sup>d</sup></td></tr><tr><td>Location</td><td>Dist_River</td><td>distance from rivers based on EU-Hydro-River Network Database: class 51 (inland water)</td><td>m</td><td>2006–2012</td><td>Copernicus<sup>d</sup></td></tr><tr><td>Location</td><td>Elev<sup>c</sup></td><td>elevation based on European digital elevation model</td><td>m</td><td>2000</td><td>Copernicus<sup>d</sup></td></tr><tr><td>Location</td><td>Slope</td><td>Slope</td><td>Value</td><td>2000</td><td>Copernicus<sup>d</sup></td></tr><tr><td>Soil</td><td>CEC<sup>c</sup></td><td>cation exchange capacity</td><td>cmol kg<sup>-1</sup></td><td>2012</td><td>Ballabio et al. (2019)</td></tr><tr><td>Soil</td><td>Clay+Silt<sup>c</sup></td><td>content of clay (%) + silt (%) (i.e., fine soil texture)</td><td>%</td><td>2009</td><td>Ballabio, Panagos, and Monatanarella (2016)</td></tr><tr><td>Soil</td><td>Cu</td><td>copper distribution in top soils</td><td>Cu mg kg<sup>-1</sup></td><td>2009</td><td>Ballabio et al. (2018) &amp; Panagos et al. (2018)</td></tr><tr><td>Soil</td><td>Hg</td><td>Hg distribution</td><td>Hg microg kg<sup>-1</sup></td><td>2009</td><td>Ballabio et al. (2021) &amp; Panagos et al. (2021)</td></tr></tbody></table>

(Continues)
<figure><img src="image_6.png" /></figure>

<table><caption>TABLE 1 (Continued)</caption><thead><tr><th>Category</th><th>Abbreviation</th><th>Description</th><th>Unit</th><th>Time scale</th><th>Reference</th></tr></thead><tbody><tr><td>Soil</td><td>Moisture</td><td>soil moisture index (SMI) (i.e., averaged daily moisture conditions based on hydrological rainfall-runoff model LISFLOOD)</td><td>Value</td><td>1995–2021</td><td>Cammalleri et al. (2017)</td></tr><tr><td>Soil</td><td>N</td><td>nitrogen content</td><td>N g kg<sup>-1</sup></td><td>2012</td><td>Ballabio et al. (2019)</td></tr><tr><td>Soil</td><td>p<sup>c</sup></td><td>phosphorus content</td><td>P mg kg<sup>-1</sup></td><td>2012</td><td>Ballabio et al. (2019)</td></tr><tr><td>Soil</td><td>pH<sup>c</sup></td><td>pH in water (H<sub>2</sub>O)</td><td>ph in H<sub>2</sub>O</td><td>2012</td><td>Ballabio et al. 2019</td></tr><tr><td>Soil</td><td>SOC</td><td>soil organic carbon</td><td>t ha<sup>-1</sup>/ kg m<sup>-2</sup></td><td>1990–2019</td><td>FAO 2002; Global Soil Organic Carbon Map-GSOCmap 1.6</td></tr><tr><td>Soil</td><td>SoilT<sup>b</sup></td><td>soil annual temperature (0–5 cm)</td><td>Degree Celsius</td><td>1979–2020</td><td>Lembrechts et al. (2021) (preprint)</td></tr></tbody></table>

(Dormann et al., 2013) (Supporting Information Appendix S5). To avoid overfitting (Vaughan & Ormerod, 2005), we identified the 10 most important predictors for the focal earthworm species (i.e., species recorded in ≥10 grid cells) from the 25 least correlated predictors. We used the MaxEnt algorithm (dismo package in R [Hijmans, Phillips et al., 2022]) to model each of the 41 species present in more than or exactly 10 grid cells (~2 × 2 km). The MaxEnt algorithm offered a good trade-off between computational time and performance (Elith et al., 2006; Hernandez et al., 2006; Valavi et al., 2021). We allowed models to be tuned individually and used the same pseudo-absence dataset as for model fitting (see below). To identify the top 10 variables determining earthworm distribution, permutation importance was calculated by permuting the values of each predictor and comparing the resulting reduction in training area under curve (AUC) values (Hijmans, Phillips et al., 2022). A large reduction in AUC (i.e., high permutation importance) indicates that the model is strongly influenced by that predictor. For species with n ≥ 100 occurrences (19 of the 41 focal species), we selected the 10 variables with the highest median permutation importance. Accordingly, the 10 predictor variables used in the SDMs were annual mean temperature, precipitation seasonality, distance to coast, proportion of area covered by agriculture, soil pH, phosphorus content, cation exchange capacity, elevation, clay + silt content, and human population density (Supporting Information Appendices S6 and S7a). We used the same criterion to identify the most important predictors for species with 10 ≤ n < 100 records, even though we did not predict their distributions (Supporting Information Appendices S6–S8).

## Model fitting

The biomod2 package (Thuiller et al., 2016) is one of the most commonly used and easily applicable tools to model species distributions with the ensemble approach (Hao et al., 2019).

We modeled the 19 earthworm species with ≥100 records to avoid overfitting (see Appendix S9 for sensitivity to number of records) and used the 10 algorithms available in biomod2 3.5-1 with tailored parameter settings to create ensemble models (Table 2) (Marmion et al., 2009; Valavi et al., 2021) (details in Supporting Information Appendix S10). We used committee averaging scores of the predictions, 10-fold cross-validation (80:20%), and the true skill statistic (TSS) (Allouche et al., 2006) to improve model performance during ensemble building. The committee averaging score is the average of the binary predictions of the individual models; it yields a prediction and a measure of uncertainty. If the prediction is close to 0 or 1, then all models agree to predict 0 and 1, respectively, whereas if the prediction is around 0.5, then half the models predict 1 and the other half 0 (Thuiller et al., 2016). During individual model building, we gave less weight to older species observations because they do not necessarily correspond to current species' occurrences. Unweighted and weighted models have a similarly good performance (Barbet-Massin et al., 2012). We assigned the weight of observations from 0.2 (sampling years 1970–1979) to 0.6 (2010–2019) and 0.7 (2020–2022). For model evaluation, we extracted Cohen's kappa (Allouche et al., 2006), the area under the receiver operating characteristic curve (AUC<sub>ROC</sub>), and TSS (all available as biomod2 output).

To identify the main environmental drivers of earthworm distributions, we estimated the relative variable importance. The function biomod2::get_variables_importance shuffles each environmental variable and computes the Pearson correlation coefficient of the reference predictions with the shuffled predictions (Thuiller et al., 2016). High variable importance values indicate a large influence of that variable on the model. In addition, we built simple linear regression models to estimate the influence of the 10 most important predictors on species richness and extracted *t* values from each variable. The higher the *t* value, the more important a given variable was for predicting earthworm species richness (i.e., number of species per

<footnote><sup>a</sup>Detailed information is in Appendix S1.</footnote>

<footnote><sup>b</sup>Predictors excluded based on Pearson correlation coefficient and variable inflation factor.</footnote>

<footnote><sup>c</sup>One of the 10 most important drivers of the investigated, highly recorded earthworm species.</footnote>

<footnote><sup>d</sup>Copernicus: European Union, Copernicus Land Monitoring Service 2018, European Environment Agency.</footnote>
<table><caption>**TABLE 2** Overview protocol for reporting the performed species distribution models according to Zurell et al. (2020) for European earthworm species.ª</caption><thead><tr><th>Element</th><th>Value</th></tr></thead><tbody><tr><td>Model objective</td><td>mapping and interpolation</td></tr><tr><td>Target output</td><td>maps of relative probability of presence, and of overall species richness</td></tr><tr><td>Focal taxon</td><td>earthworms (Annelida, Oligochaeta)</td></tr><tr><td>Location</td><td>Europe</td></tr><tr><td>Scale of analysis</td><td></td></tr><tr><td>Spatial extent</td><td>10°W, 31°E, 36°N, 70°N (xmin, xmax, ymin, ymax)</td></tr><tr><td>Spatial resolution</td><td>~5 × 5 km</td></tr><tr><td>Temporal extent</td><td>2015 (data sets from different time periods between 2001–2021) and 2055 (future climate 2041–2070)</td></tr><tr><td>Temporal resolution</td><td>1 year</td></tr><tr><td>Boundary</td><td>political</td></tr><tr><td>Biodiversity data</td><td></td></tr><tr><td>Observation type</td><td>citizen science, field survey, standardized monitoring data</td></tr><tr><td>Response data type</td><td>point occurrence, presence only</td></tr><tr><td>Predictors</td><td></td></tr><tr><td>Predictor types</td><td>climatic, edaphic, habitat, topographic</td></tr><tr><td>Hypotheses</td><td>(see “Introduction”)</td></tr><tr><td>Model assumptions</td><td>species-environment equilibrium or pseudo-equilibrium; no observation bias issues; independence of species observations; availability of all important predictors; predictors free of error; niche stability, constancy, and niche conservatism; no other extrapolation problems</td></tr><tr><td>Algorithms</td><td></td></tr><tr><td>Modeling techniques</td><td>MaxEnt; BioMod</td></tr><tr><td>Model averaging</td><td>For species with at least 100 records (i.e., present in at least 100 2 × 2 km grid cells), we performed BioMod ensemble modeling with all 10 available algorithms to estimate average probability of presence (true skill statistic [TSS] as evaluation criteria).</td></tr><tr><td>Model workflow</td><td>Prior to model building, we standardized all predictor variables and used variance inflation (VIF) analysis and Pearson correlation coefficient to avoid highly correlated variables. We included only the 10 most important (MaxEnt preanalysis) and weakly correlated (VIF ≤ 10, Pearson correlation ≤ 0.8) variables.</td></tr><tr><td>Model fitting</td><td>Model fitting was done automatically as much as possible. In addition, we generated individual background data for the models to increase their performances. Predictive model performance was assessed by TSS, area under the curve receiver operating characteristic, and kappa.</td></tr><tr><td>Software</td><td>R 4.2.0 (R Core Team 2016), RStudio 2022.02.3 (RStudio Team, 2020), ArcGIS 10.7.1 (Bajjali, 2018)</td></tr><tr><td>Code availability</td><td>code available from https://doi.org/10.25829/idiv.3524-gqvs4z; updates will be published on GitHub (https://github.com/JeMaNd-r/SoilBiodiversity)</td></tr><tr><td>Data availability</td><td>input data available from https://doi.org/10.25829/idiv.3524-gqvs4z.</td></tr></tbody></table>

5-km grid cell). All code used for data cleaning and analyses is available from the iDiv data repository (https://doi.org/10.25829/idiv.3524-gqvs4z). Updates will be published on GitHub (https://github.com/JeMaNd-r/SoilBiodiversity).

## Projections

For the 19 commonly observed earthworm species, we predicted and mapped the best-performing models (i.e., the one

<footnote>ªFull protocol with modeling details in Supporting Information Appendix S10.</footnote>
<figure><img src="image_7.png" /></figure>

<table><caption>TABLE 3 Climate projection scenarios according to Pörtner and Roberts (2022).</caption><thead><tr><th>Shared socioeconomic pathway</th><th>Explanation</th><th>Representative concentration pathway</th><th>Additional radiative forcing</th></tr></thead><tbody><tr><td>SSP1</td><td>sustainability</td><td>RCP 2.6</td><td>2.6 watt/m<sup>2</sup></td></tr><tr><td>SSP3</td><td>regional rivalry</td><td>RCP 7.0</td><td>7.0 watt/m<sup>2</sup> (middle to upper range of the bandwidth of all scenarios)</td></tr><tr><td>SSP5</td><td>fossil-fuel development</td><td>RCP 8.5</td><td>8.5 watt/m<sup>2</sup> (upper edge in the range of scenarios described in the literature)</td></tr></tbody></table>

with highest TSS) in R (Hijmans et al., 2022; Wickham et al., 2021). Maps are accessible through an R shiny app on GitHub (https://github.com/JeMaNd-r/Shiny-earthworm-distribution). We summed the discrete presence-absence predictions, which were derived from probability maps and a model-specific probability threshold that yielded the highest TSS, to get the number of species per grid cell (i.e., species richness). All maps were cropped to the area in which prediction uncertainty, averaged across the 19 SDMs, was less than 0.1 (=mean and median uncertainty); thereby, we excluded areas of large disagreement in predictions. Uncertainty was estimated as the coefficient of variation between the individual biomod2 algorithm projections of each ensemble model under current climatic conditions and varied among species (Supporting Information Appendices S11 & S12). The resulting investigated area spanned 4.3 million km² (i.e., 172,828 grid cells of 5 × 5 km). In addition, we used biomod2 (Thuiller et al., 2016) to estimate the area of potential extrapolation, given as the number of predictors per grid cell with values beyond the range of the training data. Extrapolation areas were mainly defined by one predictor being beyond the range and located at the coastal line and in Sweden or Finland; the latter matched well with the excluded high-uncertainty areas (Appendix S12).

We projected the future distribution and species richness of earthworms by using current land cover, topographic and soil variables, and future climate variables. We used future climate variables only because we were interested in the potential climate change effect rather than interactive effects between land cover and soil variables under future scenarios. In addition, many of the variables selected were not available or were difficult to predict (e.g., pH and other soil properties). Land-cover maps for 2050 were available, but were not yet in the CORINE Land Cover dataset that we used for accuracy reasons. We, therefore, did not predict potential future distributions but only potential shifts in species' distribution caused by climatic changes. Future predictions of mean annual temperature (T) and seasonality of precipitation (P) were downloaded from CHELSA 1.2 (Karger et al., 2017) for all available Intergovernmental Panel on Climate Change (IPCC) scenarios (n = 3*5) for the time period 2041–2070.

We considered 3 shared socioeconomic pathways: SSP1 (with the representative concentration pathway [RCP] 2.6), sustainability; SSP3 (RCP 7.0), regional rivalry; and SSP5 (RCP 8.5), fossil-fueled development (Pörtner & Roberts, 2022) (Table 3). Regarding the Earth system models, we considered 5: GFDL-

ESM4, IPSL-CM6A-LR, MPI-ESM1-2-HR, MRI-ESM2-0, and UKESM1-0-LL (Karger et al., 2017). Each SSP was combined with each Earth system model, resulting in a total of 15 IPCC scenarios. We centralized the future climate variables based on the scaling parameters of the current ones (e.g., mean and SD).

Temperature and precipitation are interconnected; but for modeling purposes, we tried to identify their individual effects on species distribution and vulnerability. We, therefore, predicted future earthworm distributions based on 3 adapted environmental data sets: future climate variables *T* and *P*, future *T* and current *P*, and current *T* and future *P*. This resulted in 45 future climate projections per species (i.e., 45 environmental datasets used for future projection of the fitted species' SDM). We performed one analysis of variance (lm function) that included the 10 environmental predictors to compare the 3 future projections of each scenario (factorial, *T*, *P*, *TP*) and to evaluate potential climate effects on the predicted species richness averaged across the 15 IPCC scenarios. To compare the effects of temperature and precipitation, we used the emmeans function from the corresponding package (Lenth et al., 2022) and calculated their estimated marginal means. Results are based on the *TP* scenario only and represent average predictions for each SSP across Earth system models unless otherwise indicated.

## Protection status

We used protected area networks available in the World Database of Protected Areas (UNEP-WCMC & IUCN, 2021) to estimate the area currently under protection. The 7 International Union for the Conservation of Nature (IUCN) categories of protected areas were Ia, strict nature reserve; Ib, wilderness area; II, national park; III, natural monument or feature; IV, habitat or species management area; V, protected landscape or seascape; and VI, protected area with sustainable use of natural resources (Dudley, 2008). We then calculated the area of the predicted binary range per species covered by protected areas differentiated into the IUCN categories (Dudley, 2008) in current and future climate predictions. We found that 15.8% of the investigated European area (683,577 km²) is currently protected under one of the 7 protection systems. Category V covered approximately 11% of the area; other categories covered less than 5% in total.
<figure><img src="image_8.png" /></figure>

## RESULTS

The 41 species occurring in more than 10 grid cells only occupied 503 grid cells on average (SD 1016.5). *Eiseniella tetraedra* (n = 4,611), *Aporrectodea caliginosa* (n = 2,849), and *A. rosea* (n = 2,205) occupied the most 2 × 2 km cells. For the 19 species with more than 100 occurrences, the median number of occupied grid cells was 776 (mean [SD] = 1,149 [1,190]) and the minimum was 109 (*Dendrobaena illyrica*). The SDMs for the 19 highly recorded species under the current climate showed generally good performance (mean Kappa = 0.53 [0.12], AUC = 0.93 [0.04], TSS = 0.70 [0.13]).

### Distribution of highly recorded European earthworms

The 19 focal earthworm species had broad potential distribution ranges across Europe (median [SD] = 2.1 million km² [880,943]) (Figure 2, Supporting Information Appendices S13 & S14). The species predicted to have the widest distribution were *Lumbricus terrestris* (3,4 million km²) (Figure 2c), *Dendrodrilus rubidus* (3,2 million km²), and *L. castaneus* (3,15 million km²). In contrast, *D. illyrica* showed the smallest potential distribution range (8,125 km²); it was restricted to an area close to the Tatra Mountains in Slovakia (Figure 2d). While most species were predicted to occur in western Central Europe, a few of the modeled species showed scattered occurrences in the Baltic states, Norway, and Portugal (e.g., *D. rubidus*, *D. octaedra*, and *E. fetida*). This was also reflected in earthworm species richness. Richness of the focal earthworm species (i.e., number of species in 5 × 5 km) was highest on the British Isles, 18 species per 25 km², and the lowest in eastern and northern Europe (Figure 2a). Nineteen species did not co-occur. On average, 3–5 species co-occurred within the 5 × 5 km grid cells (mean = 5.2 [5.54], median = 3), whereas areas with predicted high species richness covered nearly one-fourth of the European area (23.7% of the projected area with richness ≥10 and 10.1% with richness ≥ 15).

### Environmental drivers

Although our preanalysis showed that all 10 variables were highly influential for determining earthworm species occurrence, the climatic variables annual mean temperature and seasonality of precipitation were the most strongly related to the predicted distributions (absolute *t* > 240, adjusted R² = 0.73, *p* < 0.001; mean [SD] importance > 0.15 [0.11]) (Figure 2b). Other variables, such as area covered by agricultural land, clay + silt content, and phosphorus content, were similarly strongly correlated to earthworm species richness (absolute *t* > 200, *p* < 0.001). In contrast, cation exchange capacity, elevation, and soil pH had less importance for richness (absolute *t* < 70, *p* < 0.001) and single species distributions (mean importance < 0.09 [0.02–1.08]) (Appendix S6-S8). Individual species showed mixed variable importance; sometimes importance was

evenly distributed among environmental factors (e.g., *D. attensi* or *E. fetida*) (Supporting Information Appendix S8b).

When comparing species with high numbers of records (n ≥ 100) and those with only a few records (10 ≤n <100), the same factors strongly determined distribution of the less widely recorded species (i.e., same top 10 variables) (Supporting Information Appendix S7). For 7 of these 22 species (e.g. *Satchellius madeirensis* and *A. molleri*), climatic factors made up more than 50% of the variable importance. Sensitivity analysis of future climate scenarios revealed a significant effect of climate scenario (*P*, *T*, *TP*) on the predicted future earthworm species richness (adjusted R² = 0.711, *p* < 0.001, *F* = 1548). Temperature showed, on average, slightly higher species richness values than precipitation seasonality (estimated marginal means [SE]: *P* = 5.2, *T* = 5.5, *TP* = 5.6 [0.005]).

### Future predictions of highly recorded European earthworms

Under the future SSPs, species ranges appeared to shift from western Europe toward the east (i.e., losses in the west, gains in the east) (Figure 3, Supporting Information Appendix S14), which was also reflected in changes in local species richness. This result was confirmed across all future scenarios tested (Figure 3a), even though changes in richness appeared slightly less severe under more sustainable SSPs (i.e., SSP1 and SSP3) (Figure 3b). In large parts of the investigated areas, all 3 SSPs predicted similar changes in species richness: gain or loss of species was predicted across all 3 SSPs in 45.7 and 24.1% of the investigated areas, respectively. In only 7.1% of the area covered in this study (308,675 km²), scenarios showed mixed predictions (i.e., some showed losses and others gains). At the western coast of southern France and south of the Alps, SSP5 showed larger decreases in species richness, while fewer parts of Europe remained unchanged (SSP1 20.7% of whole investigated area remain unchanged, SSP3 11.7%, SSP5 8.2%). Especially along the coast of the Baltic Sea, SSP3 and SSP5 showed positive changes in species richness (SSP1 increase in richness in 52.3% of the area, SSP3 57.3%, SSP5 56.8%). All potentially occurring species (mean = 2.27, maximum = 9 species) were predicted to be lost under future conditions in 1.8% of the study area (78,200 km²). These scattered areas were mainly in southeastern Europe, the Pyrenees, north of Spain, and Italy. Overall, the potential earthworm species ranges declined in half of the species by up to 54%; 9 species showed comparably strong increases in range size (absolute mean [SD] = 18.9% [0.187]) (Figure 3c, Supporting Information Appendix S13). Proportional range change varied across species and tended to be more pronounced in less widely distributed species (Supporting Information Appendix S13). Accordingly, *A. caliginosa*, *A. longa*, and *D. attensi* showed the strongest decreases in their distribution ranges and little range expansion (mean > 35%) (Figure 3c; Supporting Information Appendices S13 & S14). In contrast, ranges of broadly distributed species, such as *D. rubidus*, *L. castaneus*, *L.*
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>(b) Variable importance</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) *Lumbricus terrestris*</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(d) *Dendrobaena illyrica*</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(e) Uncertainty</figcaption></figure>

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

**FIGURE 2** (a) Spatial distribution of the predicted overall European earthworm species richness (i.e., number of species, maximum 19) when the probability of species' presence is higher than the species-specific threshold with maximum true skill statistic value (dark gray, 0 species; light gray, no predictions made because of lack of environmental data and model uncertainty), (b) importance of the 10 variables used to predict earthworm presence for each species distribution model (points, average importance for specific species across 10 replicates [values outside 1.5 times the interquartile range above and below the upper and lower quartile, respectively]; variables ordered by mean importance values across all species and replicates; horizontal lines, range of values excluding potential outliers), (c) spatial distribution of the species with the largest (i.e., most predicted occupied grid cells) geographic range, (d) spatial distribution of the species with the smallest geographic range (i.e., fewest grid cells), and (e) uncertainty in predicted earthworm distributions (coefficient of variation averaged across species-specific biomod2 models).

*rubellus*, and *L. terrestris*, were predicted to change only slightly (mean < 0.15).

## Species vulnerability under nature conservation categories

The average species richness was the highest in protected areas belonging to IUCN category IV and V, both of which represent less strict protection levels (mean [SD] = 9.4 and 9.6 [5.90 and 5.16], respectively), and the lowest in Ia and Ib areas

(mean = 2.72 and 2.43 [3.5]) (Figure 4a). Unprotected areas showed on average 6-7 earthworm species per 25 km² (mean = 6.7 [5.96]) and covered 81.6% of the species range (0.04). The predicted individual species ranges were covered between 14.8% (*E. fetida*: 306,180 km²) and 31% (*D. illyrica* 2,520 km²) under one of the 7 protection systems (mean = 18.4% [0.04] (Figure 4c; Supporting Information Appendix S13). Similar to species richness, most of the protected ranges fell into IUCN category IV or V (mean > 50,000 km², 2.7 and 14.4% [0.01 and 0.04], respectively) (Figure 4c), whereas categories Ia and Ib covered only a very small percentage of the predicted species
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**FIGURE 3** (a) Agreement across the 3 shared socioeconomic pathway (SSP) scenarios of climate change (predictions, mean predictions across 5 Earth system models; gain, areas where gain of species richness is predicted in 3, 2, and 1 of the scenarios; loss, areas where loss of species richness is predicted; mixed, areas where different scenarios predict gain and loss), (b) change in number of species gained (blue, > 0) and lost (red, < 0) for each of the 3 SSPs (SSP1, most sustainable scenario; SSP5, fossil-fueled development; dark gray, areas excluded based on uncertainty > 0.1), and (c) change in range size per species ordered according to size of current distribution range (bars, proportion of area changes; error bars, SD across future scenarios; points, range sizes under current and future climate conditions; genus abbreviations: Ap., *Aporrectodea*; Al., *Allolobophora*; Ad., *Allolobophoridella*; D., *Dendrobaena*; Dd., *Dendrodrilus*; E., *Eisenia*; El., *Eiseniella*; L., *Lumbricus*; O., *Octolasion*; S. Satchellius).</figcaption></figure>

ranges (0.1–0.3% and 0–0.8%, mean < 0.1% [0.001 and 0.002], respectively).

Under future climate conditions, the average species richness was predicted to slightly increase in all protected ecosystems, except in VI areas (mean = −0.3 [2.4]) (Figure 4b). Species richness increases were largest in Ia, Ib, and IV areas (mean = −1.4, 0.7, and 0.8 [1.9, 1.1, and 2.1], respectively). The predicted range area under protection changed by an average of 11% (0.23) across IUCN categories and species (minimum = −100%, maximum = 16.5%) (Figure 4d, Supporting Information Appendix S13). Distribution ranges within IUCN category III and VI were predicted to have the greatest declines (mean = 15% [0.25 and 0.22], respectively), whereas range areas covered by category Ib and IV were predicted to increase (mean = 1.1 and 94.2% [0.4 and 3.8], respectively) (Figure 4d). In unprotected areas, species ranges decreased by 12% (0.24) on average. We found mostly agreement across SSPs with very similar patterns for the predicted changes in species richness (Figure 4c) and species ranges
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>**FIGURE 4** Status of protection of the 19 most-widespread European earthworm species under (a, b) current and (c, d) future climate conditions. (a) Average earthworm species richness in each protected area network in Europe (IUCN categories) under current climate conditions (white, unprotected areas; black dots, mean richness, error bars are standard deviations), (b) area of earthworm species ranges covered by protected area networks (unprotected area not shown; colors correspond to those in [a]), (c) change in species richness under the 3 shared socioeconomic pathway (SSP) scenarios (predictions, mean predictions across 5 Earth system models; columns [Ia-VI], protection categories; blue, average gain of species [i.e., number of species > 0]; red, average loss of species), and (d) change in protected range area under 3 SSP scenarios (+5%, future increase of 5% of the species' range size currently covered by a protected area network; dark blue, change in coverage > 100% depicted as 110%; gray, species that did not occur in certain protected area categories; bars within each cell, SSPs; genus abbreviations in legend of Figure 3).</figcaption></figure>

across protected and unprotected areas (Figure 4d). Changes in protected areas under the more sustainable SSP1 were less pronounced than under SSP3 and SSP5 (absolute mean = 20.8% vs. 33.6 and 37.3% [0.6, 1.4 and 1.5], respectively). For example, *D. illyrica* showed the more dramatic changes under SSP3 and SSP5 (mean = -99% [< 0.01]) than under SSP1 (-88% [0.06]).

DISCUSSION

In Europe, earthworm occurrence is an important indicator of ecosystems in good condition (Bardgett & van der Putten, 2014; Blouin et al., 2013; Cunha et al., 2016; Kooch & Jalilvand, 2008), but information on their conservation status is scarce (Cameron
et al., 2018; Ciobanu et al., 2019; Guerra et al., 2020; Zeiss et al., 2022). We assessed drivers of 41 European earthworm species and predicted the current and future distributions of 19 highly recorded species. Soil pH, phosphorus, and temperature were identified previously as strong drivers of earthworm distribution (Köninger et al., 2023; Phillips et al., 2019; Ruiz et al., 2021); however, we found precipitation seasonality, agricultural land area, and soil texture similarly or even more important for predicting patterns of earthworm distribution and richness (Figure 2b, Supporting Information Appendices S6–S8). The majority of the species occurred in the west of central Europe (i.e., the British Isles, Germany, the Netherlands, and Belgium) (Figure 2a).

While the maximum number of earthworm species per site in previous studies was 6 (Fourcade & Vercauteren, 2022; Phillips et al., 2019), 7 (Rutgers et al., 2016), or 8 (Decaëns, 2010), our models predicted up to 18 species that can co-occur in several European regions (140,925 km² with 18 species). Previous richness values were extrapolated based on local species richness, rather than by summing single species distributions, meaning that higher species numbers will be predicted for larger spatial scales where heterogeneity within regions or landscapes is known to promote earthworm diversity (Decaëns, 2010; Le Provost et al., 2021).

In addition, we used a more comprehensive dataset compared with previous mapping attempts and were able to fill data gaps in northern Europe, Portugal, Ukraine, Belgium, Great Britain, and Germany (Phillips et al., 2019; Rutgers et al., 2016). The areas showing high earthworm diversity also matched the areas identified to have the highest soil biodiversity in Europe (Aksoy et al., 2017). Thus, higher soil biodiversity was predicted in Ireland and southern England, south of the Alps, and along the coast of northern Germany and the Netherlands, whereas large parts of Italy and Spain showed lower potential (Aksoy et al., 2017). This pattern is in line with our findings (Figure 2a); however, the absolute species number has to be interpreted with caution due to the exclusion of earthworm species with low occurrence records and the selected spatial scale (i.e., regional at 5 × 5 km).

Our models predicted potential declines in species richness under future climatic conditions especially in western Europe and shifts in distribution ranges to the central east (Figure 3). While we did not predict potential future distributions, we estimated potential shifts in species' distribution caused by climatic changes. These shifts were mainly explained by changes in the mean annual temperature (estimated marginal means [SE] = 5.5 [0.005], variable importance *t* = 248.3). Changes under more sustainable SSPs were less dramatic than under the scenario with the highest radiative forcing (i.e., SSP5). However, areas with low variations in species richness can also be explained by shifts in species distributions leading to substitutions (i.e., one species moves into regions where another species becomes absent) (Supporting Information Appendix S14).

In line with previous findings on earthworms (Fourcade & Vercauteren, 2022; Singh et al., 2019) and vertebrates, such as mammals (Levinsky et al., 2007) and birds (Virkkala & Lehikoinen, 2017), earthworm distribution patterns are likely

to change within the next decades, leading to unpredictable consequences for ecosystem functioning. Due to their role in multiple ecosystem services (Blouin et al., 2013), earthworm range shifts could also result in alterations in the distribution of mutually related taxa, such as plant species hosting specific soil taxa (Eisenhauer et al., 2019). Although newly developed habitats could be recolonized by resilient earthworm populations, their low dispersal capabilities (Chatelain & Mathieu, 2017) and their strong dependence on migration corridors (Dupont et al., 2015, 2017) strongly limit their successful establishment into future habitats. It is, therefore, important to expand and apply knowledge on practices to mitigate climate change effects on earthworms and ecosystems in general, for example with continuous plant cover and subordinate plants that maintain high food quality for earthworms (Mariotte et al., 2016; Singh et al., 2019).

Most focal earthworm species were predicted to be widely spread across Europe (Figure 2, Appendix S14); however, we identified several endemically distributed species (Supporting Information Appendix S14), for example, *D. illyrica* was restricted to the Tatra Mountains in Slovakia. We also found *S. mammalis* was restricted to southern areas of the British Isles and to isolated areas in Romania, Bulgaria, Italy, France, and Spain (Supporting Information Appendix S14). *L. festivus* followed similar patterns as *S. mammalis* but extended toward the northeast. However, neither species showed occurrence data in Romania or Bulgaria. These distribution patterns indicate that suitable areas for *S. mammalis* and *L. festivus* in southeastern Europe could not be colonized due to geographical barriers, although more research is needed to confirm this due to scarcity of data in these regions. Our results suggest that *S. mammalis* and *L. festivus* in Romania, Bulgaria, Italy, and Spain are likely to disappear under future climatic conditions. In addition, the distribution of one-third of the less widely recorded species ($10 \le n < 100$) was mainly explained by climatic factors (i.e., summed permutation importance >50%) (Supporting Information Appendix S8), highlighting the vulnerability of earthworms to climate change (Singh et al., 2019). This finding is consistent with endemic species facing higher extinction risks because of their vulnerability to habitat changes (Urban, 2015). The existence of vulnerable earthworm species makes their protection especially important to avoid species losses and maintain regional soil biodiversity.

Our results show that the protection of earthworm species under current protected areas is clearly insufficient. Predicted species richness was lower for the more strictly protected areas (i.e., IUCN categories Ia and Ib) compared with unprotected sites; however, the highest richness was predicted in areas belonging to categories IV and V. This confirms previous findings on negligible effects of conservation areas on soil organisms (Ciobanu et al., 2019; Zeiss et al., 2022). Under future global change scenarios, regional species richness will increase in most of the protected areas, with larger increases in category Ia and Ib areas (Figure 4b). Looking at species level, earthworm species ranges were most likely to be protected by < 20% in at least one of the 7 protection categories (Ia–VI). The majority of these protected range areas belonged to conservation
<figure><img src="image_19.png" /></figure>

networks with less strict protection (i.e., IUCN categories IV and V).

We recognize that our distribution predictions represent earthworm fundamental niches rather than the actual geographic area they occupy (i.e., realized niche [Soberon & Peterson, 2005]). In addition, predicted species ranges covered by protected area networks do not reflect the real species' protection status. However, knowledge about species ranges is essential to evaluate species vulnerability under the IUCN framework. Even if models show lower accuracy in predicting the suitability of individual grid cells and fundamental niches, they are very useful surrogates for species' range size (Júnior & Nóbrega, 2018) and can be used to design protection measures. In this context, we found that distribution ranges declined for nearly half of the 19 most highly recorded earthworm species under future climate conditions with certain range shifts into protected areas (Figure 4d). Most protected areas showed increases in species richness and occurrence, indicating that current conservation measures could be sufficient in the near future as long as they specifically target earthworms (Zeiss et al., 2022; Guerra et al., 2022).

The European Green Deal set the ambitious goal of protecting at least 30% of the land and sea area by 2050. Such quantitative targets are only effective for nature conservation if cryptic, less charismatic, but functionally significant species—including those residing in the soil—are considered when designing and managing protected areas. Protected area networks are the main tool in conservation to maintain above-ground species, habitats, and ecosystems worldwide (Coetzee et al., 2014), but they do not guarantee real protection of all species and habitats (Gaston et al., 2006; Leverington et al., 2010; Barnes et al., 2018). To reach conservation goals, their design and management have to meet individual species and ecosystem needs and to consider future impacts of global change (Gaston et al., 2006; Leverington et al., 2010). Earthworms and many other soil-dwelling organisms are currently underprotected and are not explicitly targeted in nature conservation policies and regulations (Köninger et al., 2022; Zeiss et al., 2022). One solution could be to expand species-specific targets in protected areas (Loiseau et al. 2020) that can help guarantee their protection and protection of their habitats, as well as other soil taxa. Soil biodiversity data has become increasingly available over the past decades, and existing datasets already provide valuable information (cf. Cardoso et al., 2011) that can be used to identify drivers of species distribution and climate change responses so that more effective actions for soil invertebrate conservation can be taken.

Still, our aim was not to provide a comprehensive and complete spatial representation of European earthworm species richness, but to use readily available data to identify areas that are potentially at risk of losing earthworm species and in urgent need of protection and to highlight research needs. The limited number of occurrence records and the lack of true absence data for most species (i.e., repeated observations of species not being present in a location) adds uncertainty to our predictions of species occurrences and species richness distribution (Figure 2a, Supporting Information Appendices S12 & S13) (Hernandez

et al., 2006). Furthermore, we did not test for spatial autocorrelation because this may remove environmental variables that correlate with (or drive) species distributions. The effects of spatial autocorrelation may decrease transferability of the predictions into new areas and environments (Segurado et al., 2006). We, therefore, excluded areas with high uncertainty values that likely represent environmental conditions not covered by species data sampling locations.

Although the identified drivers of and predicted distributions for the targeted species aligned with the relevant literature (Rangers et al., 2016; Phillips et al., 2019; Köninger et al., 2023), evaluating the resulting potential distributions is more challenging and more prone to bias due to major geographic and taxonomic gaps (Guerra et al., 2020). More targeted and standardized sampling is needed especially in central and southern of Spain, southern of Portugal, Italy, and eastern European countries (Supporting Information Appendices S12 & S13) (Guerra et al., 2020). For example, the potential distribution of *D. attensi* was likely underestimated, possibly due to confusions with *D. octaedra*, which can be overcome by implementing harmonized sampling protocols and expert training in taxonomy (FAO et al., 2020; Guerra et al., 2021; Mathieu et al., 2022a). This possibly applies to other species across several areas in Europe. Moreover, our 19 focal species did not include several endemic species because of the low number of observations. They require a different, tailored modeling approach that allows overcoming such data limitations and an assessment at smaller spatial scales (Lomba et al., 2010; Breiner et al., 2015). Complementary, large-scale soil monitoring frameworks such as SoilBON (Guerra et al., 2021) or targeted data gathering, such as for EUdaphobase COST Action, are needed to verify and expand soil biological data, draw more robust conclusions, and take actions to limit or enhance dispersal mechanisms of earthworms among others.

We found that the most commonly observed earthworm species are broadly distributed throughout western and central Europe but there are still data gaps in northern, southern, and eastern Europe. At the regional scale, we found higher values of potential species richness than previously predicted at the local scale, with a maximum of 18 species co-occurring in 5-km grid cells. Under future climate conditions, potential species richness and individual species geographical distribution showed a shift from the west to the central east of Europe. Predicted isolated range areas of species with the most restricted distributions tended to be more vulnerable to climate change because their potential ranges are likely to shrink under future climate conditions. The existence of these vulnerable earthworm species underlines the urgency of protecting not only current above-ground species ranges, but also belowground organisms and their potential future distribution ranges, including those areas that could provide stable refugia.

## ACKNOWLEDGEMENTS

We acknowledge valuable support from K. Hohberg and D. Russell from Senckenberg Museum of Natural History Görlitz and the preparation of Figure 1 (workflow) by G. Rada from iDiv. We also acknowledge support of iDiv funded by
<figure><img src="image_20.png" /></figure>

the German Research Foundation (DFG FZT 118, 202548816), and L. Figueiredo from the iDiv Data & Code Unit for assistance with curation and archiving of the dataset. R.Z. was funded by the German Federal Environmental Foundation (DBU, 20021/752). J.M. was funded by sDiv (Grant SFW9.02) and France-Stanford Center for Interdisciplinary Studies. A.L. is supported by national funds through FCT—Fundação para a Ciência e a Tecnologia, I.P., in the context of the Transitory Norm: DL57/2016/CP1440/CT0001.

Open access funding enabled and organized by Projekt DEAL.

## ORCID

Romy Zeiss ID https://orcid.org/0000-0001-8862-9185

Maria J. I. Briones ID https://orcid.org/0000-0002-4051-8117

Jérôme Mathieu ID https://orcid.org/0000-0002-9106-6106

Angela Lomba ID https://orcid.org/0000-0003-3748-8751

Jessica Dahlke ID https://orcid.org/0000-0002-3246-3128

Laura-Fiona Heptner ID https://orcid.org/0000-0002-2081-2300

Gabriel Salako ID https://orcid.org/0000-0001-7960-8200

Nico Eisenhauer ID https://orcid.org/0000-0002-0371-6720

Carlos A. Guerra ID https://orcid.org/0000-0003-4917-2105

## REFERENCES

Aksoy, E., Louwagie, G., Gardi, C., Gregor, M., Schröder, C., & Löhnerzt, M. (2017). Assessing soil biodiversity potentials in Europe. *Science of the Total Environment, 589*, 236–249.

Allouche, O., Tsoar, A., & Kadmon, R. (2006). Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). *Journal of Applied Ecology, 43*, 1223–1232.

Aune-Lundberg, L., & Strand, G. H. (2021). The content and accuracy of the CORINE land cover dataset for Norway. *International Journal of Applied Earth Observation and Geoinformation, 96*, 102266.

Bajjali, W. (2018). *ArcGIS for Environmental and Water Issues*. Springer International Publishing.

Ballabio, C., Jiskra, M., Osterwalder, S., Borrelli, P., Montanarella, L., & Panagos, P. (2021). A spatial assessment of mercury content in the European Union topsoil. *Science of the Total Environment, 769*, 144755.

Ballabio, C., Lugato, E., Fernández-Ugalde, O., Orgiazzi, A., Jones, A., Borrelli, P., & Montanarella, L., Panagos, P. (2019). Mapping LUCAS topsoil chemical properties at European scale using Gaussian process regression. *Geoderma, 355*, 113912.

Ballabio, C., Panagos, P., Lugato, E., Huang, J. H., Orgiazzi A., Jones, A., Fernández-Ugalde, O., Borrelli, P., & Montanarella, L. (2018). Copper distribution in European topsoils: An assessment based on LUCAS soil survey. *Science of the Total Environment, 636*, 282–298.

Ballabio, C., Panagos, P., & Monatanarella, L. (2016). Mapping topsoil physical properties at European scale using the LUCAS database. *Geoderma, 261*, 110–123.

Barbet-Massin, M., Jiguet, F., Albert, C. H., & Thuiller, W. (2012). Selecting pseudo-absences for species distribution models: How, where and how many? *Methods in Ecology and Evolution, 3*, 327–338.

Bardgett, R. D., & van der Putten, W. H. (2014). Belowground biodiversity and ecosystem functioning. *Nature, 515*, 505–511.

Barnes, M. D., Glew, L., Wyborn, C., & Craigie, I. D. (2018). Prevent perverse outcomes from global protected area policy. *Nature Ecology & Evolution, 2*, 759–762.

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.

Breiner, F. T., Guisan, A., Bergamini, A., & Nobis, M. P. (2015). Overcoming limitations of modelling rare species by using ensembles of small models. *Methods in Ecology and Evolution, 6*, 1210–1218.

Buchhorn, M., Smets, B., Bertels, L., Roo, B. D., Lesiv, M., Tsendbazar, N.-E., Herold, M., & Fritz, S. (2020). Copernicus global land service: Land cover 100 m: collection 3: epoch 2018: Globe. Zenodo.

Caetano, M., Mata, F., Freire, S., & Campagnolo, M. (2006). Accuracy assessment of the Portuguese CORINE Land Cover map. *Global Developments in Environmental Earth Observation from Space, 459*, 467.

Cameron, E. K. (2018). Global gaps in soil biodiversity data. *Nature Ecology & Evolution, 2*, 1042–1043.

Cameron, E. K., Martins, I. S., Lavelle, P., Mathieu, J., Tedersoo, L., Bahram, M., Gottschall, F., Guerra, C. A., Hines, J., Patoine, G., Siebert, J., Winter, M., Cesarz, S., Ferlian, O., Kreft, H., Lovejoy, T. E., Montanarella, L., Orgiazzi, A., Pereira, H. M., ... Eisenhauer, N. (2019). Global mismatches in aboveground and belowground biodiversity. *Conservation Biology, 33*, 1187–1192.

Cammalleri, C., Vogt, J. V., Bisselink, B., & de Roo, A. (2017). Comparing soil moisture anomalies from multiple independent sources over different regions across the globe. *Hydrology and Earth System Sciences, 21*, 6329–6343.

Cardoso, P., Erwin, T. L., Borges, P. A. V., & New, T. R. (2011). The seven impediments in invertebrate conservation and how to overcome them. *Biological Conservation, 144*, 2647–2655.

Chatelain, M., & Mathieu, J. (2017). How good are epigeic earthworms at dispersing? An investigation to compare epigeic to endogeic and anecic groups. *Soil Biology and Biochemistry, 111*, 115–123.

Ciobanu, M., Eisenhauer, N., Stoica, I.-A., & Cesarz, S. (2019). Natura 2000 priority and non-priority habitats do not differ in soil nematode diversity. *Applied Soil Ecology, 135*, 166–173.

Coetzee, B., Gaston, K., & Chown, S. (2014). Local scale comparisons of biodiversity as a test for global protected area ecological performance: A meta-analysis. *PLOS ONE, 9*, e105824.

Cunha, L., Brown, G. G., Stanton, D. W. G., Da Silva, E., Hansel, F. A., Jorge, G., Mckey, D., Vidal-Torrado, P., Macedo, R. S., Velasquez, E., James, S. W., Lavelle, P., & Kille, P. (2016). Soil animals and pedogenesis: The role of earthworms in anthropogenic soils. *Soil Science, 181*, 110–125.

Decaëns, T. (2010). Macroecological patterns in soil communities. *Global Ecology and Biogeography, 19*, 287–302.

Decaëns, T., Jiménez, J. J., Gioia, C., Measey, G. J., & Lavelle, P. (2006). The values of soil animals for conservation biology. *European Journal of Soil Biology, 42*, S23–S38.

Delgado-Baquerizo, M., Oliverio, A. M., Brewer, T. E., Benavent-González, A., Eldridge, D. J., Bardgett, R. D., Maestre, F. T., Singh, B. K., & Fierer, N. (2018). A global atlas of the dominant bacteria found in soil. *Science, 359*, 320–325.

Desie, E., Van Meerbeek, K., De Wandeler, H., Bruelheide, H., Domisch, T., Jaroszewicz, B., Joly, F.-X., Vancampenhout, K., Vesterdal, L., & Muys, B. (2020). Positive feedback loop between earthworms, humus form and soil pH reinforces earthworm abundance in European forests. *Functional Ecology, 34*, 2598–2610.

Di Cecco, G. J., Barve, V., Belitz, M. W., Stucky, B. J., Guralnick, R. P., & Hurlbert, A. H. (2021). Observing the observers: How participants contribute data to inaturalist and implications for biodiversity science. *BioScience, 71*, 1179–1188.

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

Doxsey-Whitfield, E., MacManus, K., Adamo, S. B., Pistolesi, L., Squires, J., Borkovska, O., & Baptista, S. R. (2015). Taking advantage of the improved availability of census data: a first look at the gridded population of the world, version 4. *Papers in Applied Geography, 1*, 226–234.

Dudley, N. (2008). *Guidelines for Applying Protected Area Management Categories*. IUCN.
<figure><img src="image_21.png" /></figure>

Dupont, L., Grésille, Y., Richard, B., Decaëns, T., & Mathieu, J. (2015). Dispersal constraints and fine-scale spatial genetic structure in two earthworm species. *Biological Journal of the Linnean Society, 114*, 335–347.

Dupont, L., Torres-Leguizamon, M., René-Corail, P., & Mathieu, J. (2017). Landscape features impact connectivity between soil populations: A comparative study of gene flow in earthworms. *Molecular Ecology, 26*, 3128–3140.

Edaphobase.org. (2021). Edaphobase data warehouse on soil organisms.

Eisenhauer, N., Buscot, F., Heintz-Buschart, A., Jurburg, S. D., Küsel, K., Sikorski, J., Vogel, H.-J., & Guerra, C. A. (2020). The multidimensionality of soil macroecology. *Global Ecology and Biogeography, 30*, 4–10.

Eisenhauer, N., Ferlian, O., Craven, D., Hines, J., & Jochum, M. (2019). Ecosystem responses to exotic earthworm invasion in northern North American forests. *Research ideas and outcomes, 5*, e34564.

Elith, J., Graham, C. H., Anderson, R. P., Dudik, M., Ferrier, S., Guisan, A., Hijmans, R. J., Huettmann, F., Leathwick, J. R., Lehmann, A., Li, J., Lohmann, L. G., Loiselle, B. A., Manion, G., Moritz, C., Nakamura, M., Nakazawa, Y., Overton, J. M., Peterson, A. T., ... Zimmermann, N. E. (2006). Novel methods improve prediction of species' distributions from occurrence data. *Ecography, 29*(2), 129–151.

Ellis, E. C., Goldewijk, K. K., & Siebert, S. (2013). *Anthropogenic Biomes of the World, Version 2: 2000*. NASA Socioeconomic Data and Applications Center (SEDAC).

ESDAC ESDC. (2020). LUCAS: Land Use and Coverage Area frame Survey. Available from https://esdac.jrc.ec.europa.eu/projects/lucas

European Environment Agency, & European Environment Agency. (2019). CORINE Land Cover 2012 (raster 100 m), Europe, 6-yearly - version 2020_20u1, May 2020 (Version 20.01) [Data set]. European Environment Agency. https://doi.org/10.2909/A84AE124-C5C5-4577-8E10-511BFE55CC0D

FAO, ITPS, GSBI, SCBD, EC. (2020). *State of knowledge of soil biodiversity-Status, challenges and potentialities: Report 2020*. Report. Food and Agriculture Organisation.

Fourcade, Y., & Vercauteren, M. (2022). Predicted changes in the functional structure of earthworm assemblages in France driven by climate change. *Diversity and Distributions, 28*, 1050–1066.

Gaston, K. J., Charman, K., Jackson, S. F., Armsworth, P. R., Bonn, A., Briers, R. A., Callaghan, C. S. Q., Catchpole, R., Hopkins, J., Kunin, W. E., Latham, J., Opdam, P., Stoneman, R., Stroud, D. A., & Tratt, R. (2006). The ecological effectiveness of protected areas: The United Kingdom. *Biological Conservation, 132*, 76–87.

GBIF.org. (2021). Global Biodiversity Information Facility (GBIF) Home Page.

GBIF.org. (2022). GBIF occurrence download “Crassicillitellata”. *Global Biodiversity Information Facility*, https://doi.org/10.15468/dl.ghtngj (accessed March 2022)

Guerra, C. A., Heintz-Buschart, A., Sikorski, J., Chatzinotas, A., Guerrero-Ramírez, N., Cesarz, S., Beaumelle, L., Rillig, M. C., Maestre, F. T., Delgado-Baquerizo, M., Buscot, F., Overmann, J., Patoine, G., Phillips, H. R. P., Winter, M., Wubet, T., Küsel, K., Bardgett, R. D., Cameron, E. K., ... Eisenhauer, N. (2020). Blind spots in global soil biodiversity and ecosystem function research. *Nature Communications, 11*, 3870.

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, R., Briones, M. J. I., Buscot, F., Cameron, E. K., Cesarz, S., Chatzinotas, A., Cowan, D. A., Djukic, I., Van Den Hoogen, J., ... Eisenhauer, N. (2021). Tracking, targeting, and conserving soil biodiversity. *Science, 371*, 239–241. American Association for the Advancement of Science.

Guerra, C. A. (2022). Global hotspots for soil nature conservation. *Nature, 610*, 693–698.

Hao, T., Elith, J., Guillera-Arroita, G., & Lahoz-Monfort, J. J. (2019). A review of evidence about use and performance of species distribution modelling ensembles like BIOMOD. *Diversity and Distributions, 25*, 839–852.

Hernandez, P. A., Graham, C. H., Master, L. L., & Albert, D. L. (2006). The effect of sample size and species characteristics on performance of different species distribution modeling methods. *Ecography, 29*, 773–785.

Hijmans, R. J. et al. (2022a). raster: *Geographic Data Analysis and Modeling*.

Hijmans, R. J., Bivand, R., Forner, K., Ooms, J., Pebesma, E., & Sumner, M. D. (2022b). *terra: Spatial Data Analysis*.

Hijmans, R. J., Phillips, S., & Elith, J. L. (2022c). dismo: Species Distribution Modeling.

Jetz, W. (2019). Essential biodiversity variables for mapping and monitoring species populations. *Nature Ecology & Evolution, 3*, 539–551.

Júnior, P. D. M., & Nóbrega, C. C. (2018). Evaluating collinearity effects on species distribution models: An approach based on virtual species simulation. *PLOS ONE, 13*, e0202403.

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

Köninger, J., Ballabio, C., Panagos, P., Jones, A., Schmid, M. W., Orgiazzi, A., & Briones, M. J. I. (2023). Ecosystem type drives soil eukaryotic diversity and composition in Europe. *Global Change Biology, 00*, 1–14.

Köninger, J., Panagos, P., Jones, A., Briones, M. J. I., & Orgiazzi, A. (2022). In defence of soil biodiversity: Towards an inclusive protection in the European Union. *Biological Conservation, 268*, 109475.

Kooch, Y., & Jalilvand, H. (2008). Earthworms as ecosystem engineers and the most important detritivors in forest soils. *Pakistan Journal of Biological Sciences, 11*, 819–825.

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

Lavelle, P. (2011). Earthworms as ecosystem engineers. In: *Encyclopedia of agrophysics*, eds J. Glinski, J. Horabik, and J. Lipiec. 233–235. New York NY: Springer.

Le Bayon, R. C., Bullinger-Weber, G., Schomburg, A., Turberg, P., Schlaepfer, R., & Guenat, C. (2017). Earthworms as ecosystem engineers: A review. In: *Earthworms: Types, Roles and Research*. 129–178. New York NOVA Science Publishers.

Le Provost, G., Thiele, J., Westphal, C., Penone, C., Allan, E., Neyret, M., Van Der Plas, F., Ayasse, M., Bardgett, R. D., Birkhofer, K., Boch, S., Bonkowski, M., Buscot, F., Feldhaar, H., Gaulton, R., Goldmann, K., Gossner, M. M., Klaus, V. H., Kleinebecker, T., ... Manning, P. (2021). Contrasting responses of above- and belowground diversity to multiple components of land-use intensity. *Nature Communications, 12*, 3918.

Lehmitz, R., Römbke, J., Graefe, U., Beylich, A., & Krück, S. (2016). Rote Liste und Gesamtartenliste der Regenwürmer (Lumbricidae et Criodrilidae) Deutschlands. In: *Page Rote Liste der gefährdeten Tiere, Pflanzen und Pilze Deutschlands. Band 4: Wirbellose Tiere (Teil 2)*. Bonn Bundesamt für Naturschutz.

Lembrechts, J., van den Hoogen, J., Nijs, I., & Lenoir, J. (2021). Global maps of soil temperature. Preprint. EcoEvoRxiv: https://ecoevorxiv.org/repository/view/4122/

Lenth, R. V., Buerkner, P., Herve, M., Jung, M., Love, J., Miguez, F., Riebl, H., & Singmann, H. (2022). *emmeans*: Estimated Marginal Means, aka Least-Squares Means.

Leverington, F., Costa, K. L., Pavese, H., Lisle, A., & Hockings, M. (2010). A global analysis of protected area management effectiveness. *Environmental Management, 46*, 685–698.

Levinsky, I., Skov, F., Svenning, J.-C., & Rahbek, C. (2007). Potential impacts of climate change on the distributions and diversity patterns of European mammals. *Biodiversity and Conservation, 16*, 3803–3816.

Loiseau, N., Mouquet, N., Casajus, N., Grenié, M., Guéguen, M., Maitner, B., Mouillot, D., Ostling, A., Renaud, J., Tucker, C., Velez, L., Thuiller, W., & Violle, C. (2020). Global distribution and conservation status of ecologically rare mammal and bird species. *Nature Communications, 11*(1), 5071.

Lomba, A., Pellissier, L., Randin, C., Vicente, J., Moreira, F., Honrado, J., & Guisan, A. (2010). Overcoming the rare species modelling paradox: A novel hierarchical framework applied to an Iberian endemic plant. *Biological Conservation, 143*, 2647–2657.
<figure><img src="image_22.png" /></figure>

Mariotte, P., Le Bayon, R., Eisenhauer, N., Guenat, C., & Buttler, A. (2016). Subordinate plant species moderate drought effects on earthworm communities in grasslands. *Soil Biology and Biochemistry, 96*, 119-127.

Marmion, M., Parviainen, M., Luoto, M., Heikkinen, R. K., & Thuiller, W. (2009). Evaluation of consensus methods in predictive species distribution modelling. *Diversity and Distributions, 15*, 59-69.

Mathieu, J. (2022a). sOilFauna: A global synthesis effort on the drivers of soil macrofauna communities and functioning: Workshop report. *Soil Organisms, 94*, 111-126.

Mathieu, J., & Davies, T J. (2014). Glaciation as an historical filter of below-ground biodiversity. *Journal of Biogeography, 41*, 1204-1214.

Mathieu, J., Reynolds, J. W., Fragoso, C., & Hadly, E. (2022). Global worming: Massive invasion of North America by earthworms revealed. *bioRxiv, Available from https://www.biorxiv.org/content/10.1101/2022.06.27.497722v1 (accessed October 2022)

McMullin, R. T., & Allen, J L. (2022). An assessment of data accuracy and best practice recommendations for observations of lichens and other taxonomically difficult taxa on iNaturalist. *Botany, 100*, 491-497.

Mueller, K E., Eisenhauer, N., Reich, P. B., Hobbie, S. E., Chadwick, O. A., Chorover, J., Dobies, T., Hale, C. M., Jagodziński, A. M., Kałucka, I., Kasprowicz, M., Kieliszewska-Rokicka, B., Modrzyński, J., Rożen, A., Skorupski, M., Sobczyk, Ł., Stasińska, M., Trocha, L. K., Weiner, J., ... Oleksyn, J. (2016). Light, earthworms, and soil resources as predictors of diversity of 10 soil invertebrate groups across monocultures of 14 tree species. *Soil Biology and Biochemistry, 92*, 184-198.

Naimi, B. (2017). usdm: Uncertainty Analysis for Species Distribution Models.

Navarro, L. M., Fernández, N., Guerra, C., Guralnick, R., Kissling, W. D., Londoño, M. C., Muller-Karger, F., Turak, E., Balvanera, P., Costello, M. J., Delavaud, A., El Serafy, G. Y., Ferrier, S., Geijzendorffer, I., Geller, G. N., Jetz, W., Kim, E.-S., Kim, H., Martin, C. S., ... Pereira, H. M. (2017). Monitoring biodiversity change through effective global coordination. *Current Opinion in Environmental Sustainability, 29*, 158-169.

Orgiazzi, A., Panagos, P., Yigini, Y., Dunbar, M. B., Gardi, C., Montanarella, L., & Ballabio, C. (2016). A knowledge-based approach to estimating the magnitude and spatial patterns of potential threats to soil biodiversity. *Science of the Total Environment, 545-546*, 11-20.

Owojori, O. J., Reinecke, A. J., & Rozanov, A B. (2008). Effects of salinity on partitioning, uptake and toxicity of zinc in the earthworm *Eisenia fetida*. *Soil Biology and Biochemistry, 40*, 2385-2393.

Panagos, P., Ballabio, C., Lugato, E., Jones, A., Borrelli, P., Scarpa, S., Orgiazzi, A., & Montanarella, L. (2018). Potential sources of anthropogenic copper inputs to European agricultural soils. *Sustainability, 10*, 2380.

Panagos, P., Jiskra, M., Borrelli, P., Liakos, L., & Ballabio, C. (2021). Mercury in European topsoils: Anthropogenic sources, stocks and fluxes. *Environmental Research, 201*, 111556.

Phillips, H. R., Guerra, C. A., Bartz, M. L., Briones, M. J., 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., Brose, U., Decaens, T., König-Ries, B., Loreau, M., ... Nuzzo, V. (2019). Global distribution of earthworm diversity. *Science, 366*(4646), 480-485.

Phillips, H. R., Bach, E. M., Bartz, M. L., Bennett, J. M., Beugnon, R., Briones, M. J., Brown, G. G., Ferlian, O., Gongalsky, K. B., Guerra, C. A., König-Ries, B., Krebs, J., Orgiazzi, A., Ramirez, K. S., Russel, D. J., Schwarz, B., Wall, D. H., Brose, U., Decaëns, T., ... Eisenhauer, E. (2021). Global data on earthworm abundance, biomass, diversity and corresponding environmental properties. *Scientific Data, 8*(1), 1-2.

Phillips, H. R. P., Cameron, E. K., Ferlian, O., Türke, M., Winter, M., & Eisenhauer, N. (2017). Red list of a black box. *Nature Ecology & Evolution, 1*, 0103.

Pörtner, H.-O., & Roberts, D C. (2022). IPCC AR6 report-Working group II: Climate change 2022: Impacts. *Adaptation and Vulnerability, 3068*.

R Core Team. (2020). *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing.

Rodrigues, A. S. L., Pilgrim, J. D., Lamoreux, J. F., Hoffmann, M., & Brooks, T M. (2006). The value of the IUCN Red List for conservation. *Trends in Ecology & Evolution, 21*, 71-76.

RStudio Team. (2020). *RStudio: Integrated Development for R*. R Studio, PBC, Boston, MA. Available from http://www.rstudio.com/

Ruiz, S. A., Bickel, S., & Or, D. (2021). Global earthworm distribution and activity windows based on soil hydromechanical constraints. *Communications Biology, 4*, 612.

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

Segurado, P., Araújo, M. B., & Kunin, W E. (2006). Consequences of spatial autocorrelation for niche-based models. *Journal of Applied Ecology, 43*, 433-444.

Singh, J., Schädler, M., Demetrio, W., Brown, G. G., & Eisenhauer, N. (2019). Climate change effects on earthworms: A review. *Soil Organisms, 91*, 114-138.

Soberon, J., & Peterson, A T. (2005). Interpretation of models of fundamental ecological niches and species' distributional areas. *Biodiversity Informatics, 2*, 1-10.

Steen, V. A., Tingley, M. W., Paton, P. W. C., & Elphick, C. S. (2021). Spatial thinning and class balancing: Key choices lead to variation in the performance of species distribution models with citizen science data. *Methods in Ecology and Evolution, 12*, 216-226.

Tedersoo, L., Bahram, M., Pölme, S., Köljalg, U., Yorou, N. S., Wijesundera, R., Ruiz, L. V., Vasco-Palacios, A. M., Thu, P. Q., Suija, A., Smith, M. E., Sharp, C., Saluveer, E., Saitta, A., Rosas, M., Riit, T., Ratkowsky, D., Pritsch, K., Pöldmaa, K., ... Abarenkov, K. (2014). Global diversity and geography of soil fungi. *Science, 346*, 1256688.

Thuiller, W., Georges, D., Engler, R., & Breiner, F. (2016). *biomod2*: Ensemble Platform for Species Distribution Modeling. https://cran.r-project.org/web/packages/biomod2/biomod2.pdf

Torma, M., & Harma, P. (2004). Accuracy of CORINE land cover classification in Northern Finland. In IGARSS 2004. 2004 IEEE International Geoscience and Remote Sensing Symposium (Vol. 1, pp. 227-230). IEEE.

UNEP-WCMC, IUCN. (2021). Protected Planet: The World Database on Protected Areas (WDPA) and World Database on Other Effective Area-based Conservation Measures (WD-OECM) [Online]. Cambridge, UK. Available from www.protectedplanet.net

Urban, M C. (2015). Accelerating extinction risk from climate change. *Science, 348*, 571-573.

Valavi, R., Guillera-Arroita, G., Lahoz-Monfort, J. J., & Elith, J. (2021). Predictive performance of presence-only species distribution models: A benchmark study with reproducible code. *Ecological Monographs, 42*, e01 486.

van den Hoogen, J., Geisen, S., Routh, D., Ferris, H., Traunspurger, W., Wardle, D. A., De Goede, R. G. M., Adams, B. J., Ahmad, W., Andriuzzi, W. S., Bardgett, R. D., Bonkowski, M., Campos-Herrera, R., Cares, J. E., Caruso, T., De Brito Caixeta, L., Chen, X., Costa, S. R., Creamer, R., ... Crowther, T. W. (2019). Soil nematode abundance and functional group composition at a global scale. *Nature, 572*, 194-198.

van Groenigen, J. W., Lubbers, I. M., Vos, H. M. J., Brown, G. G., De Deyn, G. B., & van Groenigen, K. J. (2014). Earthworms increase plant production: A meta-analysis. *Scientific Reports, 4*, 6365.

Vaughan, I. P., & Ormerod, S J. (2005). The continuing challenges of testing species distribution models. *Journal of Applied Ecology, 42*, 720-730.

Virkkala, R., & Lehikoinen, A. (2017). Birds on the move in the face of climate change: High species turnover in northern Europe. *Ecology and Evolution, 7*, 8201-8209.

Wickham, H., Chang, W., Henry, L., Pedersen, T. L., Takahashi, K., Wilke, C., Woo, K., Yutani, H., & Dunnington, D., RStudio. (2021). *ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics*. https://ggplot2.tidyverse.org/reference/ggplot2-package.html

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.
<figure><img src="image_23.png" /></figure>

Zizka, A., Silvestro, D., Andermann, T., Azevedo, J., Duarte Ritter, C., Edler, D., Farooq, H., Herdean, A., Ariza, M., Scharn, R., Svantesson, S., Wengström, N., Zizka, V., & Antonelli, A. (2019). CoordinateCleaner: Standardized cleaning of occurrence records from biological collection databases. *Methods in Ecology and Evolution, 10*, 744–751.

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

## SUPPORTING INFORMATION

Additional supporting information can be found online in the Supporting Information section at the end of this article.

> **How to cite this article:** Zeiss, R., Briones, M. J. I., Mathieu, J., Lomba, A., Dahlke, J., Heptner, L., Salako, G., Eisenhauer, N., & Guerra, C. A. (2024). Effects of climate on the distribution and conservation of commonly observed European earthworms. *Conservation Biology, 38*, e14187. https://doi.org/10.1111/cobi.14187

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