Soil protist function varies with elevation in the Swiss Alps

Summary Protists are abundant and play key trophic functions in soil. Documenting how their trophic contributions vary across large environmental gradients is essential to understand and predict how biogeochemical cycles will be impacted by global changes. Here, using amplicon sequencing of environmental DNA in open habitat soil from 161 locations spanning 2600 m of elevation in the Swiss Alps (from 400 to 3000 m), we found that, over the whole study area, soils are dominated by consumers, followed by parasites and phototrophs. In contrast, the proportion of these groups in local communities shows large variations in relation to elevation. While there is, on average, three times more consumers than parasites at low elevation (400–1000 m), this ratio increases to 12 at high elevation (2000–3000 m). This suggests that the decrease in protist host biomass and diversity toward mountains tops impact protist functional composition. Furthermore, the taxonomic composition of protists that infect animals was related to elevation while that of protists that infect plants or of protist consumers was related to soil pH. This study provides a first step to document and understand how soil protist functions vary along the elevational gradient.

Protists are abundant and play key trophic functions in soil. Documenting how their trophic contributions vary across large environmental gradients is essential to understand and predict how biogeochemical cycles will be impacted by global changes. Here, using amplicon sequencing of environmental DNA in open habitat soil from 161 locations spanning 2600 m of elevation in the Swiss Alps (from 400 to 3000 m), we found that, over the whole study area, soils are dominated by consumers, followed by parasites and phototrophs. In contrast, the proportion of these groups in local communities shows large variations in relation to elevation. While there is, on average, three times more consumers than parasites at low elevation (400-1000 m), this ratio increases to 12 at high elevation (2000-3000 m). This suggests that the decrease in protist host biomass and diversity toward mountains tops impact protist functional composition. Furthermore, the taxonomic composition of protists that infect animals was related to elevation while that of protists that infect plants or of protist consumers was related to soil pH. This study provides a first step to document and understand how soil protist functions vary along the elevational gradient.

Background
Soil microorganisms are extremely diverse and play key roles in soil nutrient cycling and plant growth (Bardgett and Van Der Putten, 2014). Among these, protists (nonfungal microbial eukaryotes) are important contributors to nutrient cycling and soil respiration (Adl and Gupta, 2006;Margerison et al., 2020). Phototrophic protists fix atmospheric carbon and thus participate in the carbon cycle . Predatory protists ['consumers', i.e. individuals that eat other live individuals (Adl et al., 2019)] modulate populations of bacteria, archaea, fungi, algae and micro-metazoan and release nutrients from their prey, controlling soil fertility and plant growth (Bonkowski, 2004;Koller et al., 2013;Gao et al., 2019). Finally, parasitic protists control plant and animal host population, thereby having an impact on biodiversity and even human economy. For example, the oomycete Phytophthora infestans, a well-known potato pathogen, caused the Irish famine in the 19th century that forced millions of peoples to migrate from Europe to America (Large, 1940). Despite their well-documented functions in the soil ecosystem (Geisen, 2018), we still have a limited understanding of the spatial and temporal variability of the functional composition of protist soil communities (Geisen et al., 2017). In particular, there is a need to better document how functional group abundances vary across main climatic and edaphic gradients. Gathering this information is essential for making predictions about how soil organisms and associated biogeochemical cycles will be impacted by changing environmental conditions (Fierer, 2017).
Recent studies have shown that soil protist communities are mainly composed of consumers, parasites andto a lesser extentphototrophs and that their relative abundances vary widely across habitats and across the globe (Oliverio et al., 2020;Seppey et al., 2020). For example, Mahé et al. (2017) found that tropical forest soils host a high abundance of parasitic lineages while Oliverio et al. (2020) found that temperate soil habitats are enriched in consumers. Protist functional group abundance varies also across habitats at local scales: for example, between forest and adjacent agricultural field (Schulz et al., 2019) or grassland (Fiore-Donno et al., 2020). Within a given habitat type, functional group abundance also varies in relation to soil properties [e.g. soil amendment in agricultural fields, (Xiong et al., 2018)] and between seasons (Fiore-Donno et al., 2019). However, it is still unclear how protist functional group composition varies across the elevational gradient. Mountains slopes represent one of the steepest environmental gradients found on Earth and have fascinated biogeographers and ecologists for centuries (von Humbolt and Bonpland, 1805;Hoorn et al., 2018). Mountains are biodiversity hotspots and represent ideal study systems to understand the processes shaping biological diversity (Rahbek et al., 2019). While the taxonomic and functional compositions of many groups of organisms have been documented on mountain slopes (Rahbek et al., 2019), very little is known for soil protists. A recent study has shown that soil protists richness was shown to vary broadly across elevation and to be correlated to both climatic and edaphic factors  but no study, to the best of our knowledge, has documented how the composition of soil protist functional groups varies with elevation.
One way to fill this gap is to use the knowledge accumulated on the function of each taxon and to pair it with the power of next-generation sequencing to attribute the function to identified taxa (Adl et al., 2019). Such an approach was successfully applied for protists in several recent studies (Schulz et al., 2019;Oliverio et al., 2020), and can therefore also be used to document how functional group abundance and composition vary with elevation. Theory predicts that parasite abundance should decrease toward mountain tops as their specific hosts' (or vectors') biomass and/or abundance decreases (van Riper et al., 1986;Zamora-Vilchis et al., 2012). Moreover, parasite taxonomic community composition should be correlated to protist hosts' community composition while phototroph or consumer taxonomic community composition should be more correlated to the abiotic environment (Ohlmann et al., 2018).
Here, using the taxonomical information provided by the metabarcode sequencing of soil samples from 161 open habitats spanning 2600 m of elevation (from 400 m to ca. 3000 m), we ask four questions: (i) what is the relative proportion of DNA sequences of three functional groups (parasites, consumers and phototrophs) in open habitats? (ii) How do these proportions vary across the elevational gradient? (iii) Which environmental factor(s), including plant community composition, best explains changes in functional group composition along the elevational gradient? (iv) How does protist community taxonomic composition vary within functional groups along the elevational gradient?

Soil sample collection
The study region (700 km 2 ) is located in the western Swiss Alps. Open habitats soils sampling was performed in summer 2013, as described in detail in Yashiro et al. (2016). Briefly, from each plot, five soil cores were taken from the four corners and the centre of a 4 m 2 plot, pooled in a sterile bag and kept at 4 C until DNA extraction (within 36 h) and soil analyses were done.
Amplicon 18s rRNA gene sequencing DNA extractions, PCR and sequencing procedures were described in Seppey et al. (2020). Briefly, DNA was extracted from the soil samples using the MoBio PowerSoil DNA extraction kit following the manufacturer's instructions. The V4 region of the 18S rRNA gene was then amplified in triplicates using the general eukaryotic primers TAReuk454FWD1 and TAReukREV3 (Stoeck et al., 2010). PCR conditions and quantification of the amplicon were described in detail in Seppey et al. (2020). Sequencing was performed with an Illumina ® MiSeq (TruSeq Nano PCR-free Library Preparation Kit and the paired-end 2 Â 300 bp) at the University of Geneva (Molecular Systematics & Environmental Genomics Laboratory). Sequences are available on European Nucleotide Archive via the project number PRJEB30010 (ERP112373).

Sampling sites environmental characterization
Climatic and edaphic conditions were obtained for each plot. Below, we briefly describe how they were measured or derived from climatic models but more details can be found in published studies that have already used the same data (Yashiro et al., 2016(Yashiro et al., , 2018Buri et al., 2020;Seppey et al., 2020). As we were mainly interested in documenting variation of protist community composition across the elevational gradient, we intentionally included elevation as a covariate, along with edaphic and climatic factors in our models. However, the biological variation observed along the elevational gradients results from climatic and edaphic parameters co-varying with elevation. The climatic and edaphic factors were selected based on two criteria: they must (i) be known a priori to have an impact on soil microbial community composition, and (ii) be relatively uncorrelated with each other (Spearman rho <0.7). We initially selected seven edaphic variables following  and climatic variables (https://www.unil.ch/ecospat/en/home/menuguid/ecospatresources/data.html#chclim25) at 25m resolution. We then performed a principal component analysis on all 27 environmental factors to extract the axes of main variation. The first axis relates to elevation and other correlated climatic factors (temperature and precipitation) while the second and third axes represent edaphic factors, in particular total organic carbon (TOC) concentration (Additional File 1, Supp. Fig. 1). For subsequent analysis, we included elevation as it was the factor of interest in this study, mean annual temperature (bio1), TOC, soil pH, bulk soil water content (a drainage indicator) and C/N ratio (a biologically relevant summary of nutrient availability) (Additional File 1, Supp. Fig. 2). Finally, we also included plant community type (phytosociological alliances) and plant richness as they can also influence belowground communities (Delarze et al., 1998;Fierer, 2017) (Additional File 2). Plant richness was derived for 145 sites from the existing database (Dubuis et al., 2011) where all vascular plant species were exhaustively inventoried in each site.

Bioinformatic analysis
Reproducibility of the pipeline. The whole bioinformatic pipeline described below has been written in R (version 4.0.2) (R Development Core Team, 2015) with the use of the tidyverse (version 1.3) (Wickham et al., 2019) and the phyloseq packages (version 1.32) (McMurdie and Holmes, 2013). The associated R code is available on github: https://github.com/FloMazel/Protist_Altitude_alps.
Sequencing data processing and taxonomic annotation. Raw reads were demultiplexed and barcodes and primers were trimmed using cutadapt version 2.9 (with parameters e = 0, no-indels, m = 100) (Martin, 2011). Reads were then quality filtered using the filterAndTrim dada2 R function (with parameters maxEE = 6, truncQ = 2, truncLen = 0) (Callahan et al., 2016). Amplicon sequences variants (ASVs) were inferred using the dada function and reads were merged using the mer-gePairs function (with parameter minOverlap = 8). Chimaeras were removed using the removeBimeraDenovo dada2 R function. Triplicates were then pooled according to their respective samples (data are provided in Additional File 3). ASVs were chosen instead of cluster of sequences based on an arbitrary threshold (e.g. 97%) because they: (i) represent finer taxonomic units, (ii) allow comparison between studies and (iii) have become the unit of choice in environmental DNA analysis, including protist sequences (Oliverio et al., 2020). We assigned taxonomy for each ASV using the naïve Bayesian RDP classifier, as implemented in dada2 (function assignTaxonomy, parameter minBoot set to 50) with the PR 2 database (version 4.12.0 -file 'pr2_version_4.12.0_18S_dada2' provided on the dada2 website, data are provided in Additional File 4 along with bootstrap values for each assignation).
ASV filtering and functional annotation. Plant, metazoan, fungi and rare protist ASVs (defined as present in only one sample or with fewer than 100 reads overall) were removed. Each remaining ASV was considered as a protist and assigned to one of the three functional groups (Parasites, Consumers or Phototrophs) according to expert knowledge (see Additional File 5 for the assignation rules and Fig. 1 for the result). We acknowledge that these broad functional categories do not distinguish between fine-scale trophic niche differences. For example, the 'Consumers' group encapsulates protists that might feed on very different substrates (fungi, bacteria, other protists) and might either be predators or saprotrophs. Similarly, the 'Parasites' group includes protists that can infect a wide range of hosts (e.g. plants or animals). While these broad functional groups represent a gross simplification of protist trophic niches, they can be useful to depict the overall distribution of functional groups across environmental gradients and represent a first necessary step toward the functional characterization of soil protist communities (Geisen et al., 2017;Geisen, 2018;Schulz et al., 2019;Oliverio et al., 2020).
Statistical analysis. The relative abundance of the different functional groups along the environmental gradients was visualized by plotting the relative read counts of each group in each sample against environmental factors. However, due to the compositional nature of environmental DNA data, it is impossible to statistically relate these relative abundances to environmental values and ratio of read numbers between groups must be used instead (Gloor et al., 2017;Knight et al., 2018). We focused on the log-transformed ratio between consumers and parasites read counts as these two groups are the most abundant ones. To measure the univariate correlation between elevation and this log ratio, we used a Bayesian linear model using the R package MCMCglmm (version 2.29) (Hadfield, 2010) with 100 000 total iterations and 3000 iterations as burn-in. To further measure the correlation of other environmental factors on this log ratio, we built a linear model using elevation, soil pH, TOC, C/N ratio, mean annual temperature, bulk soil water content (%) and plant alliance and richness as covariates and tested the correlation of each factor using marginal sum of squares and ANOVA with F-values (i.e. testing the correlation of one factor given the correlation of all the other ones, using the drop1 R function).
Within each functional group, we measured community compositional turnover using a metric of beta-diversity (Bray-Curtis) on a standardized ASV*sample matrix (each sample counts divided by sample total sum counts). To measure the correlation of the environmental factors on beta-diversity, we used distance-based redundancy analysis as implemented in vegan (version 2.5) (Oksanen et al., 2016). Each factor effect was tested using randomisation (n = 999) and effect sizes were measured using pseudo F statistics with marginal sum of Fig. 1. Relationship between protists taxonomy and functions in alpine soils. Various taxonomic groups (rank 'Class' nested in rank 'Supergroup' in the PR 2 v 4.12 database, left bars) are connected to their predicted functions (right bars) with links. Link width is proportional to the number of reads assigned to each category in the whole dataset (total reads number = 3 520 752). squares (function adonis2 in the R vegan package, argument 'by' set to 'margins'). To be able to compare the outputs of the distance-based redundancy analysis between functional groups, we selected a common set of samples that had at least 400 read counts for each group. To test the robustness of our result to various sequencing depths between samples and functional groups, we repeated the analysis by rarefying the samples to 350 reads.

Results
A total of 3 520 752 reads corresponding to 5222 protists ASVs were recovered from 161 random-stratified open habitats soil samples in the western Swiss Alps. We then assigned each read to a functional group based on taxonomy (Fig. 1). Over the whole study area (700 km 2 ), alpine protists communities were dominated by consumers (57% of reads, 64% of ASVs) followed by parasites (25% of reads, 14% of ASVs) and phototrophs (7% of reads, 5% of ASVs), while 11% of reads could not be functionally assigned due to poor taxonomic resolution (Fig. 1).
While consumers were over the whole study area the most abundant functional group, their contribution to local (i.e. per sample) read counts varied considerably, between 17% and 89% (Fig. 2). This local variability was correlated to elevation: in general, relative read counts of both consumers and phototrophs increased while those of parasites decreased with elevation ( Fig. 3A) and this trend was observed across most consumer and parasites lineages (Supp. Fig. 1). The average ratio of consumer to parasite reads increased from ca. 3 at low elevation (0-1000 m, median ratio = 2) to ca. 12 at high elevation (2000-3000 m, median ratio = 3.6) (Fig. 4, MCMC p-value <0.001). In our data set, elevation was strongly correlated to the climatic gradients (temperature, precipitation) and, to a certain degree, to edaphic (soil pH, bulk soil water content) and biotic (plant communities, plant richness) factors as well (Additional File 1, Supp. Figs 2-3), so that these environmental factors are also related to the reads counts of the different functional groups ( Fig. 3B-D, Additional File 1 Supp. Fig. 4). For example, the soil from different plant communities [i.e. plant phytosociological alliances (Delarze et al., 1998)] host varying proportions of the three protist groups (Fig. 3D). To measure the correlation between each environmental factor and the ratio of consumer to parasite reads, while taking into account their correlation with the other environmental factors, we proceeded as follows. First, we built a model with all the environmental factors as explanatory variables and the ratio of consumer to parasite reads as the response variable (full model with 145 sites because some of the 161 sites do not have plant data). Second, we dropped each environmental factor independently and one at a time, and measured the fit of the correspondingreducedmodel. Third, we quantified whether the fit of the reduced model was significantly lower than the fit of the full model, in which case the dropped factor was considered significant in the full model. We found that, except for plant alliances, none of the environmental factors had a significant marginal correlation (Additional file 1, Supp. Table 2). In other terms, consumers disproportionally outnumber parasites at high elevations, in alkaline and dry environments and when associated to some plant alliances (e.g. the Petasition paradoxi plant alliance).
Next, we sought to measure to what extent of taxonomic variability within the functional groups across samples (i.e. taxonomic beta-diversity) also correlated with elevation and environmental factors. The taxonomic beta-diversity between samples (as quantified with the Bray-Curtis dissimilarity at the ASV level) within each functional group showed contrasted correlation to  Fig. 2. The relative abundance of soil protist functions varies widely among samples. The relative abundance (relative read counts, y-axis) of the different functional groups (colours) is represented across samples (x-axis). Samples (n = 145) are ordered by descending relative counts of consumer reads. elevation ( Fig. 5A-C for ordination plot and Fig. 5D for correlation values, Additional file 1 Supp. Fig. 5) and other environmental factors (R 2 between 0.25 and 0.35, Fig. 5E). For example, soil pH showed higher correlation with consumer (pseudo-F = 2.3, p-value <0.001, Fig. 5D) and phototrophs (pseudo-F = 2.9, p-value <0.001, Fig. 5D) than with parasite taxonomic beta-diversity (pseudo-F = 1.4, p-value = 0.04, Fig. 5D) and this was similar for the effect of plant richness (Fig. 5). The trend was the opposite for elevation (contrast pseudo-F values in Fig. 5D). Bulk soil water content showed higher correlation with consumer than with parasite or phototrophs taxonomic beta-diversity (contrast pseudo-F values in Fig. 5D). Within parasites, we distinguished between Gregarinomorphea (a clade of Apicomplexan that usually parasites animals), Oomycota (that are mostly plant  Fig. 3. Soil protist functional composition relates to the environment. All panels depict the relationship between functional group relative abundance (relative number of reads per sample) and the environment: either elevation (meters ASL, A), bulk soil water content (%, B), soil pH (C) or plant communities (ordered by descending relative counts of consumer reads, panel D). Trend lines for panels A-C represent a loess fit that is shown for illustrative purpose only.

Discussion
Soil protists play unique and multiple functions in soil ecosystems (Geisen, 2018). To understand how these functions vary along major environmental gradients, and especially along elevation, we compared their relative DNA read counts among 161 soil samples representing varying elevations and environmental conditions from open (i.e. non-forest) mountain habitats (analysis using plant data was done with 145 sites). We found that consumers dominate soil communities, which is in agreement with previous reports in other temperate soils (Oliverio et al., 2020;Seppey et al., 2020) but in contrast with tropical soils where parasites seem to predominate (Mahé et al., 2017;Schulz et al., 2019;Singer et al., 2019). The dominance of consumers in mountain open habitat soils suggests that this protist functional group could be key in the cycling and turnover of nutrients in this type of ecosystem (Foissner, 1999;Bonkowski, 2004;Adl and Gupta, 2006;Oliverio et al., 2020). While being functionally important and abundant, protist consumers are also extremely diverse (Adl and Gupta, 2006). Indeed, they are known to feed on a wide variety of substrates and to occupy different trophic niches (Geisen, 2018). For example, distinct predators have different preys, including bacteria, fungi, micro-metazoa and other protists (Glücksman et al., 2010;Seppey et al., 2017;Geisen, 2018). Given the amplicon sequencing approach taken for this study, it was not possible to finely delineate protist functions, but it could be interesting for future studies to document in more detail the trophic niches of soil protists in order to better sketch out the structure of the regional soil trophic network (Geisen, 2016). Knowledge of the nature of the trophic links between consumer and their substrates will facilitate the measurement of trophic redundancy (i.e. the number of organisms filling the same trophic niche) and vulnerability of the ecosystem to extinction cascades, which may disrupt ecosystem functions (Sanders et al., 2018). In addition, documenting local trophic web structure is believed to be key to test theories on how communities assemble and how biotic interactions vary with the environment (Guzman et al., 2019). Besides consumers, we also found that parasites account for nearly 25% of the total amplicon reads, mostly represented by Gregarinomorphea and Oomycota, to a lesser extent by Phytomyxea. These reads mostly originate from organisms that infect plant tissuesin particular plant roots for Oomycota ASVs (Thines, 2018)or soildwelling animals such as insects, other arthropods and earthworms for Gregarinomorphea ASVs [though genera infecting vertebrates are encountered as well in environmental sequences ]. While the broad functional groups used here represent a gross simplification of protist trophic niches, they can be useful to depict the overall distribution of functional groups across environmental gradients and represent a first necessary step toward the functional characterization of soil protist communities (Geisen et al., 2017;Geisen, 2018;Schulz et al., 2019;Oliverio et al., 2020). We also acknowledge that the relative amplicon read counts used in our study is not a direct measure of absolute abundance, biomass nor biological activity of protists in the field but at best describe relative abundances (Taberlet et al., 2018). Future work could include direct cell count measurements to corroborate our findings about the dominance of consumers in open habitat soils (Wilkinson and Mitchell, 2010). In addition, the relative contribution of organisms to ecosystem functioning does not necessarily fully align with their relative abundance or standing biomass but also depends on their biomass turnover rate (Odum and Barrett, 1971), and this is also the case for soil protists (Wilkinson and Mitchell, 2010). Even if we expect standing biomass and biomass turnover to be correlated (Keeling and Phillips, 2007), it would be interesting for future studies to quantify how the biomass turnover rate of soil protists varies with elevation (Wilkinson and Mitchell, 2010).
The relatively high abundance of parasite reads suggests that they could substantially regulate macroorganism population dynamics (Poulin, 2011;Mahé et al., 2017). If parasites selectively infect specific host lineages, they have in theory the potential to foster local host diversity by limiting population growth of particular host species that are locally abundant. Also, parasites can trigger host diversification and evolve in parallel through a series of Red Queen processes, i.e. the Janzen-Connell hypothesis applied to animals (Mahé et al., 2017). The degree of specificity in parasites varies across clades, with Gregarinomorphea usually assumed to be specific (Smith and Cook, 2008;Rueckert et al., 2019) while this is not as clear for Oomycota (Jiang and Tyler, 2012;Thines, 2018) and Phytomyxea (Neuhauser et al., 2014;Bass et al., 2018). While consumers dominated the regional pool of soil protists, we found large local variations: consumers represented between 17% and 89% of the sample read counts. These pronounced local variations in functional composition have also been observed in previous studies. For example, Mahé et al. (2017) found substantial variation in the protist parasites relative counts across forest soil samples but did not explore how environmental factors correlated to these variations (e.g. in soil properties). The local variability we observed in this study related to elevation and a number of correlated environmental factors (soil pH, bulk soil water content, plant communities). Several competing theories predict such variations: one hypothesis posits that parasite counts should decrease with elevation as host total biomass (and diversity) generally decreases above a certain elevation (van Riper et al., 1986;Zamora-Vilchis et al., 2012), for instance, plants (Körner, 1999) or insects (Hodkinson, 2005). Another hypothesis posits that parasites counts should increase with elevation as host organisms tend to be more physiologically stressed in a harsher environment and invest relatively more resources to survive in these environmentally stressful high alpine environments and less to defend themselves against parasites (Oppliger et al., 1998). Our data validate the first hypothesis: on average, the consumer/parasite reads ratio increases from ca. 3 at low elevation (400-1000 m) to ca. 12 at high elevation (2000-3000 m). We note that due to the compositionality nature of environmental DNA data (Gloor et al., 2017;Knight et al., 2018), the increase in consumer/parasite ratio toward mountain tops can be either interpreted as (i) a decrease in parasites (consumer remaining constant), (ii) a simultaneous increase in consumers and decrease in parasites, (iii) an increase in consumers (parasite remaining constant), or (iv) a simultaneous decrease in consumers and parasites with a more pronounced decrease for parasites. While interpretations (i), (ii) and (iv) are in line with our first hypothesis of how parasites abundance should relate to elevation, interpretation (iii) is more problematic. Indeed, this would mean an absolute decrease of consumers toward more productive lowland sites, a highly counterintuitive result. Hence, we suggest that the relative increase of consumers over parasites with elevation is the result of an absolute decrease of parasite with elevation.
The growth and reproduction of organisms are not directly influenced by elevation but depends on other ecologically meaningful factors that correlate with elevation, in particular climate and soil properties (Körner, 1999). In our study, high elevation sites tended to be colder and drier, and have calcareous and more alkaline soils. The correlation we found between protist functional composition and elevation could not be unambiguously attributed to one specific environmental factor, like soil pH, bulk soil water content or temperature, as those had confounding effects on functional composition.
To disentangle the effects of various climatic and edaphic factors, we suggest that further studies should replicate our sampling scheme along multiple independent elevation gradients and/or design a sampling scheme with a balanced set of alkaline and acidic soil at each elevation and/or make use of experimental designs targeted to answer specific questions within the larger context of climate and edaphic factors. Nonetheless, we found that some plant community composition [approximated here with plant phytosociological alliances, (Yashiro et al., 2018)] correlated with protist functional composition, independently of elevation, climatic and edaphic factors measured in this study. For example, sites colonized by plants from the Petasition paradoxi alliance (scree) tended to harbour relatively more consumers while sites colonized by plants from the Calthion alliance (temporally flooded wet meadows) tended to harbour relatively more parasites. These intriguing correlations might be the product of indirect association between specific plants species and specific soil protists [as is the case between plants and other microbes; e.g. with bacteria in the same area; (Yashiro et al., 2018)] but could also reflect the effect of other micro-edaphic factors that correlated with vegetation type and were not included in our study (Leff et al., 2018;Ohlmann et al., 2018). For instance, soils of the Calthion alliance are water-logged and anoxic, and thereby selects for plant species and bacterial communities adapted to such conditions (e.g. methanotrophs and anaerobes) (Yashiro et al., 2018). Future in situ functional studies within the protist parasites group could reveal whether parasitism is what is selected for in Calthion sites or whether other functional aptitudes that selectively promote survival in wet anoxic ecosystems are causing increasing relative abundances of the parasites functional group in our study. In addition, it will be key for future studies to document the temporal dynamic of protist functional groups abundance, as it is well established that protist communities exhibit seasonal variation, for example in water (Simon et al., 2015) and soil (Fournier et al., 2020). Even if soil protist communities seem to vary more spatially than temporally (Fournier et al., 2020), we acknowledge that some of the spatial variation reported in this study could reflect temporal dynamics as the sites were sampled at different dates. Future studies could use seasonal variation to better understand parasite abundance dynamics, for example by comparing parasites relative read counts to host biomass and diversity throughout the season.
Finally, we documented how the taxonomic identity of the protists within functional groups responds to the environmental gradient in alpine soils. In other systems, it has been shown that, while functional group composition is often strongly related to environmental gradients, trends in the taxonomic composition within functional groups appear much more difficult to detect, for example bacteria in the global ocean (Louca et al., 2016). It has been suggested that these contrasting patterns could result from functional redundancy of organisms within but not between functional groups (Louca et al., 2018). In contrast, here we show that the taxonomic composition within each of the three functional groups of protists varied strongly and predictably along the environmental gradient. We found that climatic and edaphic factors explained up to 40% of the within-group protist taxonomic composition. Interestingly, functional groups presented contrasted responses to the environmental gradients. Theory predicts that host-associated (e.g. parasite) community composition should be correlated to their host community composition while non-host-associated (e.g. phototrophs) taxonomic community composition should be more correlated to the abiotic environments (e.g. soil pH) (Ohlmann et al., 2018). Our data do not support this hypothesis: for example, the correlation between plant and protist taxonomic composition was not drastically higher for parasites than for consumers and phototrophs. Interestingly, we found that plant parasite (most of the Oomycota reads) taxonomic composition was correlated to soil pH while animal parasite (Gregarinomorphea) taxonomic composition was more correlated to elevation. The latter result can be interpreted as an indirect effect of animal host turnover along the elevation gradient because (i) Gregarinomorphea are usually considered as relatively host-specific (Mahé et al., 2017;del Campo et al., 2019;Singer et al., 2019) and (ii) soil animal communities change along elevation gradients (Hodkinson, 2005). The fact that plant parasite communities were correlated to soil pH and TOC can be explained by two non-mutually exclusive hypotheses. First, individual parasites might be sensitive to abiotic conditions (e.g. soil, climate) during the phase of their life exposed to the outside environment. This is for example the case when they disperse between hosts (i.e. at the zoospore stage) and are not living within host bodies (Thines, 2018). Second, soil pH indirectly affects the presence and/or abundance of their host which could, in turn, affect the distribution of parasite DNA found in soil. For example, plant species distributions and abundances are known to vary with soil pH  so that the distribution of protists that parasitize plants could also be related to soil pH. However, as Oomycota are relatively less hostspecific (Thines, 2018) and the DNA marker used here is not as specific as it is for animal parasites (i.e. it does not allow distinguishing closely related taxa) the second hypothesis seems less likely than the first one. We suggest future ecological studies could improve our understanding of the host-parasite association by better documenting (i) the interaction between the partners, and (ii) the host and (iii) the parasite identity. For example, contrasting the distribution of multiple lineage of parasites represents a fertile ground to better understand the ecology and evolution of protist parasites and their hosts (Poulin, 2011). Retrieving parasite DNA from the soil provides only a little evidence for a direct association with a specific host. We suggest that future studies could directly sample host tissue (e.g. roots) to provide (i) more direct evidence of the interaction as well as (ii) a better idea of the host identity. (iii) Current amplicon approaches based on one unique primer set provide limited information on rare parasites (e.g. the very low abundance of endomyxea and cryptomycota prevented us to analyse them statistically here). We suggest that focusing on multiple parasitic lineages at the same time using different sets of PCR primer targeting different lineages (e.g. as in Bass et al., 2018) represents a valuable approach to better document the breath of parasite diversity. bioinformatic pipeline code is published on github: https:// github.com/FloMazel/Protist_Altitude_alps.