The Abundance, Diversity, and Metabolic Footprint of Soil Nematodes Is Highest in High Elevation Alpine Grasslands

Nematodes are key components of soil biodiversity and represent valuable bio-indicators of soil food webs. Numerous community indices have been developed in order to track variations in soil ecosystem processes, but their use is mainly restricted to anthropogenic stresses. In this study, we propose to expand the use of nematodes’ derived ecological indices in order to shed light on variations of soil food webs in natural systems distributed along elevation gradients. For this purpose, we aimed at determining how elevation affects the community structure and the trophic diversity by studying the abundance, the composition and the functional diversity of nematode communities. Nematode communities were sampled every 200 m across five transects that span about 2000 m in elevation in the Alps. To understand the underlying ecological parameters driving these patterns we studied both abiotic factors (soil properties) and biotic factors (trophic links, relationships with plant diversity). We found that (1) nematode abundance increases with elevation of lowland forests and alpine meadows; (2) differences in nematodes communities rely on habitat-specific functional diversity (e.g. tolerance to harsh environments, “colonizer/persister” status) while most trophic groups are ubiquitous; and (3) the metabolic footprint of the complete nematode community increases with elevation. We thus conclude that the contribution of soil dwelling nematodes to belowground ecosystem processes, including carbon and energy flow, is stronger at high elevation. The resulting cascading effects on the soil food web structure are discussed from an ecosystem functioning perspective. Overall, this study highlights the importance of nematodes in soil ecosystems and brings insights in their enhanced role along ecological gradients.

Nematodes are key components of soil biodiversity and represent valuable bio-indicators of soil food webs. Numerous community indices have been developed in order to track variations in nematode-mediate soil ecosystem processes, but their use is mainly restricted to anthropogenic stresses. In this study, we propose to expand the use of nematodes' derived ecological indices in order to shed light on variations of soil food webs in natural systems distributed along elevation gradients. For this purpose, we aimed at determining how elevation affects the community structure and the trophic diversity by studying the abundance, the composition and the functional diversity of nematode communities. Nematode communities were sampled every 200 m across five transects that span about 2000 m in elevation in the Alps. To understand the underlying ecological parameters driving these patterns we studied both abiotic factors (soil properties) and biotic factors (trophic links, relationships with plant diversity). We found that (1) nematode abundance increases with elevation of lowland forests and alpine meadows; (2) differences in nematodes communities rely on habitat-specific functional diversity (e.g., tolerance to harsh environments, "colonizer/persister" status) while most trophic groups are ubiquitous; and (3) the metabolic footprint of the complete nematode community increases with elevation. We thus conclude that the contribution of soil dwelling nematodes to belowground ecosystem processes, including carbon and energy flow, is stronger at high elevation. The resulting cascading effects on the soil food web structure are discussed from an ecosystem functioning perspective. Overall, this study highlights the importance of nematodes in soil ecosystems and brings insights on their functional role along ecological gradients.
Keywords: elevation gradient, entomopathogenic nematodes, nematophagous fungi, plant-herbivore interaction, soil ecosystem functioning INTRODUCTION It has been estimated that under the earth's surface, the myriad of soil habitats shelter about 25% of the worldwide described species, thus providing crucial reservoirs of biodiversity and subsequent ecosystem functioning (Fitter et al., 2005;Decaëns, 2010; Bardgett and van der Putten, 2014). While research on soil biota continues to bear inherent challenges, the combination of traditional research with genomic tools has accelerated the exploration of soil diversity and our understanding of ecosystem dynamics (Johnson et al., 2007). In addition, numerous studies have become increasingly focused on replacing soil diversity within trophic interactions for unraveling soil ecosystem processes (Bardgett and van der Putten, 2014). Indeed, soil fauna is essential for ecosystem functioning through different processes, such as primary production and nutrient cycling of carbon, phosphorous, or nitrogen (Brussaard, 1997). The role of soil functional diversity in the decomposition of organic matter and, more importantly, in the assimilation of carbon in food webs, governs energy flows worldwide (Hunt and Wall, 2002;Krumins et al., 2013).
Several groups of soil-dwelling organisms (e.g., bacteria, fungi, protists, collembolan, enchytraeid worms or earthworms) can partition their task in order to optimize trophic interactions and energy flow. In addition, among soil inhabitants, the group of roundworms (i.e., nematodes; phylum Nematoda) is a key component of the belowground living mosaic. Indeed, nematodes, with more than 14,000 described species, are distributed in almost every habitat on Earth, and represent more than 80% of metazoan taxonomic and functional diversity in soils (Bongers and Bongers, 1998;Hodda et al., 2009). Nematodes can be assigned to basically all functional trophic guilds, and span the whole gamut of ecological adaptations, ranging from "colonizer" (r strategists) to "persister" (K strategists) along a colonizerpersister ("cp") scale (Bongers, 1990). Besides the diversity in life history traits, nematodes sustain a large range of trophic groups and eight feeding types have been described: herbivore, fungivore, bacterivore, substrate ingester, predator of animals, unicellular eukaryote feeder, parasites, and omnivore . The combination of both cp groups and feeding habits provides a wide diversity of functional guilds. In addition, nematodes occupy a central position in soil food-webs by linking microbial communities with macrofauna. Hence, nematodes are widely used as appropriate bioindicators to track changes in the environment and the resulting cascading effects on soil foodweb structure (Sochová et al., 2006;Wilson and Kakouli-Duarte, 2009). Several community and metabolic footprints indices have been developed in order to assess how nematode communities affect (or are affected by) soil quality (Bongers and Ferris, 1999;Ferris et al., 2001;Ferris, 2010), although such studies remain mostly restricted to anthropogenic systems (Salamún et al., 2014;Zhao et al., 2015). Here, we propose to expand the use of nematodes' derived ecological indices for increasing our understanding soil-driven ecosystem functioning along natural ecological gradients.
Studying the causes and consequences of species abundance and distribution along environmental clines remains crucial for providing insights into community assembly and ecosystem functioning (Gaston, 2000;Doherty et al., 2011;Oliver et al., 2015). In this context, ecological gradients act as potent environmental filters and thus provide powerful tools for dissecting biotic and abiotic factors driving species diversity and ecosystem dynamics. For instance, elevation gradients have been classically used to develop key ecological concepts such as the niche theory or the species-energy hypothesis (Grinell, 1917;Brown, 1971;Lomolino, 2001). More recently, various authors have considered elevation gradients as promising "natural experiments" to test evolutionary hypotheses in species nichebreadth or predict plant adaptation to changing environment (Körner, 2007;Alexander et al., 2015;Rasmann and Pellissier, 2015). Indeed, mountain slopes present strong variation in both biotic and abiotic factors that can to alter ecological niches, abundance in species population, or community assemblage (Hodkinson, 2005). In addition, biotic variations occur over short distances, thereby limiting the confounding effect of phylogeography when studying inter-specific interactions from a comparative ecology approach (Rasmann et al., 2014). While, numerous studies have demonstrated the ability of nematodes to colonize the harshest environments, such as the polar regions (Loof, 1971;Yeates, 2010), only few studies have been interested in studying their distribution along elevation, and to our knowledge, none of them have assessed changes in nematode communities along continuous elevation gradients (Hoschitz and Kaufmann, 2004).
With the present work, we aimed at unraveling the community ecology of soil-dwelling nematodes along steep elevation gradients, from the colline regions up to the Alpine grasslands. Specifically, we hypothesized: (1) a decrease in nematode abundance at high elevation following classic views on biodiversity gradients (McCain and Grytnes, 2010); (2) changes in the nematode communities' structures according to variations of ecological niches along the gradient; and (3) changes in the nematode-mediated metabolic footprint indices along the elevation gradient of the Alps.

Study Site
To dissect nematode food web structure along elevation gradients, between July and August 2013, we sampled soils ranging from 700 m above sea level (asl) up to 2700 m asl across five transects in the Swiss Alps ( Figure S1). The five transects were collected over 130 km in order to assess variations in nematode communities over a large scale in Alpine systems. Along each elevational transect, sampling sites of 2 × 2 m were chosen approximately separated from each other by an elevation of 200 m asl (n = 48 sites, Table S1). As we aimed to measure soil diversity in the most pristine conditions, we sampled within the climacic vegetation at each site (Delarze et al., 2015). In lowlands, soil samples were predominantly collected within Fagus sylvatica, Quercus spp., or Castanea sativa dominated forests. Sites in the mountain and the subalpine belts were mainly collected in F. sylvatica, Pinus sylvestris, Abies alba, or Picea abies dominated forests, while sampling in the alpine zone was done in Alpine grasslands found above the timberline ( Figure S1, Table S1). All along the elevation gradient we avoided cultivated, urban, or heavily grazed areas and selected sites with a similar exposition and slope for a same transect.
Sampling sites of three transects out of the five (i.e., 28 sites along the Mont d'Or, Salgesch, and Vallon de Nant transects, Table S1) were selected according to the study conducted by Pellissier et al. (2010), who described plant communities (abundance of plant taxa based on Braun-Blanquet categories) within a 40 m 2 (grasslands) or 250 m 2 (forests) quadrat at each site (Vittoz and Guisan, 2007). Hence, over three transects, each site was also described with the corresponding plant species list and plant cover estimation for each species. This allowed testing for potential spatial correlations between nematode species and the local flora (see below).

Soil and Nematode Sampling
At each site, we randomly collected 10-30 soil cores of 5 cm diameter within a 2 × 2 m area and with a maximal depth of 30 cm till we reached a total of 1.5 Kg of soil after the removal of all rock particles bigger than 2 cm in diameter. The bulk soil was homogenized before dividing it into several subsamples.
Initially, 300 g of the bulk soil were used for measuring soil traits (soil humidity, pH, conductivity and root fraction, i.e., percent root biomass). For soil humidity, we calculated the difference between soil fresh weight and soil dry weight after 7 days at 70 • C. Both pH and conductivity were measured using a 914 pH/Conductometer (Metrohn, Herisau, Switzerland) after mixing 50 g of this subsample with 100 ml of deionized water. Finally, under the microscope, we visually separated all discernible root fragments from the other soil components and weighed them to obtain the proportion in percent of root biomass, calculated as the ratio of root biomass on total soil mass.
Next, out of the initial soil bulk, a sub-sample of 200 g of fresh soil was used for extracting soil nematodes: bacterivores, fungivores, herbivores, predators, and omnivores. For this purpose, we used the sieving and Baermann funnel method (Barker, 1985). All free-living nematodes in each sample were then counted under the dissecting microscope and mounted into a slide. At least 100 nematodes in each sample were then identified under a dissecting microscope to family or genus level, and assigned to a functional guild based on their trophic group and life-histories .
Finally, in order to improve the description of nematode communities, we also performed targeted genomic approach on 200 additional grams of fresh soil to identify entomopathogenic nematodes (EPNs; i.e., nematodes parasites of invertebrates that are only free-living in the soil as the third instar infective juvenile state), and nematophagous fungi (NF) which are important in shaping nematode communities but remain virtually impossible to distinguish under a dissecting microscope. Using species-specific primers/probe and quantitative real time PCR procedures, we screened 13 EPNs species and 6 NF according to well-established methods (Atkins et al., 2005;Zhang et al., 2006;Torr et al., 2007;Campos-Herrera et al., 2011a,b, 2012, 2015Pathak et al., 2012; Table S2). The procedures for the establishment of the pure cultures of the organisms used as positive controls, and the protocols followed for the DNA extraction and the standard curves design were performed following Campos-Herrera et al. (2015). Briefly, EPNs and NF were collected using the sucrose extraction method (Jenkins, 1964), recovering the nematodes in a sieve of 25 µm mesh. The DNA from both the known quantities of the target organisms and the field samples were extracted by Power Soil R DNA Isolation Kit (MoBio laboratories, Inc., protocol for maximum yield, see Campos-Herrera et al., 2015). Details of the concentration and protocols for all the target species were described by Campos-Herrera et al. (2015). We employed a 10-fold dilution for the enumeration of nematodes in the qPCR reactions, whereas, for NF used the total DNA with no dilution. All the organisms quantified by qPCR were expressed as per 100 g of dry soil. In addition, we estimated the NF relative biomass rate by dividing the NF DNA quantity of each species by the total amount of DNA (de Rooij-van der Goes et al., 1995;Campos-Herrera et al., 2012, 2015Duncan et al., 2013).

Statistical Analyses
All statistical analyses were performed with R software, version 3.2.2 (R Core Team, 2015).

Soil Traits
The correlations between the four soil traits recorded (soil humidity, pH, conductivity, and root fraction) and the elevation were individually tested through Pearson's coefficient. Those relationships were plotted in Figure S2 using either non-linear regressions (package "nls2": Grothendieck, 2013) or mixed linear regressions (package "lme4": Bates et al., 2015) with "elevation" as fixed factor and "transects" as random factor.

Abundance and Species Diversity of Nematodes along Elevation Gradients
The total number of nematodes was expressed as number of individuals per 100 g of dry soil and the Simpson diversity index (D) was calculated as a measure of nematode diversity using the package "vegan" (Oksanen et al., 2015). Mixed linear regressions (package "lme4": Bates et al., 2015) were used to analyze the total number of nematodes, the diversity of nematode communities and the infestation rate of nematophagous fungi along elevation gradients, using "elevation" as fixed factor and "transects" as random factor.

Nematode Community Structure along Elevation Gradients
In order to describe the structure of nematode communities (composition and abundance of nematode taxa) along elevation gradients, we performed a partial least square discriminant analysis, PLS-DA (package "mixOmics"; le Cao et al., 2015). The PLS-DA is well suited for dealing with a large number of variables (47 taxa of nematodes) across a limited number of samples (43 sites). Before the analysis, a logarithmic transformation was applied to the data to improve the symmetric distribution of the variables and one outlier was removed from the original dataset based on the initial score plot. In order to plot a multivariate analysis appropriately interpretable, the final model retained was computed with the four soils traits measured and the 19 taxa capturing the most of the variations between nematode communities' structures across elevation zones (i.e., the 19 taxa with a variable importance in the projection, VIP, superior to 1).

Association between Nematode Communities and Plant Communities
First, we constructed plant community structure along the three transects where floristic inventories were conducted by performing a second PLS-DA on the 100 plants with a VIP superior to 1 (results not shown; LV1 = 40%, LV2 = 38% of intergroup variance).
Second, we assessed the correlations between nematode distribution (presence and abundance) and local flora (presence and cover percentage based on Braun-Blanquet categories) across sites using two hierarchical clusterings (package "pvclust": Suzuki and Shimodaira, 2015; the 19 taxa of nematodes and the 100 plants with a VIP superior to 1 were retained). For both clustering, a dendrogram was created using the Ward agglomerative method applied on a distance matrix computed from correlation method. Bootstrap replications were set to 10,000 and stable clusters with a significant approximativelyunbiased (AU) p-value were highlighted (significance level 0.05). Significant correlations between nematode and plant taxa across these two dendrograms were indicated by different segments (Spearman's rank correlation test, α = 0.05; p-values adjusted by the Benjamini and Hochberg method: Benjamini and Hochberg, 1995).

Variation in Nematode Trophic Function along Elevation Gradients
In order to assess the functional role of nematode-based soil food webs along elevation gradients, we calculated several indices, which are based on the abundance of functional guilds of nematodes (Bongers and Bongers, 1998;Ferris et al., 2001; Table S3).
For this purpose, first, all identified nematodes were classified into the main five trophic habits (bacterial-feeders, fungalfeeders, plant-feeders, omnivores, and predators; Yeates et al., 1993), and along the colonizer-persister (cp) scale (Bongers, 1990). The colonizer-persister (cp) scale classifies nematode families into five groups (from 1 to 5) reflecting the lifehistory characteristics similarly to the r/K scale (Bongers and Bongers, 1998). Nematodes belonging to the cp 1 group are fast-growing, bacterivore enrichment-opportunistic nematodes, which increase their population fast after soil enrichment processes; nematodes belonging to cp 2, cp 3, and cp 4 groups present progressively longer life cycles and are more sensitive to environmental perturbation. Nematodes in groups 4 and 5 are in general predators and omnivores, K-strategists, very sensitive to soil perturbation.
Next, we calculated five nematode-based ecological indicators: (1) the sigma-maturity index ( MI; Bongers, 1990), (2) the maturity index (MI; Bongers, 1990), and (3) and the plantparasitic index (PPI; Bongers, 1990). These first three indices represent the proportions of the different cp groups for the whole nematode community, the free-living nematodes, and the plant-parasitic nematodes, respectively. For all three, a higher value indicates that nematodes harboring "persiter" life history traits are predominant within each of those different nematode categories. (4) The enrichment index (EI; Ferris et al., 2001) is based on the biomass of opportunistic nematodes that respond rapidly to the increase of bacterial and fungal populations that arise from organic matter decomposition. High values indicate high soil enrichment and high fertility. Finally, (5) the channel index (CI; Ferris et al., 2001) is the ratio between the biomass of fungivore to bacterivore nematodes, and greater values indicate that fungal decomposition (the fungal "channel") predominates over bacterial decomposition for a given site. For specific calculation of each index see equations provided in Supplementary Materials (Equations 1-3).
In addition to the five nematode community indices described above, we also calculated the metabolic footprints (MF) according to the equation developed by Ferris (2010), and using the Nematode Joint Indicator Analysis tool (Sieriebriennikov et al., 2014; https://sieriebriennikov.shinyapps.io/ninja/; see the Equation 4 provided in Supplementary Materials). The MF balances the mass of carbon used by nematodes for both production (growth and egg production) and respiration (metabolism activities) components. MF can be computed either for specific functional trophic guilds or for the whole nematode community. In this later case, the so-called "composite MF" represents an indicator of the energy flow channeled by nematodes in general within soil food webs. High composite MF suggests that nematode assemblage store high amount of soil carbon (Ferris, 2010).
The effect of elevation on all nematode community indices and MF was tested using mixed linear models (package "lme4": Bates et al., 2015) with "elevation" as fixed factor and "transect" as random factor.

Abundance and Species Diversity of Nematodes along Elevation Gradients
Overall, we extracted 34.752 nematodes from the 48 soil samples collected along the five transects. Nematodes were assigned to 44 genera or three additional families when the identification Frontiers in Ecology and Evolution | www.frontiersin.org  Figure S1. The dashed lines represent the predictions of the mixed linear models. r 2 , coefficient of determination; α, regression slope; P, significance of the regression slope.
of the genus was uncertain. The total number of nematodes increased along elevation gradient ( Figure 1A). Although the specific richness was not affected by elevation (LMM, χ 2 = 0.61, df = 1, P = 0.44), the simpson index measured on nematode communities and taking into account taxa abundance was positively correlated with elevation ( Figure 1B). Finally, the increase in both nematode numbers and biodiversity up along the mountain slope was coupled with an important drop in the relative biomass of nematophagous fungi (Figure 1C).
Beyond the general increase in nematode population along mountain clines, all the trophic groups of nematodes also increased, except predators and parasites, i.e., entomopathogenic nematodes (Figure 2).

Nematode Community Structure along Elevation Gradients
The composition of nematode communities also differed along elevation gradients (Figure 3). The two axes of the PLS-DA retained for the projection explained 57 and 21% of the intergroup variance (Figure 3A). Over the two axes, we were able to discriminate clusters of nematode communities based on elevation zones. Particularly, communities observed between 1500 and 2000 m and those collected above 2000 m were clearly separated from each other ( Figure 3A). As shown in Figure 3B, nematode communities at low elevations, i.e., under 1500 m, were mainly characterized by the presence of Tripyla, Alaimus, Wilsonema and Cervidellus. Tripyla, and Alaimus genera were scarcely present between 1500 and 2000 m and completely absent at higher elevations while Wilsonema and Cervidellus remained at high elevation but in a much lower abundances. The genera Aphelenchoides, Plectus, Prodorylaimus, and Mesodorylaimus mostly dominated the nematode communities between 1500 and 2000 m. Prodorylaimus and Mesodorylaimus were not present in soil originating from other elevational levels. Among soil traits, soil humidity was driving intermediate elevation communities. Over 2000 m, four nematode genera were found in high abundance; three of them (Eudorylaimus, Paratylenchus, Teratocephalus) were also recorded at lower elevations, although in lower abundance, while Pratylenchus was only collected over 2000 m. These high elevation communities were strongly associated with high root fractions in soils.

Association between Nematode Communities and Plant Communities
We next assessed the spatial correlations between nematode and plant communities along elevation gradients. First, the dendrogram based on the distance matrix analysis of plant communities inventoried across three transects indicated a strong clustering of plants according to three elevational levels (Figure 4). Plants characterizing low (i.e., under 1500 m) and intermediate elevations (i.e., between 1500 and 2000 m) were grouped in two single clusters, while plant communities sampled over 2000 m were split in four stable clusters. Second, the dendrogram based on nematodes communities showed a broader, less structured distribution in relation to elevation. Nevertheless, the eight taxa characterizing high and low nematode communities were separated in two single clusters of the dendrogram, while the four nematode's taxa specific to intermediate altitudes were more largely spread across the classification tree.
A total of 44 significant correlations between nematode and plant taxa were observed (Figure 4). These correlations were homogenously distributed across the different trophic groups of nematodes even if 3 herbivore genera, mainly due to Pratylenchus, accounted for almost half of these correlations. Indeed, this genus characterizing alpine communities was related to 15 plant species, most of them typical to high elevations. Prodorylaimus, an omnivore genera particularly abundant in soils collected at intermediate elevations was correlated to six rather uncommon plant species in our dataset belonging to mesotrophic to eutrophic plant communities. In the same elevation zone, Mesodorylaimus was surprisingly correlated with Cuscuta epithymum, a parasite plant species. Finally, two nematode taxa associated to low elevation communities, Tripyla and Cervidellus, were related with only one plant.  Figure S1. When the mixed linear model is significant, the dashed lines represent the predictions of the models. r 2 , coefficient of determination; α, regression slope; P, significance of the regression slope.

Variation in Nematode Trophic Function along Elevation Gradients
We recorded an increase in the relative abundance of colonizer to persistent nematodes up along the transects, as showed by an increase of the sigma maturity index ( MI) with elevation ranging from 2.21 ± 0.06 under 900 m, up to 2.43 ± 0.07 for soils collected above 2000 m ( Figure S3A). This elevation pattern relied mainly on an increase of the PPI at high elevation ( Figure  S3C), while elevation did not affect the MI ( Figure S3B). These results indicate that more persistent nematodes are found at high elevation, mainly due to the high abundance of plant-parasitic nematodes.
As shown in Figures 5A,B the EI informing about the relative abundance of opportunistic nematodes was stable with elevation but the CI increased. This indicates greater relative biomass of fungivore nematodes at high elevation. The CI was multiplied by three to reach 38.11 ± 11.06 in soils located over 2000 m compared to soils under 900 m. Finally, the composite MF clearly increased with elevation, indicating that the amount of carbon entering the soil food webs from nematode food sources is progressively higher at higher elevations ( Figure 5C).

DISCUSSION
The decrease in species diversity along elevation gradient is a common feature for aboveground ecosystems. On the contrary, our results indicate that nematode communities become more abundant and richer at high elevation within the range included in this study (700-2700 m asl). In addition, a steady increase in the composite metabolic footprint (MF) of nematode communities along elevation gradient indicates that nematodes sustain greater part of the soil's energy flow at high elevation. Hence, this study brings novel insights on the role of soil fauna driving soil ecosystem functioning along elevation gradients.

Abundance and Species Diversity of Nematodes along Elevation Gradients
Species richness generally decreases with elevation for most of the organisms studied across a broad range of taxonomic groups, including soil fauna like mites (Chapin and Körner, 1995;Nagy et al., 2003;Hodkinson, 2005;McCain and Grytnes, 2010;Vittoz et al., 2010;Mumladze et al., 2015). Therefore, opposite to general predictions, we observed that both nematode abundance and nematode diversity increases at high elevation (Figures 1A,B). Because our linear models showed no sign of attenuation, it even suggests that the altitudinal threshold after which nematode abundance should decline might be located above 2700 m. Nematodes occur in every ecosystem, often providing available organic carbon sources, and this study confirms their previously reported ability to colonize harsh environments such as Antarctic or high elevation biotopes (Yeates, 2010). Nonetheless, in the Alps, above 3000 m (the alpine and nival stages), vegetation and organic soil layers becomes extremely rare, very inducing a reduction of nematode diversity.
Where vegetation is still relatively abundant (i.e., below 3000 m), different ecological factors can be proposed for understanding this elevation pattern in nematode distribution. Free water in the soil matrix is certainly one of the most important parameter controlling nematode activity and several authors have shown that water availability promotes nematode populations (e.g., Todd et al., 1999;Landesman et al., 2011). Hence, higher nematode abundance at mid-to high-elevation could be linked with the observed increase in soil moisture at high elevation (Figure 3 and Figure S2), which is related higher rainfall frequency and amount at high elevation in the Alps (Körner, 2003). High elevation nematodes are also clearly associated with denser root systems, probably shaping soil micro-habitats that offer shelters against abiotic stresses and increase water retention for free-living nematodes. Additional biotic factors like topdown pathogen pressures over the nematode community are also probably involved, since high elevation soils bare lower amounts of nematophagous fungi (Figure 1C), thereby providing enemy-free zones for nematodes to thrive.
The elevation pattern observed for nematode distribution is in line with Hoschitz and Kaufmann (2004) who recorded relatively high densities of nematodes and diversity within nematode communities collected above 1950 m in the Austrian Alps and is in line with the hypothesis that nematodes might be predominant within high elevation mesofauna because they harbor better adaptations to extreme habitats than most of the other soil-dwelling invertebrates (Procter, 1990). In this context, the ecology of nematode taxa and the functional diversity within nematode communities is expected to vary along mountain clines (see Discussion below).

Nematode and Plant Communities along Elevation Gradients
As shown in Figure 3A, three different nematode communities can be distinguished across four elevational zones. In the alpine zone, i.e., over 2000 m, the composition of nematode communities is characterized by four taxa: Eudorylaimus, Teratocephalus, Pratylenchus, and Paratylenchus ( Figure 3B). Numerous taxa of nematode have a worldwide distribution although some of them are more frequently found in specific habitats. Eudorylaimus and Teratocephalus have been previously found able to colonize arctic or high elevation soils due to their ability to cope with extreme cold temperature (Loof, 1971;Ruess et al., 1999;Hoschitz and Kaufmann, 2004). Pratylenchus, an endoparasitic herbivore, might survive harsh conditions due to its life style within root tissues, which confers appropriate protection against unfavorable environmental conditions (Jones and Fosu-Nyarko, 2014). Within the Paratylenchidae family, many infective juveniles form resistant stages to survive harsh conditions (Bongers, 1990). Consequently, the establishment and winter survival of Paratylenchus at high elevation could be promoted by the highly resistant juvenile larvae, but this needs to be further studied.
Variation in the local flora might also explain variation in nematode community clustering along elevation gradients. Indeed, the hierarchical clustering based on nematode diversity shows that nematode taxa specific to high and low elevation are distributed within two single clusters whereas nematode genera characterizing intermediate communities are spread across the whole dendrogram (Figure 4). While the hierarchical clustering indicates that practically all trophic groups are ubiquitous across the different nematode communities at all elevations, the ecology of nematodes and the functional traits conferring adaptations to elevation niches might be consequently predominant in driving patterns of community spatial variation.
Along the same lines, we could highlight a strong clustering of plant diversity along elevation gradients. At low elevation, quite pristine environments are dominated by beech forests, an habitat where 75 nematode species have been previously inventoried in Denmark (Yeates, 1972). Among nematode taxa characterizing the lowland communities, Tripyla is the most characteristic genus, mainly correlated with the dominant tree species, F. sylvatica. Cervidellus is counterintuitively correlated with the subalpine species Juniperus communis. This might be FIGURE 5 | Effect of elevation on (A) opportunistic nematodes (enrichment index), (B) fungivore to bacterivore nematodes (channel index), and (C) carbon retained by nematodes in the soil food web ( metabolic footprint). Shown are the sampling sites beloning to five mountain transects as shown in Figure S1. When the mixed linear model is significant, the dashed lines represent the predictions of the model. r 2 , coefficient of determination; α, regression slope; P, significance of the regression slope.
driven by the fact that Cervidellus is also present in the subalpine and alpine zones, although in a much lower amount compared to low elevation.
At intermediate elevation, two characteristic taxa of nematodes are correlated to several plant species that encompass a range of diverse habitats. First, Mesodorylaimus is surprisingly correlated with C. epithymum, a parasite plant species developing haustoria on the host stem. Both Mesodorylaimus and C. epithymum were collected in only one sample and we cannot exclude a biased correlation, probably indirectly driven by other, more suitable plants. Second, in the same elevation zone, Prodorylaimus is significantly correlated to six plant species that occur mainly in mesotrophic to eutrophic pastures: Plantago media, Crepis aurea, Potentilla erecta, Carum carvi, Hypericum maculatum, and Deschampsia cespitosa.
Finally, among the nematode taxa characterizing high elevation communities, the herbivore genus Pratylenchus, is correlated with 15 plant species mainly characterizing two plant communities: the high alpine calcareous grasslands with long snow cover, and the subalpine-alpine acidic grasslands or heathlands (transition between the upper subalpine forests and lower alpine grasslands). These numerous correlations confirm that Pratylenchus has developed a wide host range including high diversity of plant habitats (Jones and Fosu-Nyarko, 2014). At high elevation, the relaxation in plant defenses against herbivores could explain this wide host range (Pellissier et al., 2012;Rasmann et al., 2014), but this needs to be confirmed. Thus, considering the number of significant correlations between local flora and the abundance of nematode taxa across elevation zones, this study suggests that nematodes with broader ecological niches including more diversified vegetation might be advantaged at high elevation, i.e., in more fragmented landscapes with higher variability of the vegetation (Rasmann et al., 2014). That said, our methodological approach has the intrinsic limitation of being purely correlative, and future research should address the specificity of plant-nematode interaction across different habitats.

Variation in Nematode Trophic Function along Elevation Gradients
In recent years, numerous studies have advocated the importance of replacing taxonomical biodiversity with functional diversity for uncovering mechanisms of ecosystem functioning (Thébaud and Loreau, 2006;Reiss et al., 2009;Thompson et al., 2012;Montoya et al., 2015). Hence, a shift from studying the composition of nematode communities from a taxonomical perspective to analyzing the assemblage of these communities based on trophic functional guilds is appropriate for better understanding changes in soil ecosystem functioning along elevation gradients. Here, we analyzed several major classes of community indices in order to estimate the contribution of nematodes to soil food web structure along elevations.
First, our results show an increase of the sigma-maturity index ( MI) along mountain slopes ( Figure S3A), and, consequently, indicate that nematodes sensitive to environmental perturbations and harboring longer life cycles are more abundant at high elevation. This increase in the MI relies on both free-living and plant-parasitic nematodes. However, our study suggests that higher MI values on mountaintops are mainly driven by an increase of herbivorous nematodes with slow growing rates and longer life cycle at high elevation ( Figure S3C). This increase in "persisters" within plant-feeding nematodes at high elevation could be due to low annual soil temperature and slow turnover and nutrient cycling, which might be more suitable for nematodes with longer life cycles and low reproduction rates. Furthermore, tolerance of "persister" nematodes toward stress conditions, like those occurring at high elevation, should be more suitable for plant-feeding than for free-living nematodes (Bongers, 1990).
Second, the increase of the CI with elevation ( Figure 5B), i.e., the ratio of fungivore to bacterivore nematodes, reveals that the fungal decomposition pathways support greater nematode biomass than bacterial decomposition at high elevation. This pattern might be explained by the variation of the productivity of soils as argued by Wardle et al. (2004). At high elevation, plant traits such as slow growing and long leaf life span result in slow mineralization rates, and thus in less fertile soils (see also our results of decrease in conductivity at higher elevation, Figure  S2B). These high elevation conditions enhance fungal-based energy flows within ecosystems, consequently accompanied by slower decomposition rates compared to lowland soils with more bacterial-based pathways (Wardle and Yeates, 1993;Zhao and Neher, 2014).
The fifth community index retained in our study, the enrichment index (EI), is not affected by elevation, even if we recorded strong variability between soil locations ( Figure 5A). Hence, the ratio of nematodes indicating enrichment and basal characteristics of the food web is likely to rely on local soil conditions, independently of elevation.
In addition to the above-mentioned indices, the metabolic footprints (MF) of nematode communities can inform on how carbon assimilation in soil food web from autotrophic organisms varies with elevation (Ferris, 2010). Overall, the increase of the composite MF (i.e., the MF for the whole nematode community) along mountain slopes is similarly quite surprising. Indeed, while abundance and species diversity of most of soil invertebrates decrease with elevation (Hodkinson, 2005;Rasmann et al., unpublished), we here show that the energy flow canalized through nematodes increases with elevation. These results are in line with a previous study performed in grasslands and demonstrating the stronger role of soil mesofauna in incorporating carbon within soil foodwebs compared to macrofauna (Ostle et al., 2007). However, changes in habitats, like those occurring along elevation clines, trigger variations in soil biota, and might impact the resulting nutrient fluxes differently. For instance, the abundance and the biodiversity of nematodes and mites along grassland successional stages evolve in opposite directions, and this triggers variation in interactions between vegetation and soil biota and therefore variation in nutrient fluxes (Swift et al., 1998). In this context, our results along elevation pattern deserve further studies for better understanding to what extent nematodes replace other soil invertebrates (earthworms, collembolans, enchytraeidae, mites) for carbon cycling in Alpine soils.

CONCLUSIONS
Our results show that harsh elevation environments drive modifications in the composition of nematode communities based on ecological traits conferring local adaptations. In addition, the correlations between local flora and nematode distribution suggest that nematodes colonizing various high elevation habitats could cope with more fragmented and spatially variable vegetation on mountain tops. Further studies with higher taxonomic resolution are required to validate this hypothesis. Indeed, while genera of nematodes found at high elevation seem to harbor wider host-range, we cannot exclude that different species within these genera are specialized on different habitats and/or plant species. Finally, this study highlights the potential of nematodes' derived ecological indices in understanding ecosystem processes along ecological gradients. We could observe that the role played by nematodes in nutrient cycling increases with elevation, as they partially take over carbon assimilation in soil food web. Future accurate sampling strategies of nematodes across more specific habitats are nevertheless required to dissect how ecosystems types and ecological factors affect nematode-driven soil ecosystem processes along elevation gradient.

AUTHOR CONTRIBUTIONS
SR, RC, SS planned the experiment and collected the data. PV provided plant cover data. AK analyzed the data and wrote the manuscript. All authors reviewed and commented previous versions of the manuscript.

ACKNOWLEDGMENTS
Julia Bilat and Lea Megali provided assistance during field collection. The authors thank the editor for inviting this contribution, and the two reviewers for their insightful comments. This work was financed by Swiss National Science Foundation grants 31003A_159869 and PZ00P3_131956 to SR.