Biotic and abiotic factors shape arbuscular mycorrhizal fungal communities associated with the roots of the widespread fern Botrychium lunaria (Ophioglossaceae).

Arbuscular mycorrhizal fungi (AMF) play central roles in terrestrial ecosystems by interacting with both above and belowground communities as well as by influencing edaphic properties. The AMF communities associated with the roots of the fern Botrychium lunaria (Ophioglossaceae) were sampled in four transects at 2400 m a.s.l. in the Swiss Alps and analyzed using metabarcoding. Members of four Glomeromycota genera were identified across the 71 samples. Our analyses revealed the existence of a core microbiome composed of four abundant Glomus OTUs, as well as a low OTU turnover between samples. The AMF communities were not spatially structured, which contrasts with most studies on AMF associated with angiosperms. pH, microbial connectivity and humus cover significantly shaped AMF beta diversity but only explained a minor fraction of variation in beta diversity. AMF OTUs associations were found to be significant by both cohesion and co-occurrence analyses, suggesting a role for fungus-fungus interactions in AMF community assembly. In particular, OTU co-occurrences were more frequent between different genera than among the same genus, rising the hypothesis of functional complementarity among the AMF associated to B. lunaria. Altogether, our results provide new insights into the ecology of fern symbionts in alpine grasslands. ORIGINALITY-SIGNIFICANCE STATEMENT: This study explores for the first time and at a regional scale, the diversity and community structure of arbuscular mycorrhizal fungi (AMF) associated with the widespread fern Botrychium lunaria (Ophioglossaceae). Most community ecology studies on AMF have focused on angiosperm hosts, leaving early divergent terrestrial plants largely understudied. In addition, AMF from alpine grasslands have received less attention than those from low altitude grasslands and arable soils. Metabarcoding of AMF from four transects located at an altitude of 2400 m a.s.l. in the Swiss Alps allowed us to comprehensively assess the AMF diversity associated to B. lunaria and also to identify a core microbiome composed of Glomus OTUs. This article is protected by copyright. All rights reserved.


Introduction
Arbuscular mycorrhizal fungi (AMF) are important players in terrestrial ecosystems due to their ability to form mutualistic symbiotic associations with approximately 72% of terrestrial plant species (Brundrett and Tedersoo, 2018), including about 54% of pteridophyte species (Lehnert et al., 2017). AMF are obligate biotrophs that, as a rule, trade nutrients (phosphorus, nitrogen and others) with plant hosts for carbon (carbohydrates and lipids) (Hodge et al., 2010). In addition to this, AMF also increase plant resistance to both abiotic (e.g. drought) (Augé et al., 2015) and biotic (e.g. pathogens) stresses (Gianinazzi et al., 2010;Veresoglou and Rillig, 2012). However, there are exceptions, with several plant species actually exchanging little, or even no carbon to the fungal partner (Leake, 2005;Merckx et al., 2009;Selosse and Roy, 2009). This is the case, for instance, in Lycopodiaceae (the lycopods, or club mosses) that have a nonphotosynthetic life-stage for several years and are thus unable to exchange any carbon to their fungal partners (Winther and Friedman, 2008). Furthermore, beyond the nutritional exchange, AMF also play central roles at the ecosystem level through modifications of soil physical (e.g. aggregation) and chemical (e.g. nutrient concentrations, leaching) properties, as well as by influencing the composition of both plant and soil microbial communities (van der Heijden et al., 1998;Rillig, 2004;Powell and Rillig, 2018).
Depending on the authors, AMF are classified as belonging to the Glomeromycota phylum (Schüβler et al., 2001) or to the Glomeromycotina sub-phylum, Mucoromycota phylum (Spatafora et al., 2016). These taxa are ubiquitous in terrestrial ecosystems and several studies have investigated their diversity and distribution across the globe (Kivlin et al., 2011;Davison et al., 2015;Stürmer et al., 2018). It is assumed that AMF are more prominent in grasslands and tropical forests as compared to temperate or boreal forests (Allen et al., 1995). Treseder and Cross (2006) highlighted that a major driver of AMF abundance was fine root biomass, whatever the biome. However, most of these studies have focused on AMF associated with angiosperm hosts, leaving other plant clades such as pteridophytes poorly investigated (Berch and Kendrick, 1982;Gemma et al., 1992;Liu et al., 2004;Muthukumar and Prabha, 2013;Pressel et al., 2016). Molecular data on fern-associated AMF diversity based on high-throughput sequencing are even more scarce (Huang et al., 2018). Thus, there is currently a gap regarding the identity and the ecology of AMF associated with pteridophytes, which are a rich and widespread plant group (Kreft et al., 2010). Moreover, a large portion of these plants associate with AMF (Pressel et al., 2016).
Botrychium lunaria (Ophioglossaceae), referred to as the common moonwort, is a widespread fern distributed throughout the northern hemisphere and a few notable disjunctions to the southern hemisphere (Dauphin et al., 2014(Dauphin et al., , 2016(Dauphin et al., , 2017(Dauphin et al., , 2018Stensvold and Farrar, 2017). In Switzerland, B. lunaria is widespread in the Alps and Jura mountains Dauphin et al., 2018). It belongs to the Ophioglossaceae family, which with its sister clade Psilotaceae, belongs to an early divergent lineage that diverged from the other ferns approximately 300 million years ago (Pryer et al., 2004). Ophioglossaceae are known to be colonized by AMF, but to date most taxa from this family have not been investigated for their fungal symbionts (Lehnert et al., 2017). B. lunaria seems generalist as for its habitat, although it is considered to be a pioneer species, and highly sensitive to plant competition (Muller, 1999). As all members of the Ophioglossaceae family (Smith et al., 2008), B. lunaria is known to form subterranean, achlorophyllous and obligatory mycorrhized gametophytes (Schmid and Oberwinkler, 1994). It has also been demonstrated that both the gametophyte and the sporophyte of B. crenulatum and B. lanceolatum are colonized by AMF (Winther and Friedman, 2007).
To date, molecular studies investigating the taxonomic diversity of AMF colonizing Botrychium ferns are restricted to a cloning and Sanger sequencing approach (Kovács et al., 2007;Winther and Friedman, 2007), implying that not all the fungal diversity may have been recovered due to low sequencing depth. Additionally, biotic and/or abiotic factors structuring the AMF communities associated with Botrychium species roots have not yet been investigated. However, it is now well established in angiosperm plant communities that various factors may be involved in structuring AMF communities such as the effect of the host plant species (Sepp et al., 2019), neighbouring plants (Hausmann and Hawkes, 2009), soil pH (Dumbrell et al., 2010), AMF phylogenetic distance and spatial distance , effect of succession (Gao et al., 2019), effect of geographic distance, soil moisture, and soil temperature (Kivlin et al., 2011), as well as slope (Chai et al., 2018). Furthermore, the relative roles of these abiotic and biotic factors are known to vary depending on the spatial scale of analysis (Vályi et al., 2016;Caruso, 2018).
Our study aims at describing the AMF communities associated with the moonwort and at quantifying the factors that structure these communities. Therefore, we used high-throughput sequencing of 18S rRNA gene amplicons to investigate the taxonomic diversity of AMF associated with the roots of B. lunaria from four geographically distinct populations in the Swiss Alps. Furthermore, we investigated AMF community structure with co-occurrence and cohesion analyses. Finally, abiotic (soil physicochemical parameters, field characteristics) and biotic (floristic composition) factors, that are known to influence the community composition of AMF associated to angionsperms, were also assessed for each sampling site in order to identify and quantify the factors involved in shaping AMF communities of B. lunaria populations (Fig. 1).

Results
After quality filtering and random subsampling, 1 0 194 0 362 Glomeromycota sequences originating from 71 root samples of Botrychium lunaria were clustered into 53 OTUs. Analysis of both the rarefaction curves and the samplebased OTUs accumulation curve indicated that sequencing depth and sample size were sufficient to estimate and compare the AMF diversity of all the samples (Supporting Information Fig. S2). This was also supported by a high Good's coverage estimator of each sample (coverage > 0.99). Regarding the richness, there was on average 26 OTUs per sample and 48 OTUs per transect, indicating a low OTU turnover. No significant difference in the patterns on alpha diversity was observed between the four transects, considering either observed richness, Shannon diversity or the inverse of Simpson index (Kruskal-Wallis test, P > 0.05; Supporting Information Fig. S3). Additionally, no significant correlation was found between AMF OTU richness and any of the floristic or edaphic variable measured in the present study. However, a negative and significant relationship (r = −0.54, P < 0.001) was observed between the richness of the neighbouring plants surrounding B. lunaria sporophytes and the soil pH (Supporting Information Fig. S4).
Probabilistic modelling based on Dirichlet multinomial mixture models was used as unsupervised clustering method to partition AMF communities into community types. The optimal cluster number was two, further referred to as community types A and B. Principal coordinates analysis (PCoA) based on Canberra distances (Kuczynski et al., 2010) confirmed the existence of two community types (Fig. 3A). Additionally, AMOVA indicated that variation among community types [sum of squares (SS) = 1.14] was much lower than within community types (SS = 7.71) and that the observed separation between these two community types was statistically significant (P < 0.001). This latter result was confirmed by ANOSIM (R = 0.43, P < 0.001). Community type A included samples from the four transects, whereas community type B contains samples from three of the four transects. Then, edaphic and floristic variables were investigated within these two community types. We found that soil pH was significantly lower in community type A (mean pH value 5.34) compared to B (mean pH value 6.33) (Wilcoxon test, P < 0.001; Fig. 3B). Using compositional balances (Rivera-Pinto et al., 2018), five OTUs signatures of these two community types were identified (Fig. S5). We found that two OTUs (Glomus MOG31 and Acaulospora Acau3) were preferentially found in the community type A associated with more acidic soil pHs, whereas three OTUs (unclassified Glomus) were preferentially found in the community type B.
To further disentangle the respective contribution of geography, floristic composition, and soil properties in explaining the variation in AMF community structure, an analysis of variation partitioning was performed (Fig. 4). The largest part of the variation explained by geography (i.e. latitude and longitude) was shared with the soil properties, indicating that edaphic and spatial variables were strongly intercorrelated. Conversely, the proportion of variation in community composition explained by purely spatial structure was extremely low (adjusted R 2 = 0.006), highlighting the very weak influence of only geographic distance as part of the factor geography, and thus, indicating an absence of spatial pattern. Variation explained by floristic (adjusted R 2 = 0.017) and soil (adjusted R 2 = 0.020) variables with no spatial structure was around three time higher than the spatial variables. A significant explained variation was also shared by the interaction of the three matrices (adjusted R 2 = 0.04). Yet, a large fraction of the variation was not explained by any of these three variables (adjusted R 2 = 0.866), indicating that these variables were rather poor predictors of the AMF community structure. Co-occurrence analyses based on a presenceabsence matrix of AMF OTUs associated with B. lunaria roots across all samples revealed 9.8% of non-random co-occurrences including 94 positive and 29 negative ones and thus indicating some degree of association among AMF OTUs. Positive co-occurrences indicate that several OTUs were found associated with a greater degree than expected by chance only, while negative cooccurrences suggest exclusion between AMF OTUs. At the genus level, we observed 33 positive co-occurrences among Glomus, 17 between Glomus and Acaulospora, 14 between Glomus and Claroideoglomus, seven between Claroideoglomus and Acaulospora as well as seven among Acaulospora and four among Acaulospora (Supporting Information Table S2). On the other hand, we found nine negative co-occurrences between Glomus and Claroideoglomus, five among Glomus, four between Glomus and Acaulospora, and four between Glomus and Ambispora (Supporting Information Table S3). Overall, there were more co-occurrences between OTUs from distinct genera (49 positive and 22 negative ones) than among OTUs of the same genus (45 positive and 7 negative ones).
Because of the existence of these non-random patterns of co-occurrence, we further quantified the degree of connectivity (positive and negative cohesions) and complexity of associations among the AMF communities, using a recently developed statistical method named cohesion (Herren and McMahon, 2017). First, 10 OTUs with the lowest negative and the highest positive connectedness were investigated, representing potential keystone taxa (Herren and McMahon, 2018). OTUs with the lowest negative connectedness (nine assigned to Glomus and one to Claroideoglomus) corresponded to abundant OTUs, including the four members of the core microbiome. Interestingly, two of the three OTUs occurring preferentially in the community type B (higher soil pH) were also among these negatively connected OTUs. On the contrary, OTUs with the highest positive connectedness (eight assigned to Glomus, one to Ambispora and one to Claroideoglomus) corresponded to low abundance OTUs (relative abundance ranging from 0.05% to 4.67%).
In order to evaluate the influence of both abiotic and biotic factors on the AMF community composition, multivariate analysis of deviance of generalized linear models (GLMs) was applied. Forward selection on the explanatory variables identified positive cohesion, negative cohesion, pH and humus cover as significant variables (Table 1). Neighbouring plants, via the humus cover quantified by the Landolt's indicator H, showed a significant effect on AMF community composition, but with a lower explained deviance compared to the three other variables mentioned above (Dev = 83, P < 0.005). This indicator ranges from 1, for a bare soil with no humus cover, to 5 for raw humus or bog, with plants rooting particularly in humus. As an abiotic factor, pH showed a higher explained deviance (Dev = 573, P < 0.001) compared to the other biotic factors (Table 1). Both cohesion metrics contributed significantly to explain AMF community composition, with slightly higher explained deviance for the negative cohesion than for the positive one (Dev = 350, P < 0.001 for negative cohesion; Dev = 326, p < 0.001 for positive cohesion). As mentioned previously, soil pH was also  significantly different between the two community types identified by unsupervised clustering. A similar trend was observed for the positive cohesion (P < 0.05), but not for the negative cohesion (Supporting Information Fig. S6). Finally, the difference in humus cover was only marginally significant (P < 0.09).

Discussion
Previous studies have investigated AMF communities in the Alps, either based on spore morphology (Read and Haselwandter, 1981;Oehl et al., 2011;Oehl and Körner, 2014) or on molecular techniques such as DNA fingerprinting and cloning-sequencing (Sykorová et al., 2007;Casazza et al., 2017). With the emergence of nextgeneration sequencing, and the bypass of cloning steps, it is now possible to study the species diversity and structure of AMF communities in natural populations. However, on the regional or global scale, knowledge about the molecular diversity of AMF communities in high altitudes (> 2000 m a.s.l.) remains poor (Bueno de Mesquita et al., 2018). Here, we showed both the taxonomic diversity and the community structure of the AMF associated with the roots of the fern Botrychium lunaria from alpine grasslands at 2400 m a.s.l. Using high-throughput sequencing of 18S rRNA gene amplicons from 71 roots samples, the AMF diversity was fully covered and resulted in 53 OTUs. Despite the apparent lack of specificity of plants for their AMF symbionts (there are many more plants species than AMF species), several studies have demonstrated that old plant lineages such as liverworts (Russell and Bulman, 2004), Lycopodiaceae (Winther and Friedman, 2008), and Ophioglossaceae (Field et al., 2015) were usually associated with Glomeromycota belonging to the Glomus genus. The long co-evolution between these plant lineages and their symbiotic fungal partners may explain this apparent specificity (Brundrett, 2002). To date, only a few Botrychium species have been investigated for their AMF composition and this was done solely using a cloning sequencing approach. Winther and Friedman (2007) reported that B. crenulatum and B. lanceolatum were colonized by members of the Glomus genus. On the other hand, it was found that B. virginianum was colonized by both members of Glomus (West et al., 2009) and Scutellospora (Kovács et al., 2007). In this study, we confirmed the presence of Glomus as a dominant member of the root-associated AMF community of B. lunaria and go further in describing the taxonomic composition. We reported for the first time the presence of Acaulospora (nine OTUs), Ambispora (three OTUs), and Claroideoglomus (nine OTUs) genera. The existence in our data set of a core microbiome of four Glomus OTUs representing 62% of the reads strengthened the importance of this fungal genus among AMF associated to Botrychium and more generally with ferns (West et al., 2009). Using the online database MaarjAM (Opik et al., 2010), we also found that, in previous studies, most of the sequences associated with the Botrychium genus belonged to Glomus. Thus, our data bring new insights into the AMF diversity associated with the widespread Botrychium species and more generally to alpine AMF ecology.
Regarding OTU richness, no significant correlation was found with any of the floristic or edaphic variables measured in the present study. This contrasts with other studies in alpine grasslands (Xiang et al., 2016), orchards (Van Geel et al., 2015), and tropical montane forests (Camenzind et al., 2014) where N and P contents significantly influenced AMF OTU richness. Because AMF colonize both soil and roots (Vályi et al., 2016), they depend on both the edaphic conditions and the host plant (Vályi et al., 2015), with the latter controlling the initiation and degree of colonization by signalling compounds (Schmitz and Harrison, 2014). As a result, in order to settle into plant roots and in agreement with the modern coexistence theory (HilleRisLambers et al., 2012), AMF have to pass through two filters: an environmental (i.e. the soil) and a biotic filter (i.e. the plant host). Since we only analyzed root samples from a single fern species, the present study was not designed to test the effect of the host plant on the AMF communities. However, the presence of a core microbiome composed of four abundant OTUs strongly suggests that, at local scale, a large fraction of the AMF community is constantly associated with the host, regardless of the geographical distance between plants or of the edaphic parameters. B. lunaria might act as a biotic filter and contributes partially to the structure of its root-associated AMF community. As such, the importance of the core microbiome and the very low OTU turnover observed between samples, could also be explained by a host symbiont specificity. Indeed, such effect has been observed between plants and AMF (Becklin et al., 2012;Sepp et al., 2019). Both of these hypotheses should be tested in further experiments including several fern species and soil samples.
On the other hand, the effect of the environment on the AMF community structure was detected in different analyses. First, unsupervised clustering analysis identified two community types found in soils with significantly different pH (Fig. 3), suggesting that pH conditions might drive AMF community assembly of B. lunaria roots. Five OTUs were identified as microbial signatures of these two community types (Supporting Information Fig. S5) and among them, one OTU assigned to the Acaulospora genus was preferentially found in the more acidic samples (community type A). Interestingly, this genus has already been shown to be associated with acidic soils (Morton, 1986;Coughlan et al., 2000). Second, multivariate analysis of the deviance of generalized linear models confirmed the major influence of pH on AMF community structure, since pH was the significant factor explaining the highest deviance (Table 1). The major role of pH in structuring AMF communities has already been reported from studies with wide soil pH conditions (pH 3.72 to 8.04) (Dumbrell et al., 2010). As already shown by Lekberg et al. (2011) in grasslands, our results highlight that even within a more narrow range (pH 4.56-6.79), pH is a strong predictor of AMF community composition. It should also be pointed out that the effect of the pH gradient might reflect a spatial gradient, since in our study, the edaphic and spatial variables were intercorrelated (Fig. 4). Yet, the two community types were both present in three of the four sampled populations (on both Mase and Forclaz sampling locations, see Fig. 1), indicating that dispersal limitation was not a structuring factor at the scale of our study. This was confirmed by the variation partitioning analysis that revealed a negligible proportion of variance explained by purely spatial factors (Fig. 4). Additionally, soil pH was found to be negatively correlated with the species richness of the plants surrounding B. lunaria (Supporting Information Fig. S4), indicating that pH is a driver of both above-and belowground diversity. This pattern can be explained by the habitat hypothesis introduced by Zobel and Öpik (2014), where both plant and AMF communities follow changes in abiotic conditions. Three of the four sampled transects were grasslands at a climax state. Therefore, in such stable ecosystems at late succession, only the habitat hypothesis can explain spatial co-variation of plant and AMF communities, and this co-variation was also identified in the variation partitioning analysis (Fig. 4).
Although the influence of plant communities on AMF community composition was reported in grassland microcosms (Johnson et al., 2004;Hausmann and Hawkes, 2009), here, we did not identify a direct effect of the floristic composition on the AMF community structure. However, we found that one of the derived Landolt's ecological indicators, namely humus cover (H), significantly influenced the AMF community composition (Table 1). Humus harbours diverse soil taxa, including mycorrhizal fungi, and concentrates various biological activities (Ponge, 2003). It is also a central place and a particular environment for aboveground-belowground interactions, especially plant-soil feedbacks (Ponge, 2013). AMF have been reported to be indirectly involved in litter decomposition via the stimulation of microbial decomposers (Bunn et al., 2019). Therefore, these AMFdecomposers interactions could be under the control of the humus forms (Ponge, 2013). Alternatively, humus cover could influence AMF community composition via its impact on soil temperature and moisture. To our knowledge, this is one of the few studies reporting a significant relationship between ecological indicator values inferred from floristics and AMF community composition, highlighting again the link between above-and belowground components in the grassland ecosystem. Thus, this result also calls for a more systematic use of such ecological indicator values when studying AMF communities, especially in alpine environments. It should be underlined that, altogether spatial distance, soil physicochemical properties and floristic composition only explained a minor fraction of the AMF community composition variance and therefore, that a large proportion of variance remained unexplained (Fig. 4). This could be due to a strong selective effect of B. lunaria, as discussed above. Additionally, biotic interactions, among AMF or between AMF and other soil organisms, could impact AMF community structure.
Indeed, another factor structuring the AMF communities associated to B. lunaria was the potential interactions within the AMF communities, as demonstrated by cooccurrence and cohesion analyses. The role of non-host biotic interactions in structuring AMF communities has already been discussed (Vályi et al., 2016). For instance, the importance of interactions within the AMF community of Festuca brevipila in grasslands with environmental gradients has been suggested (Horn et al., 2017). Here, we confirm this idea with quantitative data using the cohesion metrics that measure the degree of connectivity of a microbial community (Herren and McMahon, 2017). Both positive and negative cohesions contributed significantly to explain AMF community composition, with relatively similar explained deviance values (Table 1). Together, these two biotic factors explain more deviance than the pH and the humus coverage, showing the importance of AMF connectivity compared to edaphic factors. Cohesion metrics have been used to describe community complexity (negative cohesion) and community dynamics (positive cohesion) (Danczak et al., 2018). On closer examination, we observed that OTUs with the lowest negative connectedness represented abundant OTUs including the members of the core microbiome, here supporting that negative cohesion reflects community stability and low turnover within this community. Conversely, OTUs with the highest positive connectedness represented low abundance OTUs, occurring in fewer samples and thus reflecting the fraction of the community with a higher turnover. This also highlights the importance of low abundant OTUs in the AMF community structure. The existence of these two types of OTUs could be the results of different life-history strategies of the AMF (trade-off between colonization and persistence; competitors, stress tolerators and ruderals AMF) (Hart et al., 2001;Chagnon et al., 2013), with the abundant and prevalent OTUs showing high persistence abilities, while less abundant and less prevalent OTUs would present traits promoting rapid colonization abilities. Such low abundant OTUs could only be identified in situ with an appropriate sequencing depth of the amplicons, reinforcing the need for high-throughput sequencing to properly describe the molecular diversity and community structure of AMF. Open questions remain about the importance of these low abundant AMF OTUs for ecosystem functioning (Jousset et al., 2017).
Regarding the AMF co-occurrences, we always found more co-occurrences between OTUs of different genera (49 positive and 22 negative ones) than among OTUs of the same genus (45 positive and 7 negative ones). This observation contrasts with results obtained from soil samples collected across a European transect (Bouffaud et al., 2016). Our co-occurrence pattern could support different interpretations. The high number of positive cooccurrences between OTUs from different genera (i.e. phylogenetically distant OTUs, with likely different functional characteristics) could be the result of competitive exclusion, which would prevent the co-occurrence of closely related and functionally similar OTUs. Indeed, Glomeraceae, Gigasporaceae and Acaulosporaceae present distinct functional traits, which are phylogenetically conserved Klironomos, 2007, 2012). Thus, B. lunaria could benefit from functional complementarity of OTUs belonging to different genera and/or families. The detection of negative co-occurrences could be the results of competitive interactions between AMF OTUs (Thonar et al., 2014). AMF could, for instance, compete for root space. Because we identified two community types and a significant effect of soil pH on the AMF community structure, these negative co-occurrences may also reflect environmental preferences (here, pH-related) of certain AMF OTUs leading to the exclusion of these OTUs if abiotic conditions are not suitable for their survival (i.e. non-overlapping niches between OTUs). Lastly, the higher number of positive compared to negative cooccurrences, associated with the low OTU turnover observed between our samples, support one more time the hypothesis of B. lunaria having a strong selective effect on AMF community composition.

Conclusions
We provide, at a regional scale, the first molecular inventory of the AMF associated with the widespread fern Botrychium lunaria. Community analysis revealed the existence of an abundant core microbiome and of a low OTU turnover. Contrary to what is generally observed with angiosperm hosts, we found no spatially structured pattern of the community. Soil properties (mainly pH) and aboveground vegetation (humus cover) significantly impacted the AMF community composition but overall, remained relatively poor predictors. OTU associations were not random, indicating a potential role of AMF interactions in community structure. Lastly, even if more data will be needed to confirm this hypothesis, our data suggest that B. lunaria has a strong selective effect on the composition of its AMF symbionts (i.e. host specificity).

Experimental procedures
This study was carried out in the Val d'Hérens (Valais canton, Switzerland) in June 2015. Four alpine populations (IT1, IT2, IIT1, IIT2: Fig. 1A) of Botrychium lunaria were sampled at about 2400 m above sea level. Each population consisted of six sampling plots, which were organized in two technical replicates (IT1/IT2 and IIT1/IIT2), each replicate made up of three plots (Fig. 1B-E). Finally, each of these plots contained three quadrats (0.25 m × 0.25 m) in which root samples of B. lunaria and the associated rhizospheric soil were collected (Supporting Information Fig. S1). A distance of about 0.75 m between the different sampled quadrats was respected as far as this was possible. Overall, this corresponded to 18 root samples and their associated rhizospheric soil for each population, resulting in a total of 72 samples from the four populations ( Fig. 1A-E).
A total of 72 rhizospheric soil samples were analyzed in order to measure: residual moisture (Hr %), total phosphorus (Ptot), pH, organic carbon (Corg) and total nitrogen (Ntot) contents (Supporting Information Table S1). The corresponding 72 root samples were crushed using tungsten beads and DNA extractions were carried out using the DNeasy Plant Mini Kit (Qiagen) following the manufacturer instructions (Quick-StartProtocol). The primers AMV4.5NF and AMDGR (Sato et al., 2005) were used to amplify a fragment of the 18S rRNA gene including the variable domain V4 and the 72 DNA libraries were sequenced on an Illumina 2 × 250 MiSeq platform by Fasteris SA (Plan-les-Ouates, Switzerland).
Amplicon sequence processing was performed using mothur software version 1.39.5 (Schloss et al., 2009) and by mainly following the Schloss standard operating procedure for MiSeq Illumina data (Kozich et al., 2013). To account for differences in sampling efforts, 16 822 sequences were randomly subsampled from each sample. Sequences were clustered into 97% sequence similarity OTUs using the OptiClust algorithm (Westcott and Schloss, 2017). A representative sequence of each OTU was classified using the naive Bayesian classifier with the MaarjAM database (Opik et al., 2010) and an 80% confidence threshold. Then, the OTU table was curated with the LULU algorithm to remove erroneous rare OTUs based on both sequence similarity thresholds and withinsample patterns of co-occurrence (Frøslev et al., 2017). Finally, singletons were excluded. The raw sequence data have been deposited in the NCBI Sequence Read Archive under BioProject PRJNA547499.
Rarefaction curves, Good's coverage estimator, analysis of molecular variance (AMOVA) (Schloss, 2008), ANOSIM (Clarke, 1993) and Dirichlet multinomial mixture models for unsupervised clustering (Holmes et al., 2012) were computed with mothur version 1.39.5 (Schloss et al., 2009). All the other analyses and data visualization were computed using R software version 3.4.4. (https:// www.r-project.org/). Alpha and beta diversity were visualized with the phyloseq package (McMurdie and Holmes, 2013). To quantify the contribution of spatial distance, soil properties, and floristic composition on the total variation of the AMF OTU abundances, variation partitioning (Borcard et al., 1992) was performed using the vegan R package (Oksanen et al., 2012). Regarding the spatial matrix, a set of spatial predictors from the geographical coordinates was computed by principal coordinates of neighbour matrices (PCNM). Subsequently, forward selection (Blanchet et al., 2008) of PCNM variables based on 10 5 permutation was applied with the packfor R package in order to retain only the significant variables. To identify AMF OTU signatures of the different community types, we computed compositional balances using the selbal algorithm implemented in the eponymous R package (Rivera-Pinto et al., 2018). To evaluate the cooccurrence probability between AMF OTUs, the Veech probabilistic model of species co-occurrence with a hypergeometric distribution (Veech, 2013) was applied using the cooccur R package (Griffith et al., 2016). A cohesion analysis was used to compute positive and negative cohesion metrics (Herren and McMahon, 2017) and also to identify the most connected OTUs, which have been shown to be potential keystone taxa of a community (Herren and McMahon, 2018). To determine which abiotic or biotic variables significantly explained variations in AMF community composition, a forward selection (Blanchet et al., 2008) was performed using the packfor R package. Then, the impact of abiotic and biotic variables on community structure was investigated with generalized linear models (GLMs) for multivariate abundance data using the mvabund R package (Wang et al., 2012).
Full details of all experimental procedures are given in the Supporting Information.

Supporting Information
Additional Supporting Information may be found in the online version of this article at the publisher's web-site: Appendix S1: Supplementary Information about the experimental procedures and related references. Figure S1 A. Schematic view of a sampling plot showing the dimensions and spatial organization of its three quadrats. B. Photography of a quadrat, in each of them the floristic composition was assessed, and B. lunaria root samples and associated rhizospheric soil were sampled. Figure S2 A. Rarefaction curves of the 71 AMF samples. B. Sample-based OTUs accumulation curve. Figure S3 Alpha diversity metrics (observed richness, Shannon diversity and the inverse of Simpson index) in the four transects. Figure S4 Linear relationship between the richness of the neighbouring plants surrounding B. lunaria sporophytes and the soil pH. The solid blue line shows the fit of the linear model and the shaded area corresponds to the point wise 95% confidence interval on the fitted values. R corresponds to the Pearson's product-moment correlation. Figure S5 AMF OTU signature of the different community types based on compositional balances computed with the selbal algorithm. The two groups of OTUs that form the global balance are specified at the top of the plot. The box plot represents the distribution of the balance scores for community type A and B samples. The right part of the figure contains the receiver operating characteristic (ROC) curve with its area under the curve (AUC) value (1) and the density curve for each community type. Figure S6 Values of the four significant variables impacting the AMF community composition, namely pH, humus cover (Landolt's indicator H), positive cohesion and negative cohesion, for the community type A and B. Statistical significance between the two community types was tested by the Wilcoxon signed-rank test. Table S1 Edaphic and floristic parameters measured on each sampling site.