Data collection
We collected samples within each of the 50 plots described above
underneath bare (unvegetated) soil, a grass individual (Bouteloua
eriopoda), and a shrub individual (Prosopis glandulosa) in July
2013 and 2019, resulting in three samples per plot per year, and
300 samples total (3 host types * 2 years * 5 precipitation
variation levels * 10 replications). Soil samples were collected
with soil corers and trowels to 5 cm depth from the surface within
the soil rooting zone of the focal plant. We sprayed equipment
down with 10% bleach between samples to minimize contamination.
Samples were stored on ice in a cooler in the field and
transferred to storage a -20°C freezer until August 2019, when
they were transported to the University of Georgia and stored at
-80°C for downstream applications.
DNA extraction and sequencing
To characterize soil fungal communities, we subsampled 250 mg of
soil from each sample and extracted total genomic DNA using the
QIAGEN Dneasy® PowerSoil Pro Kit®, following the manufacturer’s
protocol (QIAGEN, Hilden, Germany). Prior to DNA extraction, we
performed tissue lysis with a Spex® Genogrinder® at 1500 rpm for
10 minutes. We quantified DNA spectrophotometrically using a
Nanodrop One/One© (ND2000; NanoDrop Technologies). Genomic DNA was
stored at -20ºC prior to library preparation and sequencing.
We characterized fungi with high-throughput sequencing to detect
the variation in marker gene differences. To target fungal taxa,
we amplified the ITS2 region using primers fITS7
(5′-GTGARTCATCGAATCTTG-3′; Ihrmark et al., 2012), and ITS4
(5′-TCCTCCGCTTATTGATATGC-3′; White et al., 1990), augmented with
multiplexing barcodes. Samples were sequenced on two runs of an
Illumina MiSeq instrument (300PE V3 chemistry). Library
preparation and DNA sequencing was conducted at the Georgia
Genomics and Bioinformatics Core (GGBC) at the University of
Georgia in Athens, GA.
Bioinformatic processing
Forward and reverse reads were paired, merged, and quality
filtered with USEARCH (http://drive5.com/usearch/; Edgar, 2010,
2013). Sequences were then processessed using Mothur (Schloss et
al., 2009) to remove sequencing adapters, primers, and chimeras.
Combined reads were clustered into Operational Taxonomic Units
(OTUs) at 97% identity using the K-nearest algorithm (Wang et al.,
2007), which we classified against the UNITE (Version 8.0)
reference database (Nilsson et al., 2019). All non-fungal reads
were excluded from the analysis. The total dataset included 1308
unique fungal OTUs with 10,370,208 total sequences. Due to uneven
sequence depth, we rarefied samples to 3000 reads per sample for
all downstream analyses. Samples with less than 3000 reads were
excluded from the analysis (27 of 297) due to insufficient
sequencing depth (Figure S1).
Statistical analysis
All statistical analyses were performed in R (Version 1.2.1335, R
core team, 2017) and conducted on the rarefied OTU abundance
matrix. To evaluate the influence of increased precipitation
variability on diversity metrics, we calculated the coefficient of
variation of the amount of precipitation received for each
variability treatment from 2009 to 2013 and from 2014 to 2019,
such that precipitation variability can be analyzed as a
continuous variable (see Gherardi and Sala, 2015b). Thus, for the
remainder of this paper, we refer to increased precipitation
variability as the coefficient of variation (CV) of the amounts of
growing season precipitation received for the five precipitation
variability treatments calculated from the periods between
2009-2013 and 2014-2019, respectively. For for the indicator
species and functional guild composition anlyses, we used
increased precipitation variability as a categorical variable for
three levels of increased precipitation variability: Control,
medium (-/+50% and +/-50%) and high (-/+80% and +/-80%). We
analyzed the influence of host type, precipitation variability,
and sampling year on soil fungal diversity using a linear mixed
model (Bates et al., 2015) with Shannon’s diversity index for each
sample as a response variable, and host type, precipitation (CV),
year and their interactions as fixed effects, and plot as a random
intercept accounting for the repeated measures nature of the
experiment. We included sampling year as a categorical variable in
this and subsequent analyses to account for different ambient
conditions between sampling years 2013 and 2019 (Figure 1). We
performed post-hoc multiple pairwise comparisons using the false
discovery rate (fdr) method to correct p-values (package
'emmeans', Lenth, 2018). To further explain differences observed
in Shannon diversity, we also separately analyzed fungal richness
(number of OTUs in each sample) and evenness (distribution of
relative abundance with the Pielou’s metric of evenness) using the
same model structure and post-hoc methods as above for Shannon
diversity. To investigate soil fungal community composition, we
calculated Bray-Curtis dissimilarities among samples using the
rarefied OTU abundance matrices. Using Bray-Curtis dissimilarity
matrices, we tested the effect of host type, precipation
variability, and year on fungal community composition with a
permutational analysis of variance (PERMANOVA) (function adonis in
the 'vegan' package) (Dixon, 2003). We further identified soil
taxa that were significant indicators of treatment groups using
indicator species analysis (Dufrêne and Legendre, 1997) with the
'labdsv' package. To visualize variation in fungal community
composition, we performed non-metric multidimensional scaling
(NMDS) ordinations using Bray-Curtis dissimilarities. To assign
functional guilds, or primary lifestyles to operational taxonomic
units, we used only OTUs assigned to genus-level resolution (656
out of 1308 OTUs) and assigned genera to functional identity using
the FungalTrait database (Põlme et al., 2020).
To test whether soil fungal diversity responses were due to direct
effects of increased precipitation variability or indirect effects
via host plant induced changes, we performed stuctural equational
modelling (SEM) using the Lavaan Package (Rosseel, 2012). Our
model incorporated the direct correlation between increased
precipitation variability and soil fungal diversity, and the
direct correlation between increased precipitation variability and
percentage cover of the host type. The indirect path was modelled
as the association between increased precipitation variability and
soil fungal diversity, mediated through a change in perentage host
type cover. Models were just-identified and therefore no general
model fit was reported. We split up the OTU table by the three
host types to investigate above-ground induced changes in each
host-associated community type separately. We also modeled two
different soil fungal responses: richness and evenness. We
specifically modelled richness and evennes as separate responses —
as apposed to only modelling Shannon Diversity as a response. This
led to a total of six (host type by fungal response) models
evaluated.
