Greater change in arbuscular mycorrhizal fungal richness as a response to short-term rainfall exclusion across the North American monsoon season

As climate change alters seasonal weather patterns, interactions among organisms can become decoupled through range shifts or changing phenology (Van Der Putten et al., 2010). Arid regions around the world are predicted to experience stark reductions in seasonal precipitation while global temperatures begin to rise (Diffenbaugh and Giorgi, 2012). The North American Monsoon (NAM) season, present in Southeastern Arizona and Western New Mexico, provides roughly 40% of the regions mean annual precipitation (Mitchell et al., 2002). This seasonal effect has been shifting spatiotemporally and becoming difficult to predict under future climate conditions (Comrie and Glenn, 1999; Demaria et al., 2019). Climate models have illustrated the importance seasonal precipitation and temperature have on plants, animals, and, thus, the interactions among them (Bradie and Leung, 2017). Plants also have numerous belowground interactions with microorganisms, such as bacteria and fungi, whereby mutualistic associations provide vital access to nutrients, water, as well as resistance to pathogens (Newsham et al., 1995). For example, these associations promote plant productivity, affect plant defense strategies, and influence overall health (Van Der Heijden et al., 2008; Lau and Lennon, 2011). Yet, there remains a need to understand how changing seasonal precipitation regimes in dry, arid climates might affect plant–fungal interactions.

Arbuscular mycorrhizal (AM) fungi are a clade of plant host-associated microorganisms (Glomeromycotina) that form associations with over 70% of terrestrial plant species within the rhizosphere microhabitat (Soudzilovskaia et al., 2020). AM fungi are highly abundant in soils and have a wide distribution globally, extending from tropical forests to grasslands (Hayman, 1982). The communities of AM fungi are the key components of soil nutrient cycles, exchanging phosphorous, nitrogen, and other trace minerals with plant-derived carbohydrates and lipids (Allen, 1991). At a regional scale, AM fungal communities and the type of symbiosis (e.g., parasitism, commensalism, mutualism) are often more strongly predicted by environmental factors than plant hosts or dispersal limitation (Johnson et al., 2010; Kivlin et al., 2014). These environmental factors could be one or more of the following: soil fertility (soil organic matter content; Anderson et al., 1984), soil moisture (Miller, 2000; Gao et al., 2016), soil pH (Green et al., 1976), soil texture (e.g., silt, sand, etc.; Johnson et al., 1992), and soil nutrients (Johnson et al., 2010). The changes in soil moisture or H₂O saturation can also influence the labile nature of soil abiotic factors of the rhizosphere (Deepika and Kothamasi, 2015; Kanwal et al., 2015; Bobille et al., 2019). Therefore, it is important to investigate the relative importance perturbations in soil moisture have on AM fungal communities through both direct and indirect pathways.

Environmental stress can also have large direct consequences on AM fungal community composition and relative abundance (Hernandez et al., 2021). Stress can range from increased seasonal temperatures (Zhou et al., 2016) to changes in precipitation regimes (Waring and Hawkes, 2015). AM fungi have been shown to be exist in a wider range of annual precipitation than other mycorrhizal fungi (Tedersoo et al., 2014), and thus, AM fungi have been regarded as key regulators of plant stress brought on by drought (Birhane et al., 2012). For example, in greenhouse experiments, seasonal drought was induced through water pulsing, where irregular water regimes have a positive effect on AM fungal colonization and persistence in comparison to treatments with a continuous supply of water (Birhane et al., 2012). AM fungi confer drought tolerance in plants by increasing osmotic adjustment and regulating stomatal conductance, reducing transpiration and leading to enhanced plant growth and vigor (Wu, 2017). In natural environments, AM fungal communities may be predisposed to tolerate drought conditions based on prior rainfall history (Kuehn et al., 1991). Deveautour et al. (2020) concluded that under reduced rainfall, estimated AM fungal alpha diversity is larger in comparison to the control treatment (ambient rainfall). On the other hand, plants that experienced complete rainfall removal, or higher environmental stress, showed lower alpha diversity estimates than the control and the reduced rainfall treatment. Thus, there lacks a consensus regarding the direction and magnitude of change in composition as soil moisture varies among ecosystem types and temporal scales.

In this study, we examine the relative importance of manipulative rainfall exclusion, local soil abiotic factors, and sampling date have on rhizospheric AM fungal communities across the NAM season of 2019. We first ask, what is the strength and direction of change in AM fungal communities (i.e., composition, richness, and turnover) in response to the manipulation of seasonal rainfall in the Sonoran Desert? And secondly, what is the relative influence soil abiotic factors and sampling date have on these communities? To answer these questions, experimental rainfall exclusion treatments were set up over the course of the monsoon season (June 2019–August 2019) to remove 0% (ambient rainfall/control), 60% (moderate exclusion), and 80% (high exclusion) of seasonal rainfall. Soil abiotic factors and AM fungal communities were sampled at the beginning and end of the monsoon season. Due to the complex ecological communities found in the rhizosphere microhabitat, there are 2 predictive hypotheses that were derived from our questions. First, with reduced rainfall, AM fungal communities could show no observed net gain or loss in taxa across the monsoon season in response to rainfall exclusion, as AM fungi have been noted to have higher desiccation tolerance than other fungal clades, such as ectomycorrhizal fungi (Allen, 1991). Alternatively, AM fungal communities could differ in turnover among rainfall exclusion treatments with ambient seasonal rainfall contributing to greater community turnover than rhizosphere communities that undergo rainfall exclusion throughout the NAM season. This could be a result of the plant-host growing season contributing to AM fungal community succession (Santos-González et al., 2007; Dumbrell et al., 2011). Understanding the response of plant-host associated AM fungal communities to drier microhabitats could shed light on how AM fungi might continue to associate with hosts in more extreme, future climate conditions.

Monsoon rainfall occurs in the Sonoran Desert ecoregion of the United States around late June–mid August. This rainfall often influences how organisms interact and can help maintain population-level persistence throughout the remaining portion of the year. This season can function as a natural experiment to manipulate seasonal rainfall intensity on a stable plant population and their interactions with arbuscular mycorrhizas. Monsoon rainfall exclusion shelters were established following guidelines outlined in Beier et al. (2012). Three treatments were established at the beginning of June 2019, a control treatment that received no rainfall exclusion (n = 22), moderate rainfall exclusion (n = 22), and high rainfall exclusion (n = 22). Short-term rainfall exclusion shelters were constructed of 6-mil heavy-duty clear plastic sheeting with 60% (moderate exclusion), or 80% (high exclusion), or 100% (control treatment) of a 1-m² surface area removed (Figure 1). Rainfall exclusion shelters were supported by two 30-cm and two 64-cm wooden landscape stakes and four 1-m wooden dowels along each edge. This design raised one edge of the shelter higher than the other providing airflow and to allow rainfall water to run off plot to ensure soil moisture manipulation. Each sampled plant rhizosphere was centered in the middle of each exclusion shelter. Lupinus neomexicanus Greene and herbaceous legume were chosen as an experimental system because of an elongated flowering time allowing for consistent identification throughout the monsoon season. This allowed us to ensure that all plants sampled were within the same life stage and had similar aboveground biomass. This species is endemic to midrange elevations (1,600–2,000 m) in Southern Arizona (e.g., 31°53'50.3"N 109°16'38.5"W).

The rhizosphere soil microhabitat was sampled twice, once prior to the monsoon season (Date A; June 2019) and once at the tail end of the monsoon season (Date B; August 2019). Baseline estimates for AM fungal community structure, soil moisture, and other soil abiotic factors are from the first sampling period prior to monsoon shelter establishment. Soil was dried in a laboratory oven for 48 h at 54°C. Soil moisture content was measured as the difference in soil weight between field wet soil and oven dried soil divided by the weight of oven dried soil (Standards Association of Australia, 1977). By the end of the monsoon season, the high rainfall exclusion treatment was on average 1.66 ± 0.29 g water/g soil drier than the moderate rainfall exclusion treatment and 2.64 ± 0.17 g water/g soil drier than the control treatment. To estimate soil abiotic factors (e.g., ammonium, nitrate, and phosphorous) in the rhizosphere environment, all rhizosphere samples were subsampled twice (2 g of rhizosphere soil per subsample) and extracted with potassium sulfate and then analyzed in triplicate (K₂SO₄; Brookes et al., 1985; Vance et al., 1987; Wu et al., 1990). Since AM fungi play a pivotal role in nitrogen cycling, nitrogen was measured in 2 ways (ammonium and nitrate). Ammonium was measured using the indophenol-blue Berthelot method (Weatherburn, 1967; Verdouw et al., 1978), and a vanadium (III) chloride acid solution and Griess reagents (Mulvaney, 1996) were used to extract nitrate from the K₂SO₄ samples. Phosphorous content for each rhizosphere sample was measured by a reaction between the K₂SO₄ extractions, ammonium para-molybdate (AMP), concentrated sulfuric acid, and a malachite green solution (Lajtha et al., 1999).

AM fungal communities were sampled using sterilized, handheld soil cores and stored in Nasco Whirl-Paks (Nasco, Fort Atkinson, WI, USA). All field equipment was surface sterilized with 70% ethanol and soil particulates were removed in-between sampling. To collect AM fungal communities, the sterilized soil core was pushed approximately 5–7 cm into the rhizosphere microhabitat along the root base, where the plant stem meets the soil. All rhizosphere samples were stored on dry ice until temporary storage (−80°C) at the Southwest Research Station in Portal, AZ, USA. Samples were then shipped on dry ice to the University of Tennessee, Knoxville, and immediately stored at −80°C upon arrival. Prior to DNA isolation, all rhizosphere soil samples were thawed for 1 h at room temperature. DNA isolation proceeded with 250 mg of rhizosphere soil per sample and the DNeasy PowerSoil HTP 96 Kit (MoBio Laboratories, Inc., Carlsbad, CA, USA). The final elution volume was modified, so that 50 µL of Solution C6 was used. Additionally, the solutions sat on the DNA membrane for 5 min at room temperature prior to centrifuging the concentrated DNA down into a sterile 96 well collection plate.

AM fungal specific libraries of the 18S rRNA gene were amplified using a nested PCR protocol (Kivlin et al., 2021). The first reaction utilized the primer pair NS1/NS4, which amplified approximately a 1,200 bp product (White et al., 1990; Raja et al., 2017). The forward NS1 primer sequence was 5′-GTAGTCATATGCTTGTCTC-3′ and the reverse NS4 primer sequence was 5′-CTTCCGTCAATTCCTTTAAG-3′. Each reaction used 2.5 µL of the DNA template, 5 µL of both the forward and reverse primers, and 12.5 µL of KAPA HiFi HotStart Ready Mix (KAPA Biosystems, Wilmington, MA, USA). This reaction ran at 95°C for 3 min, followed by 30 cycles of 95°C for 30s, 40°C for 1 min, and 72°C for 1 min with a final elongation step of 72°C for 10 min. These reactions produced 25 µL of PCR product. The following reaction amplified approximately 550 bp with the mixed primer pair, NS31 and AML2, which contain the sequences 5′-TTGGAGGGCAAGTCTGGTGCC-3′ and 5′-GAACCCAAACACTTTGGTTTC-3′, respectively (Simon et al., 1992; Lee et al., 2008). These reactions required 2.5 µL of product from the first reactions, 5 µL of NS31 and AML2, and 12.5 µL of KAPA HiFi HotStart Ready Mix. They began with 95°C for 5 min followed by 40 cycles of 95°C for 45 s, 63.1°C for 1 min, and 72°C for 1.5 min with a final elongation step of 72°C for 10 min. This nested PCR produced 25 µL of AM fungal specific 18S rRNA gene product per sample. All 18S PCR products were then stored at −20°C until further use.

PCR products were combined with 20 µL of Agencourt AMPure XP magnetic beads, placed on a plate magnet, and rinsed twice with 200 µL of 80% ethanol for product purification. The remaining ethanol was allowed to evaporate before purified product was resuspended in 50 µL of Tris-HCl. Purified products were stored at −20°C until 18S rRNA communities per sample were multiplexed. Each sample was multiplexed with a forward and reverse index primer, which are specifically designed for cluster generation on an Illumina MiSeq flow cell (Illumina Corporation, San Diego, CA, USA). These primers were obtained from the i5/i7 Nextera XT Index Kit. For the indexing PCR, 5 µL of purified PCR product from the previous reaction is required, along with 5 µL of the uniquely combined forward and reverse Nextera index primers, 25 µL of KAPA HiFi HotStart Ready Mix, and 10 µL of sterile nuclease free water. These reactions ran for an initial cycle of 95ºC for 3 min followed by an 8-cycle sequence of 95°C for 3 s, 55°C for 30 s, and 72°C for 30 s and held at 72°C for 5 min. These reactions yielded 50 µL of product that was purified with 56 µL of Agencourt AMPure XP magnetic beads, and 2 washes of 220 µL freshly prepared 80% ethanol. Cleaned and purified product was then resuspended in 25 µL of Tris-HCl. All DNA libraries were then pooled by equal concentrations. The pooled solution was diluted to 4 pM with Tris and loaded onto a V3 flow cell with a 20% PhiX spike. This flow cell was run on an Illumina MiSeq instrument at the University of Tennessee Genomics Core, Knoxville, TN, USA, set to a 2 × 300 paired end cycle.

Demultiplexing was completed by the Illumina MiSeq instrument prior to all other processing. 18S rRNA gene amplicon sequencing data were processed with the open-source DADA2 pipeline version 1.8 (https://github.com/benjjneb/dada2) and R version 3.5.2. Amplification primers for reads corresponding to 18S rRNA genes were trimmed at the 5′ end according to the length of the primer (NS31: 21 and AML2: 21). For 18S, all forward and reverse reads were truncated at 275 and 270 base pairs and filtered to remove reads below Q = 30. Error rates for the forward and reverse reads were generated by performing sample inference and error estimations repeatedly until congruence was met. Next, the derepFastq function was used in the DADA2 package to dereplicate amplicons for every sample. Here, sequences with 99% similarity were grouped, and creating unique amplicon sequence variants (ASVs) and read abundance was totaled (Callahan et al., 2016). Forward and reverse sequence reads were concatenated since the amplicon length for the AM fungal specific target region within the 18S rRNA gene is too long, such that sufficient overlap in reads was negligible between pairs. Chimeras were detected across all samples and were tossed from the data set using the removeBimeraDenovo function. Across all rhizosphere samples, fungal ASVs with less than 12 reads were filtered out of the data set.

The decostand function in the vegan package v.2.5-6 (Oksanen et al., 2019) was used to proportionally standardize all AM fungal ASVs within sample (method = “total”). Next, all AM fungal ASVs were matched against the MaarjAM database, an AM fungal sequence repository for the NS31/AML2 primer pair, using BLAST, the Basic Local Alignment Search Tool (Öpik et al., 2010). The program Practical Alignment using Saté and TrAnsitivity (Mirarab et al., 2015) was used to conduct multiple sequence alignment with MAFFT software (Katoh and Standley, 2013). Then, the aligned FASTA file was used to infer a maximum likelihood phylogenetic tree using FastTree-2 (Price et al., 2010). FastTree-2 created a generalized time-reversible CAT model of rate heterogeneity to estimate 18S rRNA phylogeny. Lastly, the treePL program was utilized to correct for differences in divergence time among branch lengths using a penalized likelihood approach (Smith and O’Meara, 2012). To understand how sensitive our inferences are to differences in sampling effort, we rarefied to the lowest number of reads among samples (1,768) for 1,000 iterations, creating 1,000 rarefied AM fungal communities per sample.

To address the question, to what extent does AM fungal community composition respond to monsoon rainfall exclusion, sampling date, and soil abiotic factors, a distance-based redundancy analysis (dbRDA) was utilized (R package vegan, Oksanen et al., 2019). The quantitative Jaccard distance (also known as the Ružička distance) was used to test community dissimilarity in composition among rainfall exclusion treatments, sampling date, and soil abiotic factors. A distance matrix was calculated for each of the 1,000 rarefied communities and a mean distance matrix was calculated for incorporation into model selection. Forward stepwise model selection on the dbRDA analysis was completed using the function ordiR2step in the vegan package (permutations = 1,000). The full scope incorporated into model selection had the following variables as main effects: rainfall exclusion treatment, sampling date, log of phosphorous, log of ammonium, and the square root of nitrate. Model selection on the dbRDA using the mean quantitative Jaccard distance from 1,000 rarefied communities determined rainfall exclusion treatment and log of phosphorous as main effects. For hypothesis testing, an analysis of variation (permutations = 9,999) was conducted on this dbRDA. A constrained ordination was made from this model using the 2 most predictive axes and their associated species scores.

Hill numbers were used to examine the effective number of AM fungal ASVs (^qD) associated with each treatment and date (Hill, 1973). The effective number of AM fungal ASVs (D) is sensitive to the order of diversity (q), and the importance of rarer ASVs is reduced when calculating diversity as the value of q increases. For instance, q = 0 represents species richness, where rare and common ASVs are weighted equally. At q = 1 (exponential Shannon’s entropy), ASVs are weighted by their proportional abundance. Lastly, q = 2 is the inverse of the Simpson’s index, where rare ASVs are downweighted when calculating the effective number of AM fungal ASVs.

We fitted Bayesian linear models to assess Δ AM fungal richness and community turnover associated with rainfall exclusion treatment at q = 0, 1, and 2. This was accomplished using the brms package (Bürkner, 2017), which calls the probabilistic programming language Stan (Carpenter et al., 2017). Δ AM fungal richness was calculated as the difference in the effective number of AM fungal ASVs across the monsoon season (Date A − Date B) for each rhizosphere sample. Community turnover was calculated as 1—the pairwise beta diversity for each rhizosphere, thus providing the proportional change in community composition between Date A and Date B. Each model used 4 Markov chain Monte Carlo (MCMC) chains, each with 2,000 iterations with a starting warmup period made up of 1,000 iterations and 1,000 post-warm up draws. The default priors provided by the brms package were used for all models. Chain convergence was confirmed using the Gelman–Rubin convergence statistic (R hat = 1 for all models). Effective sample size for all models was >3,000. Inferences were based on examining the 95% highest density interval of model parameters calculated using the hdi function in the HDInterval package (v. 0.2.2) (Meredith and Kruschke, 2020). Hypotheses regarding differences in Δ AM fungal richness and community turnover among treatments were examined using the posterior joint probability distribution for the parameters using the hypothesis() function in brms, which provides the proportion of times a parameter is either greater than, or less then, another parameter across MCMC steps, in this case zero where zero represents no change in Δ AM fungal richness and community turnover. Δ AM fungal richness and AM fungal community turnover were further analyzed using standardized abundance ASV data and rarefied data to confirm that Bayesian inference was robust against the differences in reads among samples.

To examine the phylogenetic community structure among rainfall exclusion treatments, we used a phylogenetically informed Hill number approach at q = 1 for alpha diversity, which is analogous to a phylo-Horn N-assemblage overlap measure (Chao et al., 2010). This measure considers phylogenetic relationships by calculating the effective number of AM fungal ASVs (richness) given relative abundance and a phylogenetic tree. At q = 1, we used the hill_div function in the hilldiv package v.1.5.1 for proportional abundance data and all rarefied community data (Alberdi and Gilbert, 2019). We calculated the phylogenetically informed Δ AM fungal richness across the monsoon season for each sample and modeled it as a fitted Bayesian linear model as described above. The same one-sided hypothesis was used to test significant patterns in phylogenetically informed Δ AM fungal richness among treatments.

After running the DADA2 customized pipeline as described above, a total of 34,011 ASVs were generated. Based on the best practices for sequencing and analyzing AM fungal communities, ASVs with less than 10–12 reads (for this study <12 was chosen as the threshold) were omitted from the data set. This resulted in a total of 5,012 ASVs that correspond to the 18S rRNA gene region. After ASVs were confirmed to be AM fungi in the MaarjAM database, 5,011 ASVs remained in the final data set. One ASV did not match any taxa from the MaarjAM database and was removed. All other reads were assigned to the subphylum Glomeromycotina (Soudzilovskaia et al., 2020). Figure 2 illustrates the relative abundance in reads assigned to the subphylum Glomeromycotina and the distribution of these reads among orders and families within this subphylum. AM fungal communities from Date A, prior to the onset of rainfall exclusion treatment, were utilized as a baseline understanding of AM fungal diversity and relative abundance present in the L. neomexicanus rhizosphere. Of the 4 orders within Glomeromycotina; Glomerales, Diversisporales, Paraglomerales, and Archaeosporales represented 80%, 10%, 6%, and 4% of reads recovered. The genus Glomus was the most abundant genus recovered from the rhizosphere microhabitat, representing approximately 67% of all total reads at Date A. Further, the second most abundant genus recovered was Claroideoglomus with approximately 12% of total reads sampled at Date A. Lastly, the third most abundant genus was Paraglomus (Paraglomerales) comprising approximately 7% of reads in Glomeromycotina. The most abundant AM fungal variant recovered was Glomus Valera-Cervero15 BG21 representing 14% of all total reads recovered at Date A in this study.

AM fungal community composition was characterized by assigning taxonomy and documenting change in relative read abundance among treatments before and after the monsoon season (Figure 3). Across all treatments, Diversisporaceae decreased in relative abundance across the monsoon season with the largest observed difference in the high rainfall exclusion treatment. On the other hand, Figure 3 illustrates an increase in Paraglomeraceae in response to monsoon rainfall, with the most extreme observed change in the control treatment. Glomeraceae was observed to be the most common AM fungal family recovered, and relative abundance decreased over time. Claroideoglomerales decreased in relative abundance among all treatments, but most notably in the control treatment.

We used the mean quantitative Jaccard distance matrix from 1,000 rarefactions to test the strength of rainfall exclusion and the log of phosphorous on AM fungal community composition and relative abundance. This model’s constrained variables explained 3.55% of total variation in multivariate space. Additionally, axes 1 and 2 explained 47.14% and 34.04% of the 3.55% of total variation, respectively. Figure 4 shows the relative change in direction and strength of AM fungal communities given rainfall exclusion treatments relative to estimates in composition before treatments were installed. Points here represent the median centroid in multivariate space with ±1 standard error for each treatment and sampling date. A permutational analysis of variance showed that rainfall exclusion treatments were a significant predictor for AM fungal community structure and composition (P < 0.05). Based on this model, the log of phosphorous was also a significant predictor of AM fungal community composition in multivariate space (P < 0.001), contributing to the observed changes in Figure 4 from community estimates before treatments were installed. All measured soil abiotic factors were shown to not be confounded or interact with rainfall exclusion treatments (Supplemental File 1).

The differences in alpha diversity between Date A and Date B, or Δ AM fungal richness, were analyzed from diversity order q = 0–2 using linear models. From q = 0 to q = 2, the greatest observed gain in the effective number of AM fungal ASVs was in the moderate rainfall exclusion treatment (Figure 5). On the contrary, the control treatment on average lost AM fungal ASVs throughout the monsoon season. As q increases, the 95% highest density intervals decrease around the median change in AM fungal richness across the monsoon season. The control treatment was observed to have an estimated median change of −15.13, −4.85, and −2.96 in AM fungal richness across the monsoon season for q = 0, 1, and 2, respectively. We can infer that rainfall had a negative effect on Δ AM fungal richness of rhizospheres. This hypothesis was observed to be true in our posterior samples 77%, 89%, and 96% of the time for q = 0, 1, and 2 (Table 1). The median estimated Δ AM fungal richness for the moderate rainfall exclusion treatment was 4.68, 6.10, and 3.10 at q = 0, 1, and 2. For the moderate rainfall exclusion treatment, we can infer positive change in richness throughout the monsoon season for each value of q. For q = 0, 1, and 2, this hypothesis was observed 76%, 97%, and 99% of the time of the posterior samples. In the high rainfall exclusion treatment, we can infer that minimal, positive change in AM fungal richness was observed 73%, 78%, and 83% of the time of the posterior samples. The patterns here illustrate that rainfall exclusion had a positive effect on Δ AM fungal richness and that as rare AM fungal taxa are removed from our analysis (q = 2), the likelihood of this to be truer increased. These inferences in Δ AM fungal richness are consistent with inferences generated from the 1,000 rarefied samples (Supplemental File 2: Figure S1).

Phylogenetic relatedness among rainfall exclusion treatments was measured as Δ AM fungal richness across the monsoon season using the phylo-Horn N-assemblage overlap measure (q = 1 when informed by phylogeny). Median change in phylogenetic relatedness was −1.05, 1.21, and 0.41 for the control treatment, moderate rainfall exclusion treatment, and the high rainfall exclusion treatment, respectively. For the control treatment, we can infer negative change in phylogenetic relatedness 85% of the time throughout the monsoon season. The moderate and the high rainfall exclusion treatments differ in inference from the control treatment showing positive change in phylogenetic relatedness 94% and 83% of the time, respectively. A full summary of hypothesis testing for each treatment including the 95% highest density intervals and associated posterior probabilities can be found in Supplemental File 2: Table S1. Additionally, inferences per treatment for rarefied communities mirror that of inferences uncovered through this Bayesian analyses applied to proportional abundance data.

From q = 0 to q = 2, AM fungal community turnover was measured as 1—pairwise beta diversity for each sample, which represents proportional change in AM fungal community composition across the monsoon season (between Date A and Date B). Across all values of q, the control treatment had the greatest AM fungal community turnover in comparison to other treatments (Figure 6). At q = 0, 1, and 2, the median estimate for AM fungal turnover for the control treatment was estimated as 0.72, 0.57, and 0.60 and was greater than zero 100% of the time for each value of q (Table 2). This indicates that there was community turnover in the control treatment due to rainfall. Median community turnover estimates for the moderate and high rainfall exclusion treatments were more similar to each other and differed from the control treatment. From q = 0 to q = 2, rhizospheres that underwent moderate rainfall exclusion were shown to have lower estimates of turnover in comparison to the control treatment 96%, 93%, and 78% of the posterior samples were less than zero, respectively. This suggests that rainfall exclusion treatments influenced AM fungal community composition across the monsoon season (Figure 6). Further, rhizosphere microhabitats that experienced high rainfall exclusion across the monsoon season had the least AM fungal community turnover out of all treatments and differed from the control treatment where 99%, 99%, and 96% of the posterior samples were less than zero for q = 0, 1, and 2. As rare species were down weighted from calculating AM fungal community turnover, the proportion of times our hypotheses were observed in the MCMC chains decreases. Figure 6 also illustrates that as the value of q increases from 0 to 2, so does the highest density intervals around the estimated values of AM fungal community turnover. Inferences made for AM fungal community turnover using proportional abundance data are in concordance with inferences made from 1,000 rarefied samples (Supplemental File 2; Figure S2). Lastly, median point estimates for rarefied AM fungal community turnover can be found in Supplemental File 2 as Table S2. Here, one can see that the 95% highest density intervals around the median are less than the highest density intervals seen around the median point estimates calculated using proportional abundance data.

In this study, we explicitly show the importance rainfall exclusion has on AM fungal richness and community turnover in the rhizosphere microhabitats L. neomexicanus, an herbaceous plant found within the NAM region of the Sonoran Desert. Our results illustrate that variation in AM fungal community structure may be partially predicted by phosphorous and short-term rainfall exclusion. As soils get drier, AM fungal communities have been shown to have less seasonal community turnover and may be more influenced by phosphorous. Based on model selection, all other soil abiotic factors failed to explain significant variation in AM fungal community composition and richness. By the end of the monsoon season, Δ AM fungal richness at q = 0, 1, and 2 responded to rainfall exclusion treatments; however, the response differed among rainfall exclusion treatments, where the control treatment showed a loss in taxa due to rainfall and the moderate and high rainfall exclusion treatments showed increase in fungal ASVs. Rainfall exclusion also contributes to differences in phylogenetic relatedness among rhizosphere microhabitats in the Sonoran Desert, with rainfall contributing to a decrease in phylogenetic community composition. Lastly, AM fungal community turnover across the monsoon season was dependent on rainfall exclusion treatments, whereby the greatest turnover across all values of q was observed in the control treatment. Here, moderate rainfall exclusion could be a threshold that favors the coexistence of different AM fungal taxa by the end of the monsoon season (Tilman, 1977). These results support our second prediction that AM fungal composition responds to short-term rainfall exclusion treatments, and more importantly, emphasizes the relationship between rainfall perturbations and soil abiotic factors like phosphorous content in the soil.

Our results for Δ AM fungal richness across the monsoon season are consistent (Deveautour et al., 2020) as our moderate exclusion treatment showed an increase in AM fungal ASVs from q = 0 to q = 2 relative to the control treatment. Further, more extreme rainfall removal resulted in lesser changes in alpha diversity, which may be a result of taxa specific requirements for population-level persistence throughout decreased soil moisture and warmer summer temperatures (Tilman, 1977). In comparison to other ecological factors, precipitation and soil water availability more strongly influences AM fungal composition in semiarid ecosystems (Gao et al., 2016). In a study by Xiao et al. (2019), AM fungal diversity was influenced by phosphorous limitation, while total nitrogen affected the relative abundance of AM fungi. We show that in the context of naturally occurring AM fungal communities in the Sonoran Desert that phosphorous content in the rhizosphere microhabitat plays a substantial role in structuring the compositional spread of AM fungal communities and that nitrate content did not contribute to the variation in the effective number of AM fungal ASVs based on model selection. This suggests that as precipitation becomes more infrequent that AM fungal communities may be more dependent on specific, limiting soil abiotic factors in dry, arid climates.

The effects of sampling date on plant communities have been well studied and supported, providing some foundation to examine the temporal aspects of microbial communities (Smith and Huston, 1990; Guo and Brown, 1996). Temporal effects on the maintenance of AM fungal communities have been studied in a variety of ecological contexts ranging from volcanic ash fields (Higo et al., 2015), mixed hardwood forests (Davison et al., 2012), to grasslands (Deveautour et al., 2020). There is mixed evidence showing that AM fungal community composition is temporally dynamic. For example, in mixed hardwood forests, AM fungal communities have been shown to be spatially distinct but temporally stable across the growing season, perhaps due to the heterogenous plant community and a consistent AM fungal species pool over their 2-year sampling scheme (Davison et al., 2012). Yet, on the same time scale, Higo et al. (2015) state that differences in AM fungal composition are dependent on plant cover from the previous year, thus contributing to fungal community turnover year to year. Temporal effects have also been observed on much smaller time scales, such as over one season in a wet tropical forest (Waring and Hawkes, 2015). Our multivariate models show distinct clustering for the composition of AM fungal communities over the NAM season and by rainfall exclusion treatment. Additionally, we show that AM fungal community turnover responds to rainfall exclusion treatments from q = 0 to q = 2. The effect seen here may illustrate the contribution available soil moisture has on the number of ASVs throughout the monsoon season and on the relative abundance of each ASV in a dry, arid ecosystem. Our results further highlight the interconnected relationships between regional climate and local soil abiotic factors that could hinder a generalizable pattern consensus of a spatiotemporal response of AM fungal community composition.

Here, we demonstrate that rainfall exclusion affects variation in AM fungal composition and the difference in the effective number of AM fungal ASVs, where moderate rainfall exclusion led to increased estimates in alpha diversity. Additionally, we show that AM fungal community turnover is dependent on rainfall exclusion when common and rare ASVs are weighted equally and also when rare ASVs are down weighted in estimating community turnover. The work outlined in this study leaves several forward directions for AM fungal ecology. Specifically, what are the relative temporal scales needed to observe AM fungal responses to manipulative field experiments? How does temporal turnover in AM fungal communities differ between within year and among year exposure to environmental stress? AM fungi can be somewhat ubiquitous at local scales (Davison et al., 2015) and have been shown to be heterokaryotic (i.e., multinucleate with genetically distinct nuclei; Bruns et al., 2018). This could lead to lack of agreement in empirical biodiversity studies, as species richness estimates might be dependent on spatial scale, temporal scale, and the weight of rare ASVs. Therefore, while we did find that AM fungal communities respond to changes in seasonal monsoon precipitation, more studies across a wider range of ecological contexts could further uncover the direction and magnitude of changing climate conditions on plant–host AM fungal interactions.

Sequence data are available in the NCBI Sequence Read Archive under BioProject ID: PRJNA1065923.

The supplemental files for this article can be found as follows:

Supplemental File 1.xlsx

Supplemental File 2. docx

The authors would like to begin by giving thanks to the Collections Committee and the TENN Herbarium at the University of Tennessee, Knoxville (Department of Ecology and Evolutionary Biology) for recommendations regarding experimental design and funding field portions of this study. Secondly, they would like to thank the Arnold Lab at the University of Arizona, School of Plant Sciences for hosting during field collection. They additionally acknowledge Michele Lanan and the Southwest Research Station (American Museum of Natural History) for coordination and lab space while in the Chiricahua Mountains. They would like to thank the U.S. Department of Agriculture, Forest Service, and the Douglas Ranger District for permitting this study during 2019. They would also like to thank Stephanie Kivlin for providing access to lab space for 18S library preparation and soil assays during the coronavirus pandemic and for providing comments that substantial improved this work.

The Kivlin lab was funded by the NSF DEB 1936195, Texas Ecological Laboratory, and The University of Tennessee, Knoxville. This work was also supported by The National Science Foundation DEB 1638922 to JAF.

The authors declare that they have no competing or conflict of interest.

Contributed to conception and design: JRD, JAF.

Contributed to acquisition of data: JRD.

Contributed to analysis and interpretation of data: JRD, JAF.

Drafted and/or revised the manuscript: JRD, JAF.

Approved the final version for submission: JRD, JAF.

How to cite this article: Dickey, JR, Fordyce, JA. 2024. Greater change in arbuscular mycorrhizal fungal richness as a response to short-term rainfall exclusion across the North American monsoon season. Elementa: Science of the Anthropocene 12(1). DOI: https://doi.org/10.1525/elementa.2022.00116

Domain Editor-in-Chief: Steven Allison, University of California Irvine, Irvine, CA, USA

Associate Editor: Rebecca A. Lybrand, Department of Land, Air and Water Resources, University of California, Davis, CA, USA

Knowledge Domain: Ecology and Earth Systems