Environmental and biotic drivers of soil microbial β-diversity across spatial and phylogenetic scales

2144 –––––––––––––––––––––––––––––––––––––––– © 2019 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Andres Baselga Editor-in-Chief: Miguel Araújo Accepted 23 August 2019 42: 2144–2156, 2019 doi: 10.1111/ecog.04492 doi: 10.1111/ecog.04492 42 2144–2156


Introduction
It is now well-accepted that soil microbial communities play a pivotal role in the functioning of terrestrial ecosystems (Eisenhauer et al. 2013, Bradford et al. 2017, in particular through their interactions with plant communities and their influence on nutrient cycling (Bardgett et al. 2008, Van Der Heijden et al. 2008. Because soil microbial communities are diverse and composed of organisms with different evolutionary histories, ecological processes determining their structure are thus expected to be numerous and interacting, and to differ among microbial clades. Some studies have shown that, at the landscape scale, microbial communities β-diversity is mainly determined by pedoclimatic conditions (Hazard et al. 2013, Lazzaro et al. 2015; while other studies have demonstrated mutual influence between plant and microbial communities (Zinger et al. 2011, Yuan et al. 2016. Plants engineer the habitat of soil microbes (Rillig et al. 2002, Bardgett and Wardle 2003, Van Der Heijden et al. 2008. In turn, soil microbes influence plant fitness, for instance, by mobilizing and stabilizing nutrients or increasing plant host tolerance to a variety of stresses (Hacquard et al. 2015, Vandenkoornhuyse et al. 2016.
Among microbial communities, fungal communities are expected to have the strongest covariation with plant communities because of the tight mutualistic relationship between fungi and plants; in contrast, bacterial and archaeal communities are thought to be more influenced by their abiotic environment (Roy et al. 2013, Lange et al. 2014, van der Heijden et al. 2016. Nonetheless, the linkage between plant and microbial β-diversities appears modest in many studies (Barberán et al. 2015, Bahram et al. 2016, Zinger et al. 2017). These inconsistencies emphasize that despite considerable advances, there is still no general synthesis on how microbial communities change across space and whether spatial turnover is predictable from basic principles (Keddy 1992, Nemergut et al. 2013. We argue that this gap in knowledge can be addressed through the study of microbial communities across multiple spatial scales and across multiple phylogenetic scales.
Spatial scale is known to strongly influence the detection of assembly rules (Swenson et al. 2006, Chalmandrier et al. 2017). Yet spatial scale is often not appropriately considered or is ignored in analyses of soil microbial β-diversity. On the one hand, typical 'large-scale studies' focus on a large geographical area, include a broad range of environmental conditions and are based on communities defined on a large spatial grain (e.g. sampling unit area). These studies tend to detect the effects of broad environmental gradients rather than of local biotic interactions and stochastic population dynamics (Chalmandrier et al. 2013, Araújo andRozenfeld 2014). On the other hand, typical 'small spatial scale studies' focus on small geographical areas that reflect local pools of species able to survive within the same environmental conditions (de Bello et al. 2012). In such studies, communities are often defined at a small grain, where local processes, such as biotic interactions, are more prominent. A small spatial scale study therefore tends to negate the influence of broad environmental gradients on community structure and better reveal the signal of local assembly processes such as biotic interactions or ecological drift (Chalmandrier et al. 2013, Chase 2014. The integration of different spatial scales into soil microbial studies is thus essential to comprehensively analyze the multiple drivers determining the β-diversity of these communities. There is increasing evidence that evolutionary histories of taxa constrain community assemblages (Groussin et al. 2017, Graham et al. 2018. For instance, species from the same lineage may share ecological characteristics that can promote their co-occurrence in a given habitat, resulting in phylogenetic turnover across environmental gradients; conversely, these shared characteristics may hinder coexistence through increased competition, resulting in a low phylogenetic turnover across environmental gradients (Chalmandrier et al. 2015). To study evolutionary history on microbial community assembly, a promising approach is to define molecular operational taxonomic units (MOTUs) at different phylogenetic grains (or resolution, sensu Hanson et al. 2012), and then to apply classical community diversity pattern analyses. Since some specific microbial and plant traits show phylogenetic signal (Martiny et al. 2015, Valverde-Barrantes et al. 2017, one could expect microbial β-diversity patterns to be shaped by different processes depending on the phylogenetic grain at which MOTUs are defined (Graham et al. 2018). For instance, stochastic population processes are expected to generate ecological drift (Vellend 2010) but should be detectable only at fine phylogenetic grain (i.e. when lineages are defined close to the tips and should thus delineate species boundaries). Conversely, if environmental filtering favors certain lineages that exhibit adequate adaptations, we expect its imprint on community β-diversity to be more detectable at a coarser phylogenetic grain.
In this study, we propose to integrate soil organism and plant evolutionary history to investigate how abiotic processes and biotic interactions with plants shape soil microbial β-diversity at large and small spatial scale. To this end, we collected 378 soil cores in eighteen mountain grassland community plots along an elevation gradient (1856-2725 m a.s.l.). We characterized environmental dissimilarities among plots and measured the spatial distances among samples within each plot. We then used DNA metabarcoding to measure fungal, bacterial, archaeal and plant β-diversity across multiple phylogenetic grains and at two spatial scales: between community plots (hereafter referred to aslandscape scale); and between samples within community plots (hereafter, the plot scale). We then tested the following hypotheses.
(H1) At the landscape scale, we expect soil microbial β-diversities to strongly covary with environmental dissimilarities and weakly with plant β-diversity.
(H2) At the plot scale, we expect soil microbial β-diversity to strongly covary with plant β-diversity and spatial distances, reflecting the influence of interactions between microbes and vegetation and of ecological drift, respectively.
(H3) The relative importance of community assembly rules should differ between the three microbial groups: we expect plant β-diversity to covary more strongly with fungal β-diversity than with bacterial and archaeal β-diversity, irrespective of spatial scale.
(H4) As adaptative traits for soil microbes and plants may display different levels of phylogenetic signal, we expect the influence of space, environmental and biotic drivers on microbial β-diversity to vary according to the phylogenetic grain considered.

Study site and sampling
The study was conducted in the central French Alps (45.12°N, 6.40°E) in summer 2012. Nine sites within exclosures were studied along a continuous 869 m elevation gradient (1856-2725 m a.s.l.) in a pasture grazed by cattle. In this region, subalpine grasslands dominated at low elevation and alpine meadows at high elevation. All sites were on the same south-facing slope, which was composed of moderately acidic soils (pH between 4.9 and 5.7) to minimize environmental variations unrelated to elevation. In each site, we set up two 100 m 2 square plots separated from each other by a few meters. Three botanical surveys (in June, July and August) reported a total of 211 plant species across all plots with an average of 57.67 species per plot. In each plot, 21 samples from the top soil layer (0-10 cm) were collected along the two diagonals, for a total of 378 soil samples (Supplementary material Appendix 1 Fig. A1). Each soil sample contained ~50 g of soil, including both bulk soil and plant roots. The position of the samples along these transects was optimized to have the most uniform distribution of between-sample spatial distances (Supplementary material Appendix 1 Fig. A1, A2).

Environmental information
We measured an array of environmental variables to estimate the environmental dissimilarity among plots: 1) mean annual soil temperature was estimated from field meteorological stations placed in each site; 2) growing season length (GSL) and annual number of frost days were based on daily maps of snow cover and air temperature values following the methodology defined in Carlson et al. (2015); 3) top soil (0-10 cm) characteristics were determined for each plot from the average values obtained from three soil samples collected in August 2012 (C/N ratio, organic matter content and five variables describing soil nitrogen content and fluxes); 4) topographic wetness index and slope, inferred from airborne LiDAR data acquired during the year of the sampling. The environmental dissimilarity between each plot was estimated with Euclidean distances from the first five axes of a principal component analysis of all environmental variables, which were first standardized. These axes explained 90% of the total variance of the environmental dataset. Due to our sampling design, spatial distances among plots were highly correlated with elevation differences (Mantel test, Pearson's r = 0.99, p < 0.001) and environmental dissimilarities (Mantel test, Pearson's r = 0.74, p < 0.001), and were thus ignored in the analysis of microbial β-diversity at the landscape scale. Details about the sampling of environmental variables are available in Supplementary material Appendix 1.

DNA extraction, amplification and sequencing
Extracellular DNA from soil samples was extracted within 8 h after collection following a protocol adapted from Taberlet et al. (2012). Details about the extraction procedure are available in Supplementary material Appendix 1. Each sample was subdivided into two subsamples of 10 g each, which were used as DNA extraction replicates to allow for controlling amplification biases (Ficetola et al. 2015). Four primer pairs were used to amplify specifically the v8-9 region of the 16S rRNA gene in archaea (Taberlet et al. 2018), the v5-6 region of the 16S rRNA gene in bacteria (Fliegerova et al. 2014), the ITS1 in fungi (Fung02 in Pansu et al. 2015, Taberlet et al. 2018) and the P6 loop of the chloroplast trnL intron in vascular plants (Taberlet et al. 2007). For each marker, we conducted two PCRs per extraction replicate (i.e. a total of four replicates per sample).

Data processing
We used the data curation procedure detailed in Zinger et al. (2019) and Ohlmann et al. (2018). This procedure sequentially performs the following steps: assembly of paired-end reads; assignment of reads to samples; removal of low-quality reads, singletons and chimeras; and building of MOTUs (more details are available in Supplementary material Appendix 2). We chose a dissimilarity threshold of three mismatches to define MOTUs. Taxonomic assignments of MOTUs were obtained by comparing sequences against full-length references from the EMBL, SILVA and UNITE databases (Pruesse et al. 2007, Koljalg et al. 2013) and a database specific to arctico-alpine plants (Willerslev et al. 2014). Finally, we removed MOTUs with less than or equal to two reads from each sample. Basic characteristics of the sample-by-MOTUs matrices that were created for bacteria, archaea, fungi and plants (number of MOTUs, sequence length, number of reads) are available in Table 1 and Supplementary material Appendix 2 Table A2.

Analysis
Step

-Generation of MOTU phylogenies
Our approach strongly relies on the evolutionary history of all clades surveyed here. However, retrieving the evolutionary history of organisms from DNA metabarcoding data is not a trivial exercise.
First, DNA metabarcoding markers are short (i.e. between 50 and 400 bp), prone to mutational saturation, and hence yield low phylogenetic signal (Moreira andPhilippe 2000, Deagle et al. 2014). This prevents the direct use of MOTUs defined from DNA sequence dissimilarity thresholds to assess the effect of phylogenetic grain on microbial β-diversity. In addition, these markers are often polymorphic in length and cover a a large taxonomic range. Current reconstruction methods are strongly affected by these characteristics and do not produce reliable phylogenies from DNA metabarcoding markers, when these markers are used in isolation (reviewed by Coissac et al. 2012, Zinger andPhilippe 2016). Here, we applied a more rigorous approach that consists in placing the observed MOTUs on reliable reference phylogenies (i.e. based on longer DNA sequences and/or multiple loci) from their taxonomic assignment.
Second, the taxonomic resolution of such short DNA metabarcoding markers can be relatively low and/or heterogeneous across lineages, which prevents a consistent identification at the species level for every detected organism of the targeted taxon. In addition, while under constant improvement, DNA reference databases for microbes are still incomplete. Consequently, taxonomic information about MOTUs does not always allow precise location on a reference phylogeny, which can increase uncertainty in diversity estimates. To overcome these issues, we propose here a new method: based on MOTU taxonomic assignations and reference phylogenies ( Fig. 1 -blue frame) that seeks to characterize the uncertainty in phylogenetic placement. Briefly, for each marker, we produced a distribution of 50 MOTUs phylogenies to capture the uncertainty in the placement of each MOTU on the reference phylogenies.
Reference phylogenies -we used four reference phylogenies describing the evolutionary relationships of each lineage (archaea, bacteria, fungi and plants). The bacterial and archaeal phylogenies were taken from Lang et al. (2013), and the fungal phylogeny from James et al. (2006). These three phylogenies were transformed into ultrametric phylogenies using the R-function 'chronos' (R-package ape Paradis et al. 2004). The plant phylogeny was an ultrametric genus-level phylogeny of alpine plants from Roquet et al. (2013).
Taxonomic information -the taxonomic classification of each tip of the reference phylogenies and each MOTU was automatically retrieved from the NCBI taxonomy database (Sayers et al. 2009) using the R-function 'classification' (R-package taxize Chamberlain et al. 2014).
Placing MOTUs on reference phylogenies -we conceived an algorithm that 1) compares the taxonomic information of the reference phylogeny tips to that of MOTUS, and 2) optimally grafts each MOTU sequentially on the reference phylogenetic tree. The algorithm works as follows: given a clade to which the MOTU is assigned, the MOTU is stitched randomly to a tip or a branch in the relevant clade in the reference phylogeny. In some cases, the taxonomic information of the MOTU was too precise compared to the resolution of the reference phylogeny. In this case, the taxonomic information was downgraded to the most precise taxonomic level represented in the tips of the reference phylogeny. The updated phylogeny was then used for subsequent MOTU grafting until all MOTUs were placed on the reference phylogeny (Fig. 2).
Most MOTUs have an identifiable classification in the highest taxonomic ranks. Consequently, they are always associated with the same lineages when those are defined at a coarse phylogenetic grain. However, numerous MOTUs are not assigned accurately to low taxonomic ranks. Consequently, lineages defined at a fine phylogenetic grain, have more variable MOTU composition. To control for this degree of randomness, we generated 50 phylogenies for each of the four clades surveyed here (archaea, bacteria, fungi and plants).
Step 2 -Estimating β-diversity across phylogenetic grain and spatial scales Sample and plot scale MOTU relative abundance -we used DNA abundances to calculate our diversity metrics as a way to down-weight low-abundance sequences that may be false Table 1. Basic characteristics of archaeal, bacterial, fungal and plant communities. We report the total number of detected MOTUs for each clades, the median of their length and the number of reads per samples (and in brackets, 2.5% and 97.5% quantile) as well as the median of sample diversity, between-samples β-diversity within plots, plot diversity and between-plots β-diversity (and in brackets, 2.5% and 97.5% quantile). All diversity indices are based on the inverse of Simpson and are computed before exclusion of rare MOTUs (≤ 2 reads) and rarefaction (see Methods). positives (Haegeman et al. 2013, Brown et al. 2015. In other words we used read counts not as true abundances but as a way of weighting sequences that are more likely to be genuine more heavily than those sequences that are likely to be methodological artifacts (Calderón-Sanou et al. 2019). For example, we detected more plant MOTUs larger than expected from our botanical survey (543 plant MOTUs against 221 plant species), indicating that some rare DNA sequences are probably artefacts. However, because we are aware that the reliability of estimating MOTU relative abundance through read counts is a topic of debate (Deiner et al. 2017, Fonseca 2018, we also performed our analyses with presence/absence matrices and presented those in Supplementary material Appendix 4. The total number of sequences per sample was uneven across samples, a feature that may be due to variation in initial concentration of DNA in the sample or to extraction, amplification and sequencing biases. To solve this issue, we used a rarefaction procedure. We randomly sampled a fixed number of sequences per sample. This fixed number was the minimum number of sequences observed across all samples of all plots and was specific to each clade (Weiss et al. 2017). This rarefaction procedure did not induce any significant difference on β-diversity estimates (Supplementary material Appendix 2 Table A2).
We then estimated the relative abundance of each MOTU in each sample using the following relationship between the number of sequences of MOTU i in the sample j, N ij and its estimated relative abundance p ij . The log-transformation of MOTUs read abundances prior to the standardization was justified by a preliminary study on positive control samples which showed that such relative abundance estimates of plant MOTUs were better correlated to the initial DNA concentrations of their corresponding taxa before amplification and sequencing (Pearson correlation: r = 0.68, p < 0.05, Supplementary material Appendix 2). At the landscape scale, MOTU relative abundance was estimated by averaging these relative abundance estimates between samples from the same plot. Defining phylogenetic grain -for all clades, the MOTU relative abundance matrices at the landscape and at the plot scale were aggregated at 139 different phylogenetic grains (a coarse phylogenetic grain reflects ancient lineages, while a fine grain reflects recent lineages). These age values were defined to reflect the age distribution of each phylogeny nodes (see Supplementary material Appendix 2 for details). Each phylogenetic grain defines lineages of MOTUs. MOTU relative abundances were then summed within each lineage at both the landscape scale and the plot scale. We thus obtained for each clade, each spatial scale (landscape and plot) and each phylogeny, 139 matrices of MOTU lineage relative abundances (i.e. one per phylogenetic grain).
Pairwise MOTU β-diversity estimation -we estimated pairwise β-diversities between plots (landscape scale) and between samples (plot scale) across phylogenetic grains. We used the inverse of Simpson that varies between 1 (identical samples/plots) and 2 (completely dissimilar samples/ plots). This metric was used because 1) it is independent of α-diversity (Jost 2007), which is advantageous as the diversities of microbes and plants are expected to change across the elevation gradient and 2) it more heavily weights MOTUs (III -green) From this new 'MOTU lineage by sample' relative abundance table, the β-diversity distance matrix between samples is calculated. This distance matrix is then used as the response variable of the linear mixed models that aimed to measure the respective influence of plant β-diversity (a β-diversity distance matrix generated through a similar process with plants MOTUs) and spatial distance. The workflow for the landscape scale analysis is similar but uses a MOTU by plot relative abundance table obtained by summing MOTUs relative abundance across samples of each plot, and use an environmental dissimilarity matrix instead of a spatial distance matrix. that are relatively more abundant which makes the metric more robust to rare and sometimes artifactual MOTUs (Haegeman et al. 2013). Since the resulting distribution of pairwise β-diversity values was right-skewed, we transformed them with the function f ( ) ( ) β β−1 = to approach a normal distribution. Basic diversity estimates are available in Table 1.

Step 3 -Variance partitioning
We teased apart 1) at the landscape scale, the effects of plant community structure and environment and 2) at the plot scale, the effects of plant community structure and spatial distance using a variance partitioning approach. It was based on sets of linear mixed models that included a random effect factor to account for the non-independence of pairwise distances (Clarke et al. 2002, Lexer et al. 2014.
Landscape scale -we modeled the β-diversity of each microbial group as a function of the β-diversity of plants and the environmental dissimilarity between plots. All variables were standardized. We compared the linear mixed effect models with the following fixed effects: where β ij,y (MIC) is the β-diversity of the focal microbial clade at phylogenetic grain y between plots i and plot j; β ij,x (VEG) is the β-diversity of plants at phylogenetic grain x between plots i and j; Dis ij (env) is the environmental dissimilarity between plot i and plot j; and a, b, c the parameters to estimate in each model.
We then calculated the marginal variance explained by each model (Nakagawa and Schielzeth 2013) and expressed the fractions of variance explained by plant community structure alone mR xy (VEG), environmental dissimilarity alone mR xy (ENV) and the fraction where the effect of plant community structure and the effect of environmental dissimilarity cannot be distinguished mR xy (VEG∩ENV) given microbial lineage y and plant lineage age x as follows (Legendre 2008): mR 2 xy (ENV) = mR 2 xy (full model) − mR 2 xy (Biotic model) mR 2 xy (VEG) = mR 2 xy (full model) − mR 2 xy (Environmental model) mR 2 xy (VEG∩ENV) = mR 2 xy (Environmental model) + mR 2 xy (Biotic model) − mR 2 xy (full model) The analysis was repeated for 50 different couplets of unique microbial and plant phylogenies. We then visualized the median of the corresponding fractions of variance across phylogenies using heat maps (Fig. 3) and further reported the median, 2.5% and 97.5% quantiles of those fractions when exposing the results. We then compare mR 2 xy (VEG), mR 2 xy (ENV) and mR 2 xy (VEG∩ENV) to answer H1, compare their values across microbial taxa to answer H3 and across phylogenetic grains x and y to answer H4. The covariation between plant β-diversity and environmental dissimilarity was also investigated (Supplementary material Appendix 3).
Plot scale -we used the same approach as above at the plot scale, except that i and j refers to samples rather than plots. We only considered the β-diversity pairs within plots and added a random effect for each plot. The fixed effects were the plant community β-diversity (VEG) and the spatial distances (DIS). We then compare mR 2 xy (VEG), mR 2 xy (DIS) and mR 2 xy (VEG∩DIS) to answer H2, compare their values across microbial taxa to answer H3 and across phylogenetic grains x and y to answer H4. We also studied the covariation between plant β-diversity and spatial distance within plots (Supplementary material Appendix 3).
For fungi, the results suggested a covariation between fungal and plant lineages at the plot scale that was not observed for bacteria and archaea. To further explain this pattern, we conducted a supplementary analysis of the partial correlation between plant and fungal lineages relative abundances within each plot (R-package: netassoc, Morueta-Holme et al. 2016). For this, plant and fungal lineages were defined at the age values for which mR 2 xy (VEG) was maximal. We then reported the plant-fungi associations that exhibited the largest with partial correlation coefficients (in absolute values). Detailed methods and results of this analysis are available in Supplementary material Appendix 5.

Data deposition
The data and the scripts associated to this analysis are available from the Dryad Digital Repository: < http://dx.doi. org/10.5061/dryad.m905qftwt > (Chalmandrier et al. 2019).

Plot scale drivers of soil microbial β-diversity
Plot scale models of archaeal, bacterial and fungal β-diversities explained much less variance than at the landscape scale (maximum mR 2 value across taxa, phylogenetic grains and predictors was ≤ 0.15, Fig. 4). We found that spatial distances moderately explained fungal β-diversity

Influence of phylogenetic grain
At the landscape scale, the fractions of variance explained by plant β-diversity alone, environmental dissimilarity alone or their joint effects depended on the phylogenetic grain. Fungi (and to a lesser extent, bacteria and archaea) exhibited a stronger covariation with plant β-diversity when microbial or plant MOTUs were aggregated in lineages (H4). The maximum covariation between fungi and plant β-diversity was observed when fungi MOTUs were considered at a fine phylogenetic grain and plant MOTUs at a coarse phylogenetic grain corresponding roughly to plant families (Fig. 2,  3A). It is noteworthy that plant lineages do not exhibit a stronger covariation with environmental dissimilarities when defined at this phylogenetic grain compared to when they are defined at a fine phylogenetic grain (Supplementary material Appendix 3 Fig. A5).
At the plot scale, microbial β-diversity co-varied with spatial distances only when defined at fine phylogenetic grain and were only marginally influenced by the plant phylogenetic grain (Fig. 4). The effect of plant β-diversity was also dependent on both microbial and plant phylogenetic grain but changed across microbial clades: the covariation between bacterial and plant ß-diversities was maximal at a fine phylogenetic grain for both bacteria and plants; while for fungi, it was maximal at intermediate phylogenetic grains for both plants and fungi. These intermediate phylogenetic grains represented a situation where plant lineages roughly correspond to orders or families and fungal lineages to classes (Fig. 2, Supplementary material Appendix 1 Table A1). Plant lineages defined at this phylogenetic grain did not exhibit a high covariation with spatial distances (0.0384 [0.0366-0.0408], Supplementary material Appendix 3 Fig. A5).  Fig. A6). Among the strongest detected associations, we found that the relative abundances of Glomeromycetes lineages were positively associated with Fabacaeae, Euphorbioideae and Poaceae abundances (partial correlation coefficient averaged across plots are respectively of 0.0128, 0.00925, 0.0113), and were also negatively associated with Cyperoideae (−0.009). Agaricomycetes were positively associated to Cyperoideae (0.0124) and negatively to Poaceae (-0.0104). We also detected a positive association (0.00935) between Pezizomycotina and Ericaceae.

Discussion
Despite considerable advances, there is still inconsistencies on how microbial communities are expected to change across environmental gradients. Here, we show that assessing soil microbial communities at different spatial scales while accounting for their evolutionary history can inform our understanding of the processes forging their assembly.

Predominant drivers of microbial assembly
Soil microbial diversity is known to be structured by elevation gradients (King et al. 2010, Pellissier et al. 2014), but it remains unclear whether this results directly from variation in environmental conditions or whether it is mediated by biotic interactions with plant communities. Our results show an important covariation between plant and microbial communities (H1). However, plant importance differed among microbial groups (H3): direct effects were especially high for fungi, moderate for archeae and low for bacteria. For bacteria, the direct effects of environmental conditions (or joint effect of environment and plant) were more important, confirming previous studies that have consistently linked pedoclimatic variations to bacterial community turnover along elevation gradients (Zinger et al. 2011, Singh et al. 2012. At the plot scale, a large part of the turnover of soil microbial communities was left unexplained by our predictors (H2), suggesting that their assembly is more stochastic than deterministic (Stegen et al. 2012, Bahram et al. 2016 or that we did not include critical factors such as local environmental variation or microbe-microbe interactions (Darrouzet-Nardi andBowman 2011, Ohlmann et al. 2018). Nonetheless, we detected that fungal and bacterial β-diversity co-varied with that of plants with a magnitude on par with other similar studies (Pellissier et al. 2014, Barberán et al. 2015. This may indicate that microbe-plant interactions structure, to a certain extent, the local heterogeneity of those communities. Plant communities are known to interact with specific soil microbial communities in their rhizosphere, which can influence the structure and function of both plant and microbial assemblages as well as ecosystem functioning, such as nutrient fluxes (Rillig et al. 2002, Ohlmann et al. 2018. Vandenkoornhuyse et al. (2016) go as far as describing plants and their attached microbial community as 'holobionts' because of the mutualistic links between plants and soil microbes that built on, for instance, local soil modifications or root exudate production (Broeckling et al. 2008, Lange et al. 2014, Legay et al. 2016). Our study shows that these local processes likely scale up to structure microbial communities at the landscape scale and bacterial and fungal communities at the plot scale.
Fungal, and to a lesser extent bacterial, within-plots β-diversity exhibited minor distance decay in line with previous studies (Martiny et al. 2011, Sayer et al. 2013, Lear et al. 2014. Because of ecological drift, some bacterial or fungal strains may become locally dominant due to stochastic demographic processes (Vellend 2010, Nemergut et al. 2013 and would result in a patchy pattern of microbial community consistent with a local distance decay pattern (Bahram et al. 2016). Alternatively such a distance-decay pattern could be also due to the effect of unmeasured spatially structured soil variables that would generate patch dynamics (e.g. through nutrients or organic matter patches; Darrouzet-Nardi and Bowman 2011, Burns et al. 2015). Support for the ecological drift hypothesis over the latter hypothesis comes from the β-diversity pattern analysis across phylogenetic grains: the 'distance-decay' signal was strongest at the tip level (i.e. close to the species level) and broke down quickly at coarser phylogenetic grains, a result expected for ecological drift. If the pattern was generated by local environmental filtering, we would have expected the distance-decay pattern to be maintained when microbial MOTUs are considered at an intermediate phylogenetic grain because niche characteristics are expected to exhibit a degree of phylogenetic signal.
Compared to fungi and bacteria, archaeal β-diversity was not as well explained by plant β-diversity or spatial distances at the plot scale (H2, H3), confirming results from previous studies in grassland microbial communities (Zinger et al. 2011, Prober et al. 2015. However, archaeal communities were also much less diverse than bacterial and fungal communities, and exhibited reduced heterogeneity at the plot scale (Table 1). This suggests that 1) the variation of local biotic conditions may not affect archaeal communities as dramatically; and/or that 2) the phylogenetic resolution of our archaeal marker is too coarse, and lumps together archaeal species or lineages with different ecology.

An evolutionary perspective on soil microbial assembly
The main drivers of microbial β-diversity changed with the microbial and plant phylogenetic grains that were considered (H4). At both spatial scales, plant β-diversity best explained fungal and bacterial β-diversity when plant communities were considered at an intermediate phylogenetic grain. This hints that relevant innovations promoting plant-microbe associations evolved approximately at the plant family level (at the landscape scale) and at the plant order level (at the plot scale). Candidates features could be various aspects of plant ecology that exhibit phylogenetic signal. For instance, root nitrogen content is typically low within the Asterales order and high within Fabales (Valverde- Barrantes et al. 2017). Variation in proportion of those clades within a plant community could drive changes in soil chemistry and in the associated soil microbial communities (Orwin et al. 2010, Bardgett et al. 2014. Another example are plant mycorrhizal preferences: Poaceae species tend to have arbuscular mycorrhizae while Cyperaceae species have been characterized as ectomycorrhizal (Cripps and Eddington 2005, Brundrett 2009, Gao and Yang 2010. Along elevation gradients, vegetation tends to shift from Poaceae to Cyperaceae dominated grasslands (Descombes et al. 2016), which may in turn change the structure of fungal communities.
At the plot scale, fungal communities presented a singular pattern where the maximal covariation was found at intermediate phylogenetic grains that roughly correspond to fungal classes and plant orders. This can indicate that traits associated with plants-fungi relationships developed early in their evolutionary history. Our supplementary analysis showed the associations potentially responsible of this pattern (Supplementary material Appendix 5). Partial correlation coefficients had only small values which suggests that these associations are loose and/or mainly based on generalist fungi-plant relationships (Sieber and Grünig 2013). Some of the strongest observed associations point to known preferences of plant families for certain mycorrhizal associations: Cyperoidae species were positively associated with Agaricomycetes, which can be related to their preference for ectomycorrhizal associations, while the endomycorrhizal families Fabaceae and Poaceae were positively associated with Glomeromycetes; the clade of arbuscular mycorrhizal fungi (Cripps and Eddington 2005, Brundrett 2009, Gao and Yang 2010. Other associations that we detected are however, to our knowledge, not documented yet. They may be due to the limitations of our methodology or indicate gaps in the current knowledge of plant-fungi relationships.
Compared to fungi, the covariation of bacteria and archaea β-diversity with plant β-diversity decreased quickly as phylogenetic grain became coarser. This may indicate 2154 that adaptations underlying the covariation with plants have evolved recently. This conclusion is tentative, however, because the bacterial and archaeal backbone trees characterized more ancient divergences compared to plant and fungal backbone trees. Thus the absence of phylogenetically-dependent covariation between bacteria or archaea and plants may be due to the lack of characterization of bacterial and archaeal divergences contemporary to the diversification of plants.

Conclusion
Our study shows how soil microbial assembly results from multiple interacting ecological and evolutionary processes. Although the strong environmental variation along the elevation gradient was expected to exercise a substantial and direct filtering on microbial communities, vegetation was an essential factor co-varying with microbial community structure. We further observed that the relative importance of co-factors of microbial community structure varied across spatial scales and among microbial clades. Bacterial communities co-varied mainly with environmental variations and less so with vegetation at the landscape scale while fungal communities co-varied with variation in vegetation structure at both large and small spatial scales. This strengthens the idea that plant and fungal trees of life are 'interwoven' (Arnold et al. 2010) and that these eco-evolutionnary contingencies scale up to affect plant and fungal community structure. Beyond its scope, our study stresses the importance of devising more complete diversity pattern analyses of microbial communities that consider multiple spatial scales and phylogenetic grains to better disentangle interacting assembly rules and uncover new eco-evolutionary aspects of community assembly.