domain: sdv sde

text:
Cantera, I Alessia Guerrieri

Isabel Cantera
email: isa_cantera@hotmail.com

Silvio Marta

Aurélie Bonin

Alexis Carteron

Roberto Ambrosini

Marco Caccianiga

Fabien Anthelme

Roberto Sergio Azzoni

Peter Almond

| Aurélie Bonin

| Fabien Anthelme

| Roberto

Sergio Azzoni

| Pablo

Almond, P 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

Norine Khedim

| Rosa

Isela Meneses

Gwendolyn Peyre

Francesca Pittino

Angela Proietto

Antoine Rabatel

| Nurai Urseitova

Yan Yang

Vitalii Zaginaev

Andrea Zerboni

Anaïs Zimmer

Pierre Taberlet

| Guglielmina

Adele Diolaiuti

Jerome Poulenard

Diego Fontaneto

Wilfried Thuiller

| Gentile

Francesco Ficetola

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

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

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

| INTRODUC TI ON

Mountain glaciers are losing mass at an unprecedented rate due to ongoing climate change, and the pace of glacier recession is accelerating [START_REF] Hugonnet | Accelerated global glacier mass loss in the early twenty-first century[END_REF][START_REF] Rounce | Global glacier change in the 21st century: Every increase in temperature matters[END_REF]. 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; [START_REF] Pothula | Community assembly in the wake of glacial retreat: A meta-analysis[END_REF]. 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 [START_REF] Ficetola | Dynamics of ecological communities following current retreat of glaciers[END_REF][START_REF] Pothula | Community assembly in the wake of glacial retreat: A meta-analysis[END_REF][START_REF] Sytsma | Successional patterns of terrestrial wildlife following deglaciation[END_REF]. 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 [START_REF] Crowther | The global soil community and its influence on biogeochemistry[END_REF][START_REF] Ficetola | Dynamics of ecological communities following current retreat of glaciers[END_REF].

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 [START_REF] Cauvy-Fraunié | A global synthesis of biodiversity responses to glacier retreat[END_REF].

Biotic colonization after glacier retreat is generally characterized by an increase in diversity over time, at rates that differ among taxonomic groups [START_REF] Guerrieri | Metabarcoding data reveal vertical multitaxa variation in topsoil communities during the colonization of deglaciated forelands[END_REF][START_REF] Pothula | Community assembly in the wake of glacial retreat: A meta-analysis[END_REF][START_REF] Rosero | Multi-taxa colonisation along the foreland of a vanishing equatorial glacier[END_REF]. Nonetheless, it has been proposed that the colonization dynamics of deglaciated areas are modulated by local climatic conditions [START_REF] Ficetola | Dynamics of ecological communities following current retreat of glaciers[END_REF], given that the impacts of climate change on biodiversity are uneven at the global scale [START_REF] Freeman | Montane species track rising temperatures better in the tropics than in the temperate zone[END_REF]. Climatic differences in temperature and, to a lesser extent, precipitation and solar radiation can influence the evolution of ecosystems after glacier retreat [START_REF] Anthelme | Novel plant communities after glacial retreat in Colombia: (Many) losses and (few) gains[END_REF][START_REF] Franco | Response of Antarctic soil fauna to climate-driven changes since the Last Glacial Maximum[END_REF][START_REF] Kaufmann | Glacier foreland colonisation: Distinguishing between short-term and long-term effects of climate change[END_REF][START_REF] Khedim | Topsoil organic matter build-up in glacier forelands around the world[END_REF]. For instance, the accumulation through time of soil organic matter can be faster in proglacial forelands that experience relatively warm conditions [START_REF] Kaufmann | Glacier foreland colonisation: Distinguishing between short-term and long-term effects of climate change[END_REF][START_REF] Khedim | Topsoil organic matter build-up in glacier forelands around the world[END_REF]. 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 [START_REF] Ficetola | Dynamics of ecological communities following current retreat of glaciers[END_REF][START_REF] Fridley | Temperature accelerates the rate fields become forests[END_REF][START_REF] Vater | Succession of pitfall-trapped insects and arachnids on eight Norwegian glacier forelands along an altitudinal gradient: Patterns and models[END_REF]. However, broadscale 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 [START_REF] Bagousse-Pinguet | Phylogenetic, functional, and taxonomic richness have both positive and negative effects on ecosystem multifunctionality[END_REF]. 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 [START_REF] Cadotte | Beyond species: Functional diversity and the maintenance of ecological processes and services: Functional diversity in ecology and conservation[END_REF][START_REF] Loreau | Biodiversity and ecosystem functioning: Current knowledge and future challenges[END_REF][START_REF] Naeem | The functions of biological diversity in an age of extinction[END_REF], 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 [START_REF] Kuczynski | Community disassembly under global change: Evidence in favor of the stress-dominance hypothesis[END_REF].

Nematodes are the most abundant animals in soils and are excellent indicators of soil biological activity [START_REF] Neher | Role of nematodes in soil health and their use as indicators[END_REF]. Operating at all trophic levels, they play a crucial role in the development of ecosystems, as their diversity affects numerous ecosystem functions, including productivity, nutrient cycling, and decomposition of organic matter [START_REF] Crowther | The global soil community and its influence on biogeochemistry[END_REF][START_REF] Delgado-Baquerizo | Multiple elements of soil biodiversity drive ecosystem functions across biomes[END_REF] Nematodes colonized ice-free areas almost immediately. Both taxonomic and functional 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 between climatic regimes in the long term. Immediately after glacier retreat, communities 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.

| 3 of 18 GUERRIERI et al. [START_REF] Fontaneto | Long-distance passive dispersal in microscopic aquatic animals[END_REF][START_REF] Van Den Hoogen | Soil nematode abundance and functional group composition at a global scale[END_REF][START_REF] Zawierucha | A hole in the nematosphere: Tardigrades and rotifers dominate the cryoconite hole environment, whereas nematodes are missing[END_REF]. Nematodes are particularly important components of soil communities in the first decades after glacier retreat [START_REF] Devetter | Terrestrial invertebrates along a gradient of deglaciation in Svalbard: Long-term development of soil fauna communities[END_REF][START_REF] Doblas-Miranda | Changes in the community structure and diversity of soil invertebrates across the Franz Josef Glacier chronosequence[END_REF][START_REF] Rosero | Multi-taxa colonisation along the foreland of a vanishing equatorial glacier[END_REF], probably due to their small size and resistant propagules that facilitate dispersal and colonization [START_REF] Fontaneto | Long-distance passive dispersal in microscopic aquatic animals[END_REF][START_REF] Ptatscheck | The extent of wind-mediated dispersal of small metazoans, focusing nematodes[END_REF].

Previous studies on nematode colonization after glacier retreat reported a range of successional patterns, with variation among trophic groups [START_REF] Pothula | Community assembly in the wake of glacial retreat: A meta-analysis[END_REF] and study areas. For instance, in [START_REF] Lei | Soil nematode assemblages as bioindicators of primary succession along a 120-year-old chronosequence on the Hailuogou Glacier forefield, SW China[END_REF], most of richness increase occurred during the first 50 years, while colonization required much longer time in [START_REF] Devetter | Terrestrial invertebrates along a gradient of deglaciation in Svalbard: Long-term development of soil fauna communities[END_REF]. Furthermore, some studies showed linear increases in diversity over time [START_REF] Devetter | Terrestrial invertebrates along a gradient of deglaciation in Svalbard: Long-term development of soil fauna communities[END_REF], while others suggested intermediate peaks followed by plateaus or even richness declines [START_REF] Doblas-Miranda | Changes in the community structure and diversity of soil invertebrates across the Franz Josef Glacier chronosequence[END_REF][START_REF] Lei | Soil nematode assemblages as bioindicators of primary succession along a 120-year-old chronosequence on the Hailuogou Glacier forefield, SW China[END_REF][START_REF] Van Leeuwen | Soil food web assembly and vegetation development in a glacial chronosequence in Iceland[END_REF]. 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 [START_REF] Taberlet | Environmental DNA for biodiversity research and monitoring[END_REF], enable rapid, cost-effective, and consistent estimates of soil biodiversity [START_REF] Delgado-Baquerizo | Multiple elements of soil biodiversity drive ecosystem functions across biomes[END_REF][START_REF] Geisen | Integrating quantitative morphological and qualitative molecular methods to analyse soil nematode community responses to plant range expansion[END_REF] even over broad spatial scales [START_REF] White | Methods and approaches to advance soil macroecology[END_REF], which would otherwise be prohibitive. We used a standardized eDNA-based approach to sample soil nematode communities along the proglacial chronosequences [START_REF] Walker | The use of chronosequences in studies of ecological succession and soil development: Chronosequences, succession and soil development[END_REF] 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 [START_REF] Erschbamer | Glacier forelands: Lessons of plant population and community development[END_REF][START_REF] Ficetola | Dynamics of ecological communities following current retreat of glaciers[END_REF]. 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 [START_REF] Gotelli | Null models in ecology[END_REF]. 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) [START_REF] Götzenberger | Ecological assembly rules in plant communities-approaches, patterns and prospects[END_REF]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 [START_REF] Purschke | Contrasting changes in taxonomic, phylogenetic and functional diversity during a long-term succession: Insights into assembly processes[END_REF].

| MATERIAL S AND ME THODS

| 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 [START_REF] Marta | The retreat of mountain glaciers since the little ice age: A spatially explicit database[END_REF]. 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 [START_REF] Marta | The retreat of mountain glaciers since the little ice age: A spatially explicit database[END_REF].

At each site, we established an average of five regularly spaced plots (SD = 0.6) for more repeatable biodiversity estimates [START_REF] Dickie | Towards robust and repeatable sampling methods in eDNA based studies[END_REF], 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 [START_REF] Doblas-Miranda | Changes in the community structure and diversity of soil invertebrates across the Franz Josef Glacier chronosequence[END_REF], 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 [START_REF] Guerrieri | Effects of soil preservation for biodiversity monitoring using environmental DNA[END_REF]. Desiccation is a cost-effective approach allowing long-term preservation of eDNA for biodiversity assessments [START_REF] Guerrieri | Effects of soil preservation for biodiversity monitoring using environmental DNA[END_REF]. 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 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).

| 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 [START_REF] Wiesel | Determination of the optimal soil sample size to accurately characterise nematode communities in soil[END_REF]. Before extraction, we added a preliminary step where the soil was mixed with 20 mL of phosphate buffer for 15 min [START_REF] Taberlet | Soil sampling and isolation of extracellular DNA from large amount of starting material suitable for metabarcoding studies[END_REF]. 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; [START_REF] Zinger | DNA metabarcoding-Need for robust experimental designs to draw sound ecological conclusions[END_REF].

DNA amplifications were performed using the universal Euka02 F I G U R E 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.

primer [START_REF] Guardiola | Deep-sea, deep-sequencing: Metabarcoding extracellular DNA from sediments of marine canyons[END_REF][START_REF] Taberlet | Environmental DNA for biodiversity research and monitoring[END_REF], which amplifies an ~120 bp fragment of the V7 region of the 18S rDNA gene.

The short fragment length is well suited for eDNA analyses [START_REF] Ficetola | Comparison of markers for the monitoring of freshwater benthic biodiversity through DNA metabarcoding[END_REF][START_REF] Guardiola | Deep-sea, deep-sequencing: Metabarcoding extracellular DNA from sediments of marine canyons[END_REF][START_REF] Taberlet | Environmental DNA for biodiversity research and monitoring[END_REF].

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 [START_REF] Ficetola | In silico assessment of 18S rDNA metabarcoding markers for the characterization of nematode communities[END_REF]. Reverse and forward primers included 8-nucleotidelong 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 [START_REF] Taberlet | Environmental DNA for biodiversity research and monitoring[END_REF]. 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 [START_REF] Zinger | DNA metabarcoding-Need for robust experimental designs to draw sound ecological conclusions[END_REF]. PCRnegative 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 [START_REF] Guerrieri | Metabarcoding data reveal vertical multitaxa variation in topsoil communities during the colonization of deglaciated forelands[END_REF].

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 [START_REF] Ficetola | Replication levels, false presences and the estimation of the presence/absence from eDNA metabarcoding data[END_REF]. The purified amplicons were sent to Fasteris (Switzerland) for library preparation using the MetaFast protocol [START_REF] Taberlet | Environmental DNA for biodiversity research and monitoring[END_REF] and sequencing with the Illumina HiSeq 2500 platform (paired-end approach, 2 × 150 bp; about 0.5% of PhiX added in each sequencing lane). See [START_REF] Guerrieri | Metabarcoding data reveal vertical multitaxa variation in topsoil communities during the colonization of deglaciated forelands[END_REF] for additional details on molecular analyses.

We performed the bioinformatic treatment of sequence data using the OBITools software suite version 1.2.9 [START_REF] Boyer | OBITools: A unix-inspired software package for DNA metabarcoding[END_REF], following the steps described in [START_REF] Guerrieri | Metabarcoding data reveal vertical multitaxa variation in topsoil communities during the colonization of deglaciated forelands[END_REF], and then clustered sequences at a threshold of 96% sequence similarity using the sumaclust program (https:// git. metab arcod ing. org/ obito ols/ sumac lust/ wikis/ home). The 96% threshold minimizes the risk that sequences belonging to the same species are clustered in different molecular operational taxonomic units (MOTUs) [START_REF] Bonin | Optimal sequence similarity thresholds for clustering of molecular operational taxonomic units in DNA metabarcoding studies[END_REF]. 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 [START_REF] Ficetola | An in silico approach for the evaluation of DNA barcodes[END_REF]. 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 [START_REF] Boyer | OBITools: A unix-inspired software package for DNA metabarcoding[END_REF]. 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) [START_REF] Boyer | OBITools: A unix-inspired software package for DNA metabarcoding[END_REF]. 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 [START_REF] Guerrieri | Metabarcoding data reveal vertical multitaxa variation in topsoil communities during the colonization of deglaciated forelands[END_REF][START_REF] Zinger | DNA metabarcoding-Need for robust experimental designs to draw sound ecological conclusions[END_REF]. 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 [START_REF] Ficetola | Replication levels, false presences and the estimation of the presence/absence from eDNA metabarcoding data[END_REF][START_REF] Zinger | DNA metabarcoding-Need for robust experimental designs to draw sound ecological conclusions[END_REF], 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 [START_REF] Davis | Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data[END_REF].

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).

| 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 [START_REF] Fe R E N C E S Aglieri | Environmental DNA effectively captures functional diversity of coastal fish communities[END_REF][START_REF] White | Methods and approaches to advance soil macroecology[END_REF].

We extracted functional traits from the NEMAPLEX database (http:// nemap lex. ucdav is. edu/ ). We considered one life-history index (colonizer-persister index, hereafter "CP value") and eight consumption traits (parasites of animals, parasites of plants, bacterivores, fungivores, 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 [START_REF] Bongers | The maturity index: An ecological measure of environmental disturbance based on nematode species composition[END_REF][START_REF] Ferris | Form and function: Metabolic footprints of nematodes in the soil food web[END_REF][START_REF] Ferris | A framework for soil food web diagnostics: Extension of the nematode faunal analysis concept[END_REF]. 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 [START_REF] Colwell | Estimating terrestrial biodiversity through extrapolation[END_REF][START_REF] Gotelli | Estimating species richness[END_REF]. 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 [START_REF] Chazdon | Statistical methods for estimating species richness of woody regeneration in primary and secondary rain forests of NE Costa Rica[END_REF][START_REF] Colwell | Estimating terrestrial biodiversity through extrapolation[END_REF][START_REF] Gotelli | Estimating species richness[END_REF]. 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 [START_REF] Gotelli | Estimating species richness[END_REF]. 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 [START_REF] Cardoso | BAT -Biodiversity Assessment Tools, an R package for the measurement and estimation of alpha and beta taxon, phylogenetic and functional diversity[END_REF].

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 [START_REF] Taberlet | Environmental DNA for biodiversity research and monitoring[END_REF]. 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 [START_REF] Mouillot | A functional approach reveals community responses to disturbances[END_REF]. The assessment of functional diversity requires the analysis of multiple metrics because the functional attributes can impact differently ecosystem dynamics [START_REF] Grenié | funrar: An R package to characterize functional rarity[END_REF][START_REF] Bagousse-Pinguet | Phylogenetic, functional, and taxonomic richness have both positive and negative effects on ecosystem multifunctionality[END_REF]. We assessed functional diversity using functional richness and Rao's functional diversity [START_REF] Mason | Functional richness, functional evenness and functional divergence: The primary components of functional diversity[END_REF][START_REF] Petchey | Functional diversity: Back to basics and looking forward[END_REF], 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 [START_REF] Mason | Functional richness, functional evenness and functional divergence: The primary components of functional diversity[END_REF][START_REF] Petchey | Functional diversity (FD), species richness and community composition[END_REF] and was assessed using the adaptation of Faith's phylogenetic diversity developed in [START_REF] Petchey | Functional diversity (FD), species richness and community composition[END_REF]. We first estimated pairwise trait dissimilarities with Gower distances, with equal weights for each trait [START_REF] Laliberté | A distance-based framework for measuring functional diversity from multiple traits[END_REF][START_REF] Petchey | Functional diversity (FD), species richness and community composition[END_REF]. 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 [START_REF] Cardoso | BAT -Biodiversity Assessment Tools, an R package for the measurement and estimation of alpha and beta taxon, phylogenetic and functional diversity[END_REF] 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 [START_REF] Mouchet | Functional diversity measures: An overview of their redundancy and their ability to discriminate community assembly rules: Functional diversity measures[END_REF][START_REF] Petchey | Functional diversity: Back to basics and looking forward[END_REF]). 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 [START_REF] De Bello | Partitioning of functional diversity reveals the scale and extent of trait convergence and divergence[END_REF]. Using the proportion of reads can yield diversity measures that are less sensitive to parameters that are used for data filtering [START_REF] Calderón-Sanou | From environmental DNA sequences to ecological conclusions: How strong is the influence of methodological choices?[END_REF][START_REF] Mächler | Decision-making and best practices for taxonomy-free environmental DNA metabarcoding in biomonitoring using Hill numbers[END_REF]. 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).

| 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 [START_REF] Karger | Climatologies at high resolution for the earth's land surface areas[END_REF]. ChelsaClim is among the datasets that best describe climatic conditions in mountainous areas [START_REF] Karger | Climatologies at high resolution for the earth's land surface areas[END_REF] 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) [START_REF] Marta | Heterogeneous changes of soil microclimate in high mountains and glacier forelands[END_REF]. 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., [START_REF] Lembrechts | Global maps of soil temperature[END_REF][START_REF] Marta | Heterogeneous changes of soil microclimate in high mountains and glacier forelands[END_REF]. Therefore, during winter, meso-climate does not represent conditions experienced by soil organisms whereas, in summer, it provides an excellent estimate of soil conditions [START_REF] Marta | Heterogeneous changes of soil microclimate in high mountains and glacier forelands[END_REF]. 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 [START_REF] Khedim | Topsoil organic matter build-up in glacier forelands around the world[END_REF].

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 [START_REF] Gibbard | Quaternary sediments[END_REF]. TOC was measured on bulk samples after drying, using an OEA Flash200 ThermoFisher elemental analyzer (see [START_REF] Khedim | Topsoil organic matter build-up in glacier forelands around the world[END_REF]. 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 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 [START_REF] Gelman | Understanding predictive information criteria for Bayesian models[END_REF]. 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 [START_REF] Bürkner | brms: An R package for Bayesian multilevel models using Stan[END_REF]. 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 [START_REF] Smithson | A better lemon squeezer? Maximumlikelihood regression with beta-distributed dependent variables[END_REF] 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 [START_REF] Bahram | Structure and function of the global topsoil microbiome[END_REF]. 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 [START_REF] Ficetola | Dynamics of ecological communities following current retreat of glaciers[END_REF]. 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 [START_REF] Hanusch | Succession comprises a sequence of threshold-induced community assembly processes towards multidiversity[END_REF]. We thus used segmented regressions to check the occurrence of thresholds in the relationships between taxonomic richness and time since glacier retreat [START_REF] Ficetola | Ecological thresholds: An assessment of methods to identify abrupt changes in species-habitat relationships[END_REF][START_REF] Muggeo | Estimating regression models with unknown break-points[END_REF]. We used maximum likelihood to build linear mixed models with one breakpoint (segmented package in R; [START_REF] Muggeo | Estimating regression models with unknown break-points[END_REF]. 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; [START_REF] Mouillot | A functional approach reveals community responses to disturbances[END_REF]. 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 [START_REF] Borcard | Numerical ecology with R[END_REF][START_REF] Oksanen | Package 'vegan'. Community ecology package[END_REF]. 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 2 ) 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 2 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) [START_REF] Götzenberger | Ecological assembly rules in plant communities-approaches, patterns and prospects[END_REF]Weiher & Keddy, 1999), we compared the observed values of Rao's functional diversity to those obtained from random community assemblies [START_REF] Gotelli | Null models in ecology[END_REF].

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 [START_REF] Gotelli | Species co-occurrence: A meta-analysis of J. M. Diamond's assembly rules model[END_REF]. 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 [START_REF] Cardoso | BAT -Biodiversity Assessment Tools, an R package for the measurement and estimation of alpha and beta taxon, phylogenetic and functional diversity[END_REF]. This function partitions total 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-proje ct-org).

| RE SULTS

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).

| 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; 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 (Hemicycliophora_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).

| 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).

| 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).

| DISCUSS ION

Understanding how biodiversity components develop under different climates and assessing the variation of functional groups in response 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 [START_REF] Kaufmann | Glacier foreland colonisation: Distinguishing between short-term and long-term effects of climate change[END_REF][START_REF] Khedim | Topsoil organic matter build-up in glacier forelands around the world[END_REF].

F I G U R E 3
I G U R E 4 (a) (e) (f) (g) (h) (b) (c) (d)
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.

| Mesoclimate affects the development of nematode communities

Nematodes were able to colonize deglaciated terrains very soon (in 4-15 years; see also [START_REF] Lei | Soil nematode assemblages as bioindicators of primary succession along a 120-year-old chronosequence on the Hailuogou Glacier forefield, SW China[END_REF]. As nematodes are generally absent from supra-glacial environments [START_REF] Zawierucha | A hole in the nematosphere: Tardigrades and rotifers dominate the cryoconite hole environment, whereas nematodes are missing[END_REF], their rapid colonization might be related to their excellent dispersal abilities [START_REF] Fontaneto | Long-distance passive dispersal in microscopic aquatic animals[END_REF]. 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; [START_REF] Lei | Soil nematode assemblages as bioindicators of primary succession along a 120-year-old chronosequence on the Hailuogou Glacier forefield, SW China[END_REF], others suggested that nematode establishment may take one to several decades [START_REF] Devetter | Terrestrial invertebrates along a gradient of deglaciation in Svalbard: Long-term development of soil fauna communities[END_REF][START_REF] Van Leeuwen | Soil food web assembly and vegetation development in a glacial chronosequence in Iceland[END_REF]. These studies analyzed different areas of the world; slow colonization was observed in subarctic regions [START_REF] Devetter | Terrestrial invertebrates along a gradient of deglaciation in Svalbard: Long-term development of soil fauna communities[END_REF][START_REF] Van Leeuwen | Soil food web assembly and vegetation development in a glacial chronosequence in Iceland[END_REF], whereas fast colonization was observed in a warmer region [START_REF] Lei | Soil nematode assemblages as bioindicators of primary succession along a 120-year-old chronosequence on the Hailuogou Glacier forefield, SW China[END_REF]. 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 [START_REF] Purschke | Contrasting changes in taxonomic, phylogenetic and functional diversity during a long-term succession: Insights into assembly processes[END_REF]. Communities exhibiting high functional diversity often indicate limiting similarity due to niche differentiation [START_REF] Macarthur | The limiting similarity, convergence, and divergence of coexisting species[END_REF]. This can be underlined by competition for resource-rich conditions [START_REF] Bertness | Positive interactions in communities[END_REF][START_REF] Losapio | The consequences of glacier retreat are uneven between plant species[END_REF][START_REF] Pothula | Community assembly in the wake of glacial retreat: A meta-analysis[END_REF], as expected in late communities [START_REF] Ficetola | Dynamics of ecological communities following current retreat of glaciers[END_REF]. 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 [START_REF] Franco | Drought suppresses soil predators and promotes root herbivores in mesic, but not in xeric grasslands[END_REF][START_REF] Ma | Climatic warming shifts the soil nematode community in a desert steppe[END_REF][START_REF] Van Den Hoogen | Soil nematode abundance and functional group composition at a global scale[END_REF]. 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 [START_REF] Münkemüller | Dos and don'ts when inferring assembly rules from diversity patterns[END_REF].

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 [START_REF] De Bello | Functional species pool framework to test for biotic effects on community assembly[END_REF], 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 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 [START_REF] Freeman | Montane species track rising temperatures better in the tropics than in the temperate zone[END_REF] (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 [START_REF] Lei | Soil nematode assemblages as bioindicators of primary succession along a 120-year-old chronosequence on the Hailuogou Glacier forefield, SW China[END_REF][START_REF] Losapio | The consequences of glacier retreat are uneven between plant species[END_REF]. The biomass and diversity of below-ground animals, including nematodes, are tightly linked to the accumulation of organic matter [START_REF] Fan | Soil biodiversity supports the delivery of multiple ecosystem functions in urban greenspaces[END_REF][START_REF] Fierer | Global patterns in belowground communities[END_REF][START_REF] Lei | Soil nematode assemblages as bioindicators of primary succession along a 120-year-old chronosequence on the Hailuogou Glacier forefield, SW China[END_REF][START_REF] Van Den Hoogen | Soil nematode abundance and functional group composition at a global scale[END_REF], which can be particularly fast under a mild summer climate [START_REF] Khedim | Topsoil organic matter build-up in glacier forelands around the world[END_REF].

F I G U R E
Nematode colonization immediately after glacier retreat was fostered by mild climates with a pattern very similar to the one observed for soil organic matter [START_REF] Khedim | Topsoil organic matter build-up in glacier forelands around the world[END_REF], 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 is a major driver for 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 [START_REF] Doblas-Miranda | Changes in the community structure and diversity of soil invertebrates across the Franz Josef Glacier chronosequence[END_REF][START_REF] Lei | Soil nematode assemblages as bioindicators of primary succession along a 120-year-old chronosequence on the Hailuogou Glacier forefield, SW China[END_REF][START_REF] Pothula | Community assembly in the wake of glacial retreat: A meta-analysis[END_REF]. For instance, [START_REF] Doblas-Miranda | Changes in the community structure and diversity of soil invertebrates across the Franz Josef Glacier chronosequence[END_REF] 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 [START_REF] Losapio | The consequences of glacier retreat are uneven between plant species[END_REF].

| 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 and4) 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 [START_REF] Cauvy-Fraunié | A global synthesis of biodiversity responses to glacier retreat[END_REF][START_REF] Erschbamer | Glacier forelands: Lessons of plant population and community development[END_REF][START_REF] Ferris | A framework for soil food web diagnostics: Extension of the nematode faunal analysis concept[END_REF][START_REF] Hågvar | Ecosystem birth near melting glaciers: A review on the pioneer role of ground-dwelling arthropods[END_REF][START_REF] Hågvar | The role of arthropods in early colonization near melting glaciers: Contradictions between ecological assumptions and recent study results[END_REF]. 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 [START_REF] Doblas-Miranda | Changes in the community structure and diversity of soil invertebrates across the Franz Josef Glacier chronosequence[END_REF][START_REF] Pothula | Community assembly in the wake of glacial retreat: A meta-analysis[END_REF][START_REF] Van Leeuwen | Soil food web assembly and vegetation development in a glacial chronosequence in Iceland[END_REF].

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 [START_REF] Naeem | Species redundancy and ecosystem reliability[END_REF][START_REF] Pillar | Functional redundancy and stability in plant communities[END_REF]. 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 [START_REF] Bongers | The maturity index: An ecological measure of environmental disturbance based on nematode species composition[END_REF][START_REF] Mulder | Unifying the functional diversity in natural and cultivated soils using the overall body-mass distribution of nematodes[END_REF][START_REF] Sechi | Functional diversity in nematode communities across terrestrial ecosystems[END_REF].

| 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 [START_REF] Ficetola | Dynamics of ecological communities following current retreat of glaciers[END_REF][START_REF] Marta | Heterogeneous changes of soil microclimate in high mountains and glacier forelands[END_REF]. 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 [START_REF] Ficetola | In silico assessment of 18S rDNA metabarcoding markers for the characterization of nematode communities[END_REF]. 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 [START_REF] Taberlet | Environmental DNA for biodiversity research and monitoring[END_REF]. 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 [START_REF] Mulder | Nematode traits and environmental constraints in 200 soil systems: Scaling within the 60-6000 μm body size range[END_REF]. Primers amplifying longer fragments could provide a more accurate taxonomic identification and enable better differentiation of closely related species [START_REF] Ficetola | In silico assessment of 18S rDNA metabarcoding markers for the characterization of nematode communities[END_REF][START_REF] Gattoni | 18S-NemaBase: Curated 18S rRNA database of nematode sequences[END_REF]. 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 [START_REF] Doblas-Miranda | Changes in the community structure and diversity of soil invertebrates across the Franz Josef Glacier chronosequence[END_REF][START_REF] Lei | Soil nematode assemblages as bioindicators of primary succession along a 120-year-old chronosequence on the Hailuogou Glacier forefield, SW China[END_REF]. 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 underdetection, we used statistical approaches that consider the possibility of missing rare taxa [START_REF] Chazdon | Statistical methods for estimating species richness of woody regeneration in primary and secondary rain forests of NE Costa Rica[END_REF][START_REF] Colwell | Estimating terrestrial biodiversity through extrapolation[END_REF][START_REF] Gotelli | Estimating species richness[END_REF] and verified that our results remained consistent 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 [START_REF] Cantera | The importance of species addition versus replacement varies over succession in plant communities after glacial retreat[END_REF][START_REF] Treonis | Characterization of soil nematode communities in three cropping systems through morphological and DNA metabarcoding approaches[END_REF], and the observed changes in nematode communities are generally consistent with the patterns observed by both local-scale and meta-analytic studies [START_REF] Calderón-Sanou | From environmental DNA sequences to ecological conclusions: How strong is the influence of methodological choices?[END_REF][START_REF] Devetter | Terrestrial invertebrates along a gradient of deglaciation in Svalbard: Long-term development of soil fauna communities[END_REF][START_REF] Ji | Reliable, verifiable and efficient monitoring of biodiversity via metabarcoding[END_REF][START_REF] Lei | Soil nematode assemblages as bioindicators of primary succession along a 120-year-old chronosequence on the Hailuogou Glacier forefield, SW China[END_REF][START_REF] Meyer | 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[END_REF][START_REF] Pothula | Community assembly in the wake of glacial retreat: A meta-analysis[END_REF]. 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 [START_REF] Mestre | Human disturbances affect the topology of food webs[END_REF].

| 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 [START_REF] Hock | High mountain areas[END_REF]. 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 [START_REF] Cardinale | Biodiversity loss and its impact on humanity[END_REF][START_REF] Naeem | The functions of biological diversity in an age of extinction[END_REF]. 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; [START_REF] Giachello | Toward a common set of functional traits for soil protists[END_REF][START_REF] Potapov | Feeding habits and multifunctional classification of soil-associated consumers from protists to vertebrates[END_REF].

  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).

Figure

  Figure 2a,b; TableS6). Furthermore, we detected a strong interaction

F

  Changes in taxonomic and functional diversity of nematodes with age of deglaciated terrains. (a) Taxonomic richness (firstorder 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.

  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[START_REF] Martinez-Almoyna | Multi-trophic β-diversity mediates the effect of environmental gradients on the turnover of multiple ecosystem functions[END_REF]. Due to ongoing climate change, newly deglaciated terrains are expected to be increasingly important ecosystems in this century(Intergovernmental Panel on Climate Change, 2021;[START_REF] Rounce | Global glacier change in the 21st century: Every increase in temperature matters[END_REF]. 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. AUTH O R CO NTR I B UTI O N S 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: Writingreview 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. Pritam Chand: 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; writingreview 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; writingreview 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. A FFI LI ATI O N S 1 Dipartimento di Scienze e Politiche Ambientali, Università degli Studi di Milano, Milano, Italy

  13652486, 2024, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.17057 by Cochrane France, Wiley Online Library on [08/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  13652486, 2024, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.17057 by Cochrane France, Wiley Online Library on [08/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  13652486, 2024, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.17057 by Cochrane France, Wiley Online Library on [08/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  13652486, 2024, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.17057 by Cochrane France, Wiley Online Library on [08/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

	GUERRIERI et al.	| 11 of 18

  13652486, 2024, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.17057 by Cochrane France, Wiley Online Library on [08/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

13652486, 2024, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.17057 by Cochrane France, Wiley Online Library on [08/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Argaly, Bâtiment CleanSpace, Sainte-Hélène-du-Lac, France

Institute of Geosciences and Earth Resources, CNR, Pisa, Italy

Dipartimento di Bioscienze, Università degli Studi di Milano, Milano, Italy 13652486, 2024, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.17057 by Cochrane France, Wiley Online Library on [08/01/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

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

DATA AVA I L A B I L I T Y S TAT E M E N T

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

ACK N OWLED G M ENTS

CO N FLI C T O F I NTER E S T S TATEM ENT

The authors declare no competing interests.

S U PP O RTI N G I N FO R M ATI O N

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