# Environmental drivers of earthworm communities along an elevational gradient in the French Alps

Quentin Gabriac, Pierre Ganault, Isabelle Barois, Eduardo Aranda-Delgado, Elisa Cimetière, Jérôme Cortet, Gautier Montan, Mickael Hedde, Daniel F. Marchán, José Carlos Pimentel Reyes, et al.

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

► To cite this version:

Quentin Gabriac, Pierre Ganault, Isabelle Barois, Eduardo Aranda-Delgado, Elisa Cimetière, et al. Environmental drivers of earthworm communities along an elevational gradient in the French Alps. European Journal of Soil Biology, 2023, 116, pp.103477. 10.1016/j.ejsobi.2023.103477. hal-04050281

HAL Id: hal-04050281

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

Submitted on 30 Jan 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_2.png" /></figure>
1 Environmental drivers of earthworm communities

2 along an elevational gradient in the French Alps

3 Quentin Gabriac¹, Pierre Ganault¹, Isabelle Barois², Eduardo Aranda-Delgado³, Elisa

4 Cimetière ¹, Jérôme Cortet ¹, Montan Gautier ¹, Mickaël Hedde ⁴, Daniel F. Marchán ¹, José

5 Carlos Pimentel Reyes ², Alexia Stokes ⁵, Thibaud Decaëns ¹*

6 ¹ CEFE, Univ Montpellier, CNRS, EPHE, IRD, Montpellier, France

7 ² Red de Ecología Funcional, Instituto de Ecología, A.C./INECOL, Veracruz, Mexico

8 ³ Terranova Lombricultores, Coatepec, Veracruz, Mexico

9 ⁴ UMR Eco&Sols, INRAE-IRD-CIRAD-SupAgro Montpellier, Montpellier, France.

10 ⁵ AMAP, Univ Montpellier, IRD, CIRAD, CNRS, INRAE, Montferrier-sur-Lez, France

11 * Corresponding author

12 Email: [Thibaud.Decaens@cefe.cnrs.fr](mailto:Thibaud.Decaens@cefe.cnrs.fr)

13 Abstract

14 Elevational diversity gradients describe how species diversity and richness change

15 with elevation. In this study, we tested for the effect of climate, resource quality and

16 habitat heterogeneity on earthworm communities along an elevational gradient and

17 around the treeline in the French Alps. Earthworm communities and environmental

18 properties (i.e. climate, soil properties, vegetation structure and composition) were

19 sampled at six elevations ranging from 1400 to 2400 m. Results were analysed through

20 multi-table factorial analyses and structural equation modelling. We found average

21 density, biomass and species richness in the range of that usually reported in comparable
1 Abstract

2
Elevational diversity gradients describe how species diversity and richness change
3
with elevation. In this study, we tested for the effect of climate, resource quality and
4
habitat heterogeneity on earthworm communities along an elevational gradient and
5
around the treeline in the French Alps. Earthworm communities and environmental
6
properties (i.e. climate, soil properties, vegetation structure and composition) were
7
sampled at six elevations ranging from 1400 to 2400 m. Results were analysed through
8
multi-table factorial analyses and structural equation modelling. We found average
9
density, biomass and species richness in the range of that usually reported in comparable
ecosystems. We found no monotonic decrease in species richness along the elevational
10
gradient, which we explain by the species pool being dominated by taxa with high
11
environmental tolerance and dispersal capacities. Instead, we highlighted the ecotone
12
associated with the treeline as the primary driving factor of earthworm communities: at
13
1800-2000m elevation, communities were more abundant and diverse, and had a greater
14
variability in body mass. This result was largely explained by the structure and
15
composition of the vegetation, whereas soil and climate appeared to have only indirect
16
effects. Therefore, the treeline effect on earthworm communities can be explained both
17
by the effect of environmental heterogeneity and of trophic resource quality that
18
increases at the ecotone level.
19
ecosystems. We found no monotonic decrease in species richness along the elevational gradient, which we explain by the species pool being dominated by taxa with high environmental tolerance and dispersal capacities. Instead, we highlighted the ecotone associated with the treeline as the primary driving factor of earthworm communities: at 1800-2000m elevation, communities were more abundant and diverse, and had a greater variability in body mass. This result was largely explained by the structure and composition of the vegetation, whereas soil and climate appeared to have only indirect effects. Therefore, the treeline effect on earthworm communities can be explained both by the effect of environmental heterogeneity and of trophic resource quality that increases at the ecotone level.

**Key words:** earthworms, elevation gradient, treeline, community ecology
33

# 1. Introduction

34

35

36

37

38

39

40

41

42

Studying the variations in the diversity and composition of ecological communities along elevational gradients has long fascinated early naturalists and ecologists, and still represents a major topic in biodiversity research [1,2]. The study of elevational diversity gradients (EDG) has contributed to a better understanding of the links between the environment and the diversity and structure of ecological communities, and is expected to provide insights into the responses of species and ecological systems to ongoing global change [3]. Understanding EDG is therefore fundamental for establishing conservation and management strategies for mountain ecosystems that are particularly sensitive to global warming and often subject to significant anthropogenic pressures [2].

43

44

45

46

47

48

49

50

51

52

53

54

55

Many studies of EDG describe a decline of species richness with increasing elevation, but whether this decline is monotonic, or if it assumes different shapes depending on the taxa or region being studied, has been an important topic of debate [4,5]. Currently, it is broadly accepted that the most widespread pattern observed corresponds to a mid-elevation peak of species richness, and that the monotonic variations that sometimes occur correspond to results from studies using an incomplete gradient and that did not fully cover all the potential sources of environmental heterogeneity [6]. Mid-elevational peaks in species richness have been observed for different taxa, such as mammals [7,8], birds [9], arthropods [10–13], vascular plants [14] among others, but was not documented for other organisms such as soil microbial communities [15,16]. Different hypotheses have been proposed to explain EDG, most of them involving the effects of elevational gradients on habitat area, biophysical conditions, the degree of geographical isolation of montane biota and turnover among zonal communities [1,4,17]. Alternatively,
56 in some cases, the observed patterns resulted from the spurious outcome of bias in
57 sampling regimes [1].

58 Another feature of elevational gradients that can influence EDG locally is the presence
59 of a treeline, which marks the upper limit in tree distribution and the transition between
60 low elevation closed forest and upper elevation alpine prairies. The study of the dynamics
61 of this ecotone is of major importance for the conservation of mountain ecosystems,
62 especially in the context of the global change which is likely to modify both its functioning
63 and elevational position [18-20]. However, studies of ecological communities on both
64 sides of the treeline have so far focused mainly on plant communities, leaving out other
65 functionally important groups of organisms, such as soil invertebrates (but see [21,22]).
66 Among them, earthworms deserve special attention because of their capacity, as
67 ecosystem engineers, to modify soil structure and water properties, and to stimulate the
68 decomposition of organic matter and the recycling of nutrients [23]. Paradoxically, there
69 are very few studies to date on the way earthworm communities vary in structure and
70 diversity along elevational gradients and around the treeline.

71 In addition to the above factors proposed to explain variations in biological diversity
72 along elevational gradients, vegetation composition is also likely to strongly influence
73 earthworm communities. Plant species can influence soil detritivorous fauna through the
74 provision of energy and matter via living and dead plant products such as leaf and root
75 litter and rhizodeposition [24]. Deciduous trees in particular are known to produce a litter
76 of higher palatability compared to conifers [24,25], and alpine herbaceous plants also
77 produce a high quality trophic resource compared to trees and shrubs, either through
78 litter or dead roots [26,27]. At comparable elevations, earthworm abundance and
79 diversity tend to decrease with the proportion of conifers in the standing vegetation
80 [24,25], and are greater in pasture compared to forest ecosystems [28,29]. Along an
elevational gradient, changes in vegetation composition, including variations in the proportion of conifers, as well as the shift to more herbaceous vegetation above the treeline, are therefore likely to be associated with changes in earthworm biomass and diversity.

In this study, we sought to understand the factors controlling earthworm communities along an elevational gradient around the treeline in the French Alps. We tested three alternative hypotheses that could explain the variation of earthworm diversity along this elevation gradient: 1) a climate hypothesis, that assumes that earthworm biomass and diversity decrease with elevation as a result of increasingly adverse temperature conditions [1,30], and that communities are dominated by large body-sized species at higher elevations as a result of Bergmann's rule [31]; 2) the resource quality hypothesis, that predicts an increase in earthworm biomass and body mass with elevation as a result of the shift from gymnosperm stands to alpine grasslands, and the consecutive increase in the quality of litter inputs to the soil food web [32]; 3) the habitat complexity hypothesis, that assumes that taxonomic and functional community diversity peak at the ecotone level due to greater diversity in microhabitats [8].

## 2. Material and methods

### 2.1. Study sites

The study was carried out during May 2018 in the Belledonne massif in the French Alps (Chamrousse, Isère, France, N45°7'1" E5°53'35") (Fig. A.1), along a single 8 km long elevational gradient ranging from 1400 to 2400 m above sea level. The gradient transect was linear and constant in terms of its exposure. It was part of a project whose objective was to explore the relationships between biological communities (vegetation, soil invertebrates and microbial communities) and soil properties along elevational gradients
in France and Mexico (ECOPICS project, ANR-16-CE03-0009 and CONACYT-2 73659). A dataset comprising detailed information about climatic data, vegetation and soil biological, chemical and physical properties is available at [33].

The geological substrate was composed of Variscan metamorphic rocks and ophiolitic complexes, and soils were classified as Umbrisols and Cambisols [33,34]. Climatic data were estimated over 2004–2014, using the Aurelhy model [35]. Vegetation consisted of mixed and coniferous forests to alpine dwarf shrub heathland and meadows, with a transition between forest and open vegetation (i.e. ecotone) at 1800–2000 m and a treeline located between 2000 and 2100 m (Fig. A.2).

Along this elevational gradient, sampling was carried out at six elevational levels spaced every 200 m of elevation from 1400 to 2400 m. Five plots (20 x 20 m) were selected within each elevational band, with a lateral distance of at least 100 m between each plot. Plots had a south-west orientation (azimuth of 220 ± 26°) and a mean angle of 17.5 ± 2.3° from the vertical. None of them were used for grazing domestic herbivores. Areas with visible bedrock or signs of waterlogging were avoided. Plots were selected if they included three species: *Picea abies* (L.) H. Karst (Pinaceae), a tall evergreen tree; *Juniperus communis* L. (Cupressaceae), a prostrate evergreen shrub and *Vaccinium myrtillus* L. (Ericaceae), a small deciduous shrub (if these species were present at that elevational band). These three plant species were selected because of their different growth forms and occupation of different ecological niches along the elevational gradient. *Picea abies* and *V. myrtillus* are both keystone species [36,37] and although *J. communis* is not classed as a keystone species, its abundance above the treeline makes it an important species that structures plant communities. Therefore, these three species contribute to shaping the structure of the plant communities in which they are present, and we term them ‘structuring species’ [38]. Structuring species were chosen because of their different
growth forms and their ecological niche along the elevation gradient. *Picea abies* is the dominant tree below the treeline, *J. communis* is abundant in the shrub layer above the treeline, while *V. myrtillus* occurs in the shrub layer along the entire elevational gradient. Therefore, the importance of these species suggests that they could have a significant long-term influence on the biotic (i.e. plant - plant and plant - soil interactions) and abiotic (i.e. physical and chemical properties) soil properties along the elevational gradient.

The proportion of the ground covered by trees, shrubs, herbs, and bryophytes, the proportion of mineral elements and bare soil, as well as the relative proportion of the three structuring plant species within their own strata were visually estimated on each 20 x 20 m plot (Fig. A.2). *Vaccinium myrtillus* was present at all elevation levels and was more abundant below the treeline (17% vs 2% on average). The relative cover of *P. abies* among trees decreased steadily from 31% at 1400 m to 3% at 2000 m. *Juniperus communis* was present from 1800 to 2400 m with an average relative coverage of 8.5%.

## 2.2. Sampling and description of earthworm communities

Within each 20 x 20 m plot, a 25 x 25 x 15 cm soil block (or monolith) was taken from under one individual of each of the structuring plant species present. The exact location of this monolith was at the limit of the projection on the ground of the canopy of the individual plant. Earthworms were collected by hand sorting each monolith in a large tray, immediately fixed in 100% ethanol, and stored at -21° C. Total sampling varied from 10 to 15 monoliths per elevational level depending on vegetation composition.

Adult and sub-adult individuals were identified at species level on the basis of their external morphology [39]. Identification of juveniles, which often do not present the necessary morphological characters, were achieved by sequencing the DNA barcode of
the cytochrome c oxidase I (COI) gene [40] through standard protocols [41] implemented at the Canadian Centre for DNA Barcoding (CCDB), University of Guelph, Ontario, Canada. Sequences obtained for juveniles were further compared to those obtained for adult specimens, or alternatively to identified sequences already available in the Barcode of Life Data systems (BOLD; [42]). This enabled the species-level identification of all the unidentifiable specimens for which DNA barcoding was successful [43]. All sequence records, voucher deposition details, specimen data, GPS coordinates and photographs are publicly available within the DS-EWCHAM dataset in BOLD (dx.doi.org/10.5883/DS-EWCHAM).

Community functional properties were described using trait-based ecological categories and species body mass. We used data from Bottinelli et al. [44] to attribute to each species affinity scores regarding the three ecological categories defined by Bouché [39] (Table A.2). These categories can be broadly described as follow [44]: epigeics are surface dwelling species, small sized, light in weight and present pigmented body; anecic live in vertical burrows and are relatively large and pigmented species with flattened tail; endogeics are soil dwelling species deprived of pigmentation, with variable body size. Each individual was then weighed individually after being carefully dried, and body mass values were aggregated to obtain average biomass.

For each plot we calculated the community weighted mean (CWM) of ecological category scores and of body mass, and the community weighted variance (CWV) of body mass using the following formula [45-47]:

$$CWM_k = \sum_{i=0}^{S} p_{ik} \times t_i$$
$$176 \qquad CWV_k = \sum_{i=0}^{S} p_{ik} (t_i - CWM_k)^2$$

177 where S is the number of species in community k, pik is the relative abundance of species i in community k, and t_i is the mean value taken by trait t and for species i.

179 **2.3. Soil properties**

180 At each monolith location, the depth of litter and organo-mineral layers was recorded 181 on three randomly chosen sides of the pit. Bulk density was measured on a 5 cm diameter 182 and 5.9 cm depth soil core (i.e. 115.89 cm²) which was dried at 105 °C before being 183 weighed (BD expressed in g of dry soil. cm⁻³). Soil physical and chemical properties were 184 measured at the INRA-Arras (France) laboratory on a sample aliquot taken from each 185 monolith, air dried and sieved at 2 mm. Soil texture was determined on non-decarbonated 186 soil using the Robison method (NF X 31-107) [48]. Water pH was measured using the NF 187 ISO 10390 method (H₂O/KCl/CaCl₂). Different soil chemical properties, which are 188 expected to be indirectly related to the quality of litter inputs by the vegetation, were also 189 measured. NaHCO₃ extractible phosphorus was measured at a pH of 8,5 (NF ISO 11263) 190 using the Olsen et al. method [49]. The cationic exchange capacity (CEC, expressed in cmol. 191 kg⁻¹) was estimated in a cobalihexamine solution following the NF X 31-130 method [50]. 192 Total organic carbon and total nitrogen were determined after dry combustion using the 193 NF ISO 10694 and NF ISO 13878 methods, respectively [51]. NH₄⁺ and NO₃⁻ were 194 measured on soil samples conserved at -20 °C after KCl extraction, using colorimetric tests 195 as recommended by ISO 14256-2) (analyses done at US Analyses, Cirad, Montpellier).

196 **2.4. Statistical analyses**
The overall effect of elevation and vegetation (i.e. structuring plant species under which the monolith was sampled) on earthworm biomass (g of fresh mass. m⁻²), abundance (number of individuals. m⁻²), body mass CWM (g) and body mass CWV were tested using Kruskal-Wallis rank sum test with the 'kruskal.test' function of R. Pair wise differences among elevational levels and plant species were tested using pairwise Wilcoxon rank sum tests with the 'pairwise.wilcox.test' function of R 4.0.3 [52].

Local species richness was calculated as the cumulative number of species observed in the monoliths of a given study plot. To compare the richness of the species pools corresponding to each elevational level, the data from the study plots were pooled in order to plot one rarefaction curve (inflation and extrapolation) for each altitude. This was done using the 'iNEXT' package considering the number of individuals as a measure of the sampling effort. We then used the same package to calculate the asymptotic estimator of Chao for each elevational level [53–55].

The relative effect of environmental parameters on the structuring of earthworm communities was assessed through a set of factorial analyses. The data were first organised into six tables sharing 30 rows corresponding to the study plots (i.e. five plots in each of the six elevational levels). The earthworm data were thus formatted into a table containing species abundance, and a table describing community structure through a set of selected metrics (see Fig. A.4 for a list of species and metrics used in these two tables). The environmental data were organised into a *Clim* table containing climatic variables, a *Soil* table containing soil physicochemical properties, a *Veg1* table containing vegetation structure descriptors and a *Veg2* table containing vegetation compositional variables (see Fig. A.3 for a list of these variables).

The earthworm species abundance table was first analysed using a principal component analysis (PCA) after applying a Hellinger transformation to the abundance
data using the 'deconstand' function of 'Vegan' package [56]. Hellinger transformation converts species abundances from absolute to relative values (i.e. standardizes the abundances to sample totals) and then square roots them; it is recommended for clustering or ordination of species abundance data [57]. We then used the scores of the study plots on the first two axes of this PCA as two synthetic composition variables named *Com.1* and *Com.2*. These two variables were added to the community structure variables to constitute an *EW structure* table that was further analysed with another PCA. The scores of study plots on the axes of this PCA were extracted and referred to as *Str.1* and *Str.2* in further RDA and path analyses (see Fig. A.4 for a complete description of this analysis).

The main axes of environmental variability (i.e. climate, soil and vegetation) along the environmental gradient were also analysed through four individual PCAs. The scores of the study plots on the significant axes of these PCAs were used as synthetic environmental variables for further analyses. By significant axes we mean those axes of a particular PCA that when combined explained at least 70% of the total variance in the analysed table. We thus retained the first axis of the PCA of climatic variables (*Clim.1*), as well as the first two axes of the other three PCAs (denoted *Soil.1*, *Soil.2*, *Veg1.1*, *Veg1.2*, *Veg2.1* and *Veg2.2*, respectively). The synthetic interpretation of each of the retained axes is given in Table 1 (see complete interpretation in Appendix 2) All these analyses were done using the 'ade4' package for R [58].

We then ran a redundancy analysis (RDA) to figure out the relative importance of climate, soil and vegetation in explaining elevational variations in earthworm community structure. This was done using the 'rda' function of the 'Vegan' package. The community data was the *EW structure* table, and the explicative variables were the plot scores on the retained axes of the PCAs of environmental data. The model formula
247 ('EW 结构.表~Clim.1+Soil.1+Soil.2+Veg1.1+Veg1.2+Veg2.1+Veg2.2') included all
248 possible interactions among the different explanatory variables. Permutation tests (with
249 999 randomization) were used to assess the overall significance of the analysis, the
250 significance of each constrained axis, as well as the significance of each term using the
251 function 'anova.cca'. We used the 'varpart' function of 'Vegan' to partition the variation in
252 earthworm community data with respect to the main explanatory variables highlighted
253 by the RDA.

254 In order to highlight potential indirect effects of environmental variables on
255 earthworm community structure, we performed path analyses using the 'psem' function
256 of the 'piecewiseSEM' package [59]. To fix the theoretical model required for these
257 analyses, we first explored all possible correlations between both community (Str.1 and
258 Str.2) and environmental (Clim.1, Soil.1, Soil.2, Veg1.1, Veg1.2, Veg2.1, Veg2.2) synthetic
259 variables. We thus built a theoretical complete model figuring all the potential
260 interactions between community and explanatory variables, and ran the path analysis
261 using generalised linear models (glm) to build the list of structural equations.

262 3. Results

263 3.1. *Species richness and community composition*

264 We found a total of six genera and 11 species of earthworms across the whole
265 elevational gradient (Fig. 1, Table A.2). While most mean species occurrences were below
266 the treeline (2000-2100 m), seven species were distributed in both forest ecosystems and
267 alpine heathland and meadows. The species *Dendrodrilus subrubicundus* (Eisen, 1874)
268 was the most alticolous in our survey, with a mean occurrence just above 2200m.

269 The elevational variation of earthworm species density (number of species observed
270 per sample) and cumulated richness at the scale of each elevational level are shown in
Fig. 2. Species density ranged from 0 to 7 species per sample (mean ± sd: 1.4 ± 1.7; Fig. 2A), and cumulated species richness ranged from 5 to 9 species per elevational level (Fig. 2B). Species diversity generally peaked at the ecotone just below the treeline between 1800 and 2200 m of altitude. This trend was significant while considering species density, but it was partly masked when considering cumulated richness by the high inter-plot variability in species composition within each elevational level. The Chao asymptotic richness showed a similar pattern (Table 2).

The Kruskal-Wallis test did not highlight any significant effect of vegetation composition (structuring species) on earthworm species density (result not shown), while on the contrary the effect of altitude was significant (p-value < 0.01). The pairwise Wilcoxon rank sum tests show that this effect was essentially driven by the specific diversity at 1800 m, which differed significantly from that of all the other altitudes.

We found a significant pattern of variation of epigeic and endogeic scores along the elevational gradient, while no significant trend was observed for anecics (Fig. 3). The pairwise Wilcoxon rank sum tests highlight that this was mainly due to communities at 1800 and 2200 m of altitude, which had a significantly lower epigeic score and a significantly higher endogeic score compared with those at other elevational levels. No effect of vegetation was detected on ecological category scores (p-value > 0.05; result not shown).

## 3.2. Density, biomass and body mass

Community density ranged from 0 to 416 individuals per m² (mean ± sd: 51.4 ± 86.0), and biomass ranged from 0 to 135 g per m² (mean ± sd: 12.7 ± 25.2). Both density and biomass showed a similar pattern of variation along the elevational gradient, with a peak reached at the transition zone (Fig. 4A-B). The Kruskal-Wallis tests showed no significant
effect of vegetation (p-value > 0.05), while altitude significantly explained the variability in both biomass and density (p-value < 0.01 and 0.05, respectively).

No significant variation in body mass CWM was observed along the elevational gradient (Fig. 4C.), and no effect of vegetation type was also found (p-value > 0.05; result not shown). On the contrary, body mass CWV varied along the gradient in a significant way and according to a pattern comparable to that observed for biomass and density, with values peaking at 1800 and 2000m of altitude (Fig. 4D).

### 3.3. Environmental drivers of community structure

The redundancy analysis was significant after the permutation test was applied (p-value < 0.001). The first two axes explained 58.4% and 19.8% of the covariation between community and environmental tables (Fig. 5A). The proportion of constrained inertia was of 43.5% (Fig. 5D). Only the first axis appeared to be significant after the permutation test was applied (p-value < 0.001; Fig. 5E).

The first RDA axis highlighted the ecotone effect in the structuring of communities. Indeed, plots located in the transition zone just below the treeline (i.e. 1800 and 2000m elevational levels) had positive scores on this axis, and were clearly opposed to all other elevational levels (Fig. 5A). The projection of the community variables highlighted that earthworm communities in these plots were characterized by higher density, biomass and species richness, more variable body mass among co-existing species, and community composition dominated by species displaying negative scores on the community component *Com.1* (essentially *Lumbricus terrestris* L., 1758 and *Octolasion cyaneum* (Savigny, 1826); Fig. 5B, Fig. A.4B). The projection of environmental variables on RDA1 showed this ecotone effect was mainly driven by soils richer in nitrates and with a higher pH (positive projection of *Soil.2* on RDA1), and by a greater representation of the
herbaceous stratum (grasses and forbs) in the vegetation structure (negative projection of *Veg1.2* on RDA1; Fig. 5C, Fig. A.3F).

Even though axis 2 of the RDA was not significant after the permutation test (p-value = 0.18), we decided to keep it in the interpretation because it highlighted an effect of the elevational gradient that could be partially hidden by the dominant ecotone effect. Plots located below 1800m had positive averaged scores on RDA2, while plots above 2000m displayed negative scores (Fig. 5A). The projection of community variable on this axis showed that lower elevation levels hosted communities with a higher average epigeic score, mainly because of the presence of *Lumbricus rubellus* (Hoffmeister, 1843), *Dendrobaena octaedra* (Savigny, 1826) and *Dendrodilus rubidus* (Savigny, 1826), and lower body mass CWM (Figs. 5B, 5C, Fig. A.4B). This effect was mainly related to an elevational gradient in soil properties (i.e. *Soil.1*; Fig. 5C), i.e. an increase in sand content and bulk density with altitude, accompanied by a decrease in organic matter and nutrient content (Fig. A.4D). Vegetation structure was also linked, albeit less clearly, to RDA2, mainly because of a gradual replacement of *J. communis* by other shrubs with increasing altitude (i.e. *Veg2.2*; Fig. 5C, Fig. A.4H).

The variance partitioning analysis highlights that a large proportion of the variance in earthworm community structure was not explained by the monitored environmental drivers (residuals = 73.3%; Fig. A.5). Among the explained fraction of the variance, a significant proportion (i.e. 16.4% of the total variance) was explained by the interaction between the variables *Soil.2* and *Veg1.2*. Then, the main environmental predictors, in decreasing order of importance, were *Veg1.2*, *Soil.1*, *Veg2.2* and *Soil.2* (Fig. A.5).

The path analyses highlight the main direct and indirect environmental controls on earthworm communities. The ecotone effect on community structure (*Str.1*), which is described by the first axes of both the RDA (Fig. 5B) and the PCA of community structure
344 (Figs. A.4C-D), was firstly explained by an indirect effect of soil nitrate concentration and
345 pH (Soil.2) mediated by the importance of the herbaceous stratum in vegetation structure
346 (Veg1.2), and secondly by a direct effect of the composition of the shrub stratum in the
347 vegetation (Veg2.2) (Fig. 6A). The elevational effect on community structure (Str.2), which
348 is described by the second axes of both the RDA (Fig. 5B) and the PCA of community
349 structure (Figs. A.4C-D), was explained by an effect of soil properties (i.e. elevational
350 variations in sand content, bulk density, organic matter and nutrient contents; Soil.1)
351 partially determined by climate (Clim.1), and by a direct effect of shrub stratum
352 composition (Veg2.2) (Fig. 6).

353 **4. Discussion**

354 The mean values of density and biomass (51 individuals per m-2, and 12 g.m-2,
355 respectively), and species richness (1.4 species per sample, and 5 to 7 species per
356 elevational level) that we found in our study are within the range of what is generally
357 reported for comparable temperate ecosystems. Earthworm communities in montane or
358 boreal coniferous forests, for example, have abundances ranging from 0 to about 100
359 individuals. m-2, biomasses ranging from 0 to 25 g. m-2, and are generally composed of
360 fewer than 5 species [60-63]. Similarly, in mountain pastures, earthworm density can
361 vary from less than 10 to 200 individuals per m-2, biomass from 15 to 30 g. m-2, and
362 communities generally comprise 2 to 4 species [63-66]. These figures are on average
363 lower than those generally given for lowland forest and grassland ecosystems [67-70].
364 The fact that earthworm communities are less abundant and less diverse in the mountains
365 than in the lowlands can be explained by the unfavourable climatic conditions, the
366 vegetation dominated by species producing litter with a high carbon to nitrogen ratio
367 (especially conifers), as well as by the presence of shallow soils in the higher elevational
368 levels of the gradient.

369 Our study does not show a monotonic decrease in species richness along the
370 elevational gradient. Few studies have looked at the variation in species diversity of
371 earthworm assemblages along elevational gradients. Some have described a monotonic
372 decrease in species number with altitude, for example in France [39,71] or in Serbia [72].
373 It should be noted, however, that these previous studies considered a complete
374 elevational gradient (i.e. from sea level to alpine levels), which is not the case in our study.
375 In another study in Puerto Rico, Gonzales et al. [73] observed, on the contrary, an increase
376 in the number of earthworm species along an elevational gradient from 0 to 1000 m,
377 which the authors interpreted as the result of the combination of several soil and climatic
378 factors. Finally, Fontana et al. [22] described a pattern similar to ours along a gradient
379 from 1000 to 2500 m in northern Italy, with a peak in species diversity observed between
380 1800 and 2200 m altitude. Our results therefore confirm the latter study and largely
381 invalidate the hypothesis of a monotonic effect of altitude and climate. The absence of a
382 monotonic decrease in the number of species with altitude can be explained in our case
383 by the composition of the regional species pool, which is essentially dominated by species
384 with a wide distribution in Northern Europe. After the last Würm glaciation, the areas left
385 free by the retreat of glaciers were recolonised by species with a high dispersal capacity
386 and ecological plasticity [39]. This mechanism occurred simultaneously along a
387 latitudinal gradient towards northern Europe, allowing the establishment of earthworm
388 communities up to the present-day boreal forests, and in mountain areas up to the alpine
389 levels. We therefore did not observe the increase in endemism classically observed with
390 altitude for other taxa [74], nor did we observe a steady decrease in species richness along
391 our elevational gradient.
Our study highlights the importance of the ecotone associated with the treeline as the primary factor responsible for the structuring of earthworm communities. At the level of this ecotone, earthworm communities are more abundant and diverse, and show a higher CWV of body mass, notably due to the presence of large species such as *L. terrestris* and *O. cyaneum* co-occurring with smaller body-sized species such as *Aporrectodea caliginosa* (Savigny, 1826) and *Lumbricus castaneus* (Savigny, 1826). The observation of a peak in species richness in the forest-heathland and grassland transition zone at around 1800-2000 m altitude is at first sight suggestive of the mid-domain effect, which predicts a higher species diversity in this zone due to the overlapping elevational distributions of certain species [75,76]. However, this interpretation does not seem convincing in the context of our study, given the already mentioned broad ecological spectrum that characterizes the species composing the regional pool. Indeed, it is unlikely that the species present along our gradient would show a preference for a particular elevational zone, and the peak in diversity observed at the ecotone level must in our view be attributed to ecological factors rather than to a geometric overlap of species distributions.

According to our results, community structure at the ecotone is associated with high nitrate and soil pH, but this effect of soil properties is indirect and largely mediated by vegetation structure. In particular, the presence of herbaceous cover and shrub layer appears to be the main determinant of earthworm communities. This result contrasts with the findings of other works that emphasize soil properties (including pH, C:N, texture and depth) as the dominant factor in explaining the distribution of earthworm abundance, biomass or diversity [61,62,68]. However, Phillips et al. [77] demonstrated that the effect of soil properties on earthworm communities may be scale-dependent, and could be hidden by other factors such as habitat cover or annual rainfall depending on the study context. In our study, the deliberate inclusion of the ecotone in the sampling design, and
the fact that this area clearly contrasted environmentally with the other parts of the elevational gradient, probably contributed to emphasizing vegetation structure rather than soil properties, whose direct effect is relegated to axis 2 of the RDA. Moreover, the use of path analyses is not common in studies on the driving factors of earthworm communities, whereas in our study they allowed to highlight this indirect effect. The generalization of this type of approach would probably allow a better understanding of the complex pathways of environmental controls on earthworm communities.

Several studies have highlighted that transitional ecosystems represent favourable habitats for earthworms. Bernier and Ponge [78] and Grossi & Brun [79] had shown that along mountain succession gradients, transitional phases between pastures and forests host abundant and diversified earthworm communities. In lowland environments, other studies have described the same type of pattern in the early stages of post-pastoral or post-cultural successions [28,29]. This pattern can firstly be explained by a greater diversity of micro-habitats at the level of the ecotone, compared to that of the ecosystems of which it marks the transition [1,17]. For earthworm communities, the ecotone provides alternating patches of herbaceous and woody vegetation, which influence soil microclimate and in which litter accumulations vary in thickness and quality depending on the identity of the plant species (i.e. coniferous versus deciduous). This environmental variability allows species with different ecological requirements to co-exist in the same communities. This point, which is consistent with our hypothesis of habitat complexity, is further supported in our results by the observed peak in body-mass CWV at the ecotone, which suggests that earthworm communities are morphologically and therefore functionally more diverse in this area of the gradient.

The effect of the ecotone can be linked to the transition from a vegetation dominated exclusively by conifers to a mixed vegetation, i.e. with a significant proportion of
442 deciduous species. A reduced proportion of conifers in the vegetation can have an impact
443 on earthworm communities through an improvement of the quality of trophic resources.
444 Indeed, it has been shown in several studies that mixed forests can host more abundant
445 and diverse earthworm communities than coniferous forests can [61,80]. The openness
446 of the canopy at the ecotone level also allows the establishment of dense herbaceous
447 vegetation, which produces a litter layer rich in nutrients such as calcium and nitrogen
448 [24,68,81]. Herbaceous species also participate, through both the quality of their organic
449 matter input and their dense root mat, in the formation of a deep A horizon that creates
450 suitable conditions for the activity of endogeous species [24,29,67]. This phenomenon is
451 supported in our results by the average score of the endogeic ecological category peaking
452 at the ecotone. The general improvement of habitat suitability is further supported by the
453 highest earthworm abundance and biomass observed at both 1800 and 2000 m of
454 elevation. In general, our results show that, in addition to the environmental
455 heterogeneity hypothesis, the structuring of earthworm communities in response to the
456 treeline can also be explained using the resource quality hypothesis.

457 Compared to the effect of the ecotone, altitude as such is of less importance in our
458 results. However, earthworm communities at lower altitudes show a higher epigeic score
459 and a lower CWM of body mass, even though we found no evidence that Bergmann's law
460 applies to our case study. This tendency for an elevational effect is mainly related to soil
461 properties, i.e. an increase with altitude in sand content and bulk density, followed by a
462 decrease in organic matter and nutrient content, and an increase in soil C:N. A smaller
463 part of the variability in community structure is further explained by climate and
464 vegetation composition. These results confirm the already observed trend towards low-
465 abundant and epigeous-dominated earthworm communities in coniferous forests, as a
466 consequence of low soil pH and high soil C:N in this type of forests [68,69,80]. Soil texture
467 is another factor recognized for its importance in community structuring, especially in
468 mountain environments. Salomé et al. [62] showed that epigeic species were more
469 abundant in shallow, sandy soils, whereas anecics and endogeics were more abundant in
470 deeper, finer-textured soils. Behind the dominant effect of the ecotone associated with the
471 treeline, elevational variability in soil properties, rather than the effect of altitude as such,
472 can explain a moderate fraction of the variability in earthworm community structure. As
473 already suggested earlier, it is possible that this hierarchy of driving factors is scale-
474 dependent, and that soil factors might be of primary importance if the sampling design
475 had been designed to explain community structuring within each elevational level rather
476 than along the whole gradient.

477 Overall, our results largely invalidate the hypothesis of a monotonic effect of climate,
478 but fail to separate the potential effect of environmental heterogeneity and resource
479 quality, both of which can be invoked to explain the observed community pattern at the
480 ecotone level. In the context of global change, an upward shift of the treeline is expected
481 as a result of climate warming [82,83]. In contrast to what is predicted for other taxa
482 [74,84], the prospect of this shift does not appear to represent an *a priori* threat to the
483 diversity of earthworm species, which includes no montane endemics at our study site.
484 However, the functional implications of these changes could be considerable. In
485 particular, a shift of the ecotone to higher elevational levels would inevitably lead to a
486 reduction in its surface area, and consequently in the size of the associated earthworm
487 activity hotspot. In the future, it will be necessary to examine the consequences that this
488 could have on ecosystem functioning at the scale of such mountain landscapes.

489 **5. Acknowledgments**
490 The authors acknowledge the Mairie de Chamrousse, Isère, France, for access to field sites,
491 and all the participants to the ECOPICS project for their contribution to logistics during
492 field sampling. The study was founded by the Agence Nationale de la Recherche from
493 France and the Consejo Nacional de Ciencia y Tecnología from Mexico (ECOPICS project,
494 ANR-16-CE03-0009 and CONACYT-2 73659).

495 **6. References**

496 [1] M.V. Lomolino, Elevation gradients of species-density: historical and prospective
497 views, Global Ecology & Biogeography. 10 (2001) 3–13.

498 [2] C.M. McCain, J.A. Grytnes, Elevational gradients in species richness, in: Encyclopedia of Life Sciences, John Wiley & Sons, Ltd, Chichester, 2010.

500 [3] C. Rahbek, The role of spatial scale and the perception of large-scale species-richness
501 patterns, Ecology Letters. 8 (2005) 224–239. https://doi.org/10.1111/j.1461-
502 0248.2004.00701.x.

503 [4] C. Rahbek, The elevational gradient of species richness - A uniform pattern?,
504 Ecography. 18 (1995) 200–205. https://doi.org/10.1111/j.1600-0587.1995.tb00341.x.

505 [5] J.D. Gilgado, H.-P. Rusterholz, B. Braschler, S. Zimmermann, Y. Chittaro, B. Baur,
506 Six groups of ground-dwelling arthropods show different diversity responses along
507 elevational gradients in the Swiss Alps, Plos One. 17 (2022) e0271831.
508 https://doi.org/10.1371/journal.pone.0271831.

509 [6] Q. Guo, D.A. Kelt, Z. Sun, H. Liu, L. Hu, H. Ren, J. Wen, Global variation in
510 elevational diversity patterns, Scientific Reports. 3 (2013) 3007.
511 https://doi.org/10.1038/srep03007.

512 [7] J.H. Brown, Mammals on mountainsides: elevational patterns of diversity, Global
513 Ecology and Biogeography. 10 (2001) 101–109. https://doi.org/10.1046/j.1466-
514 822x.2001.00228.x.

515 [8] C.M. McCain, The mid-domain effect applied to elevational gradients: species
516 richness of small mammals in Costa Rica, Journal of Biogeography. 31 (2004) 19–31.
517 https://doi.org/10.1046/j.0305-0270.2003.00992.x.

518 [9] I. Quintero, W. Jetz, Global elevational diversity and diversification of birds, Nature.
519 555 (2018) 246+. https://doi.org/10.1038/nature25794.

520 [10] J. Beck, C.M. McCain, J.C. Axmacher, L.A. Ashton, F. Bartschi, G. Brehm, S.-W.
521 Choi, O. Cizek, R.K. Colwell, K. Fiedler, C.L. Francois, S. Highland, J.D. Holloway, J.
522 Intachat, T. Kadlec, R.L. Kitching, S.C. Maunsell, T. Merckx, A. Nakamura, E. Odell, W.
523 Sang, P.S. Toko, J. Zamecnik, Y. Zou, V. Novotny, Elevational species richness gradients in a
524 hyperdiverse insect taxon: a global meta-study on geometrid moths, Global Ecology and
525 Biogeography. 26 (2017) 412–424. https://doi.org/10.1111/geb.12548.

526 [11] G. Chaladze, S. Otto, S. Tramp, A spider diversity model for the Caucasus Ecoregion,
527 Journal of Insect Conservation. 18 (2014) 407–416. https://doi.org/10.1007/s10841-014-9649-
528 1.

529 [12] N.J. Sanders, Elevational gradients in ant species richness: area, geometry, and
530 Rapoport's rule, Ecography. 25 (2002) 25–32. https://doi.org/10.1034/j.1600-
531 0587.2002.250104.x.

532 [13] M.A. Smith, W. Hallwachs, D.H. Janzen, Diversity and phylogenetic community
533 structure of ants along a Costa Rican elevational gradient, Ecography. 37 (2014) 720–731.

534 https://doi.org/10.1111/j.1600-0587.2013.00631.x.

535 [14] J.A. Grytnes, Species-richness patterns of vascular plants along seven altitudinal

536 transects in Norway, Ecography. 26 (2003) 291–300. https://doi.org/10.1034/j.1600-

537 0587.2003.03358.x.

538 [15] J.A. Bryant, C. Lamanna, H. Morlon, A.J. Kerkhoff, B.J. Enquist, J.L. Green,

539 Microbes on mountainsides: Contrasting elevational patterns of bacterial and plant diversity,

540 Proceedings of the National Academy of Sciences of the United States of America. 105

541 (2008) 11505–11511. https://doi.org/10.1073/pnas.0801920105.

542 [16] L. Merino-Martín, D. Hernández-Cáceres, F. Reverchon, G. Angeles-Alvarez, G.

543 Zhang, D. Dunoyer de Segonzac, D. Dezette, A. Stokes, Habitat partitioning of soil microbial

544 communities along an elevation gradient: from plant root to landscape scale, Oikos. 2022

545 (2022) e09034. https://doi.org/10.1111/oik.09034.

546 [17] C.M. McCain, J.A. Grytnes, Elevational gradients in species richness, in:

547 Encyclopedia of Life Sciences, John Wiley & Sons, Ltd, Chichester, 2010.

548 [18] S. Greenwood, A.S. Jump, Consequences of treeline shifts for the diversity and

549 function of high altitude ecosystems, Arctic Antarctic and Alpine Research. 46 (2014) 829–

550 840. https://doi.org/10.1657/1938-4246-46.4.829.

551 [19] S. Dullinger, T. Dirnbock, G. Grabherr, Modelling climate change-driven treeline

552 shifts: relative effects of temperature increase, dispersal and invasibility, Journal of Ecology.

553 92 (2004) 241–252. https://doi.org/10.1111/j.0022-0477.2004.00872.x.

554 [20] M.A. Harsch, P.E. Hulme, M.S. McGlone, R.P. Duncan, Are treelines advancing? A

555 global meta-analysis of treeline response to climate warming, Ecology Letters. 12 (2009)

556 1040–1049. https://doi.org/10.1111/j.1461-0248.2009.01355.x.

557 [21] E. Guariento, K. Fiedler, Ant diversity and community composition in Alpine tree line

558 ecotones, Insects. 12 (2021) 219. https://doi.org/10.3390/insects12030219.

559 [22] V. Fontana, E. Guariento, A. Hilpold, G. Niedrist, M. Steinwandter, D. Spitale, J.

560 Nascimbene, U. Tappeiner, J. Seeber, Species richness and beta diversity patterns of multiple

561 taxa along an elevational gradient in pastured grasslands in the European Alps, Scientific

562 Reports. 10 (2020) 12516. https://doi.org/10.1038/s41598-020-69569-9.

563 [23] P. Lavelle, A. Spain, M. Blouin, G.G. Brown, T. Decaëns, M. Grimaldi, J.J. Jiménez,

564 D. McKey, J. Mathieu, E. Velasquez, A. Zangerlé, Ecosystem engineers in a self-organized

565 soil: A review of concepts and future research questions, Soil Science. 18 (2016) 91–109.

566 [24] P. Ganault, J. Nahmani, S. Hattenschwiler, L.M. Gillespie, J.-F. David, L. Henneron,

567 E. Iorio, C. Mazzia, B. Muys, A. Pasquet, L.D. Prada-Salcedo, J. Wambsganss, T. Decaens,

568 Relative importance of tree species richness, tree functional type, and microenvironment for

569 soil macrofauna communities in European forests, Oecologia. 196 (2021) 455–468.

570 https://doi.org/10.1007/s00442-021-04931-w.

571 [25] S. Scheu, D. Albers, J. Alphei, R. Buryn, U. Klages, S. Migge, C. Platner, J.-A.

572 Salamon, The soil fauna community in pure and mixed stands of beech and spruce of different

573 age: trophic structure and structuring forces, Oikos. 101 (2003) 225–238.

574 [26] F. Montané, J. Romanyà, P. Rovira, P. Casals, Aboveground litter quality changes

575 may drive soil organic carbon increase after shrub encroachment into mountain grasslands,

576 Plant and Soil. 337 (2010) 151–165.

577 [27] M. Steinwandter, J. Seeber, The buffet is open: Alpine soil macro-decomposers feed

578 on a wide range of litter types in a microcosm cafeteria experiment, Soil Biology and

579 Biochemistry. 144 (2020) 107786. https://doi.org/10.1016/j.soilbio.2020.107786.

580 [28] S. Scheu, Changes in the lumbricid coenosis during secondary succession from a

581 wheat field to a beechwood on limestone, Soil Biology & Biochemistry. 24 (1992) 1641–

582 1646. https://doi.org/10.1016/0038-0717(92)90163-R.
583 [29] T. Decaëns, T. Dutoit, D. Alard, Earthworm community characteristics during afforestation of abandoned chalk grasslands (Upper Normandy, France), European Journal of Soil Biology. 33 (1997) 1–11.

584 [30] K.J. Gaston, Global patterns in biodiversity, Nature. 405 (2000) 220–227.

587 https://doi.org/10.1038/35012228.

588 [31] C. Bergmann, Über die Verhältnisse der Wärmeökonomie der Thiere zu ihrer Grösse, Göttinger Studien. 3 (1847) 595–708.

590 [32] I. Djukic, F. Zehetner, M. Tatzber, M.H. Gerzabek, Soil organic-matter stocks and characteristics along an Alpine elevation gradient, Journal of Plant Nutrition and Soil Science. 173 (2010) 30–38. https://doi.org/10.1002/jpln.200900027.

593 [33] A. Stokes, G. Angeles, F. Anthelme, E. Aranda-Delgado, I. Barois, M. Bounous, N. Cruz-Maldonado, T. Decaens, S. Fourtier, G.T. Freschet, Q. Gabriac, D. Hernandez-Caceres, L. Jimenez, J. Ma, Z. Mao, B. Eugenia Marin-Castro, L. Merino-Martín, A. Mohamed, C. Piedallu, C. Pimentel-Reyes, H. Reijnen, F. Reverchon, H. Rey, L. Selli, C. Desiree Sieber-Grabach, K. Sieron, M. Weemstra, C. Roumet, Shifts in soil and plant functional diversity along an altitudinal gradient in the French Alps, Bmc Research Notes. 14 (2021) 54. https://doi.org/10.1186/s13104-021-05468-0.

600 [34] FAO, World Reference Base for Soil Resources 2014, update 2015, International soil classification system for naming soils and creating legends for soil maps, FAO, Italy, 2015.

601 [35] P. Bénichou, O. Le Breton, Prise en compte de la topographie pour la cartographie des champs pluviométriques statistiques, La Météorologie. 7 (1987) 23–34.

602 [36] A.E. Bjune, M. Ohlson, J. Birks, R. Bradshaw, The development and local stand-scale dynamics of a Picea abies forest in southeastern Norway, The Holocene. 19 (2009) 1073–1082.

607 [37] L. Nybakken, V. Selås, M. Ohlson, Increased growth and phenolic compounds in bilberry (Vaccinium myrtillus L.) following forest clear-cutting, Scandinavian Journal of Forest Research. 28 (2013) 319–330.

610 [38] D. Hernández-Cáceres, A. Stokes, G. Angeles-Alvarez, J. Abadie, F. Anthelme, M. Bounous, G. Freschet, C. Roumet, M. Weemstra, L. Merino-Martín, F. Reverchon, Vegetation creates microenvironments that influence soil microbial activity and functional diversity along an elevation gradient, Soil Biology and Biochemistry. 165 (2022) 108485.

611 [39] M.B. Bouché, Lombriciens de France, INRA, Paris, 1972.

615 [40] P.D.N. Hebert, A. Cywinska, S.L. Ball, J.R. DeWaard, Biological identifications through DNA barcodes, Proceedings of the Royal Society B-Biological Sciences. 270 (2003) 313–321. https://doi.org/10.1098/rspb.2002.2218.

618 [41] N.V. Ivanova, J.R. Dewaard, P.D.N. Hebert, An inexpensive, automation-friendly protocol for recovering high-quality DNA, Mol. Ecol. Notes. 6 (2006) 998–1002. https://doi.org/10.1111/j.1471-8286.2006.01428.x.

621 [42] S. Ratnasingham, P.D.N. Hebert, BOLD: The Barcode of Life Data System (www.barcodinglife.org), Mol. Ecol. Notes. 7 (2007) 355–364. https://doi.org/10.1111/j.1471-8286.2007.01678.x.

624 [43] B. Richard, T. Decaëns, R. Rougerie, S.W. James, D. Porco, P.D.N. Hebert, Re-integrating earthworm juveniles into soil biodiversity studies: species identification through DNA barcoding, Molecular Ecology Resources. 10 (2010) 606–614. https://doi.org/10.1111/j.1755-0998.2009.02822.x.

628 [44] N. Bottinelli, M. Hedde, P. Jouquet, Y. Capowiez, An explicit definition of earthworm ecological categories - Marcel Bouche's triangle revisited, Geoderma. 372 (2020) 114361. https://doi.org/10.1016/j.geoderma.2020.114361.

631 [45] E. Garnier, J. Cortez, G. Billes, M.L. Navas, C. Roumet, M. Debussche, G. Laurent, A. Blanchard, D. Aubry, A. Bellmann, C. Neill, J.P. Toussaint, Plant functional markers
633 capture ecosystem properties during secondary succession, Ecology. 85 (2004) 2630-2637.

634 https://doi.org/10.1890/03-0799.

635 [46] C. Violle, M.-L. Navas, D. Vile, E. Kazakou, C. Fortunel, I. Hummel, E. Garnier, Let

636 the concept of trait be functional!, Oikos. 116 (2007) 882-892.

637 https://doi.org/10.1111/j.2007.0030-1299.15559.x.

638 [47] G. Sonnier, B. Shipley, M.-L. Navas, Plant traits, species pools and the prediction of

639 relative abundance in plant communities: a maximum entropy approach, Journal of

640 Vegetation Science. 21 (2010) 318-331. https://doi.org/10.1111/j.1654-1103.2009.01145.x.

641 [48] M. Pansu, J. Gautheyrou, Handbook of soil analysis: mineralogical, organic and

642 inorganic methods, Springer-Verlag, Berlin Heidelberg, 2006.

643 [49] S.R. Olsen, C.V. Cole, F.S. Watanabe, Estimation of available phosphorus in soil by

644 extraction with sodium bicarbonate, US Dept. Agr. Circ., 1954.

645 [50] H. Ciesielski, T. Sterckeman, M. Santerne, J. Willery, A comparison between three

646 methods for the determination of cation exchange capacity and exchangeable cations in soils,

647 Agronomie. 17 (1997) 9-16.

648 [51] C. Girardin, A. Mariotti, Analyse isotopique du 13C en abondance naturelle dans le

649 carbone organique : un système automatique avec robot préparateur, Cahiers ORSTOM -

650 Série Pédologie. 26 (1991) 371-380.

651 [52] R Core Team, R: A language and environment for statistical computing. Version 4.0.3,

652 (2020). https://www.R-project.org/.

653 [53] T.C. Hsieh, K.H. Ma, A. Chao, iNEXT: iNterpolation and EXTrapolation for species

654 diversity. R package version 2.0.19., (2019).

655 [54] A. Chao, Species richness estimation, in: N. Balakrishnan, C.B. Read, B. Vidakovic

656 (Eds.), Encyclopedia of Statistical Sciences, Wiley, New York, 2004.

657 [55] E. Marcon, Mesures de la biodiversité, UMR Écologie des forêts de Guyane, Kourou,

658 French Guiana, 2016.

659 [56] J. Oksanen, F.G. Blanchet, M. Friendly, R. Kindt, P. Legendre, D. McGlinn, P.R.

660 Minchin, R. B. O'Hara, G.L. Simpson, P. Solymos, M.H.H. Stevens, E. Szoecs, H. Wagner,

661 Vegan: Community Ecology Package. R package version 2.5-4., (2019).

662 [57] P. Legendre, Ecologically meaningful transformations for ordination of species data,

663 Oecologia. 129 (2001) 271-280. https://doi.org/10.1007/s004420100716.

664 [58] D. Chessel, A. Dufour, J. Thioulouse, The ade4 package: One-table methods, R News.

665 4 (2004) 5-10.

666 [59] J.S. Lefcheck, PIECEWISESEM: Piecewise structural equation modelling in R for

667 ecology, evolution, and systematics, Methods in Ecology and Evolution. 7 (2016) 573-579.

668 https://doi.org/10.1111/2041-210X.12512.

669 [60] C. Sagot, J.J. Brun, J.L. Grossi, J.H. Chauchat, G. Boudin, Earthworm distribution and

670 humus forms in the development of a semi-natural alpine spruce forest, European Journal of

671 Soil Biology. 35 (1999) 163-169. https://doi.org/10.1016/S1164-5563(10)70002-9.

672 [61] S. Ammer, K. Weber, C. Abs, C. Ammer, J. Prietzel, Factors influencing the

673 distribution and abundance of earthworm communities in pure and converted Scots pine

674 stands, Applied Soil Ecology. 33 (2006) 10-21. https://doi.org/10.1016/j.apsoil.2005.09.005.

675 [62] C. Salomé, C. Guenat, G. Bullinger-Weber, J.-M. Gobat, R.-C. Le Bayon, Earthworm

676 communities in alluvial forests: Influence of altitude, vegetation stages and soil parameters,

677 Pedobiologia. 54 (2011) S89-S98. https://doi.org/10.1016/j.pedobi.2011.09.012.

678 [63] A. Rozen, R.W. Myslajek, L. Sobczyk, Altitude versus vegetation as factors

679 influencing the diversity and abundance of earthworms and other soil macrofauna in montane

680 habitat (Silesian Beskid Mts, Western Carpathians), Polish Journal of Ecology. 61 (2013)

681 145-156.

682 [64] J. Seeber, G.U.H. Seeber, W. Kossler, R. Langel, S. Scheu, E. Meyer, Abundance and
683 trophic structure of macro-decomposers on alpine pastureland (Central Alps, Tyrol): effects of abandonment of pasturing, Pedobiologia. 49 (2005) 221–228.

685 https://doi.org/10.1016/j.pedobi.2004.10.003.

686 [65] M. Steinwandter, A. Rief, S. Scheu, M. Traugott, J. Seeber, Structural and functional characteristics of high alpine soil macro-invertebrate communities, European Journal of Soil Biology. 86 (2018) 72–80. https://doi.org/10.1016/j.ejsobi.2018.03.006.

689 [66] C.G. Bueno, J.J. Jiménez, Livestock grazing activities and wild boar rooting affect alpine earthworm communities in the Central Pyrenees (Spain), Applied Soil Ecology. 83 (2014) 71–78.

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

695 [68] H. De Wandeler, R. Sousa-Silva, E. Ampoorter, H. Bruelheide, M. Carnol, S.M. Dawud, G. Danila, L. Finer, S. Hattenschwiler, M. Hermy, B. Jaroszewicz, F.-X. Joly, S.

697 Mueller, M. Pollastrini, S. Ratcliffe, K. Raulund-Rasmussen, F. Selvi, F. Valladares, K. Van Meerbeek, K. Verheyen, L. Vesterdal, B. Muys, Drivers of earthworm incidence and abundance across European forests, Soil Biology & Biochemistry. 99 (2016) 167–178.

700 https://doi.org/10.1016/j.soilbio.2016.05.003.

701 [69] S. Schelfhout, J. Mertens, K. Verheyen, L. Vesterdal, L. Baeten, B. Muys, A. De Schrijver, Tree species identity shapes earthworm communities, Forests. 8 (2017) 85.

703 https://doi.org/10.3390/f8030085.

704 [70] K. Hoeffner, M. Santonja, C. Monard, L. Barbe, M. Le Moing, D. Cluzeau, Soil properties, grassland management, and landscape diversity drive the assembly of earthworm communities in temperate grasslands, Pedosphere. 31 (2021) 375–383.

707 https://doi.org/10.1016/S1002-0160(20)60020-0.

708 [71] S. Dahmouche, Début de biogéographie des vers de terre de France, Masters thesis, University of Paris VI, 2007.

710 [72] F. Popović, M. Stojanović, S. Radosavljević, T. Trakić, J. Sekulić, Earthworm community structure along altitudinal gradients on the western slopes of Kopaonik Mountain in Serbia, Turkish Journal of Zoology. 46 (2022) 103–114.

713 [73] G. González, E. García, V. Cruz, S. Borges, M. Zalamea, M.M. Rivera, Earthworm communities along an elevation gradient in Northeastern Puerto Rico, European Journal of Soil Biology. 43 (2007) S24–S32. https://doi.org/10.1016/j.ejsobi.2007.08.044.

716 [74] M.J. Steinbauer, R. Field, J.-A. Grytnes, P. Trigas, C. Ah-Peng, F. Attorre, H.J.B. Birks, P.A.V. Borges, P. Cardoso, C.-H. Chou, M. De Sanctis, M.M. de Sequeira, M.C. Duarte, R.B. Elias, J. Maria Fernandez-Palacios, R. Gabriel, R.E. Gereau, R.G. Gillespie, J.

719 Greimler, D.E.V. Harter, T.-J. Huang, S.D.H. Irl, D. Jeanmonod, A. Jentsch, A.S. Jump, C. Kueffer, S. Nogue, R. Otto, J. Price, M.M. Romeiras, D. Strasberg, T. Stuessy, J.-C.

721 Svenning, O.R. Vetaas, C. Beierkuhnlein, Topography-driven isolation, speciation and a global increase of endemism with elevation, Global Ecology and Biogeography. 25 (2016) 1097–1107. https://doi.org/10.1111/geb.12469.

724 [75] R.K. Colwell, D.C. Lees, The mid-domain effect: geometric constraints on the geography of species richness, Trends in Ecology & Evolution. 15 (2000) 70–76.

726 https://doi.org/10.1016/S0169-5347(99)01767-X.

727 [76] R.K. Colwell, C. Rahbek, N.J. Gotelli, The mid-domain effect and species richness patterns: What have we learned so far?, American Naturalist. 163 (2004) E1–E23.

729 https://doi.org/10.1086/382056.

730 [77] H.R.P. Phillips, C.A. Guerra, M.L.C. Bartz, M.J.I. Briones, G. Brown, T.W. Crowther, O. Ferlian, K.B. Gongalsky, J. van den Hoogen, J. Krebs, A. Orgiazzi, D. Routh, B. Schwarz, E.M. Bach, J. Bennett, U. Brose, T. Decaens, B. Koenig-Ries, M. Loreau, J. Mathieu, C.
733 Mulder, W.H. van der Putten, K.S. Ramirez, M.C. Rillig, D. Russell, M. Rutgers, M.P.
734 Thakur, F.T. de Vries, D.H. Wall, D.A. Wardle, M. Arai, F.O. Ayuke, G.H. Baker, R.
735 Beausejour, J.C. Bedano, K. Birkhofer, E. Blanchart, B. Blossey, T. Bolger, R.L. Bradley,
736 M.A. Callaham, Y. Capowiez, M.E. Caulfield, A. Choi, F.V. Crotty, A. Davalos, D.J. Diaz
737 Cosin, A. Dominguez, A.E. Duhour, N. van Eekeren, C. Emmerling, L.B. Falco, R.
738 Fernandez, S.J. Fonte, C. Fragoso, A.L.C. Franco, M. Fugere, A.T. Fusilero, S. Gholami, M.J.
739 Gundale, M. Gutierrez Lopez, D.K. Hackenberger, L.M. Hernandez, T. Hishi, A.R.
740 Holdsworth, M. Holmstrup, K.N. Hopfensperger, E.H. Lwanga, V. Huhta, T.T. Hurisso, B.V.
741 Iannone, M. Iordache, M. Joschko, N. Kaneko, R. Kanianska, A.M. Keith, C.A. Kelly, M.L.
742 Kernecker, J. Klaminder, A.W. Kone, Y. Kooch, S.T. Kukkonen, H. Lalthanzara, D.R.
743 Lammel, I.M. Lebedev, Y. Li, J.B. Jesus Lidon, N.K. Lincoln, S.R. Loss, R. Marichal, R.
744 Matula, J.H. Moos, G. Moreno, A. Moron-Rios, B. Muys, J. Neirynck, L. Norgrove, M.
745 Novo, V. Nuutinen, V. Nuzzo, M.P. Rahman, J. Pansu, S. Paudel, G. Peres, L. Perez-
746 Camacho, R. Pineiro, J.-F. Ponge, M.I. Rashid, S. Rebollo, J. Rodeiro-Iglesias, M.A.
747 Rodriguez, A.M. Roth, G.X. Rousseau, A. Rozen, E. Sayad, L. van Schaik, B.C.
748 Scharenbroch, M. Schirrmann, O. Schmidt, B. Schroeder, J. Seeber, M.P. Shashkov, J. Singh,
749 S.M. Smith, M. Steinwandter, J.A. Talavera, D. Trigo, J. Tsukamoto, A.W. de Valenca, S.J.
750 Vanek, I. Virto, A.A. Wackett, M.W. Warren, N.H. Wehr, J.K. Whalen, M.B. Wironen, V.
751 Wolters, I.V. Zenkova, W. Zhang, E.K. Cameron, N. Eisenhauer, Global distribution of
752 earthworm diversity, Science. 366 (2019) 480–485. https://doi.org/10.1126/science.aax4851.

753 [78] N. Bernier, J. Ponge, Humus form fodynamics during the sylvogenetic cycle in a
754 mountain spruce forest, Soil Biology & Biochemistry. 26 (1994) 183–220.
755 https://doi.org/10.1016/0038-0717(94)90161-9.

756 [79] J.L. Grossi, J.J. Brun, Effect of climate and plant succession on Lumbricid populations
757 in the French Alps, Soil Biology & Biochemistry. 29 (1997) 329–333.
758 https://doi.org/10.1016/S0038-0717(96)00234-9.

759 [80] H. De Wandeler, H. Bruelheide, S.M. Dawud, G. Danila, T. Domisch, L. Finer, M.
760 Hermy, B. Jaroszewicz, F.-X. Joly, S. Mueller, S. Ratcliffe, K. Raulund-Rasmussen, E. Rota,
761 K. Van Meerbeek, L. Vesterdal, B. Muys, Tree identity rather than tree diversity drives
762 earthworm communities in European forests, Pedobiologia. 67 (2018) 16–25.
763 https://doi.org/10.1016/j.pedobi.2018.01.003.

764 [81] F.S. Gilliam, The ecological significance of the herbaceous layer in temperate forest
765 ecosystems, Bioscience. 57 (2007) 845–858. https://doi.org/10.1641/B571007.

766 [82] A.E. Kelly, M.L. Goulden, Rapid shifts in plant distribution with recent climate
767 change, Proceedings of the National Academy of Sciences of the United States of America.
768 105 (2008) 11823–11826. https://doi.org/10.1073/pnas.0802891105.

769 [83] I.-C. Chen, J.K. Hill, R. Ohlemueller, D.B. Roy, C.D. Thomas, Rapid range shifts of
770 species associated with high levels of climate warming, Science. 333 (2011) 1024–1026.
771 https://doi.org/10.1126/science.1206432.

772 [84] E. Marris, The escalator effect, Nature Climate Change. 1 (2007) 94–96.
773 https://doi.org/doi.org/10.1038/climate.2007.70.

774
<table><thead><tr><th>Variable code</th><th>Overall interpretation</th><th>Negative scores</th><th>Positive scores</th></tr></thead><tbody><tr><td colspan="4">Earthworm communities</td></tr><tr><td>Com.1</td><td>Epigeic relative abundance</td><td>High anecic and endogeic scores</td><td>High epigeic scores</td></tr><tr><td>Com.2</td><td>L. rubellus vs L. friendi</td><td>Presence of L. rubellus</td><td>Presence of L. friendi</td></tr><tr><td>Str.1</td><td>Abundance and diversity decrease</td><td>High density, biomass and diversity, negative score on Com.1</td><td>Low density, biomass and diversity, positive score on Com.1</td></tr><tr><td>Str.2</td><td>Epigeic and anecic relative abundance</td><td>Low epigeic and anecic scores</td><td>High epigeic and anecic scores</td></tr><tr><td colspan="4">Environmental variables</td></tr><tr><td>Clim.1</td><td>Climate gradient</td><td>Higher temperature and potential evapotranspiration</td><td>Lower temperatures and higher precipitations</td></tr><tr><td>Soil.1</td><td>Organic matter gradient</td><td>High organic matter and nutrient contents, low pH</td><td>Low organic matter and nutrient contents, high pH</td></tr><tr><td>Soil.2</td><td>Eutrophication gradient</td><td>Low nitrate concentration, low pH</td><td>High nitrate concentration, high pH</td></tr><tr><td>Veg1.1</td><td>Low ground cover by trees</td><td>High tree cover</td><td>High shrub and graminoid cover</td></tr><tr><td>Veg1.2</td><td>Lack of herbaceous cover</td><td>High forb and graminoid cover</td><td>Low forb and graminoid cover</td></tr><tr><td>Veg2.1</td><td>P. abies and V. myrtillus cover</td><td>High J. communis cover</td><td>High P. abies and V. myrtillus cover</td></tr><tr><td>Veg2.2</td><td>Ground cover by shrubs other than structuring species</td><td>Low 'other shrub' cover</td><td>High 'other shrub' cover</td></tr></tbody></table>

**Table 1.** Synthetic interpretation of the axes of the different PCA used in the analyses to describe earthworm communities and environmental variables (see Fig. A.3 for complete description of the analyses).
<table><thead><tr><th>Elevational level (m)</th><th>Sobs</th><th>S.Chao</th><th>S.Chao.se</th></tr></thead><tbody><tr><td>1400</td><td>4</td><td>6.57</td><td>3.75</td></tr><tr><td>1600</td><td>5</td><td>10.46</td><td>6.39</td></tr><tr><td>1800</td><td>8</td><td>8.00</td><td>0.10</td></tr><tr><td>2000</td><td>5</td><td>5.00</td><td>0.34</td></tr><tr><td>2200</td><td>4</td><td>6.57</td><td>3.75</td></tr><tr><td>2400</td><td>4</td><td>4.00</td><td>0.51</td></tr></tbody></table>

**Table 2.** Cumulated observed (Sobs) and estimated (SChao) earthworm species richness along the elevational gradient. SChao.se = standard error of the Chao's asymptotic estimator.

782
783 **Figure caption**

784 **Figure 1.** List of the earthworm species found in Belledonne Massif (France) and their
elevational ranges. Bars indicate the minimal and maximal altitudes where each species
785 was collected, large black dots represent the mean altitude, and small grey dots the
786 individual occurrences with a jitter-randomized position allowing them to be viewed
787 without overlap. The vertical blue line represents the elevational position of the treeline.
788
789 N= number of plots where a given species has been observed.

790 **Figure 2.** Variation of earthworm species diversity along the elevational gradient: A) Box-plots showing the variation of the average species density per elevational level; B) Sample-based rarefaction and extrapolation curves of species diversity plotted for each elevational level. In A) black dots are outliers; significance codes for Kruskal-Wallis rank sum test: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘NS’. In B) solid lines represent rarefaction (interpolated) curves, whereas dashed lines represent extrapolated curves; shaded areas are 95% confidence intervals based on a bootstrap with 200 replications.
791
792
793
794
795
796

797 **Figure 3.** Variation of the community weighted mean of ecological categories scores along the elevational gradient: A) Epigeic score, B) Anecic score, C) Endogeic score. Significance codes for Kruskal-Wallis rank sum test: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘NS’.
798
799

800 **Figure 4.** Variation of biomass (A), density (B), body mass community weighted mean (C) and community weighted variance (D) along the elevational gradient. Significance codes for Kruskal-Wallis rank sum test: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘NS’.
801
802

803 **Figure 5.** Redundancy analysis (RDA) highlighting the relative importance of environmental predictors in explaining elevational changes in earthworm community structure and composition: A) ordination of sampling plots, groups by elevational levels, on the factorial plan defined by the first two axes of the RDA (i.e. RDA1 and RDA2); B)
804
805
806
807 projection of community metrics on the first factorial plan; C) projection of synthetic
808 environmental predictors on the first factorial plan; D) relative proportion of constrained
809 (Co) and unconstrained (Un) inertia; E) eigenvalue diagram. In A) significance codes for
810 the RDA axes after permutation test: 0 '**' 0.001 '**' 0.01 '*' 0.05 'ns'. In B) S= observed
811 species richness; BM= body mass; CWM= community weighted mean; CWV= community
812 weighed variance; Endo, Epi & Ane= CWM of endogeic, epigeic and anecic scores; Com.1
813 & Com.2= plot scores on the first two axes of the PCA of earthworm community
814 composition (see Table 1). In C) synthetic environmental predictors correspond to the
815 plot scores on the first two axes of the PCAs of climatic variables (Clim.1 & Clim.2), soil
816 properties (Soil.1 & Soil.2), vegetation structure (Veg1.1 & Veg1.2) and vegetation
817 composition (Veg2.1 & Veg2.2) (see Table 1); red arrows indicate significant predictor
818 variables after permutational ANOVA (p-value < 0.05).

819 **Figure 6.** Path analyses showing the direct and indirect effects of environmental
820 predictors on earthworm community structure. Str.1 and Str.2 correspond to the plot
821 scores on the first two axes of the PCA of community structure, respectively (see Fig.
822 A.4D). Significance codes for thick arrows: 0 '**' 0.001 '**' 0.01 '*' 0.05; thin arrows
823 indicate marginally significant correlations (p-value < 0.1). See Table 1 for interpretation
824 of the synthetic environmental and community synthetic descriptors.
<caption>Figure 1</caption>

Click here to access/download;Figure;Gabriac et al Fig 1.pdf

<figure><img src="image_3.png" /></figure>
<caption>Figure 2</caption>

Click here to access/download;Figure;Gabriac et al Fig 2 richness.pdf

<figure><img src="image_4.png" /></figure>
<caption>Figure 3</caption>

A

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

*

<caption>Click here to access/download; Figure; Gabrian, et al Fig 3.pdf</caption>

B

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

C

*

<figure><img src="image_7.png" /></figure>
<caption>Figure 4</caption>

Click here to access/download;Figure; Gabrac et al Fig 4 Biomass, density and BM.pdf

B

‡

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

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

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

<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /><figcaption>Figure 5</figcaption></figure>

Click here to access/downloadFigure;Gabriac et al Fig 5 rda.pdf

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

D

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

E

<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /><figcaption>Figure 6 A</figcaption></figure>

B

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