Statistical analyses
What microbial guilds cause PSF?
To identify microbial OTUs specifically associated with positive or negative PSF for Bromus and/or Koeleria , we conducted an indicator species analysis (Dufrêne & Legendre, 1997). To this end, we first created bins of soil inocula that were associated with altered biomass relative to uninoculated plants for either of our plant species, or both. We had 9 bins, with a factorial combination of plants to be influenced (2 species: Bromus or Koeleria ) and 3 types of impact on plant growth (positive, negative or neutral) (32 = 9 combinations). To get a plant growth threshold from which plant performance would be considered increased or decreased, we computed plant growth responses to all inocula as t -scores:
\begin{equation} t=\frac{G_{\text{in}}-\ G_{\text{st}}}{\left(\frac{\text{SD}_{\text{in}}}{\sqrt{n}}\right)}\nonumber \\ \end{equation}
, where Gin represents mean growth with the inoculum, Gst is the mean growth in sterile soil only, and SDin and n are, respectively, the standard deviation in growth of inoculated pots, and the number of inoculated pots. We thus considered the soil inoculum to have a positive impact on plant growth (i.e. a positive plant-soil feedback) if growth was higher in inoculated soil relative to uninoculated soil, and vice versa for negative PSF. We identified a threshold t -score corresponding to an 80% confidence interval with our corresponding number of degrees of freedom (d.f. = 13). In other words, we considered the PSF significantly different from 0 if there was less than a 20% probability that such a growth difference inoculated and control plants had arisen by chance. As this 20% threshold was selected arbitrarily, we conducted a sensitivity analysis and show that the number of indicator OTUs retrieved is fairly robust to the threshold chosen (Fig. S2). The indicator species analyses were conducted using the R package indicspecies (De Caceres & Legendre, 2009). Despite the high number of individual tests in these analyses, we kept the α threshold at 0.05. P -value correction methods would have yielded way too stringent thresholds to consider associations as significant, given the very high number of rare microbial OTUs. We thus validated the robustness of our indicator species analysis to type I errors by running 1000 indicator species analyses on randomized datasets (microbial metacommunities shuffled by keeping number of sequence reads constant by soil sample and by microbial OTU), and found that on average, by chance, we retrieve ~ 40 and 80 indicator prokaryotes and fungi respectively (Fig. S3). However, we also compared the identity of the indicator OTUs retrieved by multiple indicator species analyses on our real dataset vs. on our randomized datasets, and found that indicator OTUs were very constant from run to run in the former, but much more variable in the latter (Fig. S4). This shows that the indicator OTUs reported here are not likely to be random subsets of our metacommunity that have arisen as type I errors.
In order to determine the contribution of different microbial guilds to PSF, we assigned indicator microbial OTUs to guilds based on the literature. For prokaryotes, we relied on published papers, manually querying scientific search engines for papers mentioning these taxa. Prokaryotic OTUs that could not be assigned to a tentative guild (21 out of 66 indicator OTUs) were left out of downstream analyses. To assign fungal OTUs to guilds, we used the open-source database FunGuild (Nguyenet al. , 2016). After guild assignment, we conducted χ2-tests to determine if certain microbial guilds were disproportionately associated with certain PSF types (e.g., pathogens overrepresented as indicator OTUs of negative PSF). We note here that with indicator species analysis, we were unlikely to assign generalist/ubiquitous taxa to any of our specific PSF “bins”, because the analysis is by essence based on variations in relative abundances. Thus, if some pathogens, for example, are very widespread in the microbial metacommunity, they are unlikely to be retrieved as indicator of negative PSF. To address this issue, we systematically assigned all fungal OTUs to a guild using FunGuild, and we assigned the most widespread prokaryotes OTUs to guilds manually. This way, we could validate that the absence of a given guild to particular PSF bins was not an artefact, due to the ubiquity of taxa from that guild.
Can we predict PSF based on microbial α or β-diversity?
We used linear mixed models (LMMs) to regress measured PSF against prokaryotic and fungal (1) richness and (2) Shannon’s diversity (exponential form; Hill, 1973), as estimators of α -diversity. We also included prokaryotic and fungal community structure as predictors using site loadings from a non-metric multidimensional scaling (using Hellinger-transformed microbial community data) constrained over two axes. We constructed separate models to predict Koeleria vsBromus growth responses, including the identity of the block from which the inoculum came from as a random factor in the models, and the state variable indicating whether the sample came, within the block, from the Bromus -invaded plot vs the native vegetation plot. We used the R package nlme (Pinheiro et al. , 2020) for LMMs construction.
Do microorganisms causing PSF respond more strongly to plants or soil?
We used redundancy analyses to regress microbial community structure against environmental data matrices. Environmental matrices combined (1) soil properties (Mehlich-III extractible P, NH4+, NO3-, pH and gravimetric moisture measured in fresh soils, as described in Chagnon et al. , 2018) and (2) plant community structure. To include the latter as a predictor of microbial community structure, we used the first 5 axes of a principal component analysis (using Hellinger-transformed plant abundance data), which collectively explained 75% of the variance in plant community structure (Table S1). This yielded RDAs with 10 explanatory variables (5 soil and 5 plant predictors, respectively). We ran RDAs separately for prokaryotes and fungi. For each, we conducted one RDA with the indicator OTUs, and one with the rest of the microbiome (i.e., the non-indicator OTUs). Then, to verify whether the indicator OTUs responded idiosyncratically to the explanatory variables (hypothesizing that they would be tied more strongly to plant community structure), we evaluated the co-inertia between the axes loadings of the explanatory variables using Procrustes rotations using the R packagevegan (Oksanen et al. , 2019). A poor correlation between these axes loadings would indicate those explanatory variables that best describe the indicator OTUs are not the same as for the rest of the microbiome.