Differences in the fungal communities nursed by two genetic groups of the alpine cushion plant, Silene acaulis

Abstract Foundation plants shape the composition of local biotic communities and abiotic environments, but the impact of a plant's intraspecific variations on these processes is poorly understood. We examined these links in the alpine cushion moss campion (Silene acaulis) on two neighboring mountain ranges in the French Alps. Genotyping of cushion plants revealed two genetic clusters matching known subspecies. The exscapa subspecies was found on both limestone and granite, while the longiscapa one was only found on limestone. Even on similar limestone bedrock, cushion soils from the two S. acaulis subspecies deeply differed in their impact on soil abiotic conditions. They further strikingly differed from each other and from the surrounding bare soils in fungal community composition. Plant genotype variations accounted for a large part of the fungal composition variability in cushion soils, even when considering geography or soil chemistry, and particularly for the dominant molecular operational taxonomic units (MOTUs). Both saprophytic and biotrophic fungal taxa were related to the MOTUs recurrently associated with a single plant genetic cluster. Moreover, the putative phytopathogens were abundant, and within the same genus (Cladosporium) or species (Pyrenopeziza brassicae), MOTUs showing specificity for each plant subspecies were found. Our study highlights the combined influences of bedrock and plant genotype on fungal recruitment into cushion soils and suggests the coexistence of two mechanisms, an indirect selection resulting from the colonization of an engineered soil by free‐living saprobes and a direct selection resulting from direct plant–fungi interactions.

Locally, the available pool of fungi then undergoes selection by soil-or plant-related parameters, including aboveground plant species composition. Fungal community composition is correlated with soil moisture (Hawkes et al., 2011) or pH (Rousk et al., 2010), but also with plant legacies, such as soil C/N ratio (Prober et al., 2015) or soil organic matter content . The soil fungal communities are further linked to the present plant communities, and fungal beta diversity (compositional dissimilarity between sites, see Anderson et al., 2011) is correlated with plant beta diversity (Geremia et al., 2016;Pellissier et al., 2014;Prober et al., 2015;Tedersoo et al., 2014). Also, different plant species give rise to diverging fungal communities from the same starting soil in controlled experiments, likely due to differences in rhizodeposition (Mouhamadou et al., 2013). Plant intraspecific genetic variation can also lead to divergent composition of dependent communities (Bangert & Whitham, 2007;Bangert et al., 2006). For endophytic fungi, this may be mediated directly by the differential production of defense compounds as shown in maize (Saunders & Kohn, 2009) or in Populus (Bailey et al., 2005;Lamit et al., 2014). However, little is known about the influence of plant genetic variation on fungal community composition in soils, which is challenging to disentangle from confounding factors in natural environments. Here, we report on a field survey providing insights into this influence.
Alpine cushion plants offer an interesting system to study fungal community assembly and to assess the role of intraspecific plant genetic variation. These foundation plants locally modify their physicochemical environment (Badano, Jones, Cavieres, & Wright, 2006) and the soil biotic community (Roy et al., 2013), both defining ecological engineering effects (see Jones, Lawton, & Shachak, 1994). Indeed, once a plant is established in cracks on bare mineral soil, the accumulating litter material below the dense green tissues allows the constitution of a de novo organic soil.
Through the buffering of microclimatic conditions (Molenda, Reid, & Lortie, 2012) and the improvement of soil organic matter and nutrients status (Roy et al., 2013), cushion plants act as nurse species for microbiota (Roy et al., 2013), soil arthropods (Maillet, Lemaitre, Chikhi, Lavenier, & Peterlongo, 2012;Molenda et al., 2012), and other plants (Butterfield et al., 2013). To our knowledge, however, the extent of the abiotic and biotic soil changes induced by various cushion plant species or subspecies has not been explicitly studied.
Silene acaulis (L.) Jacq is a common circumboreal and alpine cushion plant with slow growth and cushions that can often be older than 300 years (Morris & Doak, 1998). Its xenogamic reproduction makes that two distinct plant individuals always represent distinct genotypes. Although the evolutionary history of Silene acaulis remains to be fully unraveled, two main subspecies have been proposed in the Alps, S. acaulis exscapa and S. acaulis longiscapa, exhibiting dense-(shorter internodes) and loose-cushion morphologies, respectively (Aeschimann, Lauber, Moser, & Theurillat, 2004;Bock, 1983). They occur on adjacent territories, probably with genetic barriers, longiscapa being restricted to calcareous bedrock while exscapa is mainly, but not only, found on siliceous bedrock (Sébastien Ibanez, S. I., unpublished obs.). Subspecies determination can, however, be difficult since intermediate cushion architectures may sometimes be found ( Figure 1a).
We previously reported that the presence of a cushion affects soil chemical parameters on two bedrocks and that this shift in abiotic properties comes along with a shift in the fungal community; in particular, the fungal turnover from bare soil to cushion soil increased with environmental stress (high elevation and very low pH; Roy et al., 2013). In the present study, we first established that the two plant subspecies corresponded to distinct genotypic clusters, which then allowed us to assess the correlation between plant genetic distances and fungal beta diversity. To disentangle this genetic effect from those of the local environment (geology, elevation), and also to measure abiotic soil engineering by the plant, we sampled soils both beneath cushions and in neighboring bare soils. We finally examined the fungal functional guilds across the sampling design to address the distribution, filtering, and recruitment mechanisms.

| Plant and soil sampling
The sampling area was located in the French Hautes-Alpes along replicated ecotones from alpine to subnival environments and encompassed two neighboring mountain massifs: one mostly calcareous (Cerces) and one siliceous (Combeynot; see Roy et al., 2013 for details Silene acaulis is a gynodioecious species, but plant gender was not recorded since some individuals were not flowering. Vegetative traits are not distinguishable between plant genders, and pollinators visit both plant genders, thus disseminating the same fungal spores; therefore, we assume that gender will be only a weak source of fungal community variation compared to plant genotype or envi-

| Soil physicochemical analyses
Soil chemical parameters were quantified following the methods described in Clément et al. (2012). Briefly, soil samples were sieved at 2 mm and used to estimate soil pH in water, gravimetric soil water content (Robertson, Coleman, Bledsoe, & Sollins, 1999)

| Plant genotyping
Plants genotypes were fingerprinted using AFLP (Meudt & Clarke, 2007) by the procedure of Vos et al. (1995) with minor modifications. Total DNA was extracted from leaf tissues using the DNeasy 96 Plant Kit (Qiagen, Hilden, Germany) following the manufacturer's instructions, and plant DNA was digested for 2 hr at 37°C in a 20-µl mix using 2 U of Mse1 and 5 U of EcoRI (New England F I G U R E 2 Plant genetic groups and their geographic distribution. (a) Proportion of AFLP alleles on plant individuals along the sampling gradients. Colors and letters refer to genotypic clusters: "x," "l," and "a" for exscapa, longiscapa, and ambiguous plants, respectively, as assigned by STRUCTURE (K = 2). The gradient of origin, elevation, and morphology of the sampled plants (triangles for dense and circles for loose cushions) are indicated; letters "C" (calcareous) and "G" (granitic) indicate the bedrock type (except "Q", quartzitic). (b) NMDS plot of the Jaccard distances between AFLP fingerprints of cushion plants. Colors denote genotypic clusters, and shapes denote the bedrock of origin AFLP profiles were analyzed using Peak Scanner software (Applied Biosystems); the profiles, for which the automated process failed to detect or attribute peak size, were processed manually. The binning and scoring of the profiles were performed using the RawGeno package (Arrigo, Tuszynski, Ehrich, Gerdes, & Alvarez, 2009) in R software (R_Development_Core_Team, 2011). We kept individuals within the 5%-95% quantile confidence interval around the peak number mean and kept the peaks that corresponded to fragments longer than 100 bp to limit size homoplasy (Arrigo et al., 2009;Pompanon et al., 2005). Bins that were not reproducible at least once within a group of replicates (four groups of replicates for each pair of primers and three to four replicates per group) were systematically removed. We calculated the "bin content information" criterion (Arrigo et al., 2009;Pompanon et al., 2005) for various filtering strategies and found that our conservative strategy retained the most information. A total of 93 individual plants were reliably scored for the presence or absence of 345 reliable and polymorphic loci, which were coded in a binary matrix.
The plant population genetic structure was inferred from AFLP data using STRUCTURE 2.3.2 (Pritchard, Stephens, & Donnelly, 2000), implementing an admixture model without a priori knowledge of the geographic provenance of individuals. The likelihood of the number of genotypic clusters (K) was estimated for K ranging from 1 to 20 with five runs for each K value with 50,000 burn-in periods followed by 50,000-step Markov chain Monte Carlo (Supporting Information Figure S1). Pairwise genetic distances were assessed with the Jaccard index.

| Fungal ITS1 sequencing, clustering, and assignment
The ITS1 DNA region was amplified by PCR using DNA extracted from each soil core as previously described (Roy et al., 2013). Fungal DNA was amplified using the ITS5 (5′-GGAAGTAAAAGTCGTAACAAGG-3′) and the 5. Bioinformatic processing of sequenced amplicons used the OBITools suite (Boyer et al., 2016). Reads were quality-filtered, and ITS amplicons were reconstructed by merging the 3′ and 5′ pair-end reads using solexapairendnull. We used obigrep to discard merged reads with an alignment score ≤20 (corresponding to a perfect match on ≤5 bases for reads with the highest Phred score) together with those containing ambiguous nucleotides, or errors in the primer sequence or the sample tag. The remaining merged reads were assigned to their sample of origin by their terminal tags using ngsfilter, which trims primers and tags away from the sequence, and identical reads were counted within each sample using obiuniq. Finally, singletons were removed following Coissac, Riaz, and Puillandre (2012) (Pearson, 2000) implemented in the ecotag program, and only environmental sequences identical at a minimum of 80% to a reference sequence were assigned. Possible lifestyles were deduced from the MOTU taxonomy using FUNGuild (Nguyen et al., 2016

| Statistical analyses
MOTU read counts per sample were transformed to relative abundance. To assess the compositional changes between soil fungal communities, we summarized pairwise Bray-Curtis dissimilarities by non-metric multidimensional scaling (NMDS), using vegan (Oksanen et al., 2013). Environmental fitting was used to visualize the correla- To test for significant differences between taxonomic abundances beneath and outside cushions, we fitted two linear mixed models using lme4 (Bates, Maechler, Bolker, & Walker, 2015): a null model and a model including, as a fixed factor, the location beneath cushion or in bare soil. The two models included the sampling site as a random factor. The models were then compared in an ANOVA test, and p-values were corrected for multiple testing using the false discovery rate (Benjamini & Hochberg, 1995). Fungal diversity was quantified using the inverted Simpson index.
To identify MOTUs that were recurrently associated with plants, we grouped samples into six habitats combining the granite/limestone bedrock type, the local plant genetic cluster (see Results), and the beneath/outside cushion location. For each MOTU, we counted occurrences over a 0.1% frequency threshold in each habitat and performed a chi-square test of independence. The p-values were obtained through 10,000 permutations of the original matrix and corrected for multiple testing. We kept those MOTUs that were significantly and positively associated with at least one cushion habitat (p-value < 0.05; Pearson residual > 2) but with no bare soil habitat.

| Genetic structure of Silene acaulis populations
Two genetic groups of plants were revealed by population structure inference (Supporting Information Figure S1) and the ordination of  Table   S1; Figure 3a). By contrast, the exscapa genetic cluster was found on both calcareous and siliceous (granite and quartzite) bedrocks.
Although exscapa culminated at higher elevations than longiscapa, six populations of exscapa and three populations of longiscapa occurred within the same elevation range (2,500-2,700 m a.s.l.; Supporting Information Table S1, Figure 2a).

| Differences in soil abiotic engineering by the two plant genetic clusters
The PCA ordination of soil chemical properties revealed two axes totalizing 79% of the total variability (Figure 3a). Axis 1 was related to nutrient and water contents, and Axis 2 was related to soil pH and nitrate content. Cushions from exscapa genetic cluster induced a strong shift of soil chemical properties along Axis 1 from bare to cushion soils. This shift was smaller for longiscapa genetic cluster (Figure 3a,b). Supporting Information The contents in nitrogen, carbon, water, and ammonium all increased on both granitic and calcareous bedrocks beneath exscapa plants, and soil became less acidic on granite (Supporting Information Figure S3).
For soil beneath longiscapa cushions, the increase in nutrient contents was weaker, although significant for the carbon and ammonium contents (Supporting Information Figure S3).

| Differences in the engineering of fungal communities by the two plant genetic clusters
NMDS ordination of fungal community dissimilarities revealed a sharp difference in community composition on the two main geological substrates, which correlated with pH and nitrate contents (Figure 3c). On top of this bedrock-related divergence, fungal communities of exscapa cushions from both mountain ranges segregated away from all other habitats, which correlated with increased soil contents in carbon, total nitrogen, ammonium, and water. For longiscapa, the shift in community composition differed from the one observed for exscapa on limestone (Supporting Information Figure S2). From bare to cushion soils, fungal dissimilarities were strong for both subspecies (Figure 3d). The shifts in community composition included a consistent decrease in fungal diversity beneath the cushions compared to bare soils (Figure 4a). This diversity reduction came along with a strong MOTU turnover, which was modulated differently for the two genetic clusters and the two bedrock types (Figure 4b). F I G U R E 4 Fungal diversity and turnover across habitats. Samples were grouped into six habitats combining the plant genotype ("x" for exscapa and "l" for longiscapa), the bedrock type ("G" for granitic and "C" for calcareous), and the location ("b" for bare soil and "c" for cushion soil). From left to right, n = 40, 43, 12, 11, 22, and 20

| Links between plant genetic distances and fungal beta diversity
The variations in the composition of fungal communities among cushion soils were best explained by the combined variations in plant genotype, soil chemical composition beneath and outside the cushions, and geographic position ( Figure 5), which accounted together for 40% of the total variability. This variability was 2% more than when the geographic position was dismissed at the profit of elevation (Supporting Information Figure S4). The largest pure effect (10%) was that of the genetic distances between plants ( Figure   S4). This discrepancy suggests that each plant cluster selectively recruited a few MOTUs that rarely occur in bare soils but became highly abundant beneath cushions.

| Taxonomic identification of fungal sequences associated with the cushion habitats
Ascomycota largely dominated all fungal communities but decreased beneath cushions vs. bare soils, while Basidiomycota significantly increased ( Figure 6a). The latter included a trend for Sebacinaceae to increase (Figure 6b; 2.7% beneath vs. 0.6% outside). In both phyla, the fraction of fungi unassigned at the family level was high in bare soils and decreased strongly in cushion soils (Figure 6b). Remarkably, the abundance of Cladosporiaceae and Dermateaceae (Ascomycota) increased largely and significantly beneath the cushions, cumulating in 19.4% of the reads in cushion soils vs. 6.6% in bare soils. The nearly doubled frequency of Dermateaceae was entirely caused by MOTUs assigned to Pyrenopeziza brassicae (Supporting Information Table S2). We

| The two plant subspecies differentially engineer their soil environment
As previously reported (Roy et al., 2013), cushions engineer a soil enriched in carbon and ammonium, two key soil nutrients that are mainly derived from plant litter; here, we found that this edaphic improvement is clearly stronger for the exscapa cluster. These differences in abiotic soil engineering are likely to derive from differences in cushion architecture. A causal link between the nurse plant stature and its engineering effects has been demonstrated in a dominant coastal shrub (Crutsinger et al., 2014). The exscapa subspecies displays a denser architecture than the longiscapa one, limiting litter dispersion and certainly providing a higher capacity to limit water loss. Differences in cushion architecture between subspecies most likely have a heritable basis because they correspond to distinct genotypic groups and remain observable over a broad elevation range and on distinct bedrocks. Our sampling, however, does not encompass the whole ecological range of each subspecies and therefore does not explore their entire vegetative trait plasticity (Bonanomi et al., 2016), making it difficult to tease apart the respective influences of phenotypic plasticity and genetic differentiation.
Fungal communities strongly varied from bare soil to cushion soils, pointing to the selective filter exerted by cushions on fungi.
This filtering is highlighted first by a decrease in fungal α-diversity and by strong shifts in fungal composition. These changes were again stronger for the exscapa cluster and seem to be related to the plant subspecies rather than to the environment because they occurred with the same amplitude on the two contrasted bedrocks where exscapa plants were found. S. acaulis cushions improve edaphic properties but can also buffer soil moisture and alter the temperature regimes (Bonanomi et al., 2016;Molenda et al., 2012). This new microenvironment may indirectly favor the growth of new, especially saprophytic, fungal species. There is also a clear influence of the bedrock, as illustrated by the exscapa cushions, which remained largely different in fungal composition on granite or on limestone.
Bedrock, as seen by bare soil chemistry, contributed little as a pure part to the variability of fungal composition but much more to the fractions shared with a geographic position, owing to the spatial separation of limestone and granite sites. Ascomycota are in the red-to-yellow color range, Mucoromycota appear in blue, and Basidiomycota are in the green range. Families representing <1% of the total reads were pooled. MOTUs without assignment at the family level were pooled by phylum ("uf" prefix); "NAs" indicate MOTUs of unknown phylum. Families whose frequencies vary from bare soils to cushion soils bear a suffix ("_^", increase; "_v", decrease); the significance of the difference is indicated on sectors ("."; "*"; "**"; "***" for p < 0.1, 0.05, 0.01, and 0.001, respectively) a single host could respond to different degrees of soil engineering from the plant litter, be it in the content in organic matter, nutrient (Hanson, Allison, Bradford, Wallenstein, & Treseder, 2008;Treseder & Lennonb, 2015), or water (Crowther et al., 2014). For instance, Geomyces (Ascomycota) possibly bloomed beneath exscapa cushions because they occupy the highest elevation on limestone since this genus includes psychrophilic species able to profit from a variety of organic substrates (Hayes, 2012). Saprobes could also be host-specific epi-or endophytes turning into decomposers upon plant leaf senescence (Osono, 2002;Voriskova & Baldrian, 2013;Zhou & Hyde, 2001 (Zhang & Yao, 2015). This strongly suggests that direct interactions between the fungus and living plant tissues account for their recruitment in alpine cushion soils and is consistent with the strong costructure between phyllosphere fungal composition and intraspecific population genetic structure in beeches (Cordier, Robin, Capdevielle, Desprez-Loustau, & Vacher, 2012) as well as with the phylogenetic signal in the association of tropical trees with fungal pathogens (Gilbert & Webb, 2007). These narrow host ranges might result from a local co-evolutionary arm race (Chappell & Rausher, 2016) (Defossez et al., 2011). The cost for the plant to host endophytic fungi has been shown to be weak (Kia et al., 2016). Fungi associated with cushion plants may, for instance, offset their fitness costs because they reduce herbivore or virulent phytopathogen attacks or increase plant tolerance to abiotic stress (Rodriguez, Redman, & Henson, 2004). In a situation where nutrients are scarce, there also might be an evolutionary trade-off between the cost of hosting a biotrophic fungus and the accelerated nutrient cycling by the same fungus in the litter. Such a trade-off has been predicted for insect grazers (de Mazancourt & Loreau, 2000).

| Ecological implications
Our results suggest that the plant evolutionary dynamic of cushion plants strongly impacts the fungal composition of alpine ecosystems, possibly through diverging morphological traits. This extended phenotype sensu Whitham (Whitham et al., 2003) may have implications for plant succession. Cushion plants, including S. acaulis, help secondary plant species to settle in the landscape (Badano et al., 2006;Bonanomi et al., 2016). This settling may occur through reduction in physical stress (Bertness & Callaway, 1994), but the presence of antagonist fungi may also be important in this process. The Janzen-