Ecophylogeny of the endospheric root fungal microbiome of co-occurring Agrostis stolonifera

Plant harvesting and root sampling

Turf samples (L44 × W34 × D23 cm) were collected from a peatland in the Parc Naturel Régional des Marais du Cotentin et du Bessin, Normandy, France (49.284656–1.393090). The sampled turf vegetation was dominated by Agrostis stolonifera (∼94% coverage), a generalist plant often found in disturbed habitats, wetland margins and fields. The other co-occurring plants were Potentilla anserina, Mentha spp. and Hydrocotyle vulgaris. Turfs were placed in a growth chamber for four months at 15 °C with 16 h light per day, and 80% humidity to avoid stochastic environmental effects. Nineteen co-occurring Agrostis stolonifera plants were manually harvested from a single turf. During sampling, the turf was divided into 3 disconnected blocks but as there was no block effect in the observed variances, all data were pooled for the analyses. The sampled plants were of similar age with comparable root systems. No difference in the phenotypic traits of sampled plants was noticed. The roots were washed with tap water to remove soil and the surface roots were sterilized by washing three times with 0.1% Triton X-100 and rinsed five times with sterile distilled water (Vandenkoornhuyse et al., 2007). The roots were then transferred into RNAlater (Sigma-Aldrich) and stored at −20 °C.

Nucleic acid extraction and 18S rRNA amplicon sequencing

These sampled roots were crushed into powder with a mortar and pestle in liquid nitrogen. RNA and DNA were co-extracted from the nineteen samples using a modified RNeasy Plant Mini kit (Qiagen, Hilden, Germany) protocol. The on-column DNase digestion was skipped to keep the DNA in the total extracts. For RNA extraction half of the total nucleic extracts were treated with RQ1 DNAse (Promega, Madison, WI, USA) and the complete elimination of DNA was confirmed by 18S rRNA targeting PCR. RNA and DNA quality was checked using the RNA 6000 Pico kit or the High sensitivity DNA kit, respectively, on a 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). PCRs were performed using the primer pair SSU0817 (5′-TTAGCATGGAATAATRRAATAGGA-3′) -NS22B (5′-AATTAAGCAGACAAATCACT-3′) targeting a region of about 530 bp of the 18S rRNA gene that includes the variable regions V4 (partial) and V5 (Borneman & Hartin, 2000). These particular primers were chosen among a variety of candidate primers after an in silico analysis, using Primer Search (Rice, Longden & Bleasby, 2000). The chosen set could amplify 94% of the available fungal sequences in silico, with the exception of Microsporidia, and only 1.3% of Viridiplantae in the Silva database 115 (Quast et al., 2013). PCR amplifications were performed with fusion primers containing sequencing adaptors, multiplex identifiers (MID) and PCR primers (Table S1). The DNA was amplified by performing direct PCRs on the total nucleic acid extracts using illustra puReTaq Ready-To-Go PCR Beads (GE Healthcare). Two microliters of DNA template (at ∼1 ng µL⁻¹) were used in a final volume of 25 µL with 0.2 µM concentrations of each primer. The PCR cycling protocol consisted of 35 cycles of denaturation (95 °C for 30 s), annealing (54 °C for 30 s) and elongation (72 °C for 1 min) with an initial denaturation step (95 °C for 4 min) and a final elongation step (72 °C for 7 min). Two independent PCRs representing technical replicates were performed for each sample. The RNA was amplified by performing RT-PCRs using the Titan One Tube RT-PCR kit (Roche Molecular Systems). The reaction was carried out in a final volume of 50 µL with 0.2 µM of each primer. Reverse transcription (42 °C for 30 min) was followed by PCR amplification with an initial denaturation step at 94 °C for 3 min, 38 cycles of 30 s at 94 °C, 30 s at 54 °C and 45 s (+5 s/cycle from the 11th to the 38th cycle) at 68 °C and a final elongation at 68 °C for 7 min. Two independent RT-PCRs were performed for each sample. The quality of the PCR products was checked on High Sensitivity DNA chips (Agilent, Santa Clara, CA, USA) and purified amplicons were quantified by spectrofluorometry using the Quant-iT PicoGreen dsDNA Assay kit (Invitrogen). The libraries were pooled in equimolar amounts and purified using the AMPure XP system (Beckman-Coulter, Brea, CA, USA). Any traces of concatemerized primers were removed by subjecting the libraries to microelectrophoresis on a Caliper Labchip XT instrument (Perkin Elmer, Waltham, MA, USA). The libraries were amplified by emPCR using the GS FLX Titanium MV emPCR Kit with Lib-L chemistry and sequenced on a GS FLX+ sequencer (Roche/454, Branford, CT, USA) following the manufacturer’s instructions.

Sequence analysis using a dedicated automated pipeline

Quality trimming and filtering of amplicons, OTU identification, and taxonomic assignments were carried out with a combination of publicly available sequencing data analysis tools (Cutadapt, Mothur, Dnaclust) and in-house python scripts within a Galaxy instance at the Genouest platform (http://www.genouest.org/), as described elsewhere (Ciobanu et al., 2014; Nunes et al., 2015; Ben Maamar et al., 2015). Briefly, sequences shorter than 200 bp in length, with homopolymers longer than 8 bp or with ambiguous nucleotides, were removed from the dataset. Sequences containing errors in the MID or primer sequences were discarded. Chimeric sequences were eliminated using the chimera.uchime command of the Mothur tool suite. It is well known that both PCR and pyrosequencing can induce erroneous sequences (Shakya et al., 2013b) leading to poor diversity estimates. To improve sequence quality, two independent PCR reactions were performed for each sample and replicates were sequenced. Only sequences displaying 100% identity among these technical replicates were retained for subsequent analysis. Sequences were grouped into OTUs with a sequence identity threshold of 97%. Consequently one OTU was defined by at least two identical sequences originating from technical replicates. The taxonomic affiliations of the sequences and OTUs were determined by comparison with the SILVA database 115 (Quast et al., 2013). Two OTUs, representing 1,132 sequences, assigned to the Chloroplastida, were removed. The study accession number in the European Nucleotide Archive is PRJEB12655.

Statistical and diversity analyses

The samples were normalised to 1,288 sequences for analyses of alpha and beta diversity whereas the full dataset was used to analyze the core microbiome in co-occurring plants within a single population. All statistical analyses were performed in R (R Core Team, 2013) using the VEGAN package (Oksanen et al., 2011). To check the sequencing depth, rarefaction curves were computed using the rarefaction function. Alpha diversity and richness were estimated for each sample using Hill diversity numbers (Hill, 1973) and the Chao 1 index (Chao, 1984). Hill diversity numbers allow accurate comparison of species/OTU diversity across samples. The significance of the Hill diversity numbers depends on the value of the q parameter in the Hill formula. This parameter allows species to be weighted more or less equally. For q = 0, OTUs or species are weighted equally and the Hill diversity is equal to the OTU richness while for q = 1, the Hill diversity corresponds to the Shannon diversity and for q = 2, the Hill diversity is equal to the Simpson index. One community can be considered more diverse when all of its Hill diversity numbers are higher than those of the other communities. The Chao 1 index estimates the unseen diversity by taking rare OTUs into account. The alpha diversities and taxonomic compositions of the DNA and RNA fractions were compared using paired Student’s t-tests or Wilcoxon rank-sum tests only when the alpha diversity values did not follow a normal distribution or were heteroscedastic.

Beta diversity was studied by using non-metric multidimensional scaling (NMDS) with the Bray–Curtis dissimilarity matrix. Data were transformed using the square root and the Wisconsin double standardization implemented in the metaMDS function. Procrustes analysis was conducted on the NMDS scores to assess the concordance between the communities in the DNA and RNA fractions. Significance of the concordance was tested by permutation (10,000) using the protest function (Peres-Neto & Jackson, 2001).

We defined the core microbiome as the proportion of OTUs shared by the studied co-occurring A. stolonifera plants. To determine whether the number of shared OTUs was dependent on the sampling effort, we performed random re-sampling and increased the sampling size from one to the total number of samples for each fraction. We defined the ‘DNA core’ and the ‘RNA core’ as the proportion of OTUs present in the DNA and the RNA fraction respectively of all samples.

Phylogenetic tree construction

Representative sequences of OTUs were aligned using SINA aligner v1.2.11 (Pruesse, Peplies & Glockner, 2012), imported from the non-redundant SILVA SSU Ref ARB database (release 115). Alignments of the reference sequences and of representative OTU sequences were exported from ARB (Ludwig et al., 2004). Gaps and ambiguously aligned positions were excluded. One alignment was obtained for each main fungal phylum. The model of sequence evolution that best fitted the aligned sequences was selected using jModelTest v2.1.4 (Darriba et al., 2012). Phylogenetic trees were constructed by maximum likelihood using TREEFINDER and by bayesian inferences using MrBayes v.3.2.2 (Ronquist et al., 2012). Maximum likelihood bootstrap values were calculated from 1,000 replicates and Bayesian posterior probabilities were calculated using 100,000 generations sampled every 100 generations.

Phylogenetic structure

The community phylogenetic structure was studied by applying the Picante package (Kembel et al., 2010) to the global phylogenetic tree including all the fungal OTUs. Comparison of the phylogenetic relatedness of the OTUs in the DNA and RNA fractions of each sample with the local pool of OTUs (OTUs found in all samples) was based on two indices: the mean pairwise distance (MDP) and the mean nearest taxon distance (MNTD) described by Webb et al. (2002). The phylogenetic distances measured in the observed communities were compared with those in the null communities generated by randomization (1,000 permutations) by calculating standardized effect sizes (SES). Tip labels were randomly shuffled across the tips of the phylogeny for the null model. In this model, the community data matrix is not randomized in order to fix most of the patterns (species richness, observed occupancy rates and OTU abundance) and to allow the pattern of interest to vary (phylogenetic distance) (Swenson, 2014). However, type I error rates may be inflated if a phylogenetic signal in OTU abundance occurs (Hardy, 2008). As no phylogenetic signals in OTU abundance were found using Blomberg’s K statistic (Blomberg, Garland & Ives, 2003) the use of this null model was validated. Negative SES values and p < 0.05 indicated phylogenetic clustering and smaller than expected phylogenetic distances among co-occurring OTUs. Two analyses were performed to create the null communities, one using all the OTUs and another one using the OTUs obtained separately from the DNA and the RNA fractions. Both methods gave similar results and only results from the former analysis are therefore presented.

Phylogenetic signal

Here we measured the ‘mean expression ratio’ of each OTU, in order to access the metabolic status of the microorganisms. This metric was obtained for each OTU by dividing the mean relative abundance in the RNA fractions (mean RA_RNA) by the sum of the mean relative abundance in the DNA (mean RA_DNA) and RNA fractions among all samples. ${\displaystyle meanexpressionratio=\frac{meanR{A}_{RNA}}{meanR{A}_{RNA}+meanR{A}_{DNA}}}$ with ${\displaystyle meanRA=\frac{\sum _{i=1}^{n}RAOTUi}{n}}$ n = number of plants showing the OTUi within their DNA and/or RNA fraction The standard errors were then calculated for each OTU of the 19 samples. A zero value for this ratio would indicate an absence of the OTU from the RNA fractions while a value of one would indicate that the OTU was not detected in the DNA fractions. To determine whether the ‘mean expression ratio’ across the phylogeny was not random, the phylogenetic signal was calculated using the Kse statistic (Blomberg, Garland & Ives, 2003; Ives, Midford & Garland, 2007) for the global phylogenetic tree and for phylogenies built for each main fungal phylum (i.e., Glomeromycota, Basidiomycota, and Ascomycota). This index integrates measurement errors of the trait by calculating the standard errors of the means for each OTU. The Kse statistic compares the observed phylogenetic signal in the ‘mean expression ratio’ to the signal obtained under a Brownian motion model of trait evolution. In this model, traits evolve as a random walk along the branch length of a phylogeny. The statistical significance of Kse was evaluated by comparing the observed patterns of variance of the independent contrasts of the trait to a null model. Taxa labels were shuffled across the tips of the phylogeny for the null model. A ‘mean expression ratio’ with p < 0.05 indicated a non-random phylogenetic signal.

DNA- and RNA-based analyses provided similar alpha and beta diversity profiles but different OTU richness levels

We examined alpha and beta diversity measures of the fungal community colonizing the roots of 19 Agrostis stolonifera plants to see whether the nature of the co-extracted nucleic acids encapsulated different information. The mean Chao 1 index, that estimates ‘true OTU diversity’ by taking rare OTU into account, did not differ significantly (p = 0.11) between the DNA and RNA fractions (S.Chao1_DNA = 215, S.Chao1_RNA = 232). On average, the communities in the RNA fractions were not considered more diverse than those in the DNA fractions as not all of the Hill diversity numbers were significantly different (Fig. 1A). Nevertheless, sample-to-sample variations did exist (Fig. 1A). Interestingly, OTU richness, i.e., a special case of Hill diversity numbers when q = 0, was significantly higher in the RNA fractions than in the DNA fractions, with on average 149 OTUs in the RNA fractions and 133 OTUs in the DNA fractions (p < 0.01) (Fig. 1A). The sequencing depth was similar between both fractions (Fig. S1).The procrustes analyses indicated a significant similarity (p < 0.001) between the distances obtained from the DNA and RNA matrices. Thus, the beta diversity was not structured according to the type of nucleic acid.

Nucleic acid type impacts the perception of the taxonomic composition

The fungal communities in the roots were dominated by Pezizomycotina (Ascomycota) and Agaricomycotina (Basidiomycota) (Fig. S2). Each phylum, except for the Basidiomycota, was dominated by the same OTU in both the DNA and the RNA fractions, representing 17%–76% of the total sequence number (Fig. 1B). The most abundant OTUs in the DNA and RNA fractions were the same except for the Basidiomycota. In this phylum, the two most abundant OTUs in the DNA and RNA fractions were assigned to the Agaricomycetes. The most abundant OTUs were assigned unambiguously to the class level (Sordariomycetes) within the Ascomycota, to the order level (Agaricales) within the Basidiomycota and to the family level (Glomeraceae) within the Glomeromycota (i.e., unknown at lower taxonomic rank). The mean relative sequence abundance assigned to the Basidiomycota was significantly higher in the DNA fraction than in the RNA fraction (p < 0.05, Fig. 1B). Conversely, sequences assigned to the Ascomycota were more abundant in the RNA fraction than in the DNA fraction (p < 0.001, Fig. 1B). Glomeromycota represented 3.9% and 10.1% of the total number of sequences in the DNA and RNA fractions respectively but the differences between these fractions were not significant (Fig. 1B).

At the population level, very few fungal OTUs were shared by all co-occurring plants

Nine OTUs out of a total of 635 OTUs were found in both the DNA and RNA fractions of all of the 19 co-occurring A. stolonifera plants, accounting for 33% of the total number of sequences (Fig. 2). The number of sequences was equally divided into the RNA fraction (54%) and the DNA fraction (46%). Considering each fraction separately, 17 and 16 OTUs in the DNA and RNA fractions respectively were shared by all samples. All OTUs in the ‘DNA core’ were found in the RNA fraction of at least one sample and reciprocally. The number of OTUs shared by all samples in both fractions was negatively linked to the sampling size (Fig. S3). However, this number of shared OTUs stabilized at around 15 samples, indicating that the ‘core microbiome’ size would not be diminished by increased sampling within these co-occurring A. stolonifera (Fig. S3). The phylogenetic analysis (Fig. 2) revealed that the fungal ‘core microbiome’ was dominated by Ascomycota and the closest known sequence of the most abundant OTU was designated Ascomycete sp. (i.e., an unknown OTU) (Fig. 2). Strikingly, most of the OTUs in the fungal potential ‘core microbiome’ were only distantly related to known species (Fig. 2).

The fungal root microbiome is phylogenetically structured

We tested our working hypothesis that the host plant acts as an environmental filter on microbial communities by examining the phylogenetic structure of the fungal root microbiome. For this we measured the standardized effects size of the mean pairwise distance (SES_MPD) and of the mean nearest taxon distance (SES_MNTD) in the DNA and RNA fractions of 19 samples (i.e., 38 fractions). Negative values of SES_MPD indicated that the co-occurring fungal OTUs were more closely related than expected under a null model, and this phylogenetic pattern was significant (p < 0.05) in 20 out of the 38 DNA and RNA fractions (Table 1). This phylogenetic clustering based on the SES_MPD measures was reinforced by the SES_MNTD values (Table 1). The obtained pattern was significant for 36 DNA and RNA fractions indicating that phylogenetic clustering mainly concerned the leaves of the phylogeny (Table 1). The SES_MNTD index calculates whether the closest related OTUs (nearest taxon) tend to co-occur or not in the communities relatively to the null distribution. Significant SES_MNTD values and non-significant SES_MPD values indicated that the OTUs were phylogenetically structured near the tips of the phylogeny and randomly distributed across the tree, i.e., deeper branches contributed less to the pattern.

Phylogenetic signal

The fungal root microbiome displayed phylogenetic clustering within each sample compared to the local pool of OTUs. We investigated whether this phylogenetic relatedness was linked to ecological similarity by searching for a phylogenetic signal, using the ‘mean expression ratio’ as a proxy for a functional trait. Sample-to-sample variations were taken into account in the Kse index. We found that our trait, i.e., the ‘mean expression ratio’, displayed a significant phylogenetic signal (Table 2). This was true for the global phylogeny including all the fungal OTUs and also for the phylogenies built for each fungal phylum (i.e., Glomeromycota (Fig. S4), Basidiomycota (Fig. 3), and Ascomycota (Fig. S5). Thus, the ‘mean expression ratio’ was linked to phylogenetic relatedness. The Kse statistic was less than one for all phylogenies except the Basidiomycota phylogeny (Table 2). Kse values below one indicated that the OTUs were more divergent in their ‘mean expression ratio’ than would be expected under a Brownian motion model of trait evolution whereas a value above one indicated that the values of the ‘mean expression ratio’ were more similar between phylogenetically related OTUs than between OTUs drawn at random.

OTU richness and taxonomic composition of the root fungal microbiome differed between the DNA and RNA fractions

In previous studies, the fungal community composition in roots was essentially analyzed either by cloning of PCR products (e.g., Vandenkoornhuyse et al., 2002a; Vandenkoornhuyse et al., 2002b) or more recently, by mass sequencing of amplicons using targeted DNA (e.g., Opik et al., 2013; Shakya et al., 2013a; Bonito et al., 2014). Others studies based on RNA and DNA extractions have used the RNA/DNA ratio as a proxy to investigate microbial dormancy and to estimate the metabolically active community (Aanderud et al., 2016; Jones & Lennon, 2010). In a similar way we defined the root fungal microbiome of co-occurring A. stolonifera plants by combining the DNA- and RNA-based 18S rRNA amplicon analyses. However, the nucleic acids were co-extracted in order to limit any extraction-related bias. We found higher OTU richness levels in the RNA fractions than in the DNA fractions (Fig. 1A, q = 0), that was consistent with the study by Kuramae et al. (2013) and indicated that RNA-based approaches capture higher numbers of different OTUs in a given sample. At first glance this may seem surprising, but considering that RNA extracts contain roughly 80–90% of rRNA transcripts, the rRNA genes were much more diluted in the DNA fraction and could not be amplified. Thus, for a given quantity of nucleic acids, a higher abundance of targeted 18S rRNA genes was introduced into the RT-PCR mix possibly allowing the capture of more diverse molecules. Although not significantly different, the Chao 1 index exhibited the same trend of decreased richness in the DNA fraction. When other alpha diversity indexes were used, the differences were also not significant (Fig. 1A, q > 0.5). However, significant differences between the DNA and RNA-based analyses were detected when taxonomic composition was taken into account (Fig. 1B). The relative sequence abundance of the Basidiomycota was higher in the DNA fractions than in the RNA fractions and the opposite was true for the Ascomycota (Fig. 1B). It is tempting to conclude that the Basidiomycota was less ‘active’ than the Ascomycota phylum. However, while RNA represents expressed gene levels and is the best proxy for actively growing microorganisms at the nucleic acid level, there is not always a direct link between rRNA abundance, growth and activity in microorganisms (Blazewicz et al., 2013). The absence of a given OTU from the RNA-based analysis, but detected in the DNA, was not related to amplification bias (i.e., only OTUs found in the true replicates were kept). One possible explanation was that the OTUs detected from DNA fraction originated in part from dormant or dead cells.

These results show that it is essential to be aware of the type of nucleic acid being used for comparisons or meta-analyses and highlights the advantages and importance of taking both kinds of nucleic acids into account to gain deeper insights into microbial diversity.

Fungal root microbiome is diverse and highly variable among co-occurring A. stolonifera

In the nineteen A. stolonifera plants, harvested from a natural ecosystem, we found 635 OTUs with a dominance of Pezizomycotina (Ascomycota) followed by Agaricomycotina (Basidiomycota). The Glomeromycota represented only a small proportion of the total number of OTUs (less than 4%) in the DNA fraction, which is comparable to the fraction found in Agave species (Coleman-Derr et al., 2015). However, this proportion was higher in RNA fractions (10%), underlining the interest of RNA based-analyses. A few published studies have focused on the fungal fraction of the root microbiome, mainly in woody plant species, by mass sequencing from DNA extracts (Gottel et al., 2011; Shakya et al., 2013a; Bonito et al., 2014; Coleman-Derr et al., 2015). Although these studies involved different primer pairs and different hosts, the fungal taxonomic composition at the subphylum level was similar to our study. The fungal ‘core microbiome’ at the A. stolonifera population level (i.e., core microbiome sensu Vandenkoornhuyse et al., 2015) accounted for ∼1.5% of the total OTUs (i.e., 9/635). The core microbiome is expected to decrease in complexity whereas the level of ecological organization will increase (i.e., from individual to ecosystem) (Vandenkoornhuyse et al., 2015). Although the level used to define the ‘core microbiome’ was low (i.e., a population of co-occurring plants) and no temporal changes were taken into account, the complexity of the fungal ‘core microbiome’ was lower than expected given the fact that all the A. stolonifera sampled in this study came from the same turf. Increased sampling within this A. stolonifera population would be unlikely to modify the number of shared OTUs (Fig. S3), but might increase the total number of OTUs detected. Moreover, the OTUs shared by all the plants in both their DNA and RNA fractions were mainly abundant OTUs affiliated to Pezizomycotina (Ascomycota), which represented more than 30% of all the sequences.

An increased sequencing depth would lead to an increase of the rare, rather than abundant, OTUs. Unexpectedly, only 1 OTU affiliated to the Glomeromycota forming arbuscular mycorrhiza with A. stolonifera was shared by all individuals, while most of the other OTUs shared by all individuals were unknown Ascomycota forming unknown fungal endophytes (Fig. 2). As the different Ascomycota and Basidiomycota fungal endophytes are acknowledged to be involved in plant resistance to stresses (e.g., Soares et al., 2016; Cosme et al., 2016), our findings raise important new questions about the functions of these unknown endophytes and more widely emphasize the need for a more holistic perception and understanding of the plant holobiont (Vandenkoornhuyse et al., 2015).

A. stolonifera fungal endophytic communities exhibited significant phylogenetic clustering

We found that the fungal communities in plant roots exhibited a significant phylogenetic structure whether the analyses were based on DNA or RNA which confirmed our working hypothesis that the fungal communities in plant roots are not random assemblages (Table 1). At the host-plant scale, the phylogenetic structure of the fungal microbiome showed greater clustering than would be expected under a null model. Thus, the fungal OTUs in an individual plant were more closely related to each other than to the pool detected at the community level. It is not easy to discriminate between the different neutral or ecological processes underlying this pattern because various explanations are feasible (Losos, 2008; Revell, Harmon & Collar, 2008). Indeed, Phylogenetic clustering can be generated by dispersal limitation as well as by inter-specific interactions or environmental filtering (Bell, 2005; Cavender-Bares et al., 2009; Helmus et al., 2007). However, the spatial scale of our study was small in relation to the known dispersal capacities of fungi (Taylor et al., 2006). Dispersal limitation was therefore unlikely or, at least, would not be the main driver of the observed phylogenetic pattern. Ectomycorrhizal fungal communities in plant roots have been shown to be shaped by competitive interactions (Pickles et al., 2012) and host-plant specificity or preference has also been reported (Ishida, Nara & Hogetsu, 2007; Tedersoo et al., 2013). Similarly, assembly of the arbuscular mycorrhizal fungal community in plant roots reflects a degree of host-plant preferences (Vandenkoornhuyse et al., 2002b; Vandenkoornhuyse et al., 2003; Davison et al., 2012). In Populus deltoides roots, the endophytic fungal communities were found to be clearly distinct from the surrounding rhizospheric communities (Gottel et al., 2011; Shakya et al., 2013a). Environmental filtering by the host, which would represent a specific habitat, could thus be an important mechanism leading to phylogenetic clustering. Assuming that phylogenetically close OTUs share common phenotypic traits (Webb et al., 2002), the host plant would select for OTUs with particular biological features, adapted to a symbiotic life style. This implies that traits favouring symbiosis would be phylogenetically conserved in fungi. Environmental filtering is probably not an exclusive mechanism and other processes, such as inter-specific microbial interactions, might explain the observed phylogenetic structure. For instance, if competitive ability is assumed to be correlated with phylogenetic distance then competition can drive phylogenetic clustering (Mayfield & Levine, 2010). Mechanisms leading to an over-dispersed phylogenetic structure, especially inter-species competition, which classically leads to species exclusion between relatives (Diamond, 1975), would have a weaker effect than environmental filtering by the host. Kin selection, a strategy favouring the reproductive success of an organism’s relatives, might also be important in structuring fungal communities.

Disentangling the processes underlying the observed phylogenetic structure

The relative importance of the above-mentioned processes in explaining the observed phylogenetic clustering can be elucidated by detecting the phylogenetic signal of relevant trait(s) (Mayfield & Levine, 2010). However, phenotypic traits are especially difficult to measure in complex communities, particularly for uncultivated microbes. In this study, we used the ‘mean expression ratio’ to access the microbial metabolic status. Sample-to-sample variations were taken into account using the Kse index that allowed the power of the test to be increased and avoided an estimation bias (Hardy & Pavoine, 2012). Interestingly, a significant phylogenetic signal was found for all fungal groups (Table 2). Within the Ascomycota and Glomeromycota, the ‘mean expression ratios’ of evolutionary-related OTUs were more divergent than would be expected under a Brownian motion model of trait evolution (Table 2, Figs. S4, S5). Conversely, the OTUs in the Basidiomycota were more similar than expected, i.e., related OTUs shared similar expression ratios (Table 2 and Fig. 3). Thus, the combined information about the phylogenetic structure (Table 1) and the phylogenetic signals analyses (Table 2) suggested that the Basidiomycota assemblage was mainly governed by environmental filtering, favouring the co-existence of related and similar OTUs in their ‘mean expression ratio’. In contrast, the Ascomycota and Glomeromycota assemblages were more impacted by competitive interactions promoting the co-existence of phylogenetically related but dissimilar OTUs in their ‘mean expression ratio’. In this latter case, competitive ability would be positively correlated with phylogenetic distance or, competition would drive character displacement rather than OTUs exclusion. Other relevant functional traits now need to be used to confirm our interpretations and to improve our detection and understanding of the phylogenetic signal. One possibility would be to perform a comparative transcriptomic analysis of the fungal root microbiome, although a number of methodological padlocks would first need to be broken, notably the selective extraction of fungal RNA.

Conclusion and prospects

We describe here the root fungal microbiome associated with an A. stolonifera population using a molecular strategy combining DNA- and RNA-based approaches. We were able for the first time to draw up a comprehensive picture of the phylogenetic patterns existing in the fungal root microbiome by examining the phylogenetic structure and measuring the phylogenetic signal. We found that the fungal communities were not randomly assembled but instead appeared to be specifically filtered by their plant host. We thus provide new insights into the rules of assembly governing the root fungal microbiome community. The limited number of OTUs shared by all individuals and the clustered phylogenetic structure suggested that each plant recruits a particular microbial community to adapt to environmental conditions at a microscale.

The use of evolutionary information to describe an ecological pattern is a first step towards a full understanding of community assembly. Further experimental studies are now needed to focus on the processes underlying these patterns.