# Local climate modulates the development of soil nematode communities after glacier retreat

Alessia Guerrieri¹,² | Isabel Cantera¹ ID | Silvio Marta¹,³ ID | Aurélie Bonin² |

Alexis Carteron¹ | Roberto Ambrosini¹ | Marco Caccianiga⁴ | Fabien Anthelme⁵ |

Roberto Sergio Azzoni¹,⁶ | Peter Almond⁷ | Pablo Alviz Gazitúa⁸ |

Sophie Cauvy-Fraunié⁹ | Jorge Luis Ceballos Lievano¹⁰ | Pritam Chand¹¹ |

Milap Chand Sharma¹² | John Clague¹³ | Justiniano Alejo Cochachín Rapre¹⁴,† |

Chiara Compostella⁶ | Rolando Cruz Encarnación¹⁴ | Olivier Dangles¹⁵ | Philip Deline¹⁶ |

Andre Eger¹⁷ | Sergey Erokhin¹⁸ | Andrea Franzetti¹⁹ | Ludovic Gielly²⁰ |

Fabrizio Gili¹,²¹ | Mauro Gobbi²² | Sigmund Hågvar²³,²⁴ ID | Norine Khedim¹⁶ |

Rosa Isela Meneses²⁵,²⁶ | Gwendolyn Peyre²⁷ | Francesca Pittino¹⁹ | Angela Proietto⁶ |

Antoine Rabatel²⁸ | Nurai Urseitova¹⁸ | Yan Yang²⁹ ID | Vitalii Zaginaev³⁰ |

Andrea Zerboni⁶ | Anaïs Zimmer³¹ | Pierre Taberlet²⁰,²⁴ | Guglielmina Adele Diolaiuti¹ |

Jerome Poulenard¹⁶ | Diego Fontaneto³²,³³ | Wilfried Thuiller²⁰ |

Gentile Francesco Ficetola¹,²⁰

## Correspondence

Isabel Cantera, Dipartimento di Scienze e Politiche Ambientali, Università degli Studi di Milano, Via Celoria 10, 20133 Milano, Italy.

Email: isa_cantera@hotmail.com

## Funding information

H2020 European Research Council, Grant/Award Number: 772284; Biodiversa+; European Commission, Grant/Award Number: 101052342; MUR; ANR

## Abstract

The worldwide retreat of glaciers is causing a faster than ever increase in ice-free areas that are leading to the emergence of new ecosystems. Understanding the dynamics of these environments is critical to predicting the consequences of climate change on mountains and at high latitudes. Climatic differences between regions of the world could modulate the emergence of biodiversity and functionality after glacier retreat, yet global tests of this hypothesis are lacking. Nematodes are the most abundant soil animals, with keystone roles in ecosystem functioning, but the lack of global-scale studies limits our understanding of how the taxonomic and functional diversity of nematodes changes during the colonization of proglacial landscapes. We used environmental DNA metabarcoding to characterize nematode communities of 48 glacier forelands from five continents. We assessed how different facets of biodiversity change with the age of deglaciated terrains and tested the hypothesis that colonization patterns are different across forelands with different climatic conditions.

Alessia Guerrieri, Isabel Cantera, Silvio Marta, Aurélie Bonin, Diego Fontaneto, Wilfried Thuiller, and Gentile Francesco Ficetola contributed equally to this work.

†Died January 5, 2023.

For Affiliation refer page on 13

This is an open access article under the terms of the [Creative Commons Attribution-NonCommercial-NoDerivs License](https://creativecommons.org/licenses/by-nc-nd/4.0/), which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.
© 2023 The Authors. *Global Change Biology* published by John Wiley & Sons Ltd.
Nematodes colonized ice-free areas almost immediately. Both taxonomic and func-
tional richness quickly increased over time, but the increase in nematode diversity
was modulated by climate, so that colonization started earlier in forelands with mild
summer temperatures. Colder forelands initially hosted poor communities, but the
colonization rate then accelerated, eventually leveling biodiversity differences be-
tween climatic regimes in the long term. Immediately after glacier retreat, communi-
ties were dominated by colonizer taxa with short generation time and r-ecological
strategy but community composition shifted through time, with increased frequency
of more persister taxa with K-ecological strategy. These changes mostly occurred
through the addition of new traits instead of their replacement during succession.
The effects of local climate on nematode colonization led to heterogeneous but pre-
dictable patterns around the world that likely affect soil communities and overall eco-
system development.

KEYWORDS

climate change, colonization rates, eDNA metabarcoding, functional diversity, global scale, succession

1 | INTRODUCTION

Mountain glaciers are losing mass at an unprecedented rate due to ongoing climate change, and the pace of glacier recession is accelerating (Hugonnet et al., 2021; Rounce et al., 2023). This dramatic shrinkage will have major impacts on both the biotic and abiotic components of nearby ecosystems and the services they provide (Intergovernmental Panel on Climate Change, 2021; Pothula & Adams, 2022). Deglaciated areas undergo rapid geomorphological changes and are quickly colonized by multiple life-forms, which in turn determine the dynamics and the functioning of these emerging ecosystems (Ficetola, Marta, et al., 2021; Pothula & Adams, 2022; Sytsma et al., 2023). The loss of mountain glaciers is likely an irreversible process (Intergovernmental Panel on Climate Change, 2021); thus, we urgently need standardized, global-scale studies on how local to regional differences (e.g., differences in climate) determine the development of different biodiversity components under a changing climate (Crowther et al., 2019; Ficetola, Marta, et al., 2021). Such studies are required to identify general patterns of biotic colonization and to understand key drivers of successions in order to predict the future dynamics of these ecosystems (Cauvy-Fraunié & Dangles, 2019).

Biotic colonization after glacier retreat is generally characterized
by an increase in diversity over time, at rates that differ among taxo-
nomic groups (Guerrieri et al., 2023; Pothula & Adams, 2022; Rosero
et al., 2021). Nonetheless, it has been proposed that the colonization
dynamics of deglaciated areas are modulated by local climatic con-
ditions (Ficetola, Marta, et al., 2021), given that the impacts of cli-
mate change on biodiversity are uneven at the global scale (Freeman
et al., 2021). Climatic differences in temperature and, to a lesser
extent, precipitation and solar radiation can influence the evolution
of ecosystems after glacier retreat (Anthelme et al., 2022; Franco

et al., 2022; Kaufmann, 2002; Khedim et al., 2021). For instance, the accumulation through time of soil organic matter can be faster in proglacial forelands that experience relatively warm conditions (Kaufmann, 2002; Khedim et al., 2021). Thus, we expect that glacier forelands subjected to different climates will show differences in the development of communities, with earlier ecosystem formation and colonization in relatively warm areas (Ficetola, Marta, et al., 2021; Fridley & Wright, 2018; Vater & Matthews, 2015). However, broad-scale studies analyzing multiple proglacial landscapes located in different climates using the same approach are scarce, hampering tests of this hypothesis.

To better understand colonization dynamics due to global change, we must compare taxonomic and functional trajectories of biotic communities after glacier retreat. These biodiversity facets are complementary and influence ecosystem processes in different ways (Le Bagousse-Pinguet et al., 2019). Functional diversity encompasses the range of morphological, ecological, behavioral, and physiological traits among species within a given community. Given that species' traits are more directly connected with ecosystem processes than taxonomic diversity (Cadotte et al., 2011; Loreau et al., 2001; Naeem et al., 2012), changes in functional diversity should inform us about how changes in ecological processes might impact ecosystem functions. Analyzing functional diversity thus provides key information on the processes underlying the evolution of communities under global change, such as the relative importance of environmental filtering and competition (Kuczynski & Grenouillet, 2018).

Nematodes are the most abundant animals in soils and are excel-
lent indicators of soil biological activity (Neher, 2001). Operating at
all trophic levels, they play a crucial role in the development of eco-
systems, as their diversity affects numerous ecosystem functions,
including productivity, nutrient cycling, and decomposition of or-
ganic matter (Crowther et al., 2019; Delgado-Baquerizo et al., 2020;
<figure><img src="image_1.png" /></figure>

Fontaneto, 2019; van den Hoogen et al., 2019; Zawierucha et al., 2021). Nematodes are particularly important components of soil communities in the first decades after glacier retreat (Devetter et al., 2021; Doblas-Miranda et al., 2008; Rosero et al., 2021), probably due to their small size and resistant propagules that facilitate dispersal and colonization (Fontaneto, 2019; Ptatscheck et al., 2018). Previous studies on nematode colonization after glacier retreat reported a range of successional patterns, with variation among trophic groups (Pothula & Adams, 2022) and study areas. For instance, in Lei et al. (2015), most of richness increase occurred during the first 50 years, while colonization required much longer time in Devetter et al. (2021). Furthermore, some studies showed linear increases in diversity over time (Devetter et al., 2021), while others suggested intermediate peaks followed by plateaus or even richness declines (Doblas-Miranda et al., 2008; Lei et al., 2015; van Leeuwen et al., 2018). Nonetheless, those studies were conducted in distinct regions of the world and covered different temporal ranges, making comparisons and generalizations difficult. We thus need a global assessment using a standardized protocol to examine the successional patterns of nematode diversity and its trophic groups. This will help to define how climate modulates nematode succession.

Here, we analyzed the development of taxonomic and functional diversity of nematode communities emerging since the retreat of glaciers in five continents (Figure 1). Molecular approaches, such as environmental DNA (eDNA) metabarcoding (Taberlet et al., 2018), enable rapid, cost-effective, and consistent estimates of soil biodiversity (Delgado-Baquerizo et al., 2020; Geisen et al., 2018) even over broad spatial scales (White et al., 2020), which would otherwise be prohibitive. We used a standardized eDNA-based approach to sample soil nematode communities along the proglacial chronosequences (Walker et al., 2010) of 48 receding glaciers that experience different climatic conditions (Figure 1). First, we tested the hypothesis that local climatic conditions (i.e., meso-climate) influence the rate at which communities develop after the retreat of glaciers. Under this hypothesis, we expected that communities would develop more quickly in mild climates compared to the cold ones. We tested this hypothesis by measuring how different facets of biodiversity (i.e., taxonomic and functional diversity, as well as the richness of trophic groups) changed over time and across forelands with different climates. Second, we tested whether nematode communities show functional changes along the colonization. The colonization of deglaciated terrains is often led by colonizer species with r-ecological strategy, that may be replaced over time by more competitive species with K-ecological strategy (Erschbamer & Caccianiga, 2016; Ficetola, Marta, et al., 2021). We thus predicted a shift from r-strategy to K-strategy species. Third, we sought to unravel the coexistence mechanisms that predominantly affected the functional assembly of communities (Gotelli & Graves, 1996). Under environmental filtering, communities are expected to be more functionally homogeneous, with taxa harboring traits that allow them to tolerate specific conditions. This coexistence mechanism would lead to a functional diversity that is lower than expected based on observed taxonomic diversity. Conversely, under a limiting similarity assembly process, colonization should be primarily underlined by a few groups

that are functionally different (e.g., that use different resources) (Götzenberger et al., 2012; Weiher & Keddy, 1999). We predicted stronger environmental filtering in harsh conditions (very cold climates and/or soon after glacier retreat), whereas limiting similarity could be more important under mild conditions and/or in more developed communities, where competitive interactions are stronger (Purschke et al., 2013).

## 2 | MATERIALS AND METHODS

### 2.1 | Field activities

From 2014 to 2020, we collected 1197 soil samples in 48 glacier forelands, covering five continents from the Equator to polar regions (Figure 1). For these forelands, information on the dates of glacier retreat since the maximum extent of the Little Ice Age was available from the literature, remote sensing images and field surveys (Marta et al., 2021). For each foreland, we selected two to nine sites (mean = 5.2, SD = 1.5, Table S1) along the chronosequences of glacier retreat (that is the chronological sequence of specific geomorphological features, such as moraines or inter-moraine flat areas, along deglaciated terrains for which the date of glacier retreat is known; Figure S1). A site generally corresponded to the line representing the position of the glacier's front on a given date. We excluded sites impacted by recent surface instability and erosion. The number of sites per chronosequence depended on the number of known positions of the glacier available from the literature (Marta et al., 2021). At each site, we established an average of five regularly spaced plots (SD = 0.6) for more repeatable biodiversity estimates (Dickie et al., 2018), and, if possible, at distances of 20 m (Figure S1). Due to logistic constraints, in 9% of the sites, we sampled a slightly different number of plots. The average number of plots per site was similar between forelands in different climatic regimes and between sites with different ages since glacier retreat (linear models testing relationships between the number of plots and age or climate, both *p* ≥ .8). At each plot, we collected five soil subsamples within 1 m (Figure S1), at a depth of 0–20 cm and pooled them together to form a composite sample of ~200 g per plot. We did not include soil litter and avoided roots, leaves and other large plant organs, as we focused on soil nematodes. While the analysis of litter nematodes can provide interesting patterns (Doblas-Miranda et al., 2008), many sites completely lacked leaf litter, particularly at early successional stages, thus hampering a standardized comparison. Composite samples were homogenized; from each composite sample, we took 15 g of soil and desiccated it immediately in sterile boxes with 40 g of silica gel (Guerrieri et al., 2021). Desiccation is a cost-effective approach allowing long-term preservation of eDNA for biodiversity assessments (Guerrieri et al., 2021). Before the collection of each sample, all the sampling equipment underwent strict decontamination protocols (burned at >1000°C with a portable blow torch). Using the same approach, we also collected 150 g of soil for chemical analyses at each plot in a subset of forelands. Sampling was always
<figure><img src="image_2.png" /><figcaption>**FIGURE 1** Global distribution of the 48 studied glacier forelands. The color of the dots represents the average temperature of forelands during the warmest quarter (ranging from 2 to 17°C). The background blue dots represent the distribution of glaciers according to the GLIMS database (www.glims.org). Map lines delineate study areas and do not necessarily depict accepted national boundaries.</figcaption></figure>

performed during the warmest season (i.e., for temperate areas from late July till early September in the Northern Hemisphere and in February–March in the Southern Hemisphere; Table S1).

## 2.2 | Molecular analyses and bioinformatics treatment

Environmental DNA from the 1197 15-g soil samples was extracted in a dedicated room using the NucleoSpin Soil Mini Kit

(Macherey-Nagel). Because we collected ~five soil samples per site, the total amount of soil used for extraction was on average 75 g per site, which is close to the values proposed for robust assessments of nematode communities (Wiesel et al., 2015). Before extraction, we added a preliminary step where the soil was mixed with 20 mL of phosphate buffer for 15 min (Taberlet et al., 2012). Elution was performed in 150 µL of elution buffer. To monitor contamination in the extraction room, we also included one extraction control every ~10 samples (total: 101 extraction controls; Zinger et al., 2019). DNA amplifications were performed using the universal Euka02
<figure><img src="image_3.png" /></figure>

primer (Guardiola et al., 2015; Taberlet et al., 2018), which amplifies an ~120bp fragment of the V7 region of the 18S rDNA gene. The short fragment length is well suited for eDNA analyses (Ficetola, Boyer, et al., 2021; Guardiola et al., 2015; Taberlet et al., 2018). Bioinformatics analyses showed that Euka02 is among the markers amplifying the highest proportion of nematode species and has good taxonomic resolution, particularly at the genus level (Ficetola et al., 2023). Reverse and forward primers included 8-nucleotide-long tags on the 5' end. Each tag had at least five nucleotide differences from the others, thus allowing bioinformatic discrimination of PCR replicates after sequencing (Taberlet et al., 2018). DNA extracts were randomized in 96-well plates together with extraction controls, bioinformatic blanks, PCR negative and positive controls (total across all plates: 291 blanks, 90 negative and 53 positive controls). Positive controls consisted of 10 times diluted genomic DNA of eight bacterial and two fungal strains (ZymoBIOMICS Microbial Community DNA Standard II; Zymo Research) and were used to monitor potential cross-contamination of samples, as well as amplification and sequencing performances (Zinger et al., 2019). PCR-negative controls were used to monitor contaminations during PCR reactions and blanks allowed the identification of tag-jump issues. We used quantitative PCR (qPCR) to determine the optimal number of PCR cycles, following Guerrieri et al. (2023).

PCR amplifications of 1:10 diluted DNA samples were performed in 384-well plates and consisted of an initial step of 10 min at 95°C, followed by 45 cycles of 30 s denaturation at 95°C, 30 s annealing at 45°C, 60 s elongation at 72°C, and 7 min final elongation at 72°C. All samples and controls underwent four PCR replicates (Ficetola et al., 2015). The purified amplicons were sent to Fasteris (Switzerland) for library preparation using the MetaFast protocol (Taberlet et al., 2018) and sequencing with the Illumina HiSeq 2500 platform (paired-end approach, 2x150 bp; about 0.5% of PhiX added in each sequencing lane). See Guerrieri et al. (2023) for additional details on molecular analyses.

We performed the bioinformatic treatment of sequence data using the OBITools software suite version 1.2.9 (Boyer et al., 2016), following the steps described in Guerrieri et al. (2023), and then clustered sequences at a threshold of 96% sequence similarity using the sumaclust program (https://git.metabarcoding.org/obitools/sumaclust/wikis/home). The 96% threshold minimizes the risk that sequences belonging to the same species are clustered in different molecular operational taxonomic units (MOTUs) (Bonin et al., 2023). Finally, we performed a taxonomic assignment of cluster heads based on the EMBL reference database (version 140). The reference database was curated by carrying out an in silico PCR with the ecopcr program (Ficetola et al., 2010). We used the same primers as in the experiment, allowing three mismatches per primer and kept sequences assigned at least at the family level. Next, we assigned the detected sequences to MOTUs using the ecotag program, following the procedure described in Boyer et al. (2016). ecotag matches each sequence in the dataset against the reference database and then uses the lowest common ancestor algorithm to identify the taxonomic level of the assignation (e.g., genus, family, order) (Boyer et al., 2016). We then performed

an additional filtering in R (version 4.0) to remove contaminants and tag-jump errors based on sequences detected in controls and blanks (Guerrieri et al., 2023; Zinger et al., 2019). Specifically, we discarded MOTUs with best identity <80%, detected less than 10 times in the overall dataset and in at least one extraction or PCR negative control, as they are possible contaminants (Ficetola et al., 2015; Zinger et al., 2019), and MOTUs detected in only one sample, as they represent singletons. Alternative approaches for filtering contaminants compare the prevalence of each MOTU in "true" samples to their prevalence in negative controls, and identify the contaminants based on a threshold of occurrence on the negative controls (Davis et al., 2018). When we repeated the MOTU filtering using this approach with thresholds of 0.1 and 0.5, we removed 4 and 19 MOTUs, respectively. Estimates of taxonomic richness obtained with this filtering strategy were highly correlated (for both thresholds: r > .86, p < .0001) with the estimates obtained with our approach (22 contaminant MOTUs) and yielded models with nearly identical results (Figure S2).

## 2.3 | Assessing taxonomic and functional diversity

Due to the difficulty of relating the number of reads to the relative abundance of species in metabarcoding studies, we first measured the taxonomic diversity of PCR replicates as the total number of taxonomic units detected. Taxonomic data obtained from eDNA can be combined with available databases on functional traits to obtain estimates of functional diversity (Aglieri et al., 2020; White et al., 2020). We extracted functional traits from the NEMAPLEX database (http://nemaplex.ucdavis.edu/). We considered one life-history index (colonizer-persister index, hereafter "CP value") and eight consumption traits (parasites of animals, parasites of plants, bacterivores, fungi-vores, detritivores, omnivores, predators and feeding on unicellular eukaryotes). CP values ranged from 1 to 5, with low values representing colonizer species with short generation time and r-reproductive strategy (i.e., fast reproduction rate and release of numerous small eggs) and high values representing persister species with long generation time and K-strategy (Bongers, 1990; Ferris, 2010; Ferris et al., 2001). The NEMAPLEX database reports functional traits at genus-level resolution. Thus, we directly assigned CP values and consumption traits only for MOTUs identified at the genus level, whereas for MOTUs identified at the family level, CP values were calculated as the average of the genera within the family; consumption traits are conserved within family. The dataset resulted in 76 taxa that had sufficient taxonomic resolution (family level or better) to assign traits (Table S2). Taxa identified at coarser levels were excluded.

The number of detected taxa generally underestimates the actual diversity of a given site because some taxa often remain undetected. Therefore, several approaches have been developed to estimate the actual number of taxa based on the replicated sampling of a given site (Colwell & Coddington, 1994; Gotelli & Colwell, 2011). At each site, we estimated taxonomic richness with the first-order jackknife estimator for incidence data, which is among the best-performing approaches for evaluating
the completeness of biodiversity inventories. Analyses of datasets of completely surveyed areas and simulations suggest that this approach can provide robust estimates of the actual species richness, if a sufficiently large number of replicated samples are analyzed (Chazdon et al., 1998; Colwell & Coddington, 1994; Gotelli & Colwell, 2011). To extrapolate biodiversity values at each site, we considered the different plots and the four PCRs for each plot as replicated samples. On average, we had 20 replicates per site (~5 plots with four PCR replicates each); this replication level is appropriate for biodiversity estimation in environments with moderate diversity levels (Gotelli & Colwell, 2011). The use of estimated biodiversity allowed us to obtain diversity estimates that are comparable between sites with a different number of sampled plots. We estimated site-level taxonomic diversity using the *alpha. accum* function of the BAT R package (Cardoso et al., 2015).

Hill numbers are an alternative approach to biodiversity measurement and are particularly useful in metabarcoding studies, as they avoid issues related to the excessive weighting of rare MOTUs and are extremely robust to different filtering strategies. Hill numbers >0 take into account differences in read abundance that might be related to differences in biomass (Taberlet et al., 2018). We, therefore, repeated the analyses of taxonomic diversity using Hill numbers ($q=0$ and $q=1$, equivalent to richness and Shannon entropy, respectively) to estimate taxonomic diversity; increasing $q$ values indicates a higher weight for the MOTUs with proportionally more reads. We could not perform this analysis with Hill numbers >1 because at several sites, we detected zero nematode MOTUs. The relationships between taxonomic diversity estimated using Hill numbers, time since glacier retreat, and climatic variables (see below for details about model structure) were identical to those observed for taxonomic diversity estimated with the first-order jackknife estimators (Figure S3).

Functional diversity describes the variation of ecological traits in a community and how traits are supported by the co-occurring taxa (Mouillot et al., 2013). The assessment of functional diversity requires the analysis of multiple metrics because the functional attributes can impact differently ecosystem dynamics (Grenié et al., 2017; Le Bagousse-Pinguet et al., 2019). We assessed functional diversity using functional richness and Rao's functional diversity (Mason et al., 2005; Petchey & Gaston, 2006), which provide complementary information. We did not use approaches based on measures of the multidimensional trait space (e.g., N-dimensional hypervolumes), as they cannot be applied to communities with less than three functionally diverse taxa, and we detected <3 functionally diverse MOTUs at 61% of sites. Functional richness represents the range of traits supported by species in a community (Mason et al., 2005; Petchey & Gaston, 2002) and was assessed using the adaptation of Faith's phylogenetic diversity developed in Petchey and Gaston (2002). We first estimated pairwise trait dissimilarities with Gower distances, with equal weights for each trait (Laliberté & Legendre, 2010; Petchey & Gaston, 2002). We then built an UPGMA dendrogram using functional distances between pairs of functionally characterized taxa. Functional richness was estimated using the first-order jackknife estimator with the *alpha.accum* function (Cardoso et al., 2015) in the same way as taxonomic richness.

Rao's functional diversity (i.e., Rao's index of quadratic entropy) integrates trait dissimilarity between taxa and the functional evenness of communities by considering the relative abundance of organisms (Mouchet et al., 2010; Petchey & Gaston, 2006). Rao's functional diversity was assessed based on the relative abundance of taxa in sites (proportion of reads) and trait dissimilarities between taxa, following de Bello et al. (2009). Using the proportion of reads can yield diversity measures that are less sensitive to parameters that are used for data filtering (Calderón-Sanou et al., 2019; Mächler et al., 2021). In a preliminary analysis, we also assessed phylogenetic diversity (see Supporting Information). However, estimates of phylogenetic diversity showed an almost perfect correlation with taxonomic richness ($r=.97$). The relationships between phylogenetic diversity, time since glacier retreat, and climatic variables were identical to that observed for taxonomic richness (Figure S4).

## 2.4 | Drivers of changes in taxonomic diversity and traits

To assess whether colonization patterns are different between forelands in different climatic regimes, we defined the meso-climate of each foreland as the average conditions of the sites within the foreland during the warmest quarter of the year. We calculated mean temperature and precipitation for the period 1979–2013 based on ChelsaClim dataset at the 1-km spatial scale (Karger et al., 2017). ChelsaClim is among the datasets that best describe climatic conditions in mountainous areas (Karger et al., 2017) and provided estimates of temperature that were strongly related to field records (correlation with temperature recorded with data-loggers in the soil of 26 glacier forelands: Pearson's $r=.82$, $p<.0001$) (Marta et al., 2023). Temperatures ranged from 2 to 17°C (average 7.8°C). We used climatic conditions of the warmest quarter (hereafter: summer) because, during winter, snow cover causes a strong decoupling between air temperature and soil temperature (see e.g., Lembrechts et al., 2022; Marta et al., 2023). Therefore, during winter, meso-climate does not represent conditions experienced by soil organisms whereas, in summer, it provides an excellent estimate of soil conditions (Marta et al., 2023). In fact, previous analyses showed that summer temperature is the strongest climatic driver of soil development in proglacial environments, probably because most biological activity occurs in this season (Khedim et al., 2021).

For a subset of sites ($N=160$), we also collected additional soil samples to measure pH and total organic carbon (TOC). pH was measured with a pH-meter after having dispersed 4 g of soil into 10 mL of bi-distilled water (Gibbard, 1993). TOC was measured on bulk samples after drying, using an OEA Flash200 ThermoFisher elemental analyzer (see Khedim et al., 2021). We calculated average organic carbon and pH for each site. Sites with and without chemical features had similar ages since glacier retreat and mean temperatures (linear models; both $p>.3$).

We first used univariate Bayesian generalized linear mixed models (GLMMs) to assess how time since glacier retreat and meso-climatic conditions determine nematode diversity. We ran GLMMs with the
<figure><img src="image_4.png" /></figure>

different site-level biodiversity metrics as dependent variables (taxonomic and functional richness, Rao's functional diversity and the mean CP value of each community, calculated by averaging the CP values of all the detected MOTU in each site). As independent variables, we considered time since glacier retreat (log-transformed), and mean temperature and total precipitation (log-transformed) during the warmest season. Mixed models also included the interactions between climatic variables and time since glacier retreat. These interactions were intended to test the hypothesis that the colonization rate is different between forelands experiencing different climatic conditions. We used the widely applicable information criterion to compare models with and without interactions (Gelman et al., 2014). We included the identity of glacier forelands as a random factor. Time and climatic variables were scaled (mean=0, SD=1) to allow a comparison of their estimated effects. Models were run with three MCMC chains using 10,000 iterations and a burn-in of 5000 in the *brms* R package (Bürkner, 2017). For taxonomic and functional richness, we used square root transformed data and a normal error distribution to better meet the normality assumption of models. Rao's functional diversity was rescaled following Smithson and Verkuilen (2006) to avoid zeros and ones, and the model was run using Beta distribution. For all models, c-hat was <1.01, indicating convergence; if the 95% credible intervals of an independent variable did not overlap zero, we assumed strong evidence of a relationship between the predictor and the response variable of the model. For taxonomic and functional richness, we included all the communities (N=240), while for Rao's functional diversity and the mean CP value, we excluded communities with observed richness=0 and thus analyzed 177 communities. To confirm that the observed relationships were not affected by unaccounted biogeographical processes, we assessed the relationship between the residuals of models and latitude (Bahram et al., 2018). Furthermore, we used Moran's *I* to confirm that the residuals of our models did not show patterns of spatial autocorrelation (Figure S5). Finally, several soil characteristics, such as nutrient content, can change along chronosequences and influence soil organisms (Ficetola, Marta, et al., 2021). To confirm that our conclusions were not biased by differences in soil features, we repeated the analyses for the subset of sites for which organic carbon and pH were available (N=160).

Threshold dynamics have been suggested to underline succession of glacier forelands, with a fast increase in community richness during the first years, followed by a plateau when high richness determines strong biotic interactions (Hanusch et al., 2022). We thus used segmented regressions to check the occurrence of thresholds in the relationships between taxonomic richness and time since glacier retreat (Ficetola & Denoël, 2009; Muggeo, 2003). We used maximum likelihood to build linear mixed models with one breakpoint (segmented package in R; Muggeo, 2003). We then compared models with the breakpoint with linear mixed models based on the Bayesian information criterion. As we hypothesized that climate modulates the relationships between time since glacier retreat and diversity, we analyzed separately glaciers with summer temperature above and below the mean (7.8°C being the mean summer temperature across all the study forelands).

We next built a multivariate GLMM to assess the impact of time and climatic variables on the richness of MOTUs belonging to the different trophic groups, defined on the basis of MOTU consumption traits. Multivariate GLMMs account for the multivariate collinearity between variables while providing a better statistical power than univariate tests. For each trophic group, the MOTU richness at the site level was used in the model as a dependent variable; we used a negative binomial distribution and considered all communities with at least one detected MOTU (N=177). We used the same independent variables and random effects of univariate GLMMs.

To further summarize the variation of all the considered functional traits tested by our models, we visualized trait and community locations in the functional space (the functional identity of communities; Mouillot et al., 2013). We built a multivariate functional space by ordinating taxa based on Gower trait distances using principal coordinates analysis (PCoA) and illustrated the global functional space of nematode communities with the first two axes of the PCoA (Borcard et al., 2011; Oksanen et al., 2013). This allowed us to ordinate communities according to their functional structure, by calculating the weighted mean of the PCoA scores of all the taxa detected in each community. We then used the *envifit* function from the vegan package to fit traits, time, and climatic variables onto the PCoA. The determination coefficient (R²) was used to assess the strength of the correlations between the variables and the axes, and its significance was determined by comparing observed values with R² values obtained from 999 random permutations.

In order to evaluate the relative importance of environmental filtering (which predicts functionally homogeneous communities) versus competition (which predicts functionally different groups within the same community) (Götzenberger et al., 2012; Weiher & Keddy, 1999), we compared the observed values of Rao's functional diversity to those obtained from random community assemblies (Gotelli & Graves, 1996). We used null models to exclude the effects of taxonomic richness on functional diversity (Aros-Mualin et al., 2021) for all sites with at least two taxa and Rao's functional diversity >0 (N=142). We compared observed values of Rao's diversity to those obtained by randomizing the matrix of the 76 functionally characterized taxa 999 times, while keeping the number of detected taxa per community fixed. We then calculated standardized effect size (SES) values as the difference between the estimated values and the mean of randomly generated values of Rao's functional diversity, divided by the standard deviation of the 999 null values. Negative values of SES indicate that functional diversity is lower than expected by chance given the observed taxonomic richness, whereas positive values indicate that functional diversity is higher than expected under random assembly (Gotelli & McCabe, 2002). We then used Bayesian GLMMs to test whether SES values are related to time since glacier retreat and climatic conditions.

Finally, to characterize the processes underlying functional changes, we determined if functional dissimilarity between communities is caused by the addition of new traits or by the replacement of existing traits. We partitioned the functional β-diversity between sites within the same glacier foreland using the *beta* function in the package BAT (Cardoso et al., 2015). This function partitions total
functional β-diversity (β-Total) into the functional dissimilarity explained only by the replacement of traits (β-Replacement) and the dissimilarity explained by the loss/gain of new traits (β-Richness). We calculated functional beta diversity based on an UPGMA dendrogram built by clustering Gower functional distances between functionally characterized taxa. We only compared sites having at least one taxon (N=177) and the comparisons were performed only between sites belonging to the same chronosequence (321 comparisons). All statistical analyses were performed in the R environment (version 4.0; www.r-project-org).

## 3 | RESULTS

After data filtering and clustering, DNA metabarcoding yielded a total of 3969 MOTUs of eukaryotes, 171 of which were identified as nematodes (Table S3; see Table S4 for the number of sequences kept at each step in the analyses). The number of nematode MOTUs detected per site ranged from 0 to 22 and was unrelated to the total number of DNA reads (Spearman's correlation: r_s = .08, p > .05), suggesting that nematode richness was not affected by sequencing performance.

Comparisons between observed taxonomic diversity (number of MOTUs) and biodiversity estimated using the first-order jackknife estimator (after excluding sites with zero detections) suggested that detected MOTUs represented on average 71% of the actual taxonomic richness (range of coverage: 56%-100%). The coverage was also high for functional richness (on average, 83% of actual diversity values were detected). For both taxonomic and functional richness, detected and estimated diversity values were strongly correlated (in all cases, Pearson's r ≥ .91; Figure S6). Functional richness and Rao's functional diversity were positively correlated and correlated also with taxonomic and phylogenetic diversity (all r ≥ .4 and p < .01; Table S5).

### 3.1 | Patterns of biotic colonization under different climates

The taxonomic richness of nematodes strongly increased with time since glacier retreat and was significantly higher in forelands with mild summers (i.e., with above-the-average summer temperature: >7.8°C; Figure 2a,b; Table S6). Furthermore, we detected a strong interaction between time and temperature. In mild forelands, nematodes were able to colonize the terrains in a very short time, being already present 4 years after glacier retreat. Conversely, in forelands with cold summers (temperature ≤7.8°C), taxonomic richness was close to zero at all sites aged less than 15 years (Figure 2a). Afterward, taxonomic richness increased over time, but the increase was relatively fast in cold forelands, so that the diversity became similar across all forelands ~150 years after glacier retreat. Segmented regressions did not indicate thresholds in the increase in taxonomic richness either in warm or cold summer environments (Table S7). All results were robust to the inclusion of soil features (pH and TOC) in the models (Table S8) or the use of Shannon diversity instead of taxonomic richness (Figure S3).

Among the 76 functionally characterized taxa, 57% occurred in both early (≤65 years) and late (>65 years) communities, whereas 26 taxa were exclusively observed in late communities, and seven taxa were only found in early communities (Figure S7). Among those pioneer taxa, five were only detected in mild forelands (Hemicyclophora_01755, Neotylenchidae_02938, Oostenbrinkia_01833, Rhyssocolpus_03988, and Tobrilidae_03778), one in both mild and cold forelands (Mononchus_02143), and one (Mesodorylaimus_04367) was exclusively detected in cold forelands. Overall, most taxa (61%) were found in both cold and mild forelands, although more taxa were exclusive to mild forelands (21) than to cold forelands (9; Figure S8). Similarly, 58% of the taxa was detected in both rainy and dry forelands, with a similar number of taxa exclusively found in rainy (17) or dry (15) forelands (Figure S9).

Both functional richness and Rao's functional diversity increased over time (Figure 2c,d; Table S6). Again, the functional richness of soils exposed for <150 years was consistently higher in mild forelands. A strong interaction between time and summer temperature indicated that the effect of climate on functional richness is particularly strong immediately after glacier retreat, similar to what is observed for taxonomic richness. In contrast, this interaction was not significant for Rao's functional diversity (Table S6). All measures of taxonomic and functional diversity were unrelated to precipitation differences across forelands (Table S6).

We also evaluated whether taxonomic richness consistently increases within the different functional groups of nematodes, defined based on their consumption traits. Most consumption traits were already present soon after glacier retreat (i.e., <15 years; Figure 3). The richness of most of them (bacterivores, fungivores, omnivores, predators, detritivores, and parasites of plants) increased through time (Figure 3; Table S9). Bacterivores and detritivores showed higher richness in forelands with mild summers, and the richness of predators decreased with precipitation (Table S9). In this case, we did not find strong support for an interaction between climatic features and time since glacier retreat (Table S9).

### 3.2 | Changes in functional traits

Communities changed their mean CP values across deglaciated terrains of different ages (Figure 2e; Table S6). Soon after glacier retreat, communities were dominated by colonizer taxa with r-ecological strategy and short generation time (i.e., low CP values). At later stages, communities were dominated by more persister taxa with K-strategy and long generation times. For this trait, we found neither relationships with climate nor an interaction between climate and time (Table S6). Fitting communities in the functional space using PCoA (Figure 4) confirmed that, immediately after glacier retreat, communities were dominated by colonizer taxa, mostly feeding on bacteria and fungi or plant parasites. In later stages, communities were dominated by persister taxa and hosted a greater diversity of trophic groups, with more omnivores and predators and, to a lesser extent, detritivores and animal parasites (Figure 4; Table S10).
<figure><img src="image_5.png" /></figure>

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

FIGURE 2 Changes in taxonomic and functional diversity of nematodes with age of deglaciated terrains. (a) Taxonomic richness (first-order jackknife estimator). Blue boxes: forelands with mean temperature during the warmest quarter (summer temperature) below the average (7.8°C); yellow boxes: forelands with summer temperature above the average. (b–e) Relationships between different measures of taxonomic and functional diversity and time since glacier retreat in forelands with different summer temperature conditions: 4.3 and 11.3°C are the average ± SD. (b) Taxonomic richness. (c) Faith's functional richness. (d) Rao's functional diversity. (e) Mean colonizer–persister (CP) value of the communities, ranging from 1 (colonizer taxon with r-strategy) to 5 (persister taxon with K-strategy). (b–e) Represent conditional plots of Bayesian mixed models; partial residuals were calculated using fixed effects only; shaded areas are 95% credible intervals for forelands with different summer temperatures. For models in (b–e), conditional $R^2$: .52, .48, .36, .30; marginal $R^2$: .23, .22, .14, .07. In panels (a–c), N=240 communities. In (d, e), N=177 communities with observed richness >0.

## 3.3 | Processes underlying community assembly and functional changes

At 60% of sites, the observed functional diversity was similar to what we would expect based on taxonomic diversity (non-significant Rao's SES values; |SES|<1.75), suggesting that the accumulation of traits over time was mostly caused by the addition of new taxa. Nevertheless, SES values were negatively related to precipitation (Figure 5; Table S6), with consistently positive SES in forelands with less precipitation during summer (on average, below 292.9 mm). This suggests that, in these forelands, communities were more functionally diverse than in more humid forelands. Additionally, SES values tended to increase over time (but this effect was weak with credible intervals slightly overlapping zero), with the oldest communities having more positive SES values regardless of climate (Table S6;

Figure 5). This indicates that early communities tended to be more functionally homogeneous than late communities.

Partitioning functional dissimilarity between communities from the same foreland showed that functional changes are characterized by both the addition of new traits and the replacement of traits (functional β-Replacement=46.5% and functional β-Richness=53.5%), although the contribution of the addition of traits was significantly higher than trait replacement (Figure S10, Wilcox paired test: V=20,019, p=.012).

## 4 | DISCUSSION

Understanding how biodiversity components develop under different climates and assessing the variation of functional groups in response
<figure><img src="image_7.png" /><figcaption>FIGURE 3 Changes in the richness (number of MOTUs) of trophic groups of nematodes over time in deglaciated terrains. The blue line represents the fitted values of a negative-binomial Bayesian mixed multivariate model; shaded areas are 95% credible intervals; N = 240 communities. In (a–h), conditional R²: .5, .4, .2, .1, .2, .3, .2, .2; marginal R²: .1, .04, .06, .02, .1, .1, .05, .01. MOTU, molecular operational taxonomic unit.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 4 Results of the principal coordinate analysis (PCoA) showing changes in functional traits of nematode communities with nematode colonization. (a) Representation of detected taxa displayed according to their consumption traits within the global functional space (explained variance: 68% and 20% for axes 1 and 2, respectively). (b) Representation of communities and their trait composition based on the same PCoA. In (a, b), N = 177 communities with observed richness >0. MOTU, molecular operational taxonomic unit.</figcaption></figure>

to glacier retreat is pivotal for predicting the future of these emerging ecosystems. After glacier retreat, soil colonization by keystone animals, such as nematodes, was affected by temperature differences

among forelands (Figure 2), supporting the hypothesis that the retreat of glaciers impacts soil biodiversity in different ways depending on regional climates (Kaufmann, 2002; Khedim et al., 2021).
<figure><img src="image_9.png" /><figcaption>**FIGURE 5** Relationship between standard effect sizes of Rao's functional diversity and time since glacier retreat in forelands with different precipitation conditions: 292.9 mm is the average summer precipitation across all the forelands (gray); 539.2 (purple) and 159.2 mm (orange) are the average ± SD. The SES sign indicates if the functional diversity is lower (SES < 0) or higher (SES > 0) than expected by chance, given the observed taxonomic diversity. The dashed line indicates SES = 0. Conditional *R*²: .22; marginal *R*²: .08. *N* = 142 communities with Rao's functional diversity > 0. SES, standardized effect size.</figcaption></figure>

**4.1 | Mesoclimate affects the development of nematode communities**

Nematodes were able to colonize deglaciated terrains very soon (in 4–15 years; see also Lei et al., 2015). As nematodes are generally absent from supra-glacial environments (Zawierucha et al., 2021), their rapid colonization might be related to their excellent dispersal abilities (Fontaneto, 2019). Colonization started earlier in forelands with mild temperatures, where we found nematodes already 4 years after deglaciation. Under colder temperatures, 10–15 years were required for nematode communities to establish (Figure 2a), possibly because soil ecosystems were not developed enough. Despite these differences at the early stages, the successional trajectories observed in different climatic regimes converged and, after approximately 150 years, the coldest forelands achieved diversity values comparable to the warmest ones (Figure 2). Temporal patterns of nematode colonization in proglacial environments vary among studies: while some reported rapid colonization (within 3 years of glacier retreat; Lei et al., 2015), others suggested that nematode establishment may take one to several decades (Devetter et al., 2021; van Leeuwen et al., 2018). These studies analyzed different areas of the world; slow colonization was observed in subarctic regions (Devetter et al., 2021; van Leeuwen et al., 2018), whereas fast colonization was observed in a warmer region (Lei et al., 2015). Consistent with these findings, our results suggest that such disparities among studies can be attributed to key environmental features, such as climate.

The effect of climate was evident also on functional diversity, as temperature interacted with time to determine functional richness, and SES values were significantly related to precipitation. After accounting for taxa richness, functional diversity tended to be higher than expected under randomness (SES > 0) in late communities and

dry forelands. The effect of time on SES values was weak but followed the expected pattern (Purschke et al., 2013). Communities exhibiting high functional diversity often indicate limiting similarity due to niche differentiation (Macarthur & Levins, 1967). This can be underlined by competition for resource-rich conditions (Bertness & Callaway, 1994; Losapio et al., 2021; Pothula & Adams, 2022), as expected in late communities (Ficetola, Marta, et al., 2021). Moreover, the higher than expected functional diversity in forelands with low precipitation might arise because only a few taxa can colonize semi-arid environments, but these taxa included nematodes exhibiting diverse trophic specializations. The response of nematode communities to precipitation gradients is extremely complex and probably different across ecosystems, highlighting the importance of more studies relating nematode diversity and water availability (Franco et al., 2019; Ma et al., 2018; van den Hoogen et al., 2019). The lower than expected functional diversity obtained in early communities and rainy forelands may result from environmental filtering retaining taxa that share traits adapted to those specific conditions. However, competitive exclusion of taxa with traits associated with low competitive ability for the limiting resources can also explain those patterns (Münkemüller et al., 2020). We compared observed functional diversity to the one expected by assembling communities from a pool containing all the functionally characterized taxa detected across our sites. Therefore, we compared the functional diversity of a specific community to a global pool that included taxa capable of persisting but also those unable to persist in the conditions of the community. This makes our approach more sensitive to environmental filtering than to weaker competitor exclusion (de Bello et al., 2012), which operates among taxa able to persist given the specific abiotic constrains of a site.

Our global sampling design allowed us to identify the key role of climate, avoiding the confusion between temperature and biogeographical
<figure><img src="image_10.png" /></figure>

context. Mild temperatures occurred in forelands across a broad range of latitudes (Figure 1) because temperature also depends on foreland elevation and solar radiation. Therefore, the observed patterns are not biased by differences between tropical, temperate, and boreal areas (Freeman et al., 2021) (Figure S11), suggesting that local climate is a better proxy of the factors driving biodiversity than latitude. Forelands under different climates have thus distinct trajectories, likely determined by the balance between dispersal, environmental filtering and competition (Lei et al., 2015; Losapio et al., 2021). The biomass and diversity of below-ground animals, including nematodes, are tightly linked to the accumulation of organic matter (Fan et al., 2023; Fierer et al., 2009; Lei et al., 2015; van den Hoogen et al., 2019), which can be particularly fast under a mild summer climate (Khedim et al., 2021). Nematode colonization immediately after glacier retreat was fostered by mild climates with a pattern very similar to the one observed for soil organic matter (Khedim et al., 2021), and this effect was evident even after accounting for the heterogeneity of soil features within forelands (i.e., TOC and pH; Table S8). Overall, our results highlight that the interplay between climate and time is a major driver for the development of a broad range of features (from biotic to abiotic ones), determining ecosystem-level responses to glacier retreat.

Some studies observed nonlinear patterns in the colonization of deglaciated terrains by nematodes (Doblas-Miranda et al., 2008; Lei et al., 2015; Pothula & Adams, 2022). For instance, Doblas-Miranda et al. (2008) suggested that nematode diversity increases through time, with a peak around 5000 years after glacier retreat, and then decreases. Regardless of the climatic conditions, we did not detect analogous thresholds (Table S7). This outcome probably resulted from our shorter timespan (<400 years), which is well suited to unravel the consequences of ongoing climate change but cannot detect nonlinear dynamics that often occur over long (millennial) timescales (Losapio et al., 2021).

## 4.2 | Functional changes during colonization

Nematodes from most of the trophic groups were present since the early successional stages (Figure 3; Figure S7). The dispersal ability of nematodes, particularly the ones with r-ecological strategy and short generation time (low CP values), may allow their propagules to quickly (<15 years) reach deglaciated terrains from already structured communities in nearby sites. Indeed, early communities were mostly composed of nematodes with a r-strategy (Figures 2e and 4) that probably have high dispersal abilities and/or better fitness in these environments: a pattern observed for several organisms during the colonization of proglacial landscapes (Cauvy-Fraunié & Dangles, 2019; Erschbamer & Caccianiga, 2016; Ferris et al., 2001; Hågvar et al., 2020; Hågvar & Gobbi, 2022). With time, the richness of many trophic groups grew (Figure 3), and the functional identity of communities shifted, with an increase of functional richness through the addition of persister taxa (Figure 4). We did not find evidence for distinct successional trajectories among trophic groups, although the rates of richness increase over time differed among them (Table S9), in

accordance with previous observations (Doblas-Miranda et al., 2008; Pothula & Adams, 2022; van Leeuwen et al., 2018).

The accumulation of traits and the increase in richness of trophic groups over time might enhance the redundancy within functional groups and improve the resilience and stability of late communities (Naeem, 1998; Pillar et al., 2013). Such a pattern suggests a growing number of available niches and a potentially weaker environmental filtering over time, as also suggested by the shallow growth of SES values over time (Figure 5). Late communities are probably more sensitive to resource availability than to environmental filtering hosting taxa having better fitness in well-developed soils, but can play a larger number of functions (Bongers, 1990; Mulder & Maas, 2017; Sechi et al., 2018).

## 4.3 | Study limitations

The global scale of our analysis imposed some limitations. First, we considered average climatic conditions in each foreland, due to data availability constraints. As forelands are heterogeneous landscapes, fine-scale studies evaluating soil features and microclimatic conditions (e.g., soil humidity and below-surface temperature) are required to better understand processes acting locally (Ficetola, Marta, et al., 2021; Marta et al., 2023). Nevertheless, the inclusion of key fine-scale features of soil (pH and organic carbon content) did not affect our conclusions (Table S8), suggesting that results are robust to processes occurring at a fine scale.

We assessed nematode diversity with eDNA metabarcoding and a generalist marker targeting all eukaryotes to obtain broad-scale data that would be challenging to assemble with traditional methods. The Euka02 marker has an excellent amplification rate and a good taxonomic resolution for nematodes (Ficetola et al., 2023). This marker amplifies short DNA fragments, being particularly well suited for eDNA-based studies, where DNA is extracted from difficult substrates and is often degraded, consisting of short sequences (Taberlet et al., 2018). Nevertheless, the eDNA approach does not provide estimates of absolute biomass and often does not allow for species-level identification. This methodological limitation prevented the use of traits representing body size, as those traits can change across life stages, show phenotypic plasticity, and are imprecise when applied at the genus level (Mulder & Vonk, 2011). Primers amplifying longer fragments could provide a more accurate taxonomic identification and enable better differentiation of closely related species (Ficetola et al., 2023; Gattoni et al., 2023). However, long markers are suited for the analysis of DNA extracted directly from specimens and are not appropriate for DNA extracted from soil. We highlight that the functional database used in this study provides traits assigned at the genus or family level; thus, coarse taxonomic assignments should not strongly impact functional diversity analyses.

The number of detected taxonomic units was generally lower than what was observed using traditional approaches in glacier forefields (Doblas-Miranda et al., 2008; Lei et al., 2015). Such differences may occur due to the imperfect primer resolution, because multiple taxa are merged into the same molecular taxonomic unit or because some taxa
(e.g., the rare ones) remain undetected. To limit issues related to under-
detection, we used statistical approaches that consider the possibility
of missing rare taxa (Chazdon et al., 1998; Colwell & Coddington, 1994;
Gotelli & Colwell, 2011) and verified that our results remained con-
sistent when using metrics that account or do not account for the
relative abundance of taxa (Figure S3). Despite these limitations, soil
eDNA provides biodiversity estimates and patterns highly consistent
with those obtained with traditional approaches (Cantera et al., 2023;
Treonis et al., 2018), and the observed changes in nematode com-
munities are generally consistent with the patterns observed by both
local-scale and meta-analytic studies (Calderón-Sanou et al., 2019;
Devetter et al., 2021; Ji et al., 2013; Lei et al., 2015; Meyer et al., 2021;
Pothula & Adams, 2022). Furthermore, the large number of analyzed
chronosequences, their global distribution, and the coverage of a broad
range of conditions help to smooth the impact of these limitations and
reveal emergent signals with high generality (Mestre et al., 2022).

4.4 | Conclusion

There is growing evidence that climate change has uneven consequences across the globe, with high mountain and artic environments being among the most sensitive and impacted ecosystems (Hock et al., 2019). In these areas, biotic responses can be heterogeneous, as climatic differences across regions of the world can influence the trajectories of ecosystem evolution. We showed, on a worldwide scale, that nematodes colonizing soil after glacier retreat have successional trajectories that vary in speed among forelands. Local climate is a key driver for predicting the rate and nature of successions, and this is especially true at early stages of succession. Ecosystem functions and stability depend on local taxa and the traits they exhibit (Cardinale et al., 2012; Naeem et al., 2012). Nematodes play a key role in soil functioning, and the observed changes in their biodiversity might impact the whole dynamics of the ecosystems emerging after glacier retreat, as they can influence soil formation and the communities that follow. A key next step will be the development of multi-trophic analyses evaluating the complex relationships between ecosystem evolution and different biodiversity components (e.g., bacteria, fungi, protists, plants, and other animals; Giachello et al., 2023; Potapov et al., 2022). This will allow a better understanding of how climate change modifies habitats and trophic links in alpine environments and drives the development of new ecosystems (Martinez-Almoyna et al., 2019). Due to ongoing climate change, newly deglaciated terrains are expected to be increasingly important ecosystems in this century (Intergovernmental Panel on Climate Change, 2021; Rounce et al., 2023). Recognizing the complexity of biodiversity responses across climates will enhance our understanding of how these environments evolve and our ability to predict cascading effects at both global and local scales.

AUTHOR CONTRIBUTIONS

**Alessia Guerrieri:** Conceptualization; data curation; formal analysis; investigation; methodology; validation; visualization; writing – original draft; writing – review and editing.
**Isabel Cantera:** Conceptualization;

data curation; formal analysis; methodology; validation; visualization; writing – original draft; writing – review and editing. **Silvio Marta:** Data curation; formal analysis; investigation; methodology; visualization; writing – review and editing. **Aurélie Bonin:** Investigation; methodology; writing – review and editing. **Alexis Carteron:** Writing – review and editing. **Roberto Ambrosini:** Investigation; writing – review and editing. **Marco Caccianiga:** Investigation; writing – review and editing. **Fabien Anthelme:** Investigation; writing – review and editing. **Roberto Sergio Azzoni:** Investigation; writing – review and editing. **Peter Almond:** Investigation; writing – review and editing. **Pablo Alviz Gazitúa:** Investigation; writing – review and editing. **Sophie Cauvy-Fraunié:** Investigation; writing – review and editing. **Jorge Luis Ceballos Lievano:** Investigation; writing – review and editing. **Milap Chand Sharma:** Investigation; writing – review and editing. **John Clague:** Investigation; writing – review and editing. **Justiniano Alejo Cochachín Rapre:** Investigation. **Chiara Compostella:** Investigation; writing – review and editing. **Rolando Cruz Encarnación:** Investigation; writing – review and editing. **Olivier Dangles:** Investigation; writing – review and editing. **Philip Deline:** Investigation; writing – review and editing. **Andre Eger:** Investigation; writing – review and editing. **Sergey Erokhin:** Investigation; writing – review and editing. **Andrea Franzetti:** Investigation; writing – review and editing. **Ludovic Gielly:** Investigation; writing – review and editing. **Fabrizio Gili:** Data curation; investigation; writing – review and editing. **Mauro Gobbi:** Investigation; writing – review and editing. **Sigmund Hågvar:** Investigation; writing – review and editing. **Norine Khedim:** Investigation; writing – review and editing. **Rosa Isela Meneses:** Investigation; writing – review and editing. **Gwendolyn Peyre:** Investigation; writing – review and editing. **Francesca Pittino:** Investigation; writing – review and editing. **Angela Proietto:** Investigation; writing – review and editing. **Antoine Rabatel:** Data curation; formal analysis; investigation; writing – review and editing. **Nurai Urseitova:** Investigation; writing – review and editing. **Yan Yang:** Investigation; writing – review and editing. **Vitalii Zaginaev:** Investigation; writing – review and editing. **Andrea Zerboni:** Data curation; investigation; writing – review and editing. **Anaïs Zimmer:** Investigation; writing – review and editing. **Pierre Taberlet:** Conceptualization; methodology; writing – review and editing. **Guglielmina Adele Diolaiuti:** Investigation; writing – review and editing. **Jerome Poulenard:** Conceptualization; investigation; writing – review and editing. **Diego Fontaneto:** Data curation; formal analysis; writing – review and editing. **Wilfried Thuiller:** Conceptualization; formal analysis; writing – review and editing. **Gentile Francesco Ficetola:** Conceptualization; formal analysis; funding acquisition; investigation; methodology; project administration; resources; supervision; validation; visualization; writing – original draft; writing – review and editing.

AFFILIATIONS

<footnote n="1">¹Dipartimento di Scienze e Politiche Ambientali, Università degli Studi di Milano, Milano, Italy</footnote>

<footnote n="2">²Argaly, Bâtiment CleanSpace, Sainte-Hélène-du-Lac, France</footnote>

<footnote n="3">³Institute of Geosciences and Earth Resources, CNR, Pisa, Italy</footnote>

<footnote n="4">⁴Dipartimento di Bioscienze, Università degli Studi di Milano, Milano, Italy</footnote>
⁵Laboratory AMAP, IRD, University of Montpellier, CIRAD, CNRS, INRA, Montpellier, France

⁶Dipartimento di Scienze della Terra "Ardito Desio", Università degli Studi di Milano, Milano, Italy

⁷Department of Soil and Physical Sciences, Lincoln University, Lincoln, New Zealand

⁸Departamento de Ciencias Biológicas y Biodiversidad, Universidad de Los Lagos, Osorno, Chile

⁹INRAE, UR RIVERLY, Centre de Lyon-Villeurbanne, Villeurbanne, France

¹⁰Instituto de Hidrología, Meteorología y Estudios Ambientales IDEAM, Bogotá, Colombia

¹¹Department of Geography, School of Environment and Earth Sciences, Central University of Punjab, Bathinda, Punjab, India

¹²Centre for the Study of Regional Development - School of Social Sciences, Jawaharlal Nehru University, New Delhi, India

¹³Department of Earth Sciences, Simon Fraser University, Burnaby, British Columbia, Canada

¹⁴Área de Evaluación de Glaciares y Lfagunas, Autoridad Nacional del Agua, Huaraz, Peru

¹⁵CEFE, Univ Montpellier, CNRS, EPHE, IRD, Univ Paul Valéry Montpellier 3, Montpellier, France

¹⁶Univ. Savoie Mont Blanc, Univ. Grenoble Alpes, EDYTEM, Chambéry, France

¹⁷Mannaki Whenua - Landcare Research, Soils and Landscapes, Lincoln, New Zealand

¹⁸Institute of Water Problems and Hydro-Energy, Kyrgyz National Academy of Sciences, Bishkek, Kyrgyzstan

¹⁹Department of Earth and Environmental Sciences (DISAT), University of Milano-Bicocca, Milano, Italy

²⁰Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, LECA, Grenoble, France

²¹Department of Life Sciences and Systems Biology, University of Turin, Turin, Italy

²²Research and Museum Collections Office, Climate and Ecology Unit, MUSE-Science Museum, Corso del Lavoro e della Scienza, Trento, Italy

²³Faculty of Environmental Sciences and Natural Resource Management (INA), Norwegian University of Life Sciences, Ås, Norway

²⁴UiT - The Arctic University of Norway, Tromsø Museum, Tromsø, Norway

²⁵Herbario Nacional de Bolivia: La Paz, La Paz, Bolivia

²⁶Universidad Católica del Norte, Antofagasta, Chile

²⁷Department of Civil and Environmental Engineering, University of the Andes, Bogotá, Colombia

²⁸Univ. Grenoble Alpes, CNRS, IRD, Grenoble-INP, Institut des Géosciences de l'Environnement (IGE, UMR 5001), Grenoble, France

²⁹Institute of Mountain Hazards and Environment, Chinese Academy of Sciences, Chengdu, China

³⁰Mountain Societies Research Institute, University of Central Asia, Bishkek, Kyrgyzstan

³¹Department of Geography and the Environment, University of Texas at Austin, Austin, Texas, USA

³²CNR - Water Research Institute, Verbania, Italy

³³NBFC - National Biodiversity Future Center, Palermo, Italy

## ACKNOWLEDGMENTS

This study is funded by the European Research Council under the European Community's Horizon 2020 Programme, Grant Agreement no. 772284 (IceCommunities). This research was also funded by Biodiversa+, the European Biodiversity Partnership under the 2021-2022 BiodivProtect joint call for research proposals, co-funded by the European Commission (grant agreement no. 101052342 "PrioritIce-Vanishing habitats: conservation priorities for glacier-related biodiversity threatened by climate change") and with the funding organizations MUR and ANR. We are grateful to Rüdiger Kaufmann, Antoine Guisan, Katrin Sieron, and Marco Aurelio Morales-Martínez for help and discussions at various phases of this project. This work is dedicated to the memory of Justiniano Alejo Cochachín Rapre.

## CONFLICT OF INTEREST STATEMENT

The authors declare no competing interests.

## DATA AVAILABILITY STATEMENT

Supplementary materials including Figures S1-S11 and Tables S1-S10 underlying the main results of the study are provided in the online version of the article. Raw sequencing data are deposited in Zenodo (https://zenodo.org/record/6620359). Sequences assigned to Eukaryotes after the bioinformatics pipeline are in the "Euka02_global_cluster_filt_tag.txt" file provided in Figshare (https://figshare.com/s/839b43ed65f9ffae809c), as well as scripts for reproducing (1) bioinformatic steps, (2) taxonomic assignation, (3) MOTU filtering, (4) the measurement of estimated biodiversity using first-order jackknife estimator, (5) functional variables analyses, and (6) model computation and visualization.

## ORCID

Isabel Cantera ID https://orcid.org/0000-0003-3161-1878

Silvio Marta ID https://orcid.org/0000-0001-8850-610X

Sigmund Hågvar ID https://orcid.org/0000-0002-9020-9005

Yan Yang ID https://orcid.org/0000-0003-0858-7603

## REFERENCES

Aglieri, G., Baillie, C., Mariani, S., Cattano, C., Calò, A., Turco, G., Spatafora, D., Di Franco, A., Di Lorenzo, M., Guidetti, P., & Milazzo, M. (2020). Environmental DNA effectively captures functional diversity of coastal fish communities. *Molecular Ecology*, *30*(13), 3127–3139. https://doi.org/10.1111/mec.15661

Anthelme, F., Carrasquer, I., Ceballos, J. L., & Peyre, G. (2022). Novel plant communities after glacial retreat in Colombia: (Many) losses and (few) gains. *Alpine Botany*, *132*(2), 211–222. https://doi.org/10.1007/s00035-022-00282-1

Aros-Mualin, D., Noben, S., Karger, D. N., Carvajal-Hernández, C. I., Salazar, L., Hernández-Rojas, A., Kluge, J., Sundue, M. A., Lehnert, M., Quandt, D., & Kessler, M. (2021). Functional diversity in ferns is driven by species richness rather than by environmental constraints. *Frontiers in Plant Science*, *11*, 2128. https://doi.org/10.3389/fpls.2020.615723

Bahram, M., Hildebrand, F., Forslund, S. K., Anderson, J. L., Soudzilovskaia, N. A., Bodegom, P. M., Bengtsson-Palme, J., Anslan, S., Coelho, L. P., Harend, H., Huerta-Cepas, J., Medema, M. H., Maltz, M. R., Mundra, S., Olsson, P. A., Pent, M., Pölme, S., Sunagawa, S., Ryberg, M., ... Bork, P. (2018). Structure and function of the global topsoil microbiome. *Nature*, *560*(7717), 233–237. https://doi.org/10.1038/s41586-018-0386-6

Bertness, M. D., & Callaway, R. (1994). Positive interactions in communities. *Trends in Ecology & Evolution*, *9*(5), 191–193. https://doi.org/10.1016/0169-5347(94)90088-4

Bongers, T. (1990). The maturity index: An ecological measure of environmental disturbance based on nematode species composition. *Oecologia*, *83*(1), 14–19. https://doi.org/10.1007/BF00324627

Bonin, A., Guerrieri, A., & Ficetola, G. F. (2023). Optimal sequence similarity thresholds for clustering of molecular operational taxonomic units in DNA metabarcoding studies. *Molecular Ecology Resources*, *23*(2), 368–381. https://doi.org/10.1111/1755-0998.13709

Borcard, D., Gillet, F., & Legendre, P. (2011). *Numerical ecology with R*. Springer. http://xa.yimg.com/kq/groups/19243105/1919134110/name/Numerical.pdfscriptefunzioni:adn.biol.umontreal.ca/~numericalecology/numecolR/
<figure><img src="image_11.png" /></figure>

Boyer, F., Mercier, C., Bonin, A., Le Bras, Y., Taberlet, P., & Coissac, E. (2016). OBITools: A unix-inspired software package for DNA metabarcoding. *Molecular Ecology Resources, **16*(1), 176-182.

Bürkner, P.-C. (2017). brms: An R package for Bayesian multilevel models using Stan. *Journal of Statistical Software*, **80**(1), 1-28. https://doi.org/10.18637/jss.v080.i01

Cadotte, M. W., Carscadden, K., & Mirotchnick, N. (2011). Beyond species: Functional diversity and the maintenance of ecological processes and services: Functional diversity in ecology and conservation. *Journal of Applied Ecology*, **48**(5), 1079-1087. https://doi.org/10.1111/j.1365-2664.2011.02048.x

Calderón-Sanou, I., Münkemüller, T., Boyer, F., Zinger, L., & Thuiller, W. (2019). From environmental DNA sequences to ecological conclusions: How strong is the influence of methodological choices? *Journal of Biogeography*, **47**, 193-206. https://doi.org/10.1111/jbi.13681

Cantera, I., Carteron, A., Guerrieri, A., Marta, S., Bonin, A., Ambrosini, R., Anthelme, F., Azzoni, R., Almond, P., Alviz Gazitúa, P., Cauvy-fraunié, S., Lievano, J. C., Chand, P., Sarma, M. C., Clague, J., Rapre, J. A. C., Compostella, C., Encarnación, R. C., Dangles, O., ... Ficetola, F. (2023). The importance of species addition versus replacement varies over succession in plant communities after glacial retreat. *Research Square*. https://doi.org/10.21203/rs.3.rs-2482972/v1

Cardinale, B. J., Duffy, J. E., Gonzalez, A., Hooper, D. U., Perrings, C., Venail, P., Narwani, A., Mace, G. M., Tilman, D., Wardle, D. A., Kinzig, A. P., Daily, G. C., Loreau, M., Grace, J. B., Larigauderie, A., Srivastava, D. S., & Naeem, S. (2012). Biodiversity loss and its impact on humanity. *Nature*, **486**(7401), 59-67. https://doi.org/10.1038/nature11148

Cardoso, P., Rigal, F., & Carvalho, J. C. (2015). BAT - Biodiversity Assessment Tools, an R package for the measurement and estimation of alpha and beta taxon, phylogenetic and functional diversity. *Methods in Ecology and Evolution*, **6**(2), 232-236. https://doi.org/10.1111/2041-210X.12310

Cauvy-Fraunié, S., & Dangles, O. (2019). A global synthesis of biodiversity responses to glacier retreat. *Nature Ecology & Evolution*, **3**(12), 1675-1685. https://doi.org/10.1038/s41559-019-1042-8

Chazdon, R. L., Colwell, R. K., Denslow, J. S., & Guariguata, M. R. (1998). Statistical methods for estimating species richness of woody regeneration in primary and secondary rain forests of NE Costa Rica. In F. Dallmeier & J. A. Comiskey (Eds.), *Forest biodiversity research, monitoring and modeling: Conceptual background and Old World case studies* (pp. 285-309). Parthenon Publishing.

Colwell, R. K., & Coddington, J. A. (1994). Estimating terrestrial biodiversity through extrapolation. *Philosophical Transactions of the Royal Society B: Biological Sciences*, **345**, 101-118.

Crowther, T. W., van den Hoogen, J., Wan, J., Mayes, M. A., Keiser, A. D., Mo, L., Averill, C., & Maynard, D. S. (2019). The global soil community and its influence on biogeochemistry. *Science*, **365**, eaav0550. https://doi.org/10.1126/science.aav0550

Davis, N. M., Proctor, D. M., Holmes, S. P., Relman, D. A., & Callahan, B. J. (2018). Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data. *Microbiome*, **6**(1), 226. https://doi.org/10.1186/s40168-018-0605-2

de Bello, F., Price, J. N., Münkemüller, T., Liira, J., Zobel, M., Thuiller, W., Gerhold, P., Götzenberger, L., Lavergne, S., Lepš, J., Zobel, K., & Pärtel, M. (2012). Functional species pool framework to test for biotic effects on community assembly. *Ecology*, **93**(10), 2263-2273. https://doi.org/10.1890/11-1394.1

de Bello, F., Thuiller, W., Lepš, J., Choler, P., Clément, J.-C., Macek, P., Sebastià, M.-T., & Lavorel, S. (2009). Partitioning of functional diversity reveals the scale and extent of trait convergence and divergence. *Journal of Vegetation Science*, **20**(3), 475-486. https://doi.org/10.1111/j.1654-1103.2009.01042.x

Delgado-Baquerizo, M., Reich, P. B., Trivedi, C., Eldridge, D. J., Abades, S., Alfaro, F. D., Bastida, F., Berhe, A. A., Cutler, N. A., Gallardo, A., García-Velázquez, L., Hart, S. C., Hayes, P. E., He, J.-Z., Hseu, Z.-Y.,

Hu, H.-W., Kirchmair, M., Neuhauser, S., Pérez, C. A., ... Singh, B. K. (2020). Multiple elements of soil biodiversity drive ecosystem functions across biomes. *Nature Ecology & Evolution*, **4**(2), 210-220. https://doi.org/10.1038/s41559-019-1084-y

Devetter, M., Hanel, L., Raschmanova, N., Bryndova, M., & Schlaghamersky, J. (2021). Terrestrial invertebrates along a gradient of deglaciation in Svalbard: Long-term development of soil fauna communities. *Geoderma*, **383**, 11. https://doi.org/10.1016/j.geoderma.2020.114720

Dickie, I. A., Boyer, S., Buckley, H., Duncan, R. P., Gardner, P., Hogg, I. D., Holdaway, R. J., Lear, G., Makiola, A., Morales, S. E., Powell, J. R., & Weaver, L. (2018). Towards robust and repeatable sampling methods in eDNA based studies. *Molecular Ecology Resources*, **18**, 940-952. https://doi.org/10.1111/1755-0998.12907

Doblas-Miranda, E., Wardle, D. A., Peltzer, D. A., & Yeates, G. W. (2008). Changes in the community structure and diversity of soil invertebrates across the Franz Josef Glacier chronosequence. *Soil Biology and Biochemistry*, **40**(5), 1069-1081. https://doi.org/10.1016/j.soilbio.2007.11.026

Erschbamer, B., & Caccianiga, M. S. (2016). Glacier forelands: Lessons of plant population and community development. In F. M. Cánovas, U. Lüttge, & R. Matyssek (Eds.), *Progress in botany* (Vol. 78, pp. 259-284). Springer International Publishing. https://doi.org/10.1007/124_2016_4

Fan, K., Chu, H., Eldridge, D. J., Gaitan, J. J., Liu, Y.-R., Sokoya, B., Wang, J.-T., Hu, H.-W., He, J.-Z., Sun, W., Cui, H., Alfaro, F. D., Abades, S., Bastida, F., Díaz-López, M., Bamigboye, A. R., Berdugo, M., Blanco-Pastor, J. L., Grebenc, T., ... Delgado-Baquerizo, M. (2023). Soil biodiversity supports the delivery of multiple ecosystem functions in urban greenspaces. *Nature Ecology & Evolution*, **7**(1), 113-126. https://doi.org/10.1038/s41559-022-01935-4

Ferris, H. (2010). Form and function: Metabolic footprints of nematodes in the soil food web. *European Journal of Soil Biology*, **46**(2), 97-104. https://doi.org/10.1016/j.ejsobi.2010.01.003

Ferris, H., Bongers, T., & de Goede, R. G. M. (2001). A framework for soil food web diagnostics: Extension of the nematode faunal analysis concept. *Applied Soil Ecology*, **18**(1), 13-29. https://doi.org/10.1016/S0929-1393(01)00152-4

Ficetola, G. F., Boyer, F., Valentini, A., Bonin, A., Meyer, A., Dejean, T., Gaboriaud, C., Usseglio-Polatera, P., & Taberlet, P. (2021). Comparison of markers for the monitoring of freshwater benthic biodiversity through DNA metabarcoding. *Molecular Ecology*, **30**(13), 3189-3202. https://doi.org/10.1111/mec.15632

Ficetola, G. F., Coissac, E., Zundel, S., Riaz, T., Shehzad, W., Bessière, J., Taberlet, P., & Pompanon, F. (2010). An in silico approach for the evaluation of DNA barcodes. BMC Genomics, **11**(1), 434. https://doi.org/10.1186/1471-2164-11-434

Ficetola, G. F., & Denoël, M. (2009). Ecological thresholds: An assessment of methods to identify abrupt changes in species-habitat relationships. *Ecography*, **32**, 1075-1084.

Ficetola, G. F., Guerrieri, A., Cantera, I., & Bonin, A. (2023). In silico assessment of 18S rDNA metabarcoding markers for the characterization of nematode communities. *bioRxiv*. https://doi.org/10.1101/2023.08.01.551472

Ficetola, G. F., Marta, S., Guerrieri, A., Gobbi, M., Ambrosini, R., Fontaneto, D., Zerboni, A., Poulenard, J., Caccianiga, M., & Thuiller, W. (2021). Dynamics of ecological communities following current retreat of glaciers. *Annual Review of Ecology, Evolution, and Systematics*, **52**(1), 405-426. https://doi.org/10.1146/annurev-econsys-010521-040017

Ficetola, G. F., Pansu, J., Bonin, A., Coissac, E., Giguet-Covex, C., De Barba, M., Gielly, L., Lopes, C. M., Boyer, F., Pompanon, F., Rayé, G., & Taberlet, P. (2015). Replication levels, false presences and the estimation of the presence/absence from eDNA metabarcoding data. *Molecular Ecology Resources*, **15**(3), 543-556. https://doi.org/10.1111/1755-0998.12338
<figure><img src="image_12.png" /></figure>

Fierer, N., Strickland, M. S., Liptzin, D., Bradford, M. A., & Cleveland, C. C. (2009). Global patterns in belowground communities. *Ecology Letters*, *12*(11), 1238-1249. https://doi.org/10.1111/j.1461-0248.2009.01360.x

Fontaneto, D. (2019). Long-distance passive dispersal in microscopic aquatic animals. *Movement Ecology*, 7(1), 10. https://doi.org/10.1186/s40462-019-0155-7

Franco, A. L. C., Adams, B. J., Diaz, M. A., Lemoine, N. P., Dragone, N. B., Fierer, N., Lyons, W. B., Hogg, I., & Wall, D. H. (2022). Response of Antarctic soil fauna to climate-driven changes since the Last Glacial Maximum. *Global Change Biology*, 28(2), 644–653. https://doi.org/10.1111/gcb.15940

Franco, A. L. C., Gherardi, L. A., de Tomasel, C. M., Andriuzzi, W. S., Ankrom, K. E., Shaw, E. A., Bach, E. M., Sala, O. E., & Wall, D. H. (2019). Drought suppresses soil predators and promotes root herbivores in mesic, but not in xeric grasslands. *Proceedings of the National Academy of Sciences of the United States of America*, 116(26), 12883–12888. https://doi.org/10.1073/pnas.1900572116

Freeman, B. G., Song, Y., Feeley, K. J., & Zhu, K. (2021). Montane species track rising temperatures better in the tropics than in the temperate zone. *Ecology Letters*, 24, 1697–1708. https://doi.org/10.1111/ele.13762

Fridley, J. D., & Wright, J. P. (2018). Temperature accelerates the rate fields become forests. *Proceedings of the National Academy of Sciences of the United States of America*, 115(18), 4702–4706. https://doi.org/10.1073/pnas.1716665115

Gattoni, K., Gendron, E. M. S., Sandoval-Ruiz, R., Borgemeier, A., McQueen, J. P., Shepherd, R. M., Slos, D., Powers, T. O., & Porazinska, D. L. (2023). 18S-NemaBase: Curated 18S rRNA database of nematode sequences. *Journal of Nematology*, 55(1), 20230006. https://doi.org/10.2478/jofnem-2023-0006

Geisen, S., Snoek, L. B., ten Hooven, F. C., Duyts, H., Kostenko, O., Bloem, J., Martens, H., Quist, C. W., Helder, J. A., & van der Putten, W. H. (2018). Integrating quantitative morphological and qualitative molecular methods to analyse soil nematode community responses to plant range expansion. *Methods in Ecology and Evolution*, 9(6), 1366–1378. https://doi.org/10.1111/2041-210X.12999

Gelman, A., Hwang, J., & Vehtari, A. (2014). Understanding predictive information criteria for Bayesian models. *Statistics and Computing*, 24(6), 997–1016. https://doi.org/10.1007/s11222-013-9416-2

Giachello, S., Cantera, I., Carteron, A., Marta, S., Cipriano, C., Guerrieri, A., Bonin, A., Thuiller, W., & Ficetola, G. F. (2023). Toward a common set of functional traits for soil protists. *Soil Biology and Biochemistry*, 187, 109207. https://doi.org/10.1016/j.soilbio.2023.109207

Gibbard, P. L. (1993). Gale S. J. & Hoare P. G. 1991. Quaternary sediments. viii + 323 pp. London: 904 Belhaven Press (Pinter Publishers); New York, Toronto: John Wiley. Price £42.00 (hard 905 covers). ISBN 1 85293 017 9. *Geological Magazine*, 130(1), 126. https://doi.org/10.1017/S0016756800023852

Gotelli, N. J., & Colwell, R. K. (2011). Estimating species richness. *Biological Diversity: Frontiers in Measurement and Assessment*, 12, 39–54.

Gotelli, N. J., & Graves, G. R. (1996). *Null models in ecology*. Smithsonian Institution.

Gotelli, N. J., & McCabe, D. J. (2002). Species co-occurrence: A meta-analysis of J. M. Diamond's assembly rules model. *Ecology*, 83(8), 2091–2096. https://doi.org/10.2307/3072040

Götzenberger, L., de Bello, F., Bräthen, K. A., Davison, J., Dubuis, A., Guisan, A., Lepš, J., Lindborg, R., Moora, M., Pärtel, M., Pellissier, L., Pottier, J., Vittoz, P., Zobel, K., & Zobel, M. (2012). Ecological assembly rules in plant communities—approaches, patterns and prospects. *Biological Reviews*, 87(1), 111–127. https://doi.org/10.1111/j.1469-185X.2011.00187.x

Grenié, M., Denelle, P., Tucker, C. M., Munoz, F., & Violle, C. (2017). funrar: An R package to characterize functional rarity. *Diversity and Distributions*, 23(12), 1365–1371. https://doi.org/10.1111/ddi.12629

Guardiola, M., Uriz, M. J., Taberlet, P., Coissac, E., Wangensteen, O. S., & Turon, X. (2015). Deep-sea, deep-sequencing: Metabarcoding extracellular DNA from sediments of marine canyons. *PLoS One*, 10(10), 1–26. https://doi.org/10.1371/journal.pone.0139633

Guerrieri, A., Bonin, A., Münkemüller, T., Gielly, L., Thuiller, W., & Francesco Ficetola, G. (2021). Effects of soil preservation for biodiversity monitoring using environmental DNA. *Molecular Ecology*, 30(13), 3313–3325. https://doi.org/10.1111/mec.15674

Guerrieri, A., Carteron, A., Bonin, A., Marta, S., Ambrosini, R., Caccianiga, M., Cantera, I., Compostella, C., Doliaiuti, G., Fontaneto, D., Gielly, L., Gili, F., Gobbi, M., Poulenard, J., Taberlet, P., Zerboni, A., Thuiller, W., & Ficetola, G. F. (2023). Metabarcoding data reveal vertical multitaxa variation in topsoil communities during the colonization of deglaciated forelands. *Molecular Ecology*, 32, 6304–6319. https://doi.org/10.1111/mec.16669

Hågvar, S., & Gobbi, M. (2022). The role of arthropods in early colonization near melting glaciers: Contradictions between ecological assumptions and recent study results. *Acta Oecologica*, 114, 103820. https://doi.org/10.1016/j.actoac.2022.103820

Hågvar, S., Gobbi, M., Kaufmann, R., Ingimarsdóttir, M., Caccianiga, M., Valle, B., Pantini, P., Fanciulli, P. P., & Vater, A. (2020). Ecosystem birth near melting glaciers: A review on the pioneer role of ground-dwelling arthropods. *Insects*, 11(9), 644. https://doi.org/10.3390/insects11090644

Hanusch, M., He, X., Ruiz-Hernández, V., & Junker, R. R. (2022). Succession comprises a sequence of threshold-induced community assembly processes towards multidiversity. *Communications Biology*, 5(1), 424. https://doi.org/10.1038/s42003-022-03372-2

Hock, R., Rasul, G., Adler, C., Cáceres, B., Gruber, S., Hirabayashi, Y., Jackson, M., Kääb, A., Kang, S., Kutuzov, S., Milner, A., Molau, U., Morin, S., Orlove, B., Steltzer, H., Allen, S., Arenson, L., Baneerjee, S., Barr, I., ... Zhang, Y. (2019). High mountain areas. In H.-O. Pörtner, D. C. Roberts, V. Masson-Delmotte, P. Zhai, M. Tignor, E. Poloczanska, K. Mintenbeck, M. Nicolai, A. Okem, J. Petzold, B. Rama, & N. Weyer (Eds.), *The ocean and cryosphere in a changing climate* (pp. 131–202). Intergovernmental Panel on Climate Change.

Hugonnet, R., McNabb, R., Berthier, E., Menounos, B., Nuth, C., Girod, L., Farinotti, D., Huss, M., Dussaillant, I., Brun, F., & Kääb, A. (2021). Accelerated global glacier mass loss in the early twenty-first century. *Nature*, 592(7856), 726–731. https://doi.org/10.1038/s41586-021-03436-z

Intergovernmental Panel on Climate Change. (2021). *Climate change 2021: The physical science basis*. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press. https://www.ipcc.ch/report/ar6/wg1/#SPM

Ji, Y., Ashton, L., Pedley, S. M., Edwards, D. P., Tang, Y., Nakamura, A., Kitching, R., Dolman, P. M., Woodcock, P., Edwards, F. A., Larsen, T. H., Hsu, W. W., Benedick, S., Hamer, K. C., Wilcove, D. S., Bruce, C., Wang, X., Levi, T., Lott, M., ... Yu, D. W. (2013). Reliable, verifiable and efficient monitoring of biodiversity via metabarcoding. *Ecology Letters*, 16(10), 1245–1257. https://doi.org/10.1111/ele.12162

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

Kaufmann, R. (2002). Glacier foreland colonisation: Distinguishing between short-term and long-term effects of climate change. *Oecologia*, 130(3), 470–475. https://doi.org/10.1007/s00442-001-0815-2

Khedim, N., Cécillon, L., Poulenard, J., Barré, P., Baudin, F., Marta, S., Rabatel, A., Dentant, C., Cauvy-Fraunié, S., Anthelme, F., Gielly, L., Ambrosini, R., Franzetti, A., Azzoni, R. S., Caccianniga, M. S., Compostella, C., Clague, J., Tielidze, L., Messager, E., ... Ficetola, G. F. (2021). Topsoil organic matter build-up in glacier forelands
<figure><img src="image_13.png" /></figure>

around the world. *Global Change Biology*, 27, 1662-1677. https://doi.org/10.1111/gcb.15496

Kuczynski, L., & Grenouillet, G. (2018). Community disassembly under global change: Evidence in favor of the stress-dominance hypothesis. *Global Change Biology*, 24(9), 4417-4427. https://doi.org/10.1111/gcb.14320

Laliberté, E., & Legendre, P. (2010). A distance-based framework for measuring functional diversity from multiple traits. *Ecology*, 91(1), 299-305. https://doi.org/10.1890/08-2244.1

Le Bagousse-Pinguet, Y., Soliveres, S., Gross, N., Torices, R., Berdugo, M., & Maestre, F. T. (2019). Phylogenetic, functional, and taxonomic richness have both positive and negative effects on ecosystem multifunctionality. *Proceedings of the National Academy of Sciences of the United States of America*, 116(17), 8419-8424. https://doi.org/10.1073/pnas.1815727116

Lei, Y. B., Zhou, J., Xiao, H. F., Duan, B. L., Wu, Y. H., Korpelainen, H., & Li, C. Y. (2015). Soil nematode assemblages as bioindicators of primary succession along a 120-year-old chronosequence on the Hailuoguong Glacier forefield, SW China. *Soil Biology & Biochemistry*, 88, 362-371. https://doi.org/10.1016/j.soilbio.2015.06.013

Lembrechts, J. J., van den Hoogen, J., Aalto, J., Ashcroft, M. B., De Frenne, P., Kemppinen, J., Kopecký, M., Luoto, M., Maclean, I. M. D., Crowther, T. W., Bailey, J. J., Haesen, S., Klinges, D. H., Niittynen, P., Scheffers, B. R., Van Meerbeek, K., Aartsma, P., Abdalaze, O., Abedi, M., ... Lenoir, J. (2022). Global maps of soil temperature. *Global Change Biology*, 28(9), 3110-3144. https://doi.org/10.1111/gcb.16060

Loreau, M., Naeem, S., Inchausti, P., Bengtsson, J., Grime, J. P., Hector, A., Hooper, D. U., Huston, M. A., Raffaelli, D., Schmid, B., Tilman, D., & Wardle, D. A. (2001). Biodiversity and ecosystem functioning: Current knowledge and future challenges. *Science*, 294(5543), 804-808. https://doi.org/10.1126/science.1064088

Losapio, G., Cerabolini, B. E. L., Maffioletti, C., Tampucci, D., Gobbi, M., & Caccianiga, M. (2021). The consequences of glacier retreat are uneven between plant species. *Frontiers in Ecology and Evolution*, 8, 616562. https://doi.org/10.3389/fevo.2020.616562

Ma, Q., Yu, H., Liu, X., Xu, Z., Zhou, G., & Shi, Y. (2018). Climatic warming shifts the soil nematode community in a desert steppe. *Climatic Change*, 150(3), 243-258. https://doi.org/10.1007/s10584-018-2277-0

Macarthur, R., & Levins, R. (1967). The limiting similarity, convergence, and divergence of coexisting species. *The American Naturalist*, 101(921), 377-385. https://doi.org/10.1086/282505

Mächler, E., Walser, J.-C., & Altermatt, F. (2021). Decision-making and best practices for taxonomy-free environmental DNA metabarcoding in biomonitoring using Hill numbers. *Molecular Ecology*, 30(13), 3326-3339. https://doi.org/10.1111/mec.15725

Marta, S., Azzoni, R. S., Fugazza, D., Tielidze, L., Chand, P., Sieron, K., Almond, P., Ambrosini, R., Anthelme, F., Alviz Gazitúa, P., Bhambri, R., Bonin, A., Caccianiga, M., Cauvy-Fraunié, S., Ceballos Lievano, J. L., Clague, J., Cochachín Rapre, J. A., Dangles, O., Deline, P., ... Ficetola, G. F. (2021). The retreat of mountain glaciers since the little ice age: A spatially explicit database. *Data*, 6(10), 107. https://doi.org/10.3390/data6100107

Marta, S., Zimmer, A., Caccianiga, M., Gobbi, M., Ambrosini, R., Azzoni, R. S., Gili, F., Pittino, F., Thuiller, W., Provenzale, A., & Ficetola, G. F. (2023). Heterogeneous changes of soil microclimate in high mountains and glacier forelands. *Nature Communications*, 14(1), 6. https://doi.org/10.1038/s41467-023-41063-6

Martinez-Almoya, C., Thuiller, W., Chalmandrier, L., Ohlmann, M., Foulquier, A., Clément, J.-C., Zinger, L., & Münkemüller, T. (2019). Multi-trophic β-diversity mediates the effect of environmental gradients on the turnover of multiple ecosystem functions. *Functional Ecology*, 33(10), 2053-2064. https://doi.org/10.1111/1365-2435.13393

Mason, N. W. H., Mouillot, D., Lee, W. G., & Wilson, J. B. (2005). Functional richness, functional evenness and functional divergence: The

primary components of functional diversity. *Oikos*, 111(1), 112-118. https://doi.org/10.1111/j.0030-1299.2005.13886.x

Mestre, F., Rozenfeld, A., & Araújo, M. B. (2022). Human disturbances affect the topology of food webs. *Ecology Letters*, 25(11), 2476-2488. https://doi.org/10.1111/ele.14107

Meyer, A., Boyer, F., Valentini, A., Bonin, A., Ficetola, G. F., Beisel, J.-N., Bouquerel, J., Wagner, P., Gaboriaud, C., Leese, F., Dejean, T., Taberlet, P., & Usseglio-Polatera, P. (2021). Morphological vs. DNA metabarcoding approaches for the evaluation of stream ecological status with benthic invertebrates: Testing different combinations of markers and strategies of data filtering. *Molecular Ecology*, 30(13), 3203-3220. https://doi.org/10.1111/mec.15723

Mouchet, M. A., Villéger, S., Mason, N. W. H., & Mouillot, D. (2010). Functional diversity measures: An overview of their redundancy and their ability to discriminate community assembly rules: Functional diversity measures. *Functional Ecology*, 24(4), 867-876. https://doi.org/10.1111/j.1365-2435.2010.01695.x

Mouillot, D., Graham, N. A. J., Villéger, S., Mason, N. W. H., & Bellwood, D. R. (2013). A functional approach reveals community responses to disturbances. *Trends in Ecology & Evolution*, 28(3), 167-177. https://doi.org/10.1016/j.tree.2012.10.004

Muggeo, V. M. R. (2003). Estimating regression models with unknown break-points. *Statistics in Medicine*, 22, 3055-3071.

Mulder, C., & Maas, R. (2017). Unifying the functional diversity in natural and cultivated soils using the overall body-mass distribution of nematodes. *BMC Ecology*, 17(1), 36. https://doi.org/10.1186/s12898-017-0145-9

Mulder, C., & Vonk, J. A. (2011). Nematode traits and environmental constraints in 200 soil systems: Scaling within the 60-6000 µm body size range. *Ecology*, 92(10), 2004. https://doi.org/10.1890/11-0546.1

Münkemüller, T., Gallien, L., Pollock, L. J., Barros, C., Carboni, M., Chalmandrier, L., Mazel, F., Mokany, K., Roquet, C., Smyčka, J., Talluto, M. V., & Thuiller, W. (2020). Dos and don'ts when inferring assembly rules from diversity patterns. *Global Ecology and Biogeography*, 29(7), 1212-1229. https://doi.org/10.1111/geb.13098

Naeem, S. (1998). Species redundancy and ecosystem reliability. *Conservation Biology*, 12(1), 39-45. https://doi.org/10.1111/j.1523-1739.1998.96379.x

Naeem, S., Duffy, J. E., & Zavaleta, E. (2012). The functions of biological diversity in an age of extinction. *Science*, 336(6087), 1401-1406. https://doi.org/10.1126/science.1215855

Neher, D. A. (2001). Role of nematodes in soil health and their use as indicators. *Journal of Nematology*, 33(4), 161-168.

Oksanen, J., Blanchet, F. G., Kindt, R., Legendre, P., Minchin, P. R., O'hara, R. B., Simpson, G. L., Solymos, P., Stevens, M. H. H., & Wagner, H. (2013). *Package 'vegan'*. Community ecology package, version, 2(9).

Petchey, O. L., & Gaston, K. J. (2002). Functional diversity (FD), species richness and community composition. *Ecology Letters*, 5(3), 402-411. https://doi.org/10.1046/j.1461-0248.2002.00339.x

Petchey, O. L., & Gaston, K. J. (2006). Functional diversity: Back to basics and looking forward. *Ecology Letters*, 9(6), 741-758. https://doi.org/10.1111/j.1461-0248.2006.00924.x

Pillar, V. D., Blanco, C. C., Müller, S. C., Sosinski, E. E., Jone, F., & Duarte, L. D. S. (2013). Functional redundancy and stability in plant communities. *Journal of Vegetation Science*, 24(5), 963-974. https://doi.org/10.1111/jvs.12047

Potapov, A. M., Beaulieu, F., Birkhofer, K., Bluhm, S. L., Degtyarev, M. I., Devetter, M., Goncharov, A. A., Gongalsky, K. B., Klarner, B., Korobushkin, D. I., Liebke, D. F., Maraun, M., Mc Donnell, R. J., Pollierier, M. M., Schaefer, I., Shrubovych, J., Semenyuk, I. I., Sendra, A., Tuma, J., ... Scheu, S. (2022). Feeding habits and multifunctional classification of soil-associated consumers from protists to vertebrates. *Biological Reviews*, 97(3), 1057-1117. https://doi.org/10.1111/brv.12832

Pothula, S. K., & Adams, B. J. (2022). Community assembly in the wake of glacial retreat: A meta-analysis. *Global Change Biology*, 28(23), 6973-6991. https://doi.org/10.1111/gcb.16427
<figure><img src="image_14.png" /></figure>

Ptatscheck, C., Gansfort, B., & Traunspurger, W. (2018). The extent of wind-mediated dispersal of small metazoans, focusing nematodes. *Scientific Reports*, 8(1), 6814. https://doi.org/10.1038/s41598-018-24747-8

Purschke, O., Schmid, B. C., Sykes, M. T., Poschlod, P., Michalski, S. G., Durka, W., Kühn, I., Winter, M., & Prentice, H. C. (2013). Contrasting changes in taxonomic, phylogenetic and functional diversity during a long-term succession: Insights into assembly processes. *Journal of Ecology*, 101(4), 857–866. https://doi.org/10.1111/1365-2745.12098

Rosero, P., Crespo-Pérez, V., Espinosa, R., Andino, P., Barragán, Á., Moret, P., Gobbi, M., Ficetola, G. F., Jaramillo, R., Muriel, P., Anthelme, F., Jacobsen, D., Dangles, O., Condom, T., Gielly, L., Poulenard, J., Rabatel, A., Basantes, R., Cáceres Correa, B., & Cauvy-Fraunié, S. (2021). Multi-taxa colonisation along the foreland of a vanishing equatorial glacier. *Ecography*, 44(7), 1010–1021. https://doi.org/10.1111/ecog.05478

Rounce, D. R., Hock, R., Maussion, F., Hugonnet, R., Kochtitzky, W., Huss, M., Berthier, E., Brinkerhoff, D., Compagno, L., Copland, L., Farinotti, D., Menounos, B., & McNabb, R. W. (2023). Global glacier change in the 21st century: Every increase in temperature matters. *Science*, 379(6627), 78–83. https://doi.org/10.1126/science.abo1324

Sechi, V., Goede, R. G. M. D., Rutgers, M., Brussaard, L., & Mulder, C. (2018). Functional diversity in nematode communities across terrestrial ecosystems. *Basic and Applied Ecology*, 30, 76–86. https://doi.org/10.1016/j.baae.2018.05.004

Smithson, M., & Verkuilen, J. (2006). A better lemon squeezer? Maximum-likelihood regression with beta-distributed dependent variables. *Psychological Methods*, 11(1), 54–71. https://doi.org/10.1037/1082-989x.11.1.54

Sytsma, M. L. T., Lewis, T., Bakker, J. D., & Prugh, L. R. (2023). Successional patterns of terrestrial wildlife following deglaciation. *Journal of Animal Ecology*, 92(3), 723–737. https://doi.org/10.1111/1365-2656.13886

Taberlet, P., Bonin, A., Zinger, L., & Coissac, E. (2018). *Environmental DNA for biodiversity research and monitoring*. Oxford University Press.

Taberlet, P., Prud'homme, S. M., Campione, E., Roy, J., Miquel, C., Shehzad, W., Gielly, L., Rioux, D., Choler, P., Clement, J. C., Melodelima, C., Pompanon, F., & Coissac, E. (2012). Soil sampling and isolation of extracellular DNA from large amount of starting material suitable for metabarcoding studies. *Molecular Ecology*, 21(8), 1816–1820. https://doi.org/10.1111/j.1365-294X.2011.05317.x

Treonis, A. M., Unangst, S. K., Kepler, R. M., Buyer, J. S., Cavigelli, M. A., Mirsky, S. B., & Maul, J. E. (2018). Characterization of soil nematode communities in three cropping systems through morphological and DNA metabarcoding approaches. *Scientific Reports*, 8(1), 20366–20366-5

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

van Leeuwen, J. P., Lair, G. J., Gísladóttir, G., Sandén, T., Bloem, J., Hemerik, L., & de Ruiter, P. C. (2018). Soil food web assembly and vegetation development in a glacial chronosequence in Iceland. *Pedobiologia*, 70, 12–21. https://doi.org/10.1016/j.pedobi.2018.08.002

Vater, A. E., & Matthews, J. A. (2015). Succession of pitfall-trapped insects and arachnids on eight Norwegian glacier forelands along an altitudinal gradient: Patterns and models. *The Holocene*, 25(1), 108–129. https://doi.org/10.1177/0959683614556374

Walker, L. R., Wardle, D. A., Bardgett, R. D., & Clarkson, B. D. (2010). The use of chronosequences in studies of ecological succession and soil development: Chronosequences, succession and soil development. *Journal of Ecology*, 98(4), 725–736. https://doi.org/10.1111/j.1365-2745.2010.01664.x

Weiher, E., & Teddy, P. (Eds.). (1999). *Ecological assembly rules: Perspectives, advances, retreats*. Cambridge University Press. https://doi.org/10.1017/CBO9780511542237

White, H. J., León-Sánchez, L., Burton, V. J., Cameron, E. K., Caruso, T., Cunha, L., Dirilgen, T., Jurburg, S. D., Kelly, R., Kumaresan, D., Ochoa-Hueso, R., Ordonez, A., Phillips, H. R. P., Prieto, I., Schmidt, O., & Caplat, P. (2020). Methods and approaches to advance soil macroecology. *Global Ecology and Biogeography*, 29(10), 1674–1690. https://doi.org/10.1111/geb.13156

Wiesel, L., Daniell, T. J., King, D., & Neilson, R. (2015). Determination of the optimal soil sample size to accurately characterise nematode communities in soil. *Soil Biology and Biochemistry*, 80, 89–91. https://doi.org/10.1016/j.soilbio.2014.09.026

Zawierucha, K., Porazinska, D. L., Ficetola, G. F., Ambrosini, R., Baccolo, G., Buda, J., Ceballos, J. L., Devetter, M., Dial, R., Franzetti, A., Fuglewicz, U., Gielly, L., Łokas, E., Janko, K., Jaromerska, T., Koscinski, A., Kozłowska, A., Ono, M., Parnikoza, I., ... Takeuchi, N. (2021). A hole in the nematosphere: Tardigrades and rotifers dominate the cryoconite hole environment, whereas nematodes are missing. *Journal of Zoology*, 313, 18–36. https://doi.org/10.1111/jzo.12832

Zinger, L., Bonin, A., Alsos, I. G., Bálint, M., Bik, H., Boyer, F., Chariton, A. A., Creer, S., Coissac, E., Deagle, B. E., De Barba, M., Dickie, I. A., Dumbrell, A. J., Ficetola, G. F., Fierer, N., Fumagalli, L., Gilbert, M. T. P., Jarman, S., Jumpponen, A., ... Taberlet, P. (2019). DNA metabarcoding—Need for robust experimental designs to draw sound ecological conclusions. *Molecular Ecology*, 28(8), 1857–1862. https://doi.org/10.1111/mec.15060

## SUPPORTING INFORMATION

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

**How to cite this article:** Guerrieri, A., Cantera, I., Marta, S., Bonin, A., Carteron, A., Ambrosini, R., Caccianiga, M., Anthelme, F., Azzoni, R. S., Almond, P., Alviz Gazitúa, P., Cauvy-Fraunié, S., Ceballos Lievano, J. L., Chand, P., Chand Sharma, M., Clague, J., Cochachín Rapre, J. A., Compostella, C., Cruz Encarnación, R. ... Ficetola, G. F. (2023). Local climate modulates the development of soil nematode communities after glacier retreat. *Global Change Biology*, 30, e17057. https://doi.org/10.1111/gcb.17057