Arbuscular mycorrhizal fungi (AMF) are a globally important group of fungi that form mutualistic associations with the roots of the majority of land plants (74%; Brundrett, 2009). These mutualisms often improve plant growth and resource acquisition, and thus AMF are recognized as drivers of plant community structure and function and biogeochemical cycling (van der Heijden et al., 1998; Eom et al., 2000). Despite their ecological importance, our understanding of AMF diversity lags far behind that of other groups of fungi (Sanders and Rodriguez, 2016), largely due to the cryptic diversity of phenotypically similar species, late adoption of modern molecular methods, and a study bias toward grassland ecosystems. This limits our understanding of how AMF diversity feeds back to influence plant species' distribution, productivity, diversity, and community assembly (Davison et al., 2012; Bainerd et al., 2014), as well as how AMF–plant interactions vary in response to abiotic drivers such as temperature and precipitation regimes (Egerton-Warburton et al., 2007; Camenzind et al., 2014) and edaphic stresses (e.g., eutrophication, metal-contaminated soils; Cabello, 1997; Bunn et al., 2009; Hassan et al., 2011).

Recent advances in molecular genetic approaches, notably amplicon-targeted high-throughput sequencing technologies such as 454 pyrosequencing, have largely revolutionized the characterization of AMF diversity, phylogeny, and biogeography (Öpik et al., 2009; Schüßler and Walker, 2010) and have revealed high levels of species diversity and complex relationships between AMF and their host plants. For example, these technologies have led to the identification of more than 350 well-characterized molecular taxa (Öpik et al., 2010); shown fine-scale spatial and temporal structuring of AMF communities (Lekberg et al., 2007; Dumbrell et al., 2010, 2011; Bainerd et al., 2014), especially to edaphic constraints (Wang et al., 2016; Wilson et al., 2016); and enabled a sweeping systematic revision of the Glomeromycota (Schüßler and Walker, 2010; Redecker et al., 2013). Although these technologies specifically target AMF species, there is a need to achieve deeper sequencing depth in environmental samples to more completely characterize levels of AMF diversity that are relevant to ecosystem function (Smith and Peay, 2014; Sanders and Rodriguez, 2016).

Ultra-high-throughput amplicon sequencing capacity using the Illumina MiSeq platform far surpasses 454 technology. It follows that the use of barcoded AMF primers that take advantage of the deep sequencing capacities of the Illumina MiSeq platform could be a promising way to address these needs (Caporaso et al., 2010, 2011). Our objective was to adapt an AMF-specific primer pair (NS31/AML2) with barcodes and to develop and apply a broadly applicable protocol for AMF community amplicon sequencing on the Illumina MiSeq platform. Here, we report on these primers, detail the sequencing protocol and bioinformatics pipeline, and demonstrate the efficacy of our approach by sequencing AMF communities in complex environmental samples from tropical dry forests.

MATERIALS AND METHODS

Sampling sites and edaphic factors—Primers and protocols were tested on environmental samples collected in two seasonally dry tropical forest sites located in the Yucatán Peninsula, Mexico: Reserva Ecológico El Eden (21.195°N, 87.167°W) and Rancho La Higuera (20.445°N, 87.352°W). These sites are privately operated conservation areas (>60 yr after the last disturbance). The local landscape consists of porous Cenozoic limestone (<150 m a.s.l.) overlain by shallow (10–20 cm deep) soils with high organic matter content (20–39% carbon by combustion). The plant community comprises secondary successional forests and is dominated by a high diversity of woody tree species. The climate of the region is warm, subhumid (mean annual temperature 26°C), and strongly seasonal with a wet season (May–October) followed by a marked dry season (November– April). Precipitation ranges from 1100 to 1600 mm/yr, with ~80% of the total precipitation received during the wet season (López-Martinez et al., 2013).

Single soil samples comprising mixed soil and roots (to 20 cm deep) were collected at the drip line of mature Brosimum alicastrum Sw., Vachellia cornigera (L.) Seigler & Ebinger, Metopium brownei (Jacq.) Urb., Ceiba pentandra (L.) Gaertn., Bursera simaruba (L.) Sarg., and Manilkara zapota (L.) P. Royen trees in February 2013 (dry season; N = 48). Soil samples were imported into the United States under a U.S. Department of Agriculture (USDA) Animal and Plant Health Inspection Service (APHIS) Permit to Receive Soil (P330-11-00358), and all further handling and processing of these samples were conducted in the APHIS authorized containment facility at the Chicago Botanic Garden. Samples were stored frozen (−20°C) until analysis.

Each soil sample was extracted using a ratio of 1 : 10 soil : deionized water, after which the extracts were analyzed for levels of ammonium (NH₄; Weatherburn, 1967), nitrate (NO₃; Doane and Horwath, 2003), phosphate (P; Baykov et al., 1988), and pH and electrical conductivity (EC; Hach Instruments, Loveland, Colorado, USA). In general, soils in La Higuera were of higher pH and contained significantly higher levels of P and NH₄ than those at El Eden (Appendix 1). These differences resulted in a significantly higher soil N : P supply ratio in La Higuera.

Table 1.

Individual components and complete sequences of Illumina MiSeq-compatible custom sequencing primers used to sequence amplicon libraries in this study and their estimated melting temperatures.

Genomic DNA extraction and amplicon library generation—Genomic DNA was extracted from each sample using 0.25 g soil/plant material with the PowerSoil DNA Isolation kit (QIAGEN-MO BIO, Carlsbad, California, USA), following manufacturer's protocols.

Libraries of small subunit ribosomal RNA (hereafter referred to as 18S) fragment amplicons were prepared for each sample with a single-step PCR reaction using primers NS31 (Simon et al., 1992) and AML2 (Lee et al., 2008) that we modified for use with Illumina sequencing platforms following the protocol of Caporaso et al. (2011). These modifications include the addition to both primers of technical adapter sequences for annealing to Illumina flow cells, a standard “pad” sequence, and a novel two-base linker sequence (Table 1); the pad sequence and the two-base linker sequence were both designed to reduce secondary structure formation. Reverse primer constructs were also modified to include a 12-base Golay error-correcting barcode (or index) to enable demultiplexing during data processing (Appendix 2). Primer Prospector (Walters et al., 2011) was used with the MaarjAM database of AMF 18S sequences (Öpik et al., 2010) to optimize linker sequences and to test for secondary structure formation in all barcoded primer constructs. Complete sequences of primer constructs NS31f-il and AML2r-il and of all barcodes checked for secondary structure formation are provided in Table 1 and Appendix 2, respectively.

For each sample, PCR was carried out using 10 µL of 5PRIME HotMasterMix, 0.5 µL of NS31f-il, 13 µL of molecular biology–grade water (Fisher Scientific Bio-Reagents, Fair Lawn, New Jersey, USA), 0.5 µL of uniquely barcoded AML2r-il, and 1 µL of genomic DNA extract. The PCR reaction was run using the following thermal cycler program: initial denaturation for 3 min at 94°C; followed by 35 cycles of 45 s at 94°C (denaturation), 60 s at 63.1°C (annealing), and 90 s at 72°C (extension); followed by a final extension step of 10 min at 72°C. PCR reactions were carried out in triplicate for each sample and pooled prior to final sequencing library preparation. The final sequencing library and sequencing reactions were performed at Argonne National Laboratories (Lemont, Illinois, USA). All individual sample amplicon libraries were quantified fluorometrically with a Quant-iT Pico-Green dsDNA assay (Invitrogen Molecular Probes, Eugene, Oregon, USA). An equimolar sequencing library was produced, cleaned using a MO BIO UltraClean PCR Clean-up kit (QIAGEN-MO BIO), and sequenced on an Illumina MiSeq using version 2, 2 × 250-bp paired-end chemistry (Illumina, San Diego, California, USA).

Data processing and bioinformatics—Sequence read processing and bioinformatic analyses (Appendix 3) were performed using QIIME version 1.9.1 (Caporaso et al., 2010), BLAST (Altschul et al., 1990), and the vegan package (Oksanen et al., 2015) in the R statistical environment (R Core Team, 2014). Sequence reads were included in the analyses only if the index read matched a barcode sequence used in this study with two or fewer errors.

During data processing, we found that a large proportion of sequenced amplicons were too long to allow for overlap with the Illumina MiSeq version 2, 2 × 250-bp sequencing technology, and thus could not be aligned and assembled (Appendix 4). The average amplicon length was ~530 bp based on alignment of paired forward and reverse reads to full-length amplicon regions in the MaarjAM Virtual Taxon (VT) database of 348 well-resolved AMF species-level sequence clusters ( http://maarjam.botany.ut.ee). Instead, we conducted all downstream analyses using only the 250-bp forward reads based on evidence provided by Davison et al. (2012). More specifically, Davison et al. (2012) found that artificially truncating NS31/AML2 reads from 400 to 170 bp resulted in a near identical capacity to capture AMF diversity because the majority of taxonomically informative characters occurred in the 5′-most 170 bases.

Table 2.

Total number of operational taxonomic units (OTUs) clustered at three similarity thresholds showing the number of OTUs retained after each filtering step during data processing, and used in each of the nine analyzed data sets.

Raw reads were demultiplexed and quality filtered with default parameters and a quality threshold of 20. Sequencing reads were truncated after four consecutive base calls with quality scores less than 20 (99% confidence interval). Truncated reads that were less than 75% of their original length were removed from further analyses. Because input reads were 250 bp, the resulting data set contained reads with a minimum 187 bp.

Fig. 1.

Alpha rarefaction curves for 90% similar (A), 95% similar (B), and 97% similar (C) operational taxonomic units (OTUs) in samples from El Eden (red) and La Higuera (blue); Venn diagrams illustrating the number of 90% similar (D), 95% similar (E), and 97% similar (F) OTUs unique to El Eden (EE) or La Higuera (LH), or present at both sites. Data are shown for the most inclusive, all OTUs data sets.

Operational taxonomic units (OTUs) were clustered using reads from all 48 AMF community libraries using an open reference strategy and the MaarjAM VT database. We clustered OTUs at 97% similarity, which is conventionally used as a species-level threshold, and also at 95% and 90% similarities, as numerous reports have documented intraspecific, and even intraindividual, variation in AMF ribosomal DNA that exceeds the 3% dissimilarity threshold (e.g., Clapp et al., 1999; Rodriguez et al., 2004).

Probable artifacts and sequences from nontarget taxa were excluded by BLAST against the complete MaarjAM database of 5934 AMF 18S sequences. OTUs were excluded if they did not hit a MaarjAM database sequence with an E-value below 10⁻⁵⁰ and alignment of at least 90% of the read length with sequence identity equal to or greater than the OTU clustering similarity threshold (e.g., 90% sequence identity for 90% OTUs). We also used a second BLAST search against a database of primer constructs to screen for OTUs that may have incorporated any sequence originating from the primer (E-value threshold of 10⁻¹⁰; identity threshold of 50%), but no further OTUs were removed in this step. Chimeric OTUs were identified using the USearch v6 implementation of UCHIME (Edgar et al., 2011) with default parameters and a database of all 5934 sequences from the MaarjAM database and their reverse complements. Any OTU with a UCHIME score greater than 1 was removed from further analyses. Finally, all OTU-type sequences were BLASTed against the National Center for Biotechnology Information (NCBI) nonredundant nucleotide database, and any OTU with a best hit that was not identified as an AMF in the GenBank record was removed from further analyses. Taxonomy was assigned to remaining OTUs using the QIIME implementation of RDP Classifier (Wang et al., 2007) retrained using the complete MaarjAM database, and with a minimum confidence threshold of 0.8. Finally, to reduce sampling imbalance between sites, samples with fewer than 1000 reads remaining after all quality-filtering steps were removed. This approach resulted in a total of 33 samples, representing 14 samples from El Eden and 19 samples from La Higuera.

For each collection of OTUs (El Eden, La Higuera), we explored the effects of rare taxa on patterns of diversity by conducting three separate sets of analyses at both sites that included all OTUs, excluded singletons, or excluded all OTUs with fewer than 10 constituent sequences to capture only core diversity (Smith and Peay, 2014). The removal of rare OTUs did not result in read numbers falling below the 1000-reads/sample threshold. Table 2 demonstrates the number of OTUs retained in each clustering threshold following the quality-filtering steps.

We used QIIME to generate Bray–Curtis dissimilarity matrices for each data set using rarefied data (1000 reads/sample) and tested for significant differences in AMF community structure between sites and in relation to soil chemical variables (log-transformed except pH) using principal coordinates analysis (PCoA) and PERMANOVA tests with the capscale and ADONIS functions in the vegan R package. Next, we tested whether there were significant differences in AMF community dispersion among sites using the betadisper function in vegan. We used QIIME to generate Chao₁, Shannon–Wiener, and Simpson (1-D) alpha diversity metrics for each sample (not rarefied). Differences in OTU richness and diversity and soil chemical variables between sites were analyzed using ANOVA and Tukey's honest significant difference posthoc tests in R. Finally, we used the vegdist function in vegan to generate a Euclidean distance matrix based on soil chemical variables (NO₃, NH₄, P, N : P supply ratio, EC, pH), and used a Mantel test with Pearson and Spearman coefficients to test for correlation between the soil chemical and the AMF community matrices.

RESULTS

The 5,977,389 quality-filtered, demultiplexed sequence reads used in this study have been uploaded to the NCBI Sequence Read Archive, and are associated with BioProject PRJNA329250.

General primer performance—After quality filtering and exclusion of reads belonging to artifacts, nontarget taxa, and rare OTUs (if any), each of the nine data sets contained a total number of reads ranging from 2,325,440 (97%, ≥ 10-ton set) to 2,776,849 (90%, all OTUs set). Mean reads per sample ranged from 70,468 ± 116,616 (mean ± SD) to 84,146 ± 133,491, while the median number of reads for samples in a data set ranged from 22,932 to 28,368. Alpha rarefaction indicated that ∼20,000 to >60,000 reads per sample may be necessary to adequately sample these AMF communities, with more reads required at higher similarity thresholds (Fig. 1A–C). This was particularly notable in the El Eden data sets (Fig. 2) and in more rare-OTUinclusive data ( Appendix S1 (apps.1700017_s1.pdf)).

AMF species identification—The total number of OTUs increased markedly with increasing clustering similarity threshold and rare OTU inclusivity in both sites (Fig. 1D–F), and ranged from 196 (90% similar ≥ 10-ton OTUs) to 6416 (97% similar all OTUs; Table 2).

RDP Classifier assigned 91.5% (90% similarity), 99.2% (95% similarity), and 99.8% (97% similarity) of OTUs to genus-level or species-level accessions in the VT database with at least 80% confidence. Overall, these OTUs accounted for >99.5% of the sequencing reads in each data set.

AMF communities were characterized by a small group of highly abundant AMF species and numerous rare species. In each clustering threshold, AMF communities were dominated by 23–25 OTUs that, in total, represented ∼78–87% of reads. The remaining reads were distributed among the numerous subordinate (rare) OTUs ( Appendix S2 (apps.1700017_s2.xlsx)). This pattern is consistent with the lognormal model distribution of AMF species noted in previous studies (Dumbrell et al., 2010) and suggests that a large number of specialist and/or narrowly endemic AMF fungi occurred in these forests.

Nine genus-level assignments were made to OTUs and were consistent with well-characterized AMF genera from community studies (Redecker et al., 2013). As in other studies of AMF communities (Eom et al., 2000; Egerton-Warburton et al., 2007), Glomus sensu lato (s.l.) (taxa formerly classed as Glomus Group A) dominated the AMF community and accounted for 76–97% of OTUs and >90% of all sequence reads in every data set (Table 3; Fig. 3;  Appendix S3 (apps.1700017_s3.pdf)).

Fig. 2.

Averaged alpha rarefaction curves of observed operational taxonomic unit (OTU) richness at 90%, 95%, and 97% similarity clustering thresholds in El Eden (EE) and La Higuera (LH). Vertical bars represent the standard deviation of the mean.

RDP Classifier assigned 2–8% of OTUs in each clustering threshold to a species-level VT. The majority of species-level assignments were to Glomus VT; this genus comprised 54%, 70%, and 94% of AMF communities in 90%, 95%, and 97% clustering thresholds, respectively. In each clustering threshold, at least one species-level assignment was made in each genus except Gigaspora (Table 3). Complete RDP Classifier assignments for all OTUs at each clustering similarity threshold are provided in  Appendices S4 (apps.1700017_s4.txt),  S5 (apps.1700017_s5.txt), and  S6 (apps.1700017_s6.txt).

AMF community diversity, structure, and composition—Our analyses revealed strong positive effects of both increasing OTU clustering similarity threshold and increasing rare OTU inclusivity on observed and estimated taxonomic richness (Tables 2, 4). Specifically, increasing the OTU clustering similarity threshold significantly inflated the observed OTU richness owing to the clustering of rare OTUs (singletons, <10 tons) in the most read-rich samples in El Eden (Table 4;  Appendices S7 (apps.1700017_s7.xls),  S8 (apps.1700017_s8.xls), and  S9 (apps.1700017_s9.xls)). Even so, there were consistent patterns of OTU richness and diversity across all data sets. For example, OTU richness (observed, Chao₁ estimates) was significantly higher in El Eden than La Higuera (P < 0.02) while Shannon–Wiener and Simpson index values did not differ significantly between sites (Table 4). In addition, OTU richness was negatively correlated with soil pH (P < 0.020) and NH₄ levels (P < 0.028) in all data sets.

Table 3.

Number of operational taxonomic units assigned to each of nine arbuscular mycorrhizal fungal genera using three clustering similarity thresholds.

Fig. 3.

The percentage of operational taxonomic units (OTUs) (A, B) and percentage of sequence reads (C, D) assigned to arbuscular mycorrhizal fungal genera at El Eden (EE) and La Higuera (LH) study sites. Data are shown for 97% similar OTUs at two inclusivity levels: all clustered OTUs (A, C) and ≥ 10-ton OTUs (B, D). Colors indicate genus assigned by RDP Classifier: Glomus (green); Ambispora, Archaeospora, and Gigaspora (white); Scutellospora (orange); Acaulospora (black); Paraglomus (aqua); Claroideoglomus (yellow); Diversispora (blue); and no genus-level assignment (red).

AMF community composition was also impacted by OTU clustering threshold and rare OTU inclusivity. However, some consistent trends were apparent. For example, AMF communities were dominated by species of Glomus s.l., and those in La Higuera contained a greater abundance and diversity of Diversispora species and lower levels of Claroideoglomus than El Eden (Fig. 3). Removing the rare OTUs from analyses increased the similarity in AMF community composition between sites with respect to the proportion of OTUs assigned to each genus (Fig. 3A, B;  Appendix S3 (apps.1700017_s3.pdf)) but had little effect on read abundances (Fig. 3C, D). These results support the presence of site-specific AMF taxa and suggest that the primary difference in proportion of OTUs assigned to each genus among sites was due to relatively large numbers of rare taxa (i.e., rare Diversispora OTUs at La Higuera), while differences in read abundance assigned to genera between sites were due to a few highly abundant OTUs.

Table 4.

Observed operational taxonomic unit (OTU) richness per individual sample and site, number of OTUs unique to each site, and indices of diversity (Chao₁, Simpson, Shannon–Wiener) for arbuscular mycorrhizal fungal communities at both study sites. Values represent the site mean with standard deviation in parentheses; El Eden, N = 14 samples; La Higuera, N = 19 samples.

These results were further supported by significant differences in AMF community structure (P < 0.0047, PERMANOVA; Fig. 4), but no significant difference in dispersion between sites (P > 0.28 for all data sets, PERMDISP2). AMF community structure was significantly influenced by pH (P < 0.0035, PERMANOVA) and NH₄ (P < 0.0062, PERMANOVA; Fig. 4) because OTU richness was negatively influenced by increasing soil pH or NH₄ level (Fig. 5). These patterns were consistent across clustering thresholds ( Appendix S10 (apps.1700017_s10.pdf)). Comparisons between PCoA ordinations for environmental variables (Euclidean distance; Fig. 6A) and OTUs at all clustering thresholds and inclusivity levels (Bray–Curtis distance; 90% similar, ≥2-ton data shown in Fig. 6B) support a highly significant correlation between AMF community composition and soil properties using either Pearson or Spearman coefficients (Mantel test: P < 0.004, R² > 0.20).

Fig. 4.

Principal coordinates analysis (PCoA) ordination plots of arbuscular mycorrhizal fungal (AMF) communities sampled at El Eden (EE, red circles) and La Higuera (LH, blue triangles) using Bray–Curtis dissimilarities based on all operational taxonomic units (OTUs) clustered at 90% (A), 95% (B), and 97% (C) similarity thresholds. Percentage values on the axes represent the variation in AMF community dissimilarity explained by each axis. Ellipses represent the central tendency of communities at each site. Vectors denote the magnitude and direction of statistically significant effects of soil properties on AMF community dissimilarity.

Data sets using the three OTU clustering similarity thresholds also showed that AMF community structure was responsive to soil P levels ( Appendix S10 (apps.1700017_s10.pdf)). Negative relationships between taxonomic richness and soil P levels were only supported (P < 0.05, ANOVA) in 95% and 97% similar data sets, while AMF community structure responses to soil P were only supported in 90% and 97% similar data sets (P < 0.030, PERMANOVA). Rare OTU inclusivity had little effect on these patterns.

Fig. 5.

Relationship between operational taxonomic unit (OTU) richness and soil NH₄ (A) and pH (B) in the 97% ≥10-ton data set, which is representative of patterns observed across all data sets. Red dashed lines show the best fit of linear regression models with P < 0.003.

DISCUSSION

The overarching goal of our study was to determine whether AMF diversity and variations in AMF communities could be adequately captured on the Illumina MiSeq platform, and to determine its potential utility in large-scale surveys of AMF communities. To address this goal, we modified existing 18S primers for barcoding, applied robust protocols with which to undertake 18S amplicon analysis on the Illumina MiSeq platform, and developed bioinformatics pipelines for data processing.

Our results clearly demonstrate that the application of barcoded NS31/AML2 primers improves the level of resolution in AMF species identification, diversity, and community composition in complex environmental samples. This primer pair effectively amplified a wide diversity of AMF genera and species, and did not appear to exclude taxa that have been omitted previously due to primer bias (e.g., Archaeosporaceae and Paraglomeraceae; Lee et al., 2008). Along with the deeper sequence coverage provide by the Illumina MiSeq, this approach also revealed one of the highest levels of AMF species richness recorded to date (Öpik et al., 2010), ranging from 196 OTUs in the most conservative data set (≥10-ton, 90% threshold) to more than 6000 at the highest levels of OTU clustering similarity and rare OTU inclusivity (all samples, 97% threshold). A large percentage of AMF taxa were also present in extremely low abundance, thereby supporting the capacity of our protocols to capture rare taxa. In contrast, previous studies using 454 pyrosequencing indicated that AMF communities hosted, on average, 70 AMF taxa (e.g., Dumbrell et al., 2010), while estimates using morphological methods indicated ~45 AMF taxa within a community (Eom et al., 2000; Egerton-Warburton et al., 2007).

We also captured biogeographic differences in AMF communities between the two study sites. Across all data sets, there were consistent and significant differences in OTU richness abundance. For example, La Higuera contained a greater abundance and diversity of Diversispora species and lower levels of Claroideoglomus than El Eden. In addition, the significantly higher pH and levels of NH₄ and P at La Higuera appeared to drive the observed site effect on AMF community composition and structure. These results are in general agreement with spore-based studies of AMF communities in other systems (Egerton-Warburton et al., 2007). Our study was not designed to examine the mechanistic basis of these results. Based on studies elsewhere, however, it is possible that the dry season constrains the AMF community to taxa that are physically or physiologically tolerant of low soil moisture (e.g., Glomus and Diversispora; Augé, 2001) or to high soil P fertility or pH (Wang et al., 2016). Alternatively, these shifts may reflect differences in host plant requirements during the dry season for AMF that increase host drought tolerance (e.g., stomatal control, cytokinin production; Augé, 2001) or increase the acquisition of limiting nutrients from carbonate substrates (N, Fe, Zn; Labidi et al., 2012). Irrespective, our results indicate a high potential to use our AMF protocol in large-scale sequencing projects to address AMF diversity with sufficient taxonomic precision, to determine the extent to which AMF assemblages vary across host plants and ecosystems, and to resolve AMF species' responses to edaphic stressors, such as anthropogenic N deposition, in complex environmental samples.

Our study also highlighted a number of technical considerations. First, a key finding was OTU inflation, particularly at the 97% clustering threshold, which is the level applied in most mycorrhizal fungal studies. Between the 97% and 90% clustering thresholds, there was a twofold increase in OTU richness at the ≥10-ton level, an 8-fold increase when considering ≥2-ton OTUs, and a 17-fold increase when considering all clustered OTUs (Table 2). In our study, this result was primarily due to the recovery of numerous rare or unique AMF taxa in the most readrich samples (see Table 4;  Appendix S1 (apps.1700017_s1.pdf)), rather than to issues such as uneven number of sequences among samples. To avoid overestimation of AMF community diversity, excluding all OTUs with fewer than 10 constituent sequences (regardless of clustering threshold) will result in levels of taxonomic (OTU) richness consistent with current estimates of AMF species richness (Öpik et al., 2010; Schüßler and Walker, 2010).

Fig. 6.

Principal coordinates analysis (PCoA) ordination plot of soil properties (A) and arbuscular mycorrhizal fungal (AMF) communities (B) in El Eden (EE) and La Higuera (LH). Percentage values on the axes represent the variation in soil properties (A) and AMF operational taxonomic unit (OTU) read abundance (B) explained by each axis. Ellipses represent the central tendency of communities, and vectors denote the magnitude and direction of the effects of significant soil nutrients on AMF communities.

Second, the majority of OTUs could not be assigned to any species-level accession in the MaarjAM database. While this result indicates that novel AMF species likely occur in the Yucatán as they do in other tropical systems (Chaiyasen et al., 2014), it raises questions about our ability to identify AMF species and catalog their diversity in environmental samples. One possibility is that OTU matching was hampered by the limited availability of well-characterized AMF taxa from tropical forests. Alternatively, the clustering of sequences to generate AMF VT (Öpik et al., 2010) may have reduced the potential for OTU matching if sequences of multiple species were lumped into the same OTU (Bruns and Taylor, 2016) or if sequences from a single species were assigned to multiple OTUs (House et al., 2016). Either scenario could mask phylogenetic diversity (and inferences of functionality) or AMF species with large amounts of intraspecific variation. Thus, a more comprehensive reference database of AMF sequences is now needed to improve our capacity to identify AMF to species level and to address within-species sequence variation (House et al., 2016). In particular, more direct assessments that use well-characterized sequences from individual spores or single-spore cultures are needed.

Finally, there are limitations to using MiSeq version 2, 2 × 250-bp chemistry with AMF-barcoded samples. Improvements in the MiSeq chemistry with version 3 (2 × 300 bp) is expected to further improve the differentiation of OTUs and the taxonomic resolution of AMF species by allowing consistent assembly of forward and reverse reads into a ∼530-bp sequence fragment. Preliminary analyses of recent 2 × 300-bp sequencing data suggest that OTUs clustered using assembled reads are assigned to species level at approximately two times the rate of OTUs clustered from forward reads alone using similar screening stringency levels (Morgan and Egerton-Warburton, unpublished data).

Our results support the continuing development and use of highthroughput sequencing approaches to address the AMF “black box.” As a first step, the tools detailed herein allow the detection of ecologically relevant levels of AMF diversity that shape plant community composition, diversity, and nutrient acquisition in natural and restored communities, including rare and unique species (Sanders and Rodriguez, 2016). While this approach is broadly applicable to most ecosystems, it is especially important in those where major gaps remain in our understanding of AMF species richness and their role in plant community composition and function.