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.