Archaea and bacteria mediate the effects of native species root loss on fungi during plant invasion

15 Pages·2017·1.74 MB·English

by  Steven D Mamet

TheISMEJournal(2017)11,1261–1275
©2017InternationalSocietyforMicrobialEcology Allrightsreserved 1751-7362/17
www.nature.com/ismej
ORIGINAL ARTICLE
Archaea and bacteria mediate the effects of native
species root loss on fungi during plant invasion
Steven D Mamet1, Eric G Lamb2, Candace L Piper2, Tristrom Winsley1
and Steven D Siciliano1
1Department of Soil Science, University of Saskatchewan, Saskatoon, Saskatchewan, Canada and
2Department of Plant Science, University of Saskatchewan, Saskatoon, Saskatchewan, Canada
Althoughinvasiveplantscandriveecosystemchange,littleisknownaboutthedirectionalnatureof
belowgroundinteractionsbetweeninvasiveplants,nativeroots,bacteria,archaeaandfungi.Weused
detailed bioinformatics and a recently developed root assay on soils collected in fescue grassland
along a gradient of smooth brome (Bromus inermis Leyss) invasion to examine the links between
smooth brome shoot litter and root, archaea, bacteria and fungal communities. We examined
(1) aboveground versus belowground influences of smooth brome on soil microbial communities,
(2)theimportanceofdirectversusmicrobe-mediatedimpactsofplantsonsoilfungalcommunities,
and (3) the web of roots, shoots, archaea, bacteria and fungi interactions across the A and B soil
horizonsininvadedandnon-invadedsites.Archaeaandbacteriainfluencedfungalcomposition,but
not vice versa, as indicated by redundancyanalyses. Co-inertia analyses suggested that bacterial–
fungal variance was driven primarily by 12 bacterial operational taxonomic units (OTUs). Brome
increased bacterial diversity via smooth brome litter in the A horizon and roots in the B horizon,
whichthenreducedfungaldiversity.ArchaeaincreasedabundanceofseveralbacterialOTUs,andthe
key bacterial OTUs mediated changes in the fungi’s response to invasion. Overall, native root
diversitylossandbacterialmediationweremoreimportantdriversoffungalcompositionthanwere
thedirecteffectsofincreasesinsmoothbrome.Critically,nativeplantspeciesdisplacementandroot
lossappearedtobethemostimportantdriveroffungalcompositionduringinvasion.Thiscausalweb
likelygivesrisetotheplant–fungifeedbacks,whichareanessentialfactordeterminingplantdiversity
ininvadedgrasslandecosystems.
TheISME Journal (2017) 11,1261–1275; doi:10.1038/ismej.2016.205; published online31January 2017
Introduction cycling (Elgersma et al., 2011; Shannon-Firestone
et al., 2015).
Invasive plant species are important drivers of
Several recent studies suggest that plant invasion
ecosystem change through their influence on below-
successislinkedtosuppressionand/orpromotionof
groundmicrobialcommunities(vanderPuttenetal.,
microbial species (Bever et al., 2012; van der Putten
2001). Soil microbial interactions are critical deter-
et al., 2013). For example, plant invasion may
minants of aboveground diversity (Klironomos,
selectively suppress dominant bacterial species,
2002), yet little is known about how invasion alters
allowing for rare and/orammonia-oxidizing bacteria
the directional nature of interactions among micro-
(AOB) to increase in relative abundance, which in
bial groups (although see Schröter et al., 2004;
turn has an impact on nitrogen cycling (Ehrenfeld,
Scherber et al., 2010). For example, the switch from
2010; Piper et al., 2015a). However, plant effects are
negative plant–soil microbial feedbacks in native
notconstantamongplantspecies.Inparticular,some
habitats to positive feedbacks in invaded habitats
invasive species may decrease fungal abundance
may be crucial for invasive plant success (Callaway
(Vogelsang and Bever, 2009), thus increasing the
et al., 2004; Diez et al., 2010). How this switch
bacteria to fungi ratio in the soil; other species may
occurs is still open to debate; however, it is likely
increase particular fungal groups (Kourtev et al.,
linkedtoabovegroundcarbonsupply(Bradleyetal.,
2002). Temporal lags between invasion and accu-
2006)orchangesinthecompositionandstructureof
mulation of species-specific pathogens in the soil
soil microbial communities and altered nutrient
may also allow invasive plants to flourish at the
expense of native species (Klironomos, 2002),
Correspondence: SD Siciliano, Department of Soil Science, although in some instances pathogenic fungi may
University of Saskatchewan, 51 Campus Drive, Saskatoon, prevent invasion (Reinhart et al., 2005). Despite the
Saskatchewan,S7N5A8,Canada. clear importance of plant–microbial interactions
E-mail:[email protected]
during invasion, the causal chain linking changes
Received 4 October 2016; revised 16 November 2016; accepted
15December2016;publishedonline31January2017 in aboveground to belowground diversity and the
Bacteriamediatefungiduringplantinvasion
SDMametetal
1262
evolution of plant–microbial feedbacks is not clear. issomeevidencethatmoregenerallinksbetweenthe
Some authors suggest that plant effects are dam- fungal and bacterial communities may have a major
pened as they percolate through soil trophic levels role in structuring microbial communities (for
(Scherber et al., 2010), whereas others suggest that example, Siciliano et al., 2014).
plant effects are limited to specialist microbial Here we present a model to study interkingdom
groups (van der Putten et al., 2007). interactions and responses to plant invasion in a
Fungal communities are critical in any considera- smooth brome-invaded fescue grassland. Using
tion of plant–soil feedbacks, as fungi appear to samples from the smooth brome invasion gradient,
determine vegetation composition through structur- we assessed how the structure of the aboveground
ing plant–plant interactions (Bennett and Cahill, plant and its root communities interact with below-
2016) and plant abundance (Reinhart et al., 2005). ground archaea, bacteria and fungi to influence
Howeveritisnotclearhowdirectionalrelationships soil ecosystem composition. We hypothesized that
among fungi, bacteria and archaea evolve during smoothbrome-structuredfungiviaroots(cf.Johnson
invasion. Fungi may metabolize root exudates, and etal.,2004)andthat,inturn,thefungiaffectedother
those metabolites then drive bacterial community microbes. To address this hypothesis, we asked a
composition (Millard and Singh, 2009) or, alterna- series of linked questions: (1) does the influence
tively, bacteria can suppress fungal pathogens of invasive smooth brome on soil microbial commu-
(Kobayashi and Crouch, 2009). Relative to bacteria, nities differ between aboveground litter and below-
muchlessisknownabouttherelationshipofarchaea ground root effects? (2) What is the relative
withothermicrobes.Weknowarchaeaarecommon, importance of direct versus microbe-mediated
although not major, inhabitants of roots (Buée impacts of invasive smooth brome on soil fungal
et al., 2009, although see Simon et al., 2000 for an communities? (3) Can we determine the directional
exception), and in some instances fungi may reduce nature of the web of roots, litter, archaea, bacteria
ammonia-oxidizingarchaeaandAOBabundancevia and fungi interactions across A and B horizons of
plant root pathways (Chen et al., 2013), although soils in invaded and non-invaded sites?
often there is no relationship between archaeal and
fungal diversity at broad spatial scales (Prober et al.,
2015). This is likely because archaea, unlike fungi, Materials and methods
typically respond more to soil and environmental
properties than to plant cover (Nielsen et al., 2010). Study area
Thus, from a plant–microbe perspective we may The field site is a 14.6ha native fescue prairie
conceptualize microbial community responses to undergoing invasion by smooth brome, ~120km
plants at a field scale as a gradient, where fungi southofSaskatoon,Saskatchewan,Canada(51°12′N,
respond mainly to vegetation and only slightly to 107°17′W). The site is near the border of the moist
soil,bacteriarespondstronglytobothvegetationand mixed and mixed prairie ecoregions and within
soil, and archaea respond mainly to soil and not to the Orthic Dark Brown Chernozemic soil order (Soil
vegetation. However, there is still considerable Classification Working Group, 1998). The grasses
uncertainty in how the structure of intermicrobe Festuca hallii and several species of Hesperostipa,
interactions changes in communities undergoing ElymusandPascopyrumdominatethenativeprairie
plant invasion. at this site. Smooth brome is invading the site from
Smooth brome (Bromus inermis Leyss) invasion disturbededges(primarilyroads),andmanyinvaded
lowers plant species diversity in both mixed-grass patchescanbefoundintheinteriorofthesite(Piper
prairie (Fink and Wilson, 2011) and fescue grass- et al., 2015a, b). At the time of sampling there were
lands(Wilson andPinno,2013).Thisinvasion alters no cattle grazing on the site or other management
soil conditions (that is, increased nitrogen availabil- practices applied. The mean annual temperature is
ity; Piper et al., 2015b), thereby initiating a cascade 3.5°C, and 376.9mm of precipitation falls each year
of soil community changes, including the suppres- (Rock Point weather station, ~7km from the field
sion of dominant members of the bacterial commu- site, Environment Canada, 2016).
nity (Piper et al., 2015a). However, no one has
conclusively identified an above- or belowground
mechanism linking brome invasion to bacterial Sample collection
diversity. Rather, structural equation models (SEM) We collected plant and soil samples from a broad
have indicated that an unknown direct pathway range of smooth brome cover classes (Piper et al.,
driven by smooth brome abundance is responsible 2015a, b). A total of sixty 50×50cm plots were
for the changes in bacterial diversity (Piper et al., establishedinastratifiedrandommanner.Toensure
2015a).Onehypothesisisthattheunknownmechan- even sampling across the full range of the brome
ism may result from the interactions between invasion (0, 40–50, 51–85 and 485%), each plot
bacterial and fungal communities (Frey-Klett et al., had14–16 samples ineach offourclasses of smooth
2011). These bacterial–fungal interactions could be brome shoot abundance. Sampling locations were
specific pathogen interactions (for example, as out- determined using the random point generator in
lined in van der Putten et al., 2007); however, there ArcMap (ESRI, Redlands, CA, USA). Within each
TheISMEJournal
Bacteriamediatefungiduringplantinvasion
SDMametetal
1263
plot, we assessed plant species cover and gathered 2009) and A344F and A589R primers for archaea
grass, forb, shrub and litter biomass. Smooth brome (Teske and Sørensen, 2007). These amplicon
biomass was collected separately from other grass libraries were sequenced from the forward primer
species. on the 454 FLX titanium platform (Roche, Branford,
Biomass samples were dried for 2 days at 60°C CT, USA) at Genome Quebec. We note here,
and weighed. We measured the depth of the A however, that the ITS1F and ITS4 primers may not
horizon (defined as mineral soil with less than 17% effectively amplify arbuscular mycorrhizal fungi
organic carbon) using changes in soil color and (Stockingeretal.,2010).Thus,wesuggestourresults
texture (Soil Classification Working Group, 1998). be interpreted in terms of the overall microbial
UsinganAMSsoilcorer(AMS,Inc.,AmericanFalls, community, and not specifically the arbuscular
ID), we extracted two 5-cm-diameter soil cores from mycorrhizal fungi community—which is important
the top 5cmof each of the A andB horizons at each but beyond the scope of this study.
plot. To minimize cross contamination, A and B
horizon samples were collected and analyzed sepa-
rately because (1) soil horizons are created by Bioinformatics processing
The 16Sgeneisoftenusedasaspeciesdiagnosticfor
fundamentally different soil-forming processes, and
bacteria and archaea (Stackebrandt and Goebel, 1994;
(2) ammonia-oxidizing microbes respond differently
to smooth brome abundance in each horizon (Soil Bakeretal.,2003),althoughithasfewerhypervariable
domainsinfungi.Amongtheregionsoftheribosomal
Classification Working Group, 1998; Piper et al.,
cistron,theinternaltranscribedspacer(ITS)regionhas
2015b). The two soil cores from each horizon were
then combined, resulting in one composite soil core the highest probability of successful identification
for the widest range of fungi (Schoch et al., 2012).
from each horizon at each plot. These samples were
frozen at –20°C until they could be analyzed. Therefore, we used the 16S and ITS sequences to
identify and quantify bacterial and archaeal composi-
To characterize roots, we used a recently devel-
tion and fungal composition, respectively.
oped assay capable of identifying and quantifying
plant roots from DNA extracts of mixed-species root
samples (Hiiesalu et al., 2012; Lamb et al., 2016; Bacterial and archaeal 16S. Data from the Ion
see Appendix S3 in Supplementary Information). torrent and 454 Platforms were processed for
Briefly, we modified the trnL primers used pre- bacterial and archaeal 16S rRNA gene analysis,
viously by other researchers (Taberlet et al., 1991; respectively (see Appendix S2 in Supplementary
2007)andusedanimprovedbioinformaticspipeline Information). Seed sequences for each bacterial and
that referenced site-specific taxon lists to identify archaeal operational taxonomic unit (OTU) were
plants to a species level (for detailed information on taxonomicallyidentifiedusingtheGreenGenesdata-
the approach, see Lamb et al., 2016). Bacterial, base (August 2013 release) using a naive Bayesian
archaeal and fungal community compositions were classifier algorithm implemented in mothur (Wang
assessed by high-throughput amplicon library et al., 2007; Schloss et al., 2009). We calculated two
sequencing, as described in Piper et al. (2015a) and measures of alpha diversity: community evenness
Siciliano et al. (2014). Briefly, DNA was extracted using Evar (Smith and Wilson, 1996) and observed
from 0.5g of 2-mm sieved soil using the Ultraclean richnessasthetotalnumberofOTUsineachsample
Soil DNA Extraction Kit (MoBio, Carlsbad, CA, normalizedtoaspecificnumberofreadspersample.
USA). DNA concentration was determined using a In addition to these measures of alpha diversity, we
ultraviolet–visible spectrophotometer (Nanodrop calculated phylogenetic diversity for archaea and
2000, ThermoScientific, Wilmington, DE, USA). For bacteria by considering the sequence data and the
bacteria,weusedtheuniversal16SrRNAprimerset evolutionary relationships in the communities by
515F/806R,whichamplifiesa291-bpfragmentnear implementing the UniFrac algorithm (Lozupone and
the bacterial v4 region (Caporaso et al., 2011). Knight, 2005; Lozupone et al., 2011). The UniFrac
Primers contained an Ion Torrent adapter and a approach incorporates evolutionary distance at the
unique barcode sequence for sample pooling. Sam- sequence level and provides a more robust esti-
ples were amplified in triplicate using a 25-μl mate of community composition than OTU-based
reaction mix containing 18μl Platinum Blue Super- approaches. For fungal community composition, we
mix (Invitrogen), and sample replicates were pooled were unable to use the UniFrac approach because
and purified using a QiaQuick gel extraction kit the differing read lengths of the intragenic spacer
(Qiagen Inc., Toronto, ON, Canada). Following region would bias UniFrac’s sequence-based analy-
pooling, the samples were sequenced according to sis. Thus, we used principal component analysis to
the Ion PGM 200 Sequencing Kit v2 (Life Technol- characterize the fungal OTU matrix and used the
ogies). Similarly, fungal and archaeal diversity were scores for principal component one as a metric of
assessed via triplicate reactions using ITS1F and phylogenetic diversity.
ITS4 primers (Gardes and Bruns, 1993; Manter and
Vivanco, 2007) for fungi (acceptable for differentiat- Fungal ITS. ITS sequences were dereplicated and
ing fungi at the genera level and coarser taxonomic denoised, and chimeric sequences were removed.
resolutions, cf. Krüger et al., 2009; Stockinger et al., Remaining sequences (4200bp) were further
TheISMEJournal
Bacteriamediatefungiduringplantinvasion
SDMametetal
1264
dereplicated, such that reads identical with longer matrix (Legendre and Legendre, 2012). From one
reads were added to the longer read count (cluster point of view, researchers may consider rare species
fragments function in mothur). Next, we used the as potential indicators of special environmental
OTU picking function in QIIME (pick_otus.py with conditions.However,theprimarypurposeofordina-
usearch61 and a 0.90 sequence similarity threshold) tion is to display the main axes of variation of the
to generate a list of OTUs and pick_rep_set.py to data, not to deal with exceptions. On the other
establish a file of seed sequences for the selected hand, researchers may consider the occurrence of
OTUs (Caporaso et al., 2010). These OTUs were rare species as a chance event that should not be
classified against the UNITE database (mothur for- heavily favored in the analysis. To mitigate the
mat version –v6_sh_99) to provide a list of taxo- bias of zero-inflated data, we used a stepwise
nomic identifications (Kõljalg et al., 2005). OTUs elimination procedure to remove these species,
with identical classification were aggregated into which progressively downweighted and removed
single OTUs and their abundance summed using rare species (Figures 1a and b; see Supplementary
Primer6 (Clarke, 1993). Appendix S3 in Supplementary Information;
Legendre and Legendre, 2012).
We acknowledge that the topic of rare species
Statistical analyses removal is a contentious and unresolved issue (for
Overview. We used redundancy analysis (RDA) to example, Marchant, 1999; Cao et al., 2001) and
identifysignificantcovariatesamongmicrobialcom- recognize both the potential ecological significance
munities, followed by co-inertia analysis (CoIA) to of some rare species and statistical quagmire of rare
identify keystone OTUs (addressing questions 1 and species present due to chance (Lynch and Neufeld,
2 in the introduction), and then confirmed these 2015). Numerous methods exist for downweighting
causal assumptions using SEM (question 3 and see rare species (for example, Poos and Jackson, 2012),
SEM of plant, root and soil microbial relationships although many lack strong ecological foundations
below). (seeadiscussiononthistopicinLynchandNeufeld,
2015). Therefore, we ran our RDA analyses on both
Rare species. Species that occur in low quantities thefullanddownweighteddatamatricestoevaluate
across a study site (few occurrences and low the effect of rare species removal on variation
abundance) generate a large number of zeros in data in community composition. We found statistical
matrices. Because zeros have high leverage, they relationships among microbial composition and
contribute heavily to the total inertia of the data environmentalcovariateswerestronger(asindicated
Figure1 Sequenceofstatisticalanalysesusedtoassessfunctionalrelationshipsamongsmoothbromeandsoilmicrobialcommunities.
(a)Distinctdatamatriceswereusedforarchaea,bacteria,fungi,andexplanatorycovariates(bottomleft)alongabromeinvasiongradient.
(b)OTUswithsmalloccurrencevaluesonlyappearinafractionofthesamplingplotsandgeneratealotofzerosinthedatamatrices.
Becausezeroshavehighleverageandcontributeheavilytothetotalinertiaoftheordinations,wechosetoeliminaterarespeciesusinga
stepwiseeliminationprocedure(LegendreandLegendre,2012;seeMaterialsandmethods).(c)RDAwasusedtosummarizesignificant
linear relationships among microbial communities, and between communities and above- and belowground explanatory covariates.
(d)Principalcomponentsanalysis(PCA)wasrunoneachmicrobialdatamatrixtofindthemaximalinertiaprojectedontwoaxes.The
reductionindimensionsofeachdatasetresultedinindependentstructuresforeachcommunity.CoIAtestsforagreementbetweentwo
structures(microbialand/orrootcommunities)bycalculatingaxesthatmaximizethecovariancebetweenthefactorialscoresofsamples.
(e)SEMusestwoormorestructuralequationstomodelmultivariaterelationships.HereweusetheresultsoftheSEManalyses,basedon
priorexperienceandsubstantiveknowledgegainedfromthepreviousstatisticalanalyses,toreflectcausaldependenciesbetweensmooth
bromeandsoilmicrobialcommunitiesinanaturalgrasslandunderinvasionbysmoothbrome.
TheISMEJournal
Bacteriamediatefungiduringplantinvasion
SDMametetal
1265
Figure 2 Conceptual model of aboveground (green arrows), belowground (brown and gray arrows), and interkingdom influences
(large light blue arrows) on soil microbial communities along cover classes of smooth brome. Conceptual arrows representing
aboveground,belowground,andinterkingdominfluencesarebasedonsignificantexplanatoryvariablesasdeterminedthrougheachRDA.
BiplotsofRDAsforeachkingdom((a)bacteria,(b)fungi,and(c)archaea)arepresentedwithsignificantexplanatoryvariables(boldblue
arrows within biplots; P⩽0.05). Box color in each RDA represents smooth brome cover classes and point shape represents A and B
horizons(seelegendattopoffigure).
Table1 ResultsofRDAofsoilmicrobialabundance
Kingdom Covariate Varianceexplained F VIF P-value
Archaea Litternitrogen 0.015 1.627 1.000 0.001***
Finalmodel 0.015 1.627 — 0.001***
Bacteria Litternitrogen 0.024 12.803 1.046 0.001***
Rootbiomass 0.012 6.189 1.582 0.001***
SOC 0.008 4.376 1.605 0.003**
AOA 0.006 3.225 1.057 0.006**
Finalmodel 0.207 6.649 — 0.001***
Fungi Rootbiomass 0.022 3.265 1.622 0.001***
Bacterialrichness 0.014 2.094 1.166 0.003**
SOC 0.013 1.966 1.847 0.002**
AOB 0.012 1.762 1.427 0.006**
Litternitrogen 0.011 1.640 1.178 0.017*
Archaearichness 0.010 1.444 1.136 0.044*
Finalmodel 0.111 2.028 — 0.001***
Abbreviations:AOA,ammonia-oxidizingarchaea;AOB,ammonia-oxidizingbacteria;ANOVA,analysisofvariance;RDA,redundancyanalysis;
SOC,soilorganiccarbon;VIF,varianceinflationfactor.
Belowgroundexplanatoryvariables:soilpH,SOCandtotalsoilnitrogen.Abovegroundvariables:litternitrogen,rootbiomass,rootevenness,and
smoothbromecover.Interkingdomvariables:AOA,AOBandrichness,evennessanddiversityforarchaea,bacteria,andfungi.Shownare
explainedvariance,F-values,VIFs,andP-values(ANOVA).Thefinalmodelincludesalllistedcovariates.*Po0.05,**Po0.01and***Po0.001
significance.
TheISMEJournal
Bacteriamediatefungiduringplantinvasion
SDMametetal
1266
Table2 ResultsofPERMANOVAtestingtheeffectsofsmoothbromecoverclass(0,0–50,50–85,85–100%),soilhorizon(A,B),andtheir
interactionsonmicrobialcommunitycompositionalsimilaritywithinasmoothbrome-invadedfescuegrassland(Hellingerdistance)
Kingdom Effect d.f. SumofSqs MeanSqs F Varianceexplained Pr(4F)
Archaea Coverclass 1 1.053 1.053 2.147 0.020 0.001***
Soilhorizon 1 0.182 0.182 0.371 0.004 1.000
Cover×horizon 1 0.124 0.124 0.253 0.002 1.000
Residuals 103 50.509 0.490 — 0.974 —
Total 106 51.868 — — 1.000 —
Bacteria Coverclass 1 0.134 0.134 2.692 0.022 0.025*
Soilhorizon 1 0.600 0.600 12.072 0.101 0.001***
Cover×horizon 1 0.109 0.109 2.197 0.018 0.074
Residuals 103 5.118 0.050 — 0.859 —
Total 106 5.960 — — 1.000 —
Fungi Coverclass 1 0.460 0.460 1.901 0.017 0.006**
Soilhorizon 1 1.581 1.581 6.534 0.060 0.001***
Cover×horizon 1 0.214 0.214 0.883 0.008 0.609
Residuals 100 24.195 0.242 — 0.915 —
Total 106 24.450 — — 1.000 —
Abbreviations:d.f.,degreeoffreedom;sq.,square.
Communitydatamatriceswerestatisticallydownweighted(rarespeciesremoved)beforethePERMANOVAtest.*Po0.05,**Po0.01,and
***Po0.001significance.
by higher F-statistics) after rare species removal isolate larger-abundance influential consortia, and
(Figure 2, Supplementary Figure S4; compare present this Poisson–RDA method as an alternative
F-values in Tables 1 and 2 of the downweighted approach to identify influential conditionally rare
data set to the corresponding Supplementary Tables taxa such as pathogens (see Appendix S4 in
S1 and S2, which are of the full data set). We Supplementary Information and SEM results).
interpretedthestrongerpatterninthedownweighted
dataasanincreaseinoursignal-to-noiseratio,anda Soil microbial composition and keystone species.
winnowingofthefullmicrobialdatamatricestoonly We identified significant covariates among below-
the most influential abundant species. Therefore, ground, aboveground, and interkingdom influences
we focused all subsequent analyses on the down- on microbial composition using RDA on Hellinger-
weighted data sets. transformed reduced data (Figure 1c). RDA is a
As another check of rare species influence, we method of summarizing linear relationships in a set
assessed whether individual bacterial OTU abun- of dependent variables influenced by a set of
dancewasassociatedwithsmooth bromecover,and independent variables, using a blend of multiple
whether these OTUs were associated with environ- linear regression and principal components analysis
mental covariates. To do this, we first assessed the (LegendreandLegendre,2012).RDAisaconstrained
relationship between OTU abundance and smooth ordinationmethodthatcomputesaxesthatarelinear
brome cover via zero-inflated Poisson regression combinations of the explanatory variables (Borcard
models (Zuur, 2009) followed by Bonferroni correc- et al., 2011). Thus, the axes produced using RDA
tion (Shaffer, 1995), and then ran an RDA using model the variation of the dependent variable,
these correlated OTUs as covariates on the environ- allowing one to test a hypothesis (H ) of absence of
0
mental data matrix. We found that 88 OTUs were a linear relationship between microbial composition
stronglycorrelatedwithsmoothbromecover(n=88) and explanatory independent variables (Figure 1c).
using a Po0.001 threshold. Of these 88 OTUs, Thefollowingfourgroupsofexplanatoryvariables
only eight were correlated with the environmental wereusedtomodelmicrobialcomposition(bacteria,
matrix (Po0.05). These OTUs were from the phyla archaea, and fungi): (1) soil influences (soil pH,
Acidobacteria, Actinobacteria, and Planctomycetes soil organic carbon (SOC), and soil total nitrogen),
and represented a mean relative abundance of (2) aboveground plant influences (litter nitrogen,
0.29±0.03%, compared with 1.27±0.04% for the smooth brome shoot abundance), (3) belowground
downweighted data OTUs. We included these eight plant influences (root biomass, root evenness,
OTUs as endogenous variables in our final SEM and root abundance of smooth brome), and (4)
models with direct links between brome roots and microbial community influences (ammonia-oxidizing
shoots, and fungal evenness (see below). We found archaea and AOB abundance; bacterial, fungal, and
only two significant relationships: (1) smooth brome archaealevenness,richness,anddiversity[H′];Piper
roots increased the abundance of a bacterial OTU of et al., 2015a, b). Explanatory variables selected for
class Actinobacteria and (2) pH reduced the abun- retention were based on permutation tests (Borcard
dance of a bacterial OTU of class Kueneniae, which et al., 2011) and we used a permutational multi-
itselfreducedfungalevenness.Thus,wepresentour variate analysis of variance (PERMANOVA; Anderson,
downweighting–RDA–CoIA approach as a means to 2001)totestforsignificantdifferencesamongbrome
TheISMEJournal
Bacteriamediatefungiduringplantinvasion
SDMametetal
1267
cover categories and soil horizon. We assumed a relationships (Figure 1e). Whereas many multivari-
direction of influence based on the nature of ate methods are largely descriptive and more appro-
explanatory variables. For example, we found bac- priate for exploratory analyses, SEM is capable of
terial richness was a significant predictor of fungal testing a network of causal hypotheses and is
composition,butfungaldiversitywasnotapredictor recommended for the evaluation of multivariate
of bacterial composition. We interpreted this as a hypotheses (Grace, 2006; Grace et al., 2012). Speci-
primarily unidirectional influence of bacteria on fically, we used SEM because it allows the evalua-
fungi. Though there may be exceptions to this, our tion of simultaneous influences (for example, brome
statistical framework supported this directional may influence fungal composition both directly and
nature and thus we maintained this direction of through altering bacterial diversity) rather than
influence through subsequent analyses. individual causes (for example, brome influences
We implemented CoIA (Dolédec and Chessel, fungal composition only directly). The method is
1994) using the package ‘ade4’ version 1.7-2 in R thus appropriate for establishing probable causality
(Dray and Dufour, 2007) to examine important atthesystem(forexample,brome–archaea–bacteria–
correlations between kingdoms as revealed by the fungi) rather than the individual level (for example,
RDAs(Figure1d).InCoIA,multivariateanalyseslike brome–fungi).
principal components analysis are processed on SEM rely on specification by researchers of a
each community matrix separately to find the network of a priori causal assumptions based on a
maximal inertia (that is, variance) projected onto scientific body of evidence (for example, that pH is
two axes. The reduction in dimensions of each animportantdriverofbacterialrichnessandnotthat
data set results in independent structures for each bacterial richness determines soil pH), and then
community. CoIA tests for agreement between these testing whether that causal network is consistent
independent structures (microbial and/or root com- with empirical data (Grace, 2006; Pearl, 2012;
munities) by calculating axes that maximize the Shipley, 2016). The postulated causal network gives
covariance between the factorial scores of samples. risetoaseriesoflinearequations,whichinturngive
Here we first reduced the dimensionality of the rise to an expected covariance matrix. The expected
archaeal,bacterial,andfungalmatricesusingPrinci- covariance matrix is then compared with the
pal components analysis, and then selected the observed covariance matrix arising from the data.
dominant components (axes). The CoIA then gener- A statistically acceptable congruence between the
atednewaxesbyrotationinmultidimensionalspace expected (causal model implied) and observed
soastomaximizethecovariancebetweentheaxesin covariance matrices is thus an empirical validation
the two data sets (Dolédec and Chessel, 1994). of the causal assumptions used. In other words, an
Two eigenvectors were selected from each kingdom SEM is a ecological theory describing a particular
matrix for inclusion into the CoIA. RV coefficients system, and if congruent with the data, the theory is
(a multivariate generalization of the Pearson correla- supported. SEM is a well-established tool widely
tion coefficient) were calculated for each CoIA to used in the natural and social sciences for testing
measure global correlations between matrices. The causal inferences with correlative data sets; it is
significance of the associations between kingdom howevercriticallyimportantthatthecausalassump-
matrices was tested with a Monte Carlo permutation tions made by the researcher be well grounded in
test using 9999 permutations. prior studies, scientific knowledge, logical argu-
Preliminary analyses revealed significant correla- ments, and other evidence. An SEM model fitting
tionsbetweensitedistancesfromtheoriginoftheco- thedatadoesnotprovethecausalassumptionsused,
inertia plane and smooth brome cover for the soil A butreplicationofagivenmodelacrossmanysystems
horizonandnottheBhorizon;therefore,wefocused representsaverystrongtestoftheunderlyingtheory
ouranalysesontheAhorizon.WithintheAhorizon, (Grace et al., 2012).
we used t-tests followed by Bonferroni correction The initial (theoretical) SEM used smooth brome
(Shaffer,1995)todeterminewhichplotscontributed shoot biomass and root abundance as measures of
most to the common variance (total co-inertia) the above- and belowground severity of smooth
within the two matrices (that is, sites furthest from brome invasion; it also used bacterial richness,
theplotorigin),aswellasthemostimportantspecies evenness, and phylogenetic diversity as measures
drivingtheco-inertia(Figure1d).Wethenvisualized of bacterial diversity. Using RDA, we identified a
the configuration of the matrix pairs with joint directional link between archaeal diversity and both
site plots, in which arrows represent the distances bacterialandfungaldiversity(forexample,Figure2);
between matrices in n-dimensional space and thus, we included a directional pathway from
points represent the arrow midpoint (Legendre and archaea to bacteria and fungi in our SEMs. As pH
Legendre, 2012). is known to influence microbial community compo-
sition (Fierer and Jackson, 2006), we included
SEM of plant, root, and soil microbial relationships. direct relationships from pH to microbial diversity.
We used SEM (Grace, 2006; Lamb et al., 2011) Both competitive and mutualistic interactions exist
to investigate a range of potential mechanisms between bacteria and fungi (de Boer et al., 2005).
underlying smooth brome–microbial community However, our RDA analyses suggested the influence
TheISMEJournal
Bacteriamediatefungiduringplantinvasion
SDMametetal
1268
ofbacteriaonfungiisgreaterthanthereverseinthis Approximation, and Standardized Root Mean Square
system, so we included a directional pathway from Residuals), we evaluated this model structure for
bacteria to fungi in our model. eachbacterialOTUofinterestidentifiedthroughthe
We used multigroup SEMs with observed vari- CoIA(n=12),aswellasbacterialrichness,evenness,
ables, as we collected two subsets of data (A and B and phylogenetic diversity (n=3) for each fungal
horizon) from the same sample points (Grace, 2006; diversity metric (n=3) for a total of 108 models
Piperetal.,2015a).InamultigroupSEM,modelsare (12×3×3). As an additional check on the direction
initially constrained to ensure path coefficients are of causality of the diversity measures, we inverted
equal between groups (Grace, 2006; Lamb et al., the causal link between fungi and bacteria, that is,
2011). Constraints are then progressively released to that fungal richness, evenness, and phylogenetic
improve model fit. Different path coefficients diversity drove bacteria. For all three metrics,
between horizons suggest significant differences in evenness, richness, and phylogenetic diversity, the
the biological processes represented by that path χ2 was greater when bacteria were dependent on
betweenhorizons.Modelswerefitusingthe‘lavaan’ fungi. For example, for evenness, the χ2 increased
package v. 0.5–20 in R (Rosseel, 2012). In order to from 5.4 to 8.9, for richness the χ2 increased from
determine differences in biological processes 14.7 to 21.3 and for phylogenetic diversity the
between groups in the multigroup SEM, we initially χ2 increased from 6.0 to 6.1. The SEM analyses
constrainedthemodelssothatpathcoefficientswere suggested three classes of relationships among
equal between horizons. Constraints were progres- smooth brome and microbial diversity based on
sively released to improve model fit (cf. Piper frequency of microbial interactions. Specifically,
et al., 2015a). Path coefficient differences between if we found significant relationships between
horizons indicated a significant distinction in the endogenous variables in more than one SEM, we
biological processes represented by that path designated this as a ‘class’. For brevity, we present a
between horizons. Before fitting the SEM, we used summary of the three model classes here and an
generallinearmodelswithaquadratictermtocheck exampleofeachclassinthesupplementarymaterial
for nonlinear relationships among exogenous and (Supplementary Tables S9–S17).
endogenous variables; no nonlinearities were found.
To equalize variances we standardized variables by
dividing raw values by their group maximum.
Results
We constructed our SEM model stepwise. Our
initial multigroup model assumed that the effects of Aboveground, belowground, and interkingdom
smooth brome shoots and roots differed between predictors of microbial diversity
horizons. We separated roots and shoots as they Removing rare OTUs from each microbial matrix
likely influence nutrient availability differently resulted in a total of 82 archaeal (from 1339), 29
between soil horizons (Soil Classification Working bacterial (from 6747), and 119 fungal (from 4014)
Group, 1998). We then tested for horizon-specific OTUs used in subsequent analyses (Supplementary
effects of (1) pH on richness, evenness and phyloge- Figures S1–S3). The range of relative abundances
netic diversity of archaea, bacteria, and fungi, followingdownweightingchangedasfollows:archaea:
(2) archaea on bacteria and fungi, and (3) bacteria 0.001–1.170%to0.418–2.666%,bacteria:0.001–5.151%
on fungi. Improved model fit was assessed based on to 1.438–11.326%, and fungi: o0.001–2.607% to
reduced χ2 and Akaike Information Criterion. 0.067–3.424%.
There were no horizon-specific effects for pH on Using RDA, we found compositional differences
bacterial or fungal composition (Supplementary within fungi and bacteria between the A and B
Tables S6–S8). Finally, we fit alternative models for horizons; this variance appeared driven by
each bacterial OTU that exerted a significant influ- root biomass and SOC (Figures 2a and b,
ence on fungal diversity with relaxed parameter Supplementary Table S2). Archaea composition
constraints identical to those mentioned above. did not vary between soil horizons (Figure 2c).
InourfirstgroupofSEMs,wehadthreemetricsof We identified significant (Po0.05) belowground,
fungal diversity (richness, evenness, and phyloge- aboveground, and interkingdom correlations on
netic diversity), and evaluated the following four bacterialcomposition(Figure2a,Table1),represent-
separate causal assumptions: (1) no direct links ing 20% of the total variance in the bacterial matrix.
between smooth brome roots or litter to fungal These influences included litter nitrogen, root
diversity, (2) direct links between brome litter and biomass, SOC, and ammonia-oxidizing archaea
fungal diversity, (3) direct links between smooth abundance. Bacteria and archaea both influenced
bromerootsandfungaldiversity,and(4)directlinks fungal composition (although not vice versa as
between both smooth brome litter, and roots and inferred from the RDAs), along with several below-
fungal diversity. Thus, weran 12 separate models to ground and aboveground influences (archaeal and
test each of these potential causal assumptions (that bacterial richness, AOB abundance, root biomass,
is, 3×4; Supplementary Tables S6–S8). Once we SOC, and litter nitrogen); the constrained variance
found the best model configuration (as determined represented 11% of the total variance in fungal
byComparativeFitIndex,RootMeanSquareErrorof composition (Figure 2b). Of the potential explanatory
TheISMEJournal
Bacteriamediatefungiduringplantinvasion
SDMametetal
1269
variables, litter nitrogen was the only significant biplot origin (Figure 3a; Supplementary Table S3),
predictor of archaeal composition (Figure 2c; although this relationship appeared independent of
Table 1), although it represented only 2% of the smooth brome cover (4% explained variance,
variance. Microbial community compositions of P=0.509). The correlation between archaea and
archaea, bacteria, and fungi were similar within bacteria community composition was significant
microbial groups within a given brome cover class (RV=0.357, P=0.036) for the first and second axes
(PERMANOVA, brome cover class explained var- of the co-inertia plane (Figure 3b). The first axis
iance: 2% for each congener, Po0.023; Table 2). represented 30% of the variance and co-varied with
Bacterial and fungal communities were similar 19 archaeal OTUs. Nine of these OTUs were
within horizons (explained variance: 10% for bac- significantly distant (Po0.05) from the biplot origin
teria and 6% for fungi, Po0.001). (Supplementary Table S4) and all from the Candi-
datus phylum. There was not a clear separation
between invaded and uninvaded plots in terms of
smooth brome (2% explained variance, P=0.855).
Interkingdom soil microbial relationships
The association between archaea and fungi was
The association between bacteria and fungi was
significant (RV=0.610, P=0.030): the first axis
significant as indicated by the CoIA (RV=0.365,
represented12%ofthevariance,althoughsiteswere
Po0.001); the first axis represented 47% of the
differentiatedalongthesecondaxisbasedonsmooth
variance. Bacterial–fungal variance in the co-inertia
bromecover(22%explainedvariance,P=0.042)and
planewasdrivenmainlyby12bacterialOTUs.Sixof
theseOTUsinthephylaAcidobacteria,Actinobacteria,
andVerrucomicrobiawereparticularlyinfluentialas
indicated by significant distances (Po0.05) from the
Figure3 Co-inertiaanalysesofAhorizoninterkingdominterac-
tionsbetween(a)bacteriaandfungi,(b)archaeaandbacteria,and
(c)archaeaandfungialongasmoothbromeinvasiongradient.The
largeboxescontainjointsiteplotscoloredbysmoothbromecover, Figure4 Co-inertiaanalysesofAhorizonroot-kingdominterac-
where arrowslink kingdom diversity. Top left, clockwise: arrow tionsalongasmoothbromeinvasiongradient.Thelargepanelsare
tails represent bacteria and the heads represent fungi; tails jointsiteplotscoloredbysmoothbromecover,wherearrowslink
represent archaea, heads represent bacteria; tails represent root-kingdom diversity. (a) Arrow tails represent roots, and the
archaea, heads represent fungi. Arrow length indicates the headsrepresentbacteria.(b)Arrowtailsrepresentroots,andheads
strengthoftheassociationbetweenmatrices:shortarrowsindicate represent fungi. Arrow length indicates the strength of the
strongconcordance;longarrows,weakconcordance.Smallboxes: association between matrices: short arrows indicate strong con-
projection of the kingdom of interest onto the co-inertia plane. cordance;longarrows,weakconcordance.Smallboxes:projection
BacterialandarchaealOTUswithsignificantloadings(P⩽0.05)on of roots, fungi, or bacteria of interest onto the co-inertia plane.
the smooth brome gradient are indicated by bold red arrows Boldredarrowsindicaterootspecieswithconsiderable(P⩽0.05)
(SupplementaryTablesS3andS4). loadingsonthesmoothbromegradient(SupplementaryTableS5).
TheISMEJournal
Bacteriamediatefungiduringplantinvasion
SDMametetal
1270
29archaealOTUs(ninesignificantlydistantfromthe Root influence on soil microbial composition
biplot origin (Po0.05)), and represented 8% of the Roots were significantly linked to both bacteria and
variance (Figure 3c; Supplementary Table S4). fungiviabiomassandSOC(Figures2aandb).Inthe
root–bacteriaco-inertiaanalyses,thefirstandsecond
axes together represented 47% of the variance, with
clear separation between smooth brome-invaded
and non-invaded sites (20% explained variance,
P=0.022; Figure 4a); however, the concordance
was insignificant (RV=0.273, P=0.105). Plots were
differentiated into invaded and non-invaded by
roots from 13 plant species (Figure 4a). There
was significant concordance (RV=0.472, P=0.049)
between roots and fungi; the first axis represented
15% of the variance and 10 root species differen-
tiated sites (Figure 4b; Supplementary Table S5).
Smooth brome clearly separated invaded and non-
invaded sites (29% explained variance, P=0.010).
Smooth brome reduced the abundance of several
nativespecies(Bouteloua gracilis,Carexduriuscula,
C. filifolia, Elymus sp., and Galium boreale), which
appeared to be linked to bacterial and fungal
community composition (Figure 4b; Supplementary
Table S5).
Effects of smooth brome abundance on soil microbial
structure
Although including direct pathways between
smooth brome and fungi diversity resulted in better
fit for the SEMs, as indicated by reduced χ2-values
Figure5 MultigroupSEMsofsmoothbromeandpHinfluences
on bacterial and fungal diversity (richness, evenness, and
phylogenetic diversity) and interkingdom interactions. (a–c)
Influential abundant taxa identified through downweighting,
RDA, and CoIA. (d, e) Influential rare taxa identified through
Poisson regression and RDA. (a) Increased archaeal diversity
increases abundance of four bacterial OTUs. (b) Smooth brome
shoots reduced abundance of two bacterial OTUs. (c) Smooth
brome shoots and roots increase bacterial diversity resulting in
decreased fungal diversity. (d) Smooth brome roots increased
abundanceofOTU5113.(e)pHhadaninverseeffectonOTU2472
(that is, alkaline pH, reduced abundance; acidic pH, increased
abundance)andthatOTUreducedfungaldiversity.Solidarrows
represent significant relationships (Po0.05, dashed lines are
nonsignificant), the thickness of the arrow indicates the strength
of the relationship and the color indicates the direction of the
relationship (green=positive, red=negative). Standardized path
coefficientsareshownnexttopathways.Pathcoefficientscanbe
interpreted as this: if, in a for example, brome shoot abundance
increasesbyones.d.fromthemean,thenfungaldiversitywould
increaseby0.027s.d.'sfromitsownmean.Thus,thegreaterthe
path coefficient, the larger the effect. Where parameters differ
betweenhorizons,twopathcoefficientsareshown:thefirstisthe
coefficient for the A horizon and the second for the B horizon.
All models achieved adequate fit (Supplementary Table S20).
Significant path coefficients between smooth brome roots and
shoots, bacterial parameters on fungal diversity and archaeal
diversity on bacterial parameters are illustrated in the inset bar
graph(SupplementaryTablesS9–S17).Barswithaheightof0.025
indicatenonsignificantpaths.Fordetailsoneachmodelclass,see
Supplementary Information (fungal evenness: Supplementary
Table S10, fungal richness: Supplementary Table S9, fungal
phylogenetic diversity: Supplementary Table S11, bacterial
OTU 5113: Supplementary Table S18, bacterial OTU 2472:
SupplementaryTableS19).
TheISMEJournal