The relative importance of ecological drivers of arbuscular mycorrhizal fungal distribution varies with taxon phylogenetic resolution

The phylogenetic depth at which arbuscular mycorrhizal (AM) fungi harbor a coherent ecological niche is unknown, which has consequences for operational taxonomic unit (OTU) delineation from sequence data and the study of their biogeography. We tested how changes in AM fungi community composition across habitats (beta diversity) vary with OTU phylogenetic resolution. We inferred exact sequence variants (ESVs) to resolve phylotypes at resolutions ﬁner than provided by traditional sequence clustering and analyzed beta diversity proﬁles up to order-level sequence clusters. (cid:1) At the ESV level, we detected the environmental predictors revealed with traditional OTUs or at higher genetic distances. However, the correlation between environmental predictors and community turnover steeply increased at a genetic distance of c . 0.03 substitutions per site. Furthermore, we observed a turnover of either closely or distantly related taxa (respec-tively at or above 0.03 substitutions per site) along different environmental gradients. (cid:1) This study suggests that different axes of AM fungal ecological niche are conserved at different phylogenetic depths. Delineating AM fungal phylotypes using DNA sequences should screen different phylogenetic resolutions to better elucidate the factors that shape communities and predict the fate of AM symbioses in a changing environment.


Introduction
Arbuscular mycorrhizal (AM) fungi are a widespread group of plant symbionts (Brundrett & Tedersoo, 2018) belonging to the subphylum Glomeromycotina (Spatafora et al., 2016), and play a key role in ecosystem functioning (Powell & Rillig, 2018). To understand how organisms interact with their environments, it is necessary to adequately group individuals into a biological entity that harbors a coherent ecological niche or function. Descriptions of fungal communities are now largely based on DNA sequencing data, where the sequences obtained are usually clustered into operational taxonomic units (OTUs; Lindahl et al., 2013); that is, groups of similar sequences at a given percentage of sequence similarity. However, the appropriate sequence similarity level (i.e. the phylogenetic resolution) at which these OTUs should be defined remains debated, despite renewed efforts to appropriately delineate evolutionary and/or ecologically coherent units in AM fungi (Powell et al., 2011;Lekberg et al., 2014). This is because, in AM fungi, we have still little understanding of the ecological variation (and underlying trait variation) at different levels of phylogenetic divergence (e.g. intra-vs interspecific genetic variation; Hazard & Johnson, 2018).
There is now compelling evidence that the community composition of AM fungi is interactively influenced by different abiotic and biotic drivers (e.g. soil pH and ecosystem type; Davison et al., 2015). According to classic community ecology theory (Vellend, 2016), these drivers filter organisms from a regional pool based on their functional traits (e.g. capacity to associate with a given group of plants, tolerance to acidic soils). However, trait information is often unavailable for AM fungi, for which we often only have access at best to the organisms' identities (e.g. OTUs) and their distribution. In this case, the phylogenetic depth at which these ecological traits are conserved is pivotal for elucidating aspects of the environment that drive organism distribution (Martiny et al., 2015). For instance, consider two closely related phylotypes (a1 and a2) adapted to colonize grasslands and two other closely related phylotypes (b1 and b2) adapted to colonize forests. Phylotypes a1 and a2 (and b1 and b2) are found in the same broad habitat but occur in different microhabitats: phylotypes a1 and b1 are present in acidic soils, whereas phylotypes a2 and b2 are found in alkaline soils. Because the habitat preference (grasslands/forests) is conserved at the level of broad clades, its correlation with organism distribution will be strong at this broad phylogenetic resolution but detectable to a lesser extent at the phylotype level. By contrast, lumping the two fine-scale phylotypes into the broader clade will obscure soil pH effects. In summary, studying community composition at one phylogenetic resolution (typically at the species or one OTU level) most probably gives only partial information on the relative importance of the drivers of organism distribution (Hanson et al., 2012). These drivers could be identified by screening the turnover of AM fungal communities at different OTU phylogenetic resolutions (Groussin et al., 2017).
The beta diversity patterns (i.e. similarity in community composition across spatial and/or environmental gradients) of AM fungi are likely to change with OTU phylogenetic resolution because their functional traits are conserved at different phylogenetic depths. Some traits are similar among closely related species and others are not. For example, spore size and amount produced, as well as investment in extra-or intraradical mycelium, appear conserved at the family level (Hart & Reader, 2002;Maherali & Klironomos, 2007;Powell et al., 2009;Chagnon et al., 2013). Likewise, family-or order-level community composition is well predicted by spatial variation in soil pH, soil phosphorus (P) and nitrogen (N) content, or soil depth (Camenzind et al., 2014;Rodr ıguez-Echeverr ıa et al., 2017;Roy et al., 2017;Sosa-Hern andez et al., 2018a;St€ urmer et al., 2018a;Treseder et al., 2018), supporting a phylogenetic conservatism of niche and traits at a coarse phylogenetic resolution. On the other hand, AM fungi exhibit extensive genetic variation within morphologically (spore-based) defined species, which can affect compatibility with their host plant (Angelard et al., 2014) and traits such as extraradical hyphal density (Munkvold et al., 2004;Koch et al., 2006;Mensah et al., 2015): the variability in these traits and the consequences for host performance were higher within than between AM fungal species. Accordingly, though AM fungal community composition has been shown to vary consistently across land-use management or vegetation types at a regional scale regardless of OTU phylogenetic resolution, the association with a particular plant community type at the local scale was better resolved at a fine OTU phylogenetic resolution (Lekberg et al., 2014;Powell and Sikes 2014). These observations suggest that screening multiple phylogenetic resolutions, ideally both at the inter-and intraspecific levels (Johnson et al., 2012;Sanders & Rodriguez, 2016), would help to better understand the ecology of AM fungi. Ecological variation at the intraspecific level is challenging to infer from DNA sequence data owing to the difficulties in differentiating genuine sequences from PCR and sequencing artifacts. These latter are often corrected by clustering sequences at a given similarity threshold (Lindahl et al., 2013;Hart et al., 2015), which is often a compromise between the rate of errors produced by molecular techniques and biological variation, hence resulting in the loss of intraspecific variability. New bioinformatics methods now allow better discriminating genuine sequences from artifactual ones. These genuine sequences are termed exact sequence variants (ESVs, also known as amplicon sequence variants) and may diverge from one another by differences as subtle as one nucleotide (Eren et al., 2015;Callahan et al., 2016;Edgar, 2016). Inferring ESVs has advantages over classical OTUs. The first advantage is because ESVs represent a fundamental unit whereas traditional OTUs are nota genuine DNA sequence that is comparable across studies (Callahan et al., 2017). The second advantage is because ESVs keep valuable genetic information that is lost in traditional OTUs (Selosse et al., 2016), thereby increasing the phylogenetic resolution of observation. However, we still do not know whether the finest phylogenetic resolution (i.e. ESVs) is useful for inferring ecological processes in AM fungal communities, or if it only increases noise-to-signal ratio.
We here inferred ESVs to resolve phylotypes at resolutions finer than provided by traditional sequence clustering, and we implemented a phylogenetic decomposition of beta diversity to test whether and how beta diversity patterns change across phylogenetic resolutions from the ESV level to order-level sequence clusters. We used two datasets of AM fungal large-subunit (LSU) ribosomal DNA (rDNA) sequences generated with an Illumina MiSeq platform: a chronosequence of 52 yr of agricultural recultivation after open-cast mining (Roy et al., 2017; hereafter 'Chronosequence') and a soil depth diversity survey in an agricultural field (Sosa-Hern andez et al., 2018a; hereafter 'Soil depth'), both in western Germany. Each dataset spanned variations in several aspects of the environment that likely impose selection on traits that are conserved at different phylogenetic depths. Given the variable degree of phylogenetic trait conservatism in AM fungi, we hypothesize that, depending on the environmental gradient considered, we will observe a turnover of either closely or distantly related taxa.

Dataset description
The Chronosequence dataset consists of 10 fields spanning a recultivation gradient of 52 yr since open-cast mining near J€ ulich, Germany. Briefly, the recultivation process starts from a mix of old agricultural soil and the overburden loess layer from where lignite is extracted. It then consists of three phases: (1) 3 yr of soil restoration with alfalfa, (2) 2 yr of conversion to conventional agriculture, and (3) implementation of conventional agriculture by local farmers. In each field, five replicated soil cores were collected. The Soil depth dataset was obtained from an agricultural field near Bonn, Germany. It consists of soil samples at depths of 10-30 cm and 60-75 cm, and distributed in three soil compartments: the rhizosphere (soil directly influenced by roots), drilosphere (soil directly influenced by earthworms), and bulk soil, replicated three times. The same protocol (including primers, PCR conditions, and sequencing platform) was used to characterize the composition of the AM fungal communities in both datasetssee Roy et al. (2017) and Sosa-Hern andez et al. (2018a) for the detailed protocols. Amplicons of the LSU ribosomal (rRNA) gene were sequenced on an Illumina MiSeq platform using the 2 9 300 paired-end chemistry. Raw reads and sample metadata are freely available at https://doi.org/10.5061/dryad. t096s (Chronosequence dataset) and https://doi.org/10.6084/m9.f igshare.8850119 (Soil depth dataset).

ESV inference
We used DADA2 (Callahan et al., 2016) in R (R Core Team 2017) to obtain denoised, chimera-free, nonsingleton ESVs. Primers were removed and the forward and reverse reads were trimmed at 270 bp and 220 bp, respectively, due to the base-call quality degradation in the 3 0 end of the reads. Sequence trimming at these positions permits a minimal overlap of 20 bp between the forward and reverse paired-end reads, which is necessary for their assembly, and an additional 70 bp of sequence length variation, which corresponds to that observed in the reference sequences of Kr€ uger et al. (2012) for our amplified LSU region (length of reference sequences is 349-419 bp). Sequences with ambiguous bases were excluded, and the maximum expected number of errors was set to two and five for forward and reverse reads, respectively. We chose this threshold, instead of one, the default parameter of the USEARCH OTU pipeline (Roy et al., 2017), because DADA2 infers genuine sequences using an error-rate model. Setting a higher threshold of expected errors allowed keeping a substantial sequencing depth equivalent to that of the USEARCH procedure. ESVs were inferred on a sample basis. Singleton ESVs were considered as artifactual and removed. ESVs for which the sequence corresponded to subsequences of two more abundant sequences were considered chimeras and removed. Taxonomic identity of ESVs was inferred using BLAST against reference Glomeromycotina sequences (Kr€ uger et al., 2012). Taxonomic assignment was done for queries having at least 90% of sequence similarity, 90% of sequence coverage, and a minimum e-value of 10 À50 against their closest reference. The sequences not fulfilling these criteria were considered nonGlomeromycotina and were excluded from the analysis.

Additional filtering of putative PCR errors
Despite our conservative sequence data curation and the low proportion of error sequences retrieved with DADA2 (Callahan et al., 2016), our dataset could still contain PCR and sequencing errors. Indeed, DADA2 assumes independent errors in the model, but some errors tend to be repeatable and can co-occur with the true sequences from which they originate (Coissac et al., 2012). The co-occurrence of errors with their genuine sequence, besides inflating diversity estimate, can hence artificially inflate the correlation between phylogenetic beta diversity and environmental gradients. To correct for this potential bias, we used LULU in R (Frøslev et al., 2017). We used the default parameters based on the co-occurrence and relative abundance relationships between closely related ESVs. We considered as a PCR error the sequences that are consistently at a lower abundance than their true sequence relative across all samples. Putative PCR error sequences were removed without summing their read counts to the inferred true sequence relative. LULU curation served two purposes: (1) to examine the robustness of the results to putative repeatable sequence errors; (2) to estimate lower and higher boundaries of phylogenetic resolution at which biological units correlate with the environment in case true diversity is removed.

Phylogenetic tree reconstruction
ESVs were aligned to the reference alignment of partial small subunit (SSU)-internal transcribed spacer (ITS)-LSU rDNA Glomeromycotina sequences (Kr€ uger et al., 2012), which we used as a phylogenetic backbone. The alignments were performed using MAFFT (Katoh et al., 2002) with the -addfragments option to align short sequences to an alignment of long sequences. A phylogenetic tree was built with RAXML (Stamatakis, 2014) from these alignments, containing both ESVs and reference sequences, by conducting a bootstrap analysis (100 bootstraps) and searching for the best-scoring maximum-likelihood tree under a GTRGAMMA model of nucleotide substitution. The nucleotide sequence divergence expressed in mean substitution per site (subs/site) is hereafter referred to as genetic distance. The phylogenetic tree was analyzed and visualized in R using APE (Paradis et al., 2004) and GGTREE (Yu et al., 2017).

Decomposition of beta diversity across phylogenetic resolutions
Beta diversity through time In order to assess beta diversity patterns at different OTU phylogenetic resolutions, we used the beta diversity through time (BDTT) approach (Groussin et al., 2017) and clustered ESVs into OTUs by steps of nucleotide sequence divergence (hereafter referred to as OTU BDTT ; ESVs are OTU BDTT0s ) across the phylogenetic tree. The set of samples in which an OTU BDTT occurs is defined as the union of the samples where each ESV belonging to this OTU BDTT is found. The abundance of the OTU BDTT is the sum of the count of its respective ESVs. For each genetic distance, a new OTU contingency table was created. The function was adapted to cluster tree edges in nonultrametric trees to maintain cluster monophyly and, as such, uses directly the estimated evolutionary distances. This method will tell us at which phylogenetic resolution the OTUs best correlate with the different environmental gradients by testing whether the genetic distance used to define phylotypes is too narrow (i.e. there is noise, and it is therefore necessary to lump finer OTUs into broader OTUs to detect a signal) or has insufficient resolution to detect a correlation with the environmental (i.e. it is necessary to split broader OTUs into finer OTUs to detect a signal).

UPARSE OTU clustering
We used also UPARSE OTU clustering, a common approach in AM fungal ecology Roy et al., 2017;Sosa-Hern andez et al., 2018a), to test the consistency of beta diversity patterns between different methods. Sequence processing followed Roy et al. (2017). Briefly, pairedend reads were merged and quality filtered in USEARCH (Edgar 2010 (Schloss et al., 2009;Kozich et al., 2013;Schloss, 2013). Sequences were clustered into OTUs using UPARSE (hereafter referred to as OTU UPARSE ; Edgar, 2013) at sequence similarity ranging from 99% to 70% to test the robustness of our phylogenetic decomposition of beta diversity to OTU clustering methods. LULU curation was run, as described earlier, to discard putative PCR errors.
Mean nearest taxon distance and mean pairwise distance Both the BDTT and the UPARSE clustering rely on fixed genetic distance thresholds, and hence do not take into account the different rates of trait evolution across clades, which are nonnegligible (e.g. € Opik et al., 2010; Powell et al., 2011;Lekberg et al., 2014). To consider this aspect, we additionally used beta diversity metrics that measure the shared evolutionary history between communities (measured in branch length) to further estimate the depth of phylogenetic turnover across different environmental gradients without an OTU clustering step. To this end, we calculated the mean nearest taxon distance (MNTD) and the mean pairwise distance (MPD) using PICANTE in R (Kembel et al., 2010). MNTD and MPD are averaged measures of the phylogenetic distance between pairs of taxa (ESVs here) drawn from two communities at fine (MNTD) or coarse (MPD) phylogenetic resolution.
Taxonomic resolution To test whether taxonomic annotations reflect phylogenetic OTU clustering and phylogenetic beta diversity measures, we used the taxonomic annotations of ESVs to cluster them (and to sum their abundances) at taxonomic ranks from the species to order level.

Statistical analyses
Estimating the correlation between beta diversity at different OTU phylogenetic resolutions and the environmental predictors To assess how the beta diversity patterns at different OTU phylogenetic resolutions correlated to environmental gradients, we used nonparametric multivariate ANOVA (PERMANOVA; Anderson, 2001) in VEGAN (Oksanen et al., 2016) based on Bray-Curtis dissimilarity. This method partitions the variance in beta diversity among the different environmental predictors. For the Chronosequence dataset, 24 soil variables were measured for the 50 samples considered here. We summarized these soil variables into three orthogonal ecological gradients, corresponding to the first three principal components of a principal components analysis. Principal component 1 (PC1; 44.6% of variance) corresponds to a gradient of soil extractable P and total N contents that increase from young (1-to 3-yr-old fields and 4-to 5-yr-old fields) to old fields (> 10 yr old; Supporting Information Fig. S1). Principal component 2 (PC2; 17.7%) corresponds to a gradient of inorganic N, which peaks in the 5-to 10-yr-old fields. Principal component 3 (PC3; 9.7%) differentiated fields from one another across time with a peak in magnesium in the 10-yr-old field. From now and on, PC1, PC2, and PC3 refer to the environmental predictors used for the Chronosequence dataset. For the Soil depth dataset, we used the soil layer (topsoil vs subsoil) and compartment (rhizosphere, drilosphere, bulk soil) of 18 samples from one agricultural field (Sosa-Hern andez et al., 2018a). The significance of the environmental predictors was assessed with 999 Monte Carlo permutations. Beta diversity patterns at different phylogenetic resolutions and for each beta diversity metric were visualized using a principal coordinates analysis (PCoA).
To ensure that our results were due to a genuine phylogenetic signal (i.e. departing from random phylogenetic relationships), we randomized the tip labels of the phylogeny while keeping constant the community composition as in Groussin et al. (2017) 100 times. Variance partitioning was conducted for each permutation step to generate a null distribution of the explained variance (PERMANOVA R 2 ) of each variable. R 2 values were considered not due to chance when departing from the 5-95% quantiles of the null distribution. A departure from the null distribution further indicates that an increase in explained variance at coarser phylogenetic resolutions is not solely due to a reduction in OTU number. We also took into account the uncertainty in phylogenetic tree topology in our analysis. Variance partitioning was repeated for each of the 100 bootstrapped RAXML trees, and predictors were deemed significant if the P values were < 0.01 for at least 90% of the trees.
Decomposition of MNTD and MPD phylogenetic beta diversity measures into phylogenetic distance classes To estimate the phylogenetic similarity between communities across the environmental gradients, we further decomposed the correlation of MNTD and MPD phylogenetic beta diversity with the environmental gradients into small (shallow phylogenetic turnover) and large (deep phylogenetic turnover) phylogenetic distance classes using Mantel correlogram analysis in ECODIST (Goslee & Urban, 2007). The sign of the correlation indicates whether communities within a phylogenetic distance class are found in similar (positive correlation) or dissimilar (negative correlation) environments compared with communities in the other distance classes. The significance was assessed with 999 Monte Carlo permutations of the phylogenetic distance matrix.
Phylogenetic signal of ESV ecological traits To complement the beta diversity analyses, we further tested whether closely related ESVs tend to share similar ecological traits (i.e. phylogenetic signal) and estimated the phylogenetic depth of this correlation. To this end, we conducted phylogenetic correlograms analyses using the PHYLOSIGNAL R package (Keck et al., 2016) to test a phylogenetic signal in ESV co-occurrence (i.e. whether closely related ESVs tend to occur in the same samples) and in the optimum of each of the environmental predictors (as a measure of ecological traits). Two closely related ESVs may not cooccur yet be found in samples of similar environmental values. The co-occurrence between any pair of ESVs was measured as the Euclidean distance of ESV presence/absence in each sample. We inferred the ecological trait of each ESV by calculating their optimum for each of the environmental predictors (i.e. PC1, PC2, PC3, soil depth, and compartment). These optima corresponded to the abundance weighted mean of each of the environmental predictors for each ESV. Phylogenetic correlograms were calculated on a continuous basis by computing Moran's I index www.newphytologist.com (Moran, 1948) in the case of PCs, soil depth, and compartment class optima, or with the Mantel statistic in the case of co-occurrence. A matrix of phylogenetic weights was computed with default parameters to account for the nonuniform distribution of tips within the phylogeny (e.g. caused by different rates of diversification between lineages or phylogenetic sampling extent). Significance was tested using 100 nonparametric bootstrap resamplings of the phylogenetic tree to estimate the 95% confidence interval (CI) of the correlation. CI curves not overlapping with zero are significant.
Identification of clades that covary with the environmental predictors To further identify clades at different phylogenetic depths that covary in abundance with each of the environmental predictors, the correlation of node abundance (sum of the count of a node's respective ESVs) with each environmental predictor, for all nodes along the phylogenetic tree, was further tested with a linear model and corrected for multiple testing using the false discovery rate (Benjamini & Hochberg, 1995).

Sensitivity analysis
To assess to what extent the approaches and OTU phylogenetic resolution mentioned earlier yielded similar results in terms of beta diversity, we conducted a sensitivity analysis by comparing beta diversity patterns with Mantel tests (Pearson's R correlation coefficient) using VEGAN (Oksanen et al., 2016). It revealed that beta diversity patterns varied according to the phylogenetic resolution at which they were assessed and not the method used (Figs S2-S5). Therefore, in the following, we discuss the BDTT analysis using the LULU curated dataset.

Number of OTUs across phylogenetic depths
In total, we obtained 574 and 233 Glomeromycotina ESVs in the Chronosequence and Soil depth datasets, respectively, of which 240 (42%) and 75 (33%), respectively, were kept after LULU curation, accounting for 84% and 77% of the reads. In both the Chronosequence and Soil depth datasets, we observed a drop in the number of OTUs from 0.01 to 0.03 subs/site (Figs 1-3a

Environmental predictors of AM fungal beta diversity across OTU phylogenetic resolutions
For the Chronosequence dataset, PC1 was the strongest predictor of AM fungal beta diversity and was significant at all OTU phylogenetic resolutions, but the variance explained varied strongly (Fig. 3b): PC1 explained 11% of the variance at the ESV level, but rose strongly to 22% at 0.03 subs/site and rose again to 30% at 0.1 subs/site. In addition, PC2 also explained 8% of the variance at the ESV level, rose to 12% at 0.03 subs/site, where it peaked, and then dropped. PC3 explained little but significant variance in community composition at the ESV level and was not significant at higher genetic distances. These results are in agreement with our PCoA analyses (Fig. S3): at 0.2 subs/site, samples clustered into two groups corresponding to fields recultivated for less (low PC1 value) and for more (high PC1 value) than 10 yr. At 0.03 subs/site, communities from the three first years of recultivation (low PC2 value) and the 4-and 5-yr old fields (high PC2 value) clustered apart from each other and formed two new clusters with low variability within each cluster. At the ESV level, fields from each of these clusters were further separated according to year of recultivation. Decomposition of MNTD and MPD beta diversity into classes of low to high phylogenetic dissimilarity revealed estimates of community phylogenetic turnover across environmental gradients similar to the BDTT analysis. Communities within an MNTD of 0.03 subs/site and an MPD of 0.3-0.6 subs/site were positively correlated with PC1 (i.e. these communities tend to occur in samples of similar PC1 values), whereas communities with higher phylogenetic dissimilarity (MNTD > 0.03 subs/site and MPD > 0.03-06 subs/site) were negatively correlated with PC1 (i.e. tend to occur in samples of dissimilar PC1 values; Fig. 4a,b). The correlation with PC2 was difficult to decompose because the relationship with MNTD and MPD was nonmonotonic, indicating low phylogenetic turnover for high PC2 distances (Fig. S8); yet, the positive correlation of MPD to PC2 was observed for lower genetic distance classes than PC1; and at high distance classes, MPD negatively correlated to PC2 while it was still positively correlated to PC1 (Fig. 4b).
In the Soil depth dataset, soil layer was the unique predictor of beta diversity and was significant at all genetic distances (Figs 3, S5), but, again, the explained variance varied strongly. Soil layer explained 28% of the variance at the ESV level and rose to 45% to reach a plateau at 0.03 subs/site. It then increased slightly from 0.15 to 0.16 subs/site. Soil compartment effects were never significant. MNTD was strongly correlated with soil layer within the first genetic distance class of 0.03 subs/site (Fig. 4c). MPD was positively and negatively correlated at intermediate and high distance classes, respectively (Fig. 4d). Neither MNTD nor MPD correlated with soil compartment.
In summary, our results show that clustering at c. 0.03 subs/ site has lumped ESVs into phylogenetically and ecologically coherent units that strongly correlate with the environment. They also show a turnover of distantly related AM fungal clades (> 0.03 subs/site) across PC1 and soil layer, and a turnover of relatively closely related AM fungal clades (c. 0.03 subs/site) across PC2. PERMANOVA P values did not change across the 100 bootstrapped trees, supporting that these results were robust to phylogenetic uncertainty. PERMANOVA R 2 were not comprised within the 5-95% distribution of null distribution of R 2 from 100 phylogenetic randomizations, further indicating that these New Phytologist (2019)

Research
New Phytologist results derived from a genuine phylogenetic signal and that none of the changes in the correlation strength across phylogenetic resolutions could be solely attributed to change in OTU number.

Phylogenetic signal of ESV ecological traits
In the Chronosequence dataset, we observed a positive autocorrelation (Mantel's R = 0.05) in ESV co-occurrence up to deep phylogenetic distances classes (Fig. 5a). However, PC1 optimum had the highest positive autocorrelation (Moran's I = 0.3) and the signal was significantly positive up to 0.4 subs/site patristic distance and significantly negative at largest genetic distances (Fig. 5b). PC2 optimum showed positive autocorrelation (Moran's I = 0.09) up to 0.2 subs/site patristic distance (Fig. 5c). We observed similar results in the nonLULU curated dataset, but the phylogenetic signal for PC2 was deeper (Fig. S9). In the Soil depth dataset, a positive phylogenetic signal was significant only for co-occurrences (Fig. 5d), but in the noncurated LULU dataset www.newphytologist.com the soil layer optimum, and soil compartment to a lesser extent, was positively autocorrelated for ESVs as distant as 0.6 subs/site patristic distance (Moran's I = 0.2) and negatively autocorrelated after 1 sub/site patristic distance (Fig. S9). These results are consistent with the beta diversity analyses. They indicate that closely related ESVs tend to be found in the same samples; and if not, then in samples with similar soil conditions. PC1 and soil layer optimum had the strongest phylogenetic signal, whereas PC2 optimum was conserved at a shallow phylogenetic depth.

Identity of clades that covary with the environmental predictors
In the Chronosequence dataset, multiple nodes at all phylogenetic depths showed significant correlation with PC1: notably, Fig. 2 Maximum-likelihood phylogenetic tree of arbuscular mycorrhizal fungal exact sequence variants (ESVs) for the Soil depth dataset. The phylogenetic clustering of ESVs at a genetic distance of 0.03 mean substitutions per site is highlighted in blue and gray rectangles. Significant correlation between the abundance of clades (nodes) across the phylogeny and soil layer (blue) or soil compartment (orange) is depicted at each node. Bar, 0.1 mean substitutions per site. The heatmap on the right of the tree shows the centered and scaled soil layer and soil compartment optima. (2019)

Research
New Phytologist members of Rhizophagus, Funneliformis, and Claroideoglomus, including the Glomerales node and to a lesser extent Diversispora, are mostly found in young fields (low PC1), whereas Archaeospora members are found in old fields (high PC1; Figs 1a, S10). By contrast, only nodes mostly at or contained within 0.03 subs/site genetic distance correlated with PC2, including members of Funneliformis, Diversispora, and Claroideoglomus, found at either low or high PC2 values. In the Soil depth dataset, the Claroideoglomus node and its members were found in subsoil, whereas Diversispora and the node Diversispora/Scutellospora (Diversisporales) were found in topsoil (Figs 1b, S11). The significant increase in (relative) abundance of the Glomerales node with soil depth is probably driven by a strong increase in the abundance of Claroideoglomeraceae and does not correspond with an increase in abundance of other members of this node (i.e. Funneliformis, Glomus, and Rhizophagus).

Discussion
In this study, we tested whether and how beta diversity patterns change with OTU phylogenetic resolution from the ESV to order-level sequence clusters. At the ESV level, we were able to detect changes in community composition similar to traditional OTU clustering methods or at coarser phylogenetic resolutions. However, we showed a stronger correlation between the environmental predictors and beta diversity assessed at coarser phylogenetic resolutions, especially at c. 0.03 subs/site. Furthermore, we observed a turnover in either distantly or closely related clades along different environmental gradients, such that the variance explained by different environmental predictors varied with OTU phylogenetic resolution. This suggests that accounting for the evolutionary history of AM fungi could help better understanding of the drivers of their spatial distribution. At the ESV level, we statistically detected the environmental predictors revealed with traditional OTUs or at higher genetic distances, consistent with previous results for bacteria (Thompson et al., 2017) and fungi using the ITS marker (Glassman & Martiny, 2018). This was because closely related ESVs tend to occur in the same samples. However, we showed a strong increase in the correlation between most environmental predictors and beta diversity at a genetic distance of c. 0.03 subs/site, and even

Research
New Phytologist higher for certain environmental predictors. This shows that clustering ESVs into OTUs at 0.02-0.04 subs/site (and 0.03-0.05 subs/site for the nonLULU-curated dataset), or higher, reduced the noise in the biological signal by lumping closely related ESVs that occur in otherwise similar environments but in different communities. These results support the idea that fungal community assembly, and in particular AM fungal communities, are governed by environmental selection of organisms exhibiting strong phylogenetic niche conservatism (Lekberg et al., 2014;Powell & Sikes, 2014;Botnen et al., 2018). In this case, given that AM fungi exhibit a significant degree of phylogenetic niche conservatism, clustering ESVs into OTUs at c. 0.03 subs/site or higher may serve a better purpose for niche modeling by increasing the power to detect correlations with environmental predictors (Powell & Sikes, 2014).
Interestingly, for both the Chronosequence and Soil depth datasets, the beta diversity approaches based on a phylogenetic tree (BDTT, MNTD, and MPD) converged towards c. 0.03 subs/site as a level of evolutionary divergence carrying a strong ecological specialization signal. This threshold to delineate AM fungal OTUs with LSU will certainly be lower for SSU and higher for ITS (Thi ery et al., 2016); this may change when studying other organisms and, as shown here, with other environmental gradients. Interestingly, the number of OTUs plateaued after 0.03 subs/site (both after and before LULU curation). These characteristics indicate that OTUs at 0.03 subs/site can be defined as ecologically coherent biological units, which could be species (Bruns et al., 2018). More broadly, biological units that are simultaneously evolutionarily and ecologically coherent can be functional groups amenable to quantitative PCR monitoring.
Our results suggest an additional ecological structure in beta diversity at the ESV level: in the Chronosequence dataset, we observed a weak but significant correlation between beta diversity and PC3, and a split among each fields recultivated from 1 to 5 yr, which was not the case at coarser phylogenetic resolutions, including when considering OTU UPARSE . This result is supported in our sensitivity analysis by the slight clustering of beta diversity patterns at the ESV level apart from the other beta diversity patterns inferred at relatively fine phylogenetic resolutions. What biological unit an ESV captures remains unknown, however, and it would ultimately depend on the rate of evolution of the targeted gene and clades. The observed diversity of LSU variants within a range of sequence divergence from 0 to 0.03 subs/ site is consistent with estimates of intraspecific diversity obtained from single isolates and spores within Rhizophagus irregularis and Gigaspora margarita (Thi ery et al., 2016), suggesting that successional population processes were captured. Through ESV inference, we could separate the process of sequence denoising from the biological unit delineation, with the potential to infer the ecological variation at intra-vs interspecific levels.
Our results show that the degree of phylogenetic turnover among AM fungal communities differs along different environmental gradients. We observed a turnover of distantly related clades across PC1 and soil depth (clades diverging at least by 0.2 subs/site), whereas clades at or within a genetic distance of c. 0.03 subs/site varied along PC2. These results suggest that ecological specializations to different aspects of the environment (e.g. soil P and N content, soil depth) are conserved at different phylogenetic depths in AM fungi. This could reflect the primary factors imposing pressure on the AM symbiosis. For instance, fungi reported to have a ruderal strategy, such as Funneliformis, Rhizophagus, and Claroideoglomus within the Glomeraceae (Hart & Reader, 2002;Chagnon et al., 2013), were abundant in young fields, whereas Archaeospora only occurred in soils recultivated for > 10 yr. PC1 represented a strong temporal gradient, with steep variations in plant-available P content. The occurrence and abundance of AM fungi across this gradient are probably based on their ability to disperse, colonize, and persist in soil and roots, and maintain an association with plants in soils of different P content. The strong mutual exclusion of Diversispora/ Scutellospora found in topsoil and of Claroideoglomus found in subsoil (Soil depth dataset) and in young fields where the subsoil : topsoil mix has been recently deposited (Chronosequence dataset) indicates that some Claroideoglomus fungi cannot persist in topsoil and/or have specialized to subsoil (Sosa-Hern andez et al., 2018b), a habitat with low plant-carbon supply probably selecting for stress-tolerant fungi. Our results are consistent with colonization-persistence or competitors-stress tolerators-ruderals life-history strategies probably being conserved at the family level (Hart & Reader, 2002;Chagnon et al., 2013). They are also consistent with P uptake being the primary function of AM fungi (Smith & Read, 2008). PC2 represented a strong increase in inorganic N availability but also in N : P in fields recently converted to agriculture. The turnover of closely related AM fungi along PC2 contradicts the idea that AM fungal association with N availability (plant-available N in our case) is conserved at family or genus level (Treseder et al., 2018). This discrepancy could be due to the difference in the study ecosystems (natural or agricultural ecosystems), to fungal physiological requirement (e.g. stoichiometry) but also to the difficulty in disentangling confounding abiotic and/or biotic factors. Along with investigations on isolates, further research conducted at multiple phylogenetic resolutions will certainly help test whether soil P, more than soil N content, could have driven the early diversification of AM fungi.
Inferring life-history strategies and functions in fungi remains difficult (Powell & Rillig, 2018;Treseder et al., 2018). We show that one could infer such traits from distribution data based on DNA sequencing when coupled with contextual parameters, which is a promising result given that DNA-based assessments of biodiversity are now accumulating. Increasing spatial extent and decreasing spatial grain will increase environmental heterogeneity and capture environmental gradients of varying steepness and with additional importance for the AM symbiosis (e.g. vegetation characteristics, climate) in addition to an increased breadth of AM phylogenetic diversity. It will further reveal the ecological variation at different levels of phylogenetic divergence in AM fungi and narrow the estimated depth of phylogenetic conservatism of these ecological traits. For example, traits variation at very fine phylogenetic resolution could be relevant for subtle environmental variation (Powell & Sikes, 2014)  instance, between extra-or intraradical mycelium (if different portions of the mycelium express functionally different ribosomes; Gilbert, 2011;Filipovska & Rackham, 2013), when considering the selection of different nucleotypes at the arbuscule scale (Limpens & Geurts, 2014), or when other deterministic processes occur (Verbruggen and Kiers 2010). At the other extreme of phylogenetic depth, differences between biomes or ecoregions could be clear at broad phylogenetic resolutions (St€ urmer et al., 2018b). Alternatively, a biogeographic structure might be revealed at the ESV level in the case of dispersal limitation of cryptic taxa and stochastic population-level processes (e.g. genetic drift; Bruns & Taylor, 2015;Savary et al., 2018). A phylogenetic decomposition of beta diversity may prove difficult using ITS, a genetic marker used by most molecular fungal ecologists (Schoch et al., 2012), due to limitations in phylogeny reconstruction. Clustering ESVs based on (morphology-based) taxonomy (Cline et al., 2017) or using traditional clustering methods (Lekberg et al., 2018) can alternatively be used to reveal significant changes of community composition across environmental gradients. Our sensitivity analysis showed both methods can be used to assess change in community composition at different phylogenetic resolutions. Finally, phylogenetic placement methods, together with the recent release of improved phylogenies for the entire fungal kingdom (e.g. Tedersoo et al., 2018) will certainly anchor fungal molecular ecology using ITS into a phylogenetic framework.

Conclusion
The establishment of an organism in a community is determined by a suite of processes, including ecological drift, dispersal, and abiotic and biotic selection. The interaction of these processes imprints on the phylogenetic structure of communities at multiple phylogenetic resolutions. We used this information to identify at which phylogenetic resolution AM fungal communities covary with different environmental factors. We showed that beta diversity patterns of AM fungi change with OTU phylogenetic resolution, so that one phylogenetic resolution gives only partial information on the factors shaping their spatial distribution. We suggest that assessing beta diversity patterns at varying OTU phylogenetic resolutions offers a deeper understanding of the ecology and evolution of fungi. Identifying ecological variation among taxa of AM fungi at different phylogenetic depths will ultimately help define functional groups and serve as a guide for the engineering of AM fungal communities (Rillig et al., 2016).

Supporting Information
Additional Supporting Information may be found online in the Supporting Information section at the end of the article.