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

Belowground biodiversity distribution does not necessarily reflect aboveground biodiversity 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 estimated 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 occurrences 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.


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 biogeographical interest, species-level information allows for more effective conservation actions and the evaluation of the conservation 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 conditions 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 Appen-dices 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érome 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.,  Table 1]; n, number of occurrences at 4-km 2 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).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)   (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 ROC ), 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

Model workflow
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.

Model fitting
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.Software R 4.2.0 (R Core Team 2016), RStudio 2022.02.3 (RStudio Team, 2020), ArcGIS 10.7.1 (Bajjali, 2018)

Projections
For the 19 commonly observed earthworm species, we predicted and mapped the best-performing models (i.e., the one 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-earthwormdistribution).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 2 (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).Landcover 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.
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 2 ) 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.

Distribution of highly recorded European earthworms
The 19 focal earthworm species had broad potential distribution ranges across Europe (median [SD] = 2.1 million km 2 [880,943]) (Figure 2, Supporting Information Appendices S13 & S14).The species predicted to have the widest distribution were Lumbricus terrestris (3,4 million km 2 ) (Figure 2c), Dendrodrilus rubidus (3,2 million km 2 ), and L. castaneus (3,15 million km 2 ).In contrast, D. illyrica showed the smallest potential distribution range (8,125 km 2 ); 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 2 , and the lowest in eastern and northern Europe (Figure 2a).Nineteen species did not co-occur.On average, 3-5 species cooccurred 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 2 = 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. attemsi 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 2 = 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 2 ), 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 2 ).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. attemsi 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 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).
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%, max-imum = 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 ), (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).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 2 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(Dupont et al., , 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 ≤ 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 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 aboveground 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 (Rutgers 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. attemsi 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 aboveground species ranges, but also belowground organisms and their potential future distribution ranges, including those areas that could provide stable refugia.

FIGURE 1
FIGURE 1 Workflow of study of earthworm distribution and protection in Europe under current and future climatic conditions (small squares, predictor variables [sources inTable 1]; n, number of occurrences at 4-km 2 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).

FIGURE 4
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 Figure3).

TABLE 1
Environmental predictors used for species distribution modeling (SDM) of European earthworm species.

TABLE 1 (
Continued) Soil Moisture soil moisture index (SMI) (i.e., averaged daily moisture conditions based on hydrological rainfall-runoff model LISFLOOD) a Detailed information is in Appendix S1. b Predictors excluded based on Pearson correlation coefficient and variable inflation factor.c One of the 10 most important drivers of the investigated, highly recorded earthworm species.d Copernicus: European Union, Copernicus Land Monitoring Service 2018, European Environment Agency.

TABLE 2
Zurell et al. (2020)r reporting the performed species distribution models according toZurell et al. (2020)for European earthworm species.a.
Model averagingFor 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).
2 (upper edge in the range of scenarios described in the literature)