Arbuscular mycorrhizal fungi (AMF) facilitate ecosystem functioning through provision of plant hosts with phosphorus (P), especially where soil P is limiting. Changes in soil nutrient regimes are expected to impact AMF, but the direction of the impact may depend on context. We predicted that nitrogen (N)-only enrichment promotes plant invasions and exacerbates their P limitation, increasing the utility of AMF and promoting AMF diversity. We expected that enrichment with N, P and other nutrients similarly promotes plant invasions, but decreases the benefit and diversity of AMF because P is readily available for both native and exotic plants. We tested these hypotheses in eucalypt woodlands of south-western Australia, that occur on soils naturally low in P. We evaluated AMF communities within three modified ground-layer states representing different types of nutrient enrichment and associated plant invasions. We compared these modified states to near-natural reference woodlands. AMF richness varied across ground-layer states. The moderately invaded/N-enriched state showed the highest AMF richness, while the highly invaded/NP-enriched state showed the lowest AMF richness. The reference state and the weakly invaded/enriched state were intermediate. AMF richness and colonisation were higher in roots of exotic than native plant species. AMF community composition differed among ground-layer states, with the highly invaded/NP-enriched state being most distinct. Distinctions among states were often driven by family-level patterns. Reference and moderately invaded/N-enriched states each supported distinct groups of zero-radius operational taxonomic units (zOTUs) in Acaulosporaceae, Gigasporaceae and Glomeraceae, whereas Gigasporaceae and Glomeraceae were nearly absent from the highly invaded/NP-enriched state. Further, Diversisporaceae and Glomeraceae were most diverse in the moderately invaded/N-enriched state.

 Synthesis. Both the nature of soil nutrient enrichment and plant provenance matter for AMF. N-only enrichment of low-P soils increased AMF richness, likely due to introduction of AMF-dependent exotic plant species and exacerbation of their P-limitation. In contrast, multi-nutrient enrichment, decreased AMF richness potentially due to a decrease in host dependence on AMF, regardless of host provenance. The changes in AMF community composition with nutrient enrichment and plant invasion warrants further research into predicting the functional implications of these changes.

Site selection

Sampling was centred on two nature reserves with areas of York gum (Eucalyptus loxophleba subsp. loxophleba) – Jam (Acacia acuminata) woodlands (hereafter York gum woodlands) in close to reference condition. The reserves are located in the central wheatbelt, Western Australia: Mount Caroline (31°45'25.3"S, 117°38'38.3"E) and Namelkatchem Nature Reserves (31°10’47.9”S, 117°11’18.1”E), and are ~ 70 km apart (Fig. 1). These reserves have a history of minimal livestock grazing and plant invasion, permitting the persistence of large areas of diverse understoreys dominated by native perennial and annual forbs and grasses (Prober and Wiehl 2012). The climate is Mediterranean type, with long-term (1990–2022) mean annual temperature and rainfall, respectively, of 17.8 °C and 321 mm at Mt Caroline and 17.6 °C and 333 mm at Namelkatchem (BoM, 2023).

Our experimental design involved four ground-layer states representing different levels of nutrient enrichment and plant invasions as described above: a reference (control) state in near natural conditions with naturally low soil P and N, and three degraded states that we sought from different parts of the landscape: a weakly invaded/enriched state, a moderately invaded/N-enriched, and a highly invaded/NP-enriched state (Table 1; Fig. 1). We chose the term ‘ground-layer state’ to be consistent with previous research of this ecosystem (Prober et al. 2012). ‘Ground-layer’ is used rather than ‘understorey’ to distinguish the herbaceous ground-layer from a shrub layer that can occur in the woodland understorey.

Two areas within each of the two nature reserves with few-to-no exotic plant species, were selected to sample the ‘reference’ state. Each reserve included localised patches invaded by exotic annuals, that likely arose historically due to localised disturbances (e.g. by introduced rabbits). These patches were used to represent the ‘moderately invaded/N-enriched state’, given that other studies have demonstrated that such exotic-invaded states of York gum woodlands are typically N-enriched (Prober and Wiehl, 2012).

Because P-enriched areas rarely occur inside the nature reserves, to represent the ‘highly invaded/NP-enriched’ state we sampled four fertilised 2 m × 2 m plots from a pre-existing long-term nutrient addition experiment located ~1.5 km from Mount Caroline Nature reserve (31°46'56.43"S, 117°36'41.61"E; Fig. 1). The experiment was established in a grazed York gum woodland remnant in 2009, as part of the global Nutrient Network (NutNet) experiment (Borer et al. 2014). We also sampled four unfertilised plots from the same experiment, representing the ‘weakly invaded/enriched’ state that arose through historical sheep grazing, resulting in some plant invasion and N and P enrichment. To match the NutNet design, plots of 2 m × 2 m were established at both nature reserves for the reference and moderately invaded/N-enriched states. Plots were at least 20 m apart, with four replicate plots of each ground-layer state. Remnants of York gum woodlands are rare in the landscape; hence we chose the two closest nature reserves near the NutNet experiment. This was done to incorporate as much of the natural variation within Reference sites as possible. We note that a previous study found little variation in AMF communities among remnants of York gum woodlands across 200 km distance (Prober et al. 2015).

The experimental plots from NutNet were arranged in a randomised complete block design with four blocks. Nutrients were added annually in autumn to plots representing the highly invaded/NP-enriched state, as 10 g N per m²·yr⁻¹ of timed-release urea, 10 g P per m²·yr⁻¹ as triple superphosphate, and 10 g K per m²·y⁻¹ as potassium sulphate. These plots also received a once-off addition of other macro- and micronutrients in 2009: 100 g per m⁻² of a mix containing iron (15%), sulphur (14%), magnesium (1.5%), manganese (2.5%), copper (1%), zinc (1%), boron (0.2%), and molybdenum (0.05%). Weakly invaded/enriched plots were dominated by native plant species and highly invaded/NP-enriched plots were dominated by exotic plant species. The experimental plots had been open to livestock grazing since European colonisation (1860s) until 2015.

Sample collection and processing      

Sampling occurred in August 2021 during the growing season. In each plot, plant community composition and cover were recorded. Then, rhizosphere soil and roots were collected from two plants of the six most abundant plant species. Samples from the two plants per species were pooled for a total of six samples per plot. Soils were thoroughly mixed in a sealable bag, collecting a subsample for DNA analyses. Roots were immediately stored in 98% ethanol pending processing. Soil and root samples were stored in 15 ml tubes and placed in dry ice during sampling and transport to the laboratory. Samples were stored at -80°C thereafter. Root samples were thoroughly cleaned with deionised water, and fine roots of < 2 mm were retained. Clean root samples were split in two: one part for DNA analyses and one for root colonisation assessment. Two plots had only five plant species, resulting in 94 soil and 94 root samples (Table S1). No permit was needed to access and sample plots from the NutNet experiment. Permits to access and sample the two nature reserves were granted by the Department of Biodiversity, Conservation and Attractions of Western Australia (Licence: FT61000839; Regulation 4: CE006388).

Soil chemical analyses

Soil samples were sent to CSBP Laboratories (Bibra Lake, Western Australia) for nutrient analyses. Plant-available P and K were measured using the Colwell test (Colwell, 1963; Rayment & Higginson, 1992). Organic carbon (OC) was determined according to Walkley & Black (1934). Ammonium-N, nitrate-N, and total N were measured as per Searle (1984). Soil pH was measured in CaCl₂ in a solution ratio of 1:5 (Rayment & Lyons, 2012).

Root colonisation

Root subsamples allocated for measurement of root colonisation were cleared in 1 M KOH and stained with ink in vinegar (5% v/v) as described by Vierheilig et al. (1998). Colonization by AMF, including Glomeromycotina-AMF and Mucoromycotina-AMF, was scored using the line intercept method (McGonigle et al. 1990). One hundred intercept points were scored for each sample, and the percentage of root length colonized by AMF was calculated.

DNA extraction and sequencing

Root samples allocated for DNA analyses were cut into 5 mm pieces, homogenised, and ground with beads. DNA was extracted from 20 mg of root and 250 mg of soil material using the DNeasy Plant Pro kit and DNease PowerSoil Pro kit, respectively, (Qiagen, Carlsbad, USA). PCR amplification and sequencing was performed by the Australian Genome Research Facility. For each sample, 15 ng DNA were used to amplify the 18S rRNA gene using the AMF primer set AMV4.5NF and AMDGR (Sato et al., 2005). These primers accurately retrieve a wide range of AMF taxa, including both Glomeromycotina and Mucoromycotina subphyla (Orchard et al., 2017; Albornoz et al., 2022). Thermocycling was completed with an Applied Biosystem 384 Veriti and using Platinum SuperFi II mastermix (Life Technologies, Australia) for the primary PCR. Thermocycling consisted of an initial denaturation at 98℃ for 30 s followed by 30 cycles of 98℃ for 10 s, 60℃ for 10 s and 72℃ for 30 s. The final extension was at 72℃ for 5 min. The first stage PCR was cleaned using magnetic beads, and samples were visualised on 2% Sybr Egel (Thermo-Fisher). A secondary PCR to index the amplicons was performed with the same conditions with Platinum SuperFi II mastermix (Life Technologies, Australia). The resulting amplicons were cleaned using magnetic beads, quantified by fluorometry (Promega Quantifluor) and normalised. The eqimolar pool was cleaned a final time using magnetic beads to concentrate the pool and then measured using a High-Sensitivity D1000 Tape on an Agilent 2200 TapeStation. The pool was diluted to 5nM and molarity was confirmed using a Qubit High Sensitivity dsDNA assay (ThermoFisher). This was followed by sequencing on an Illumina MiSeq (San Diego, CA, USA) with a V3, 500 cycle kit (2 x 250 base pairs paired-end).

Bioinformatics

Following sequencing, adapters were trimmed using ‘cutadapt’ (Martin, 2011), retaining only sequences that contained primers. After trimming, further quality checks and sequence processing were done in VSEARCH v2.14.1 with default parameters (Rognes et al., 2016). Trimmed paired-end sequences were merged with a minimum of 10 bp overlap. Merged sequences were filtered using a maximum error rate of 0.1 and a minimum length of 200 bp. Filtered sequences were dereplicated at 100% identity and singletons were discarded. Unique sequences were clustered into zero-radius operational taxonomic units (zOTUs), chimeras were detected (denovo), and a zOTU abundance table was produced using the UNOISE3 algorithm in USEARCH v11 (Edgar, 2016). Finally, zOTUs were queried against the SILVA SSUref v138.1 database (Quast et al., 2013). Taxonomy was assigned to zOTUs with a threshold of > 95% match and query cover of > 90%. Sequences matching Glomeromycotina or Mucoromycotina-AMF were classified as AMF.

Statistical analyses

Two samples failed to amplify DNA, resulting in 186 samples for statistical analyses. The initial denoised zOTU abundance table was rarefied to the smallest sequencing depth (4,927 sequences) to avoid sequencing depth bias (Dickie, 2010) (Fig. S1).

All data were analysed, and figures were created, in R (R Core Team, 2016). To visualise variation in plant and AMF communities among ground-layer state, non-metric multidimensional scaling with Bray Curtis dissimilarity of the log-transformed matrix was used. To visualise species contributions to the ordination, two-way tables were constructed with the ‘inkspot’ function from the rioja package (Juggins and Juggins, 2020). Replicate plots were ordered on the x-axis by ground-layer states, while each individual zOTU was placed on the y-axis. Two-way tables summarise and present the raw community data as a powerful way of visualising the distribution of taxonomic units along environmental, spatial, or temporal gradients.

To test for differences in community compositions among ground-layer states and source material (i.e., soil vs roots), permutational multivariate analysis of variance were performed using ‘adonis2’ within the vegan package (Oksanen et al., 2017). To test which vegetation and soil variables correlated with plant and AMF communities among ground-layer states, vectors of maximum correlation were calculated with the vector-fitting procedure using ‘envfit’ (9,999 permutations) within the vegan package (Oksanen et al., 2017). A biplot was drawn on the ordination to display the relationships between explanatory variables and ordination axes. Because floristic variables strongly covaried (Fig. S2), principal component analyses (PCA) were used to distil covariates to PCA axes (Fig. S3). An ‘exoticPCA1’ (i.e., first axis of a PCA between exotic plant cover and richness) and a ‘nativePCA1’ (i.e., first axis of a PCA between native plant cover and richness) indices were created to represent the difference in cover and richness for exotic and native plants, respectively.

Differences in plant cover, plant richness, rarefied zOTU AMF richness (hereafter ‘AMF richness’), and soil chemistry among ground-layer states were evaluated. These differences were analysed using linear mixed effect models using ‘block’ as random effect. To test how environmental attributes associated with transition among ground-layer states related to AMF richness and community composition, structural equation models (SEM) were built using the lavaan R package (Rosseel, 2012). To meet lavaan’s requirements, ground-layer state was transformed to an ordinal variable based on exotic plant abundance and levels of nutrient enrichment (reference = 0; weakly invaded/enriched = 1; moderately invaded/N-enriched = 2; highly invaded/NP-enriched = 3). To include a proxy of community composition in the SEMs, the first axis of a PCA of the log-transformed AMF community matrix was used. Due to high covariation among most soil variables (Fig. S2), models were simplified to only include variables with a priori knowledge of being involved in state transitions in eucalypt woodlands. These variables were OC, nitrate-N, ammonium-N, and available P, which are hypothesised to promote exotic annuals and AMF richness (Prober and Wiehl, 2012; Prober et al., 2014; Albornoz et al., 2023). Neither nitrate-N nor ‘nativePCA1’ were selected in the best model.

We acknowledge that feedbacks between soil chemistry and plants are likely. However, in our SEMs, we chose soil chemistry as the driver for exoticPCA1 for two reasons: first, there is a priori knowledge that N enrichment promotes plant invasions in these ecosystems (Prober et al. 2012). Second, the highly invaded/NP-enriched state was created by deliberate nutrient addition, meaning it was demonstrated experimentally that nutrient addition promoted the observed plant invasions, not the other way around.

Comparisons of AMF richness and communities between native and exotic plant species across treatments (i.e., Plant provenance × Ground-layer state) were not possible due to the absence of sufficient native plant species in the highly invaded/NP-enriched state, and exotic plant species from the reference state. Further, for most samples, no root material was left after prioritising the subsample for DNA analyses, leaving 31 (32%) samples for scoring colonisation. Hence, AMF richness and root colonisation comparisons between native and exotic plants are presented using ground-layer state as a random effect.