Exotic plant invasions alter native plant productivity and ecosystem function. Mechanisms of plant invasions can include altered disturbance regimes or altered plant-soil feedbacks. Some important hypotheses regarding the invasion of cheatgrass in North America include pathogen spillover and mutualism disruption. Since an important mutualistic interaction in ecosystems is that of mycorrhizal fungi, here we test the mutualism disruption hypothesis by closely examining mycorrhizal allocation and community composition. We conducted a greenhouse study where we grew a native perennial grass in association with invaded and uninvaded soils. We found altered mycorrhizal communities and mycorrhizal allocation patterns associated with cheatgrass-invaded soils supporting the mutualism disruption hypothesis.

Sources of plants, soil, and inoculum 

Seeds and soil were collected from the west side of the Kaibab Plateau (Coconino County, Arizona, USA) at an elevation of 2,064 m with approximately 43 cm of precipitation annually (PRISM Climate Group). The site was a semi-arid piñon-juniper woodland that burned in 1996 and had a mix of intact native grass communities as well as areas invaded by cheatgrass. Invaded patches were relatively small but were generally monocultures of cheatgrass whereas uninvaded patches were biodiverse with blue gramma as a dominant species.

Blue gramma seed was collected from the uninvaded portions of the site using the Seeds of Success protocol (http://www.nps.gov/planTs/sos/protocol/index.htm). Live soil inoculum was collected every 10 m from the rooting zone of blue gramma along three 100 m transects established from a random origin (azimuths of 0˚, 90˚ and 270˚) in the uninvaded portion of the site and from the rooting zone of cheatgrass (using similar transects) in the invaded portion of the site. Soil subsamples within each collection were pooled together and mixed. We justify homogenizing inoculum from each sampling location because we were interested in seedling responses to average soil biotic conditions across sites, rather than within a single site or extrapolating to a broader geography than our sampling sites (a “type C” design; Gundale et al. 2017, 2019). Inoculum soil was refrigerated for two weeks until its use in the experiment. The abundance of different soil organisms in the two inoculum soils was determined using phospholipid fatty acid (PLFA) and neutral lipid fatty acid (NLFA) analysis. Lipids were extracted from 5 g of freeze-dried inoculum soil by vortex mixing in a one-phase mixture of citrate buffer, methanol, and chloroform (0.8:2:1: v/v/v, pH 4.0). The biomass of AM fungi was estimated from the PLFA 16:1 w5, 20:1 w9, and 22:1 w13, biomass of other fungi was estimated from 18:2 w9:12, and bacterial biomass was estimated from signature PLFAs for gram-positive and gram-negative bacteria (Olsson et al., 1995). This analysis indicated that the initial soil inoculum from the uninvaded and invaded sites had similar overall abundances of AM fungi, other fungi, and bacteria (Supporting Information Table S1).    

Experimental design

Mesocosms were constructed from 21 L plastic containers (43 cm x 28 cm x 18 cm) with six 0.3 cm diameter holes drilled into the bottom for drainage. To remove the effects of any variation in soil physical and chemical characteristics in the cheatgrass-invaded and uninvaded soils, we created a sterilized common soil using a 1:1 mixture of the two soils that was steam-sterilized at 125°C for 48 hours. Our experimental design matches type C in Gundale et al. (2017) because unique and variable sub-populations of plant subjects (a random draw of seeds collected from a site) are confronted with one of two soil biota assemblages. This design is preferred when the goal is to detect differences among two or more groups of subjects, and when within-site or regional spatial variation is not a focus (Gundale et al., 2017, 2019). Each mesocosm was filled with approximately 15 liters of sterilized soil and topped with a 1 cm thick band of either live or sterilized (dead) inoculum soil that had been steam sterilized twice at 125°C for 24 hours. Mesocosms were prepared with three inoculum treatments: live soil from uninvaded areas, live soil from cheatgrass-invaded areas and sterilized soil. Each inoculum treatment was replicated 10 times, resulting in 30 mesocosms. Blue gramma seed was sprinkled onto the inoculum soil at a rate of 60 seeds per mesocosm and later thinned to 10 seedlings per mesocosm. Mesocosms were placed in fully randomized locations (to account for microclimatic variation within the glasshouse) and watered three times per week.

Analysis of soil biota assemblages

The composition of soil biota in cheatgrass invaded and uninvaded inoculum treatments were compared after the experiment. Samples of soil were collected from the 30 mesocosms and DNA was extracted from 0.5 g of soil using a PowerSoil DNA Extraction Kit (MO BIO Laboratories, Inc., Carlsbad, CA, USA). Genomic DNA was normalized to 2 ng/mL, diluted 10-fold, and amplified in triplicate PCR using the universal eukaryotic primer WANDA and the AM fungal specific primer AML2 for the small subunit (SSU) rRNA gene (Lee et al., 2008; Dumbrell et al., 2011). Purified products were quantified with PicoGreen fluorescence. Indexing PCR was completed using 8 bp dual indexed WANDA and AML2 primers. Indexed PCR products were purified using a 1,1 carboxylated magnetic bead solution, quantified, and combined into a final sample library. The library was purified, concentrated, and quantified using quantitative PCR against Illumina DNA standards on an Illumina MiSeq System (Illumina, Inc., San Diego, CA) running in paired end 2 x 300 bp mode. Forward reads were trimmed to 250 bp to remove low quality tails and demultiplexing was carried out using a minimum quality threshold of q20 and default parameters in QIIME 2.0 (Caporaso et al., 2010). Taxonomy was assigned to sequences using BLAST with 90% similarity and an E-value less than 10^-4, against the online MaarjAM database (http,//maarjam.botany.ut.ee; accessed 06/20/2023, Öpik et al., 2010). Bioinformatics for 16S and ITS sequences were performed using QIIME 2 version 2021.11 (Bolyen et al. 2019). Forward and reverse 16S reads were joined after trimming to a length of 240. Due to low quality reverse reads, only forward ITS reads were used, which were trimmed to a length of 222. Denoising was performed using DADA2 (Callahan et al. 2016). Taxonomic assignment for 16S and ITS was performed using Emu version 3.4.3 (https://gitlab.com/treangenlab/emu; Curry et al.2022) using NCBI and UNITE, respectively, as taxonomic references. Taxa that made up less than 1% of relative abundance were labeled as ‘other’, otherwise species were recorded to the genus level for community comparisons. Many species remained unidentified or classified only to order or family.

Plant and fungal variables

After eight months, all aboveground biomass was clipped, dried at 60°C for 24 hours, and weighed. Root biomass was sampled by taking four soil cores (5 cm diameter and 18 cm deep). Roots were cleaned, dried, and weighed and the weight of roots per volume of core was used to estimate root biomass in the total volume of the mesocosm. A 10 g subsample of fresh root material was refrigerated until it was cleared with 5% KOH and stained with ink in vinegar (Vierheilig et al., 1998). Colonization by AM fungi and other root endophytes was determined using the gridline intersect method at 200 × magnification (McGonigle et al., 1990). Mycorrhizal root colonization was distinguished as arbuscules, vesicles and hyphae. Colonization by non-AM fungi, including dark septate endophytes (DSEs) was also quantified. External hyphae of AM fungi were extracted from the soil cores after root removal, using the methods of Sylvia (1992), and quantified using a gridded eyepiece graticule in an inverse compound microscope at 250 × magnification.  At points where hyphae intersected gridlines, hyphae were counted, and counts were converted to length of hyphae per g of soil. Hyphae of AM fungi were distinguished from other fungal hyphae based on their morphology and color.

Statistical analysis

The soil biota effect was calculated using a pooled standard deviation for the treatments, following (Kulmatiski et al.2008)to quantify plant biomass responses to AM fungi and other soil organisms relative to plants grown in the absence of living inoculum. The total biomass of blue gramma in each mesocosm inoculated with living cheatgrass invaded or uninvaded inoculum was compared to the average total blue gramma biomass of mesocosms with sterile inoculum.

Soil biota effect = (X_treatment - M_sterile) / SD_sterile

Where X is the value of an individual experimental mesocosm, M is the mean of the sterile treatment group and SD is the standard deviation of the sterile treatment group.  

One-way ANOVA was used to compare the effect of the inoculum origin on final plant biomass, soil biota effect, density of external AM hyphae, and percent root length colonized by AM fungi and DSEs. Sterile controls were excluded from all ANOVA models. Linear regression was used to determine relationships between soil biota effect and density of external AM hyphae, and percent root length colonized by different AM fungal structures and DSEs. Model assumptions were checked using the Shapiro-Wilk test of normality and the Levene’s test of heterogeneity of variance. All statistics were conducted in R (version 3.3.1). A permutation analysis of variance (PERMANOVA) was performed on each microbial community using Bray-Curtis dissimilarities to determine if cheatgrass invasion resulted in significantly different assemblages of bacteria, AM fungi and non-AM fungi. We used the permanova function from skbio (version 0.5.6; scikit-bio.org) and the pairwise_distances from sklearn (version 1.3.0; Pedregosa et al., 2011) to determine differences. Indicator species in mesocosms inoculated with cheatgrass invaded or uninvaded soils were identified using a random forest classification using the RandomForestClassifier function from sklearn (Pedregosa et al.2011). We averaged three training and test data splits to ensure there was no variation due to the randomness of splitting the data. We used GridSearchCV from sklearn.model_selection to perform a grid search of the hyperparameter space and identify the optimal parameters for our random forest model with 3-fold cross-validation to avoid overfitting. The grid search explored combinations of the following hyperparameters: the number of trees in the forest (n_estimators : [50, 100, 200]), the estimation of maximum features method (max_features: ['sqrt', 'log2']), the maximum depth of the trees (max_depth: [10, 20, 30, None], the minimum number of samples required to split an internal node (min_samples_split: [2, 5, 10]), and the minimum number of samples required to be at a leaf node (min_samples_leaf: [1, 2, 4]), and whether to use bootstrapping (bootstrap: [True, False]).

Changes after Dec 1, 2023: Uploaded the correct data sheet that contains the relevant abundances of organisms and deleted the other data files that contained redundent information.