More bang for the buck? Can arbuscular mycorrhizal fungal communities be characterized adequately alongside other fungi using general fungal primers?

Abstract

High throughput sequencing has revolutionized microbial ecology. We now have the ability to characterize both prokaryotic and eukaryotic microbial communities more thoroughly and disentangle relationships among broad organismal groups (Bödeker et al., 2016; Morriën et al., 2017; Toju et al., 2017). Arbuscular mycorrhizal fungi (AMF; subphylum Glomeromycotina, Spatafora et al., 2016) are no exception, and high throughput sequencing has revealed greater richness within communities compared with Sanger sequencing (Öpik et al., 2009), as well as an increased understanding about global distribution patterns and their potential underlying drivers (Davison et al., 2015). AMF are apparently asexual (but see Corradi & Brachmann, 2017), obligate symbionts that exchange nutrients and other services for plant carbon (Smith & Read, 2008). There are c. 300 morphologically-defined (Öpik & Davison, 2016) and 1000 or more molecularly-defined species of AMF (Kivlin et al., 2011; Öpik et al., 2014). This is low compared to other mycorrhizal fungi, and surprising given the prevalence and potential ability by AMF to colonize > 200 000 plant species (Brundrett, 2009). As molecular techniques have evolved, different marker regions have been used to characterize AMF communities, each with their strengths and weaknesses (Hart et al., 2015). The small subunit rRNA (SSU rRNA) gene is the region most often targeted because primer pairs exist that amplify most known AMF families. Also, sequences can be aligned against virtual taxa (VTs) using the curated database MaarjAM (Öpik et al., 2010), which facilitates comparisons across studies. The internal transcribed spacer (ITS) region is the universal barcode for fungi (Schoch et al., 2012) and AMF specific primers exist for ITS (Krüger et al., 2009), but this region is less frequently used in AMF community surveys. This may be due to poor amplification of AMF with general fungal ITS primers (Tedersoo et al., 2015, 2018), or high variation within species and even individuals that hampers clade recognition (Thiéry et al., 2016). Yet with the exponential increase in sequencing output with high throughput sequencing, AMF sequences are now detected and discussed in studies using general fungal primers that target all soil fungi. For example, Tedersoo et al. (2014) used ITS to show that AMF sequences make up a small proportion of the total fungal abundance and that AMF richness is positively correlated with soil pH and potential evapotranspiration. Likewise, Leff et al. (2015) used a shift in total number of AMF sequences relative to other fungal groups to discuss potential changes in AMF abundance with fertilizer additions. Targeting the ITS rather than the SSU is appealing as it allows fungal ecologists to make broad comparisons that include all fungi. It is unclear, however, if general fungal ITS primers detect similar AMF communities and show comparable responses to shifts in environmental conditions as AMF-specific SSU primers do. If they do, targeting the ITS rather than the SSU will provide broader comparisons using the same effort and resources without the loss of precision in detecting shifts in AMF communities. To answer this question, we took advantage of a recent glasshouse experiment and a season-long field survey where the same samples were amplified using both AMF-specific SSU primers and general fungal ITS primers. Detailed methods are provided in Supporting Information Methods S1, but briefly, the glasshouse experiment involved spotted knapweed (Centaurea stoebe) and mountain brome (Bromus marginatus) grown for 3 months in spotted knapweed rhizosphere soil under either high (c. 30%) or low (c. 10%) soil moisture before harvesting roots for analyses. The field survey consisted of a seasonal (April, July, September, 2016) sampling of showy milkweed (Asclepias speciosa) roots from an intermountain grassland in Montana, USA. Here we use significant shifts observed in the SSU datasets due to soil moisture (R. Bunn, unpublished) and season (B. G. Larkin, unpublished) to ask if the AMF amplified using general fungal ITS primers and the same DNA templates report similar shifts in composition. Portions of the SSU and the ITS regions were amplified using the primer pairs WANDA-AML2 (Lee et al., 2008; Dumbrell et al., 2011) and fITS7:fITS7o/ITS4 (Ihrmark et al., 2012; Kohout et al., 2014), respectively, where the fITS7o is a modification of fITS7 with improved detection of all AMF families. We used a semi-nested PCR for both primer pairs where barcodes and tags were added in the second PCR step. PCR products were cleaned, quantified and pooled within studies at equimolar concentrations before Illumina sequencing using 2 × 300 bp MiSeq at the University of Idaho (www.ibest.uidaho.edu). Sequences were quality filtered, and thereafter identified against the MaarjAM SSU and ITS reference databases (Öpik et al., 2010) by BLAST using scripts of Vasar et al. (2017). SSU data in MaarjAM database are pre-assigned to species proxies (VT). For comparability, ITS data in MaarjAM was clustered into operational taxonomic units (OTUs), and used in the same manner as VT for new read identification. Full details of bioinformatical procedures are provided in Methods S1. Sample-taxa matrices were generated for each study and target region by pooling sequences identified as the same VT or OTU. Any sequences that did not have at least 97% sequence similarity with a reference AMF sequence across a minimum of 95% of the length of alignment were removed from further analyses, as were VTs and OTUs represented by single sequences. While this closed reference identity-assignment approach resulted in the removal of some putative AMF sequences, the majority of AMF sequences could be identified against MaarjAM databases, showing that AMF diversity in the study systems is well represented therein (Table S1). Comparisons between SSU and ITS were based on the same samples and AMF sequence numbers (200 sequences per sample in the glasshouse and 1000 sequences per sample in the field survey). This choice of data standardization represents a trade-off between maximal taxon capture and losing as few samples as possible (Figs S1, S2). We used three complementary approaches to evaluate the results. First, we used Analysis of Similarities (ANOSIM, 999 permutations; Clarke, 1993), which is a multivariate approach used for categorical data, to compare treatment effects between ITS and SSU outputs in the two studies. We then used Procrustes analysis to statistically examine the concordance between the SSU- and ITS-based communities in principal coordinate analysis (PCoA) space, as well as Mantel tests. ANOSIM was run within the open source Quantitative Insights into Microbial Ecology software package (QIIME, vs 1.6.0; Caporaso et al., 2010), and Procrustes and Mantel analyses were performed using the procrustes and mantel functions in the vegan R package (Oksanen et al., 2017). We constructed phylogenetic trees using representative sequences of VT and OTUs and calculated extrapolated richness. We also compared the sequence abundance within AMF families between the two marker regions for both studies using a Wilcoxon signed rank sum test on paired sequence data (stats package; R Core Team, 2017). For this last comparison, we chose family-level rather than genus- or species-level contrasts because some AMF genera, as taxonomically recognized, are difficult to delimit when including environmental sequences where more diversity is apparent than has been possible to characterize by traditional taxonomy (Öpik et al., 2014). Also, different taxon delimitation has been used in the two markers (phylogenetic curator-based VT delimitation of SSU sequences vs automatic clustering algorithm-based OTU-delimitation of the highly variable ITS2 sequences). Until the nonoverlapping markers can be analyzed as a single amplicon and synonymies established, these aspects do not allow direct comparison of SSU- and ITS-based AMF species proxies (Öpik & Davison, 2016). Finally, to assess the sequence diversity independent of database coverage, we used DADA2 to identify unique sequences (Callahan et al., 2016) and clustered those at 97% similarity using VSEARCH (Rognes et al., 2016) in QIIME2 (version 2017.10; https://qiime2.org/; Caporaso et al., 2010). All samples contained AMF sequences, but the proportion of AMF differed substantially between the two markers. More than 90% of quality filtered sequences in both SSU datasets matched with MaarjAM VTs, which is not too surprising given that we used AMF-specific primers. The proportion of AMF sequences was substantially lower in the two ITS datasets obtained with general fungal primers (Table S1). We lost 33 ITS samples and two SSU samples (out of 97 samples in total) from the glasshouse study, and seven ITS samples in the field survey (out of 59 in total) due to insufficient numbers of AMF sequences in those samples. Thus, comparisons between ITS and SSU were based on 64 glasshouse samples and 52 field samples. Of those ITS samples, 11% of total sequences (range 3.5–30.5%) matched with MaarjAM ITS sequences in the glasshouse experiment compared with 43% (range 7.5–73.5% across samples) in the field survey. The two ITS datasets highlight the fact that the proportion of AMF amplified by general fungal primers can vary dramatically among samples, and may differ among studies and host species. It is also possible that proportions differ depending on the particular ITS primers used (Ihrmark et al., 2012; Kohout et al., 2014) and whether DNA is extracted from soil or roots. Indeed, studies that have extracted DNA from soil rather than roots report AMF sequence proportions of < 5% (Orgiazzi et al., 2012; Tedersoo et al., 2014; Leff et al., 2015), and direct comparisons by Berruti et al. (2017) show lower proportions in soil than roots. Therefore, we recommend preliminary sequencing of selected samples to ensure sufficient coverage of AMF (i.e. where taxon accumulation curves are approaching saturation) among samples before any large-scale efforts targeting ITS. It is worth noting, however, that the 500–5000 AMF sequences we detected per sample in the two ITS datasets are an order of magnitude greater than many sequence numbers in previous studies using AMF-specific primers, cloning and Sanger sequencing (Helgason et al., 1998; Kohout et al., 2014). Thus, even though AMF sequences often constitute a relatively small proportion of the total sequence data, general fungal ITS primers and high throughput sequencing can characterize AMF with greater precision than older methodologies that targeted AMF exclusively. Both markers reported highly significant (P = 0.001 in all cases) responses of AMF community composition to soil moisture (RSSU = 0.36; RITS = 0.20) and season (RSSU = 0.20; RITS = 0.31). This was supported by significant Procrustes analyses (P < 0.001 for both studies) and Mantel tests (Mantel RGlasshouse = 0.66, P = 0.001; RField = 0.77, P = 0.001), with visual separation of AMF communities between moisture treatments and season (Fig. 1). However, the high m2 values (sum of squared deviations between sample pairs, where lower values mean a better fit between matrices) in the Procrustes analyses indicate that the two target regions can capture different communities. This was reflected in the phylogenetic trees where the two markers display different shares of taxa among families (Figs 2, S3–S6). Typical of many AMF communities (Oehl et al., 2005; Davison et al., 2015), both target regions reported a dominance by Glomeraceae richness and sequence abundance (Fig. 2; Table S2). However, ITS consistently reported higher Glomeraceae sequence numbers than SSU, whereas some subordinate families were only detected with SSU and not with ITS2 approach (Table S2). ITS2 marker length is similar across families (Table S3) and thus, the apparent bias in family-level amplification is not due to high throughput sequencing favoring shorter amplicons (Lindahl et al., 2013). Nor is it a database issue because all families are represented in the MaarjAM ITS database, and the great majority of putative AMF excluded from analyses were either Glomeraceae or Claroideoglomeraceae (Table S4). Future comparisons using mock communities and in silico simulations would allow testing of trends highlighted in these datasets, help avoid potential discrimination against specific clades of Glomeromycotina, and ensure that critical functional groups of AMF are not systematically undetected. Regardless of the discrepancies we found in phylogenetic community composition, there was no consistent difference in observed or extrapolated richness between ITS and SSU approaches in the two datasets. When clustering unique sequence variants at 97% similarity, however, SSU richness was remarkably lower than ITS, possibly due to pooling of multiple VTs into individual clusters given that the sequence similarity within VTs range from 97% to 99% (Öpik et al., 2014; Table S5). In summary, the two molecular approaches generated slightly different AMF community compositions but reported similar responses to environmental shifts. This may be because both target regions amplified the abundant Glomeraceae family in a similar manner as indicated by significant positive correlations in SSU and ITS Glomeraceae sequence abundances across samples in both studies (RGlasshouse = 0.73, P < 0.001; RField = 0.76, P < 0.001). While the comparison made here is valid only for the primer pairs and conditions employed in these studies, our results agree remarkably well with a recent comparison using the same marker regions but different primers and fewer samples (Berruti et al., 2017). Like us, Berruti et al. (2017) found that even though general fungal ITS2 primers preferentially amplified Glomeraceae and detected fewer families than AMF-specific SSU primers, comparable estimates of AMF community structure and responses to environmental variables were obtained. All current approaches used to characterize AMF communities have trade-offs that should be considered when choosing the specific marker. For example, AMF-specific SSU primers may be able to amplify more families and provide a broader view of the AMF community than general fungal ITS primers. However, if the main goal of a study is to assess treatment responses of mainly dominant taxa, then our results indicate that general fungal ITS primers could be used, so long as a sufficient number of Glomeromycotinan sequences are obtained in all samples to detect adequate completeness of the existing diversity. Using general fungal ITS primers is advantageous as they can reveal community-level response and interactions across the kingdom Fungi using the same effort and resources as when targeting AMF alone. The authors are grateful to Peter Kennedy, Björn Lindahl and three anonymous reviewers for constructive comments on an earlier version of this manuscript. Y.L., L.S.B. and B.G.L. thank MPG Ranch for funding. M.Ö., S-K.S. and M.V. are supported by the Estonian Ministry of Education and Research (grant IUT20-28), the European Union through the European Regional Development Fund (Center of Excellence EcolChange) and ERA-NET Cofund BiodivERsA3 (Project SoilMan). P.M.A. is funded by the Natural Sciences and Engineering Research Council of Canada, and R.B. acknowledges WWU and WWU Biology Glasshouses. Y.L. and M.Ö. conceived of the idea for this comparison. B.G.L., L.A.B. and Y.L. conducted the field survey included, and R.B., P.M.A. and Y.L. conducted the glasshouse experiment. Y.L., L.S.B., M.V. and S-K.S. performed the analyses presented here and all co-authors contributed to their interpretation. Y.L. wrote the first draft with input from all co-authors. Please note: Wiley Blackwell are not responsible for the content or functionality of any Supporting Information supplied by the authors. Any queries (other than missing material) should be directed to the New Phytologist Central Office. Fig. S1 VT and OTU accumulation across sequence depth per sample in the glasshouse experiment. Fig. S2 VT and OTU accumulation across sequence depth per sample in the field survey. Fig. S3 Detailed phylogenetic tree depicting family distribution of representative ITS sequences in the glasshouse experiment. Fig. S4 Detailed phylogenetic tree depicting family distribution of representative SSU sequences in the glasshouse experiment. Fig. S5 Detailed phylogenetic tree depicting family distribution of representative ITS sequences in the field survey. Fig. S6 Detailed phylogenetic tree depicting family distribution of representative SSU sequences in the field survey. Table S1 Total sequence numbers matching with databases, potential AMF and non-AMF among SSU and ITS2 sequences in the glasshouse and field study Table S2 Total sequence abundance, range and median per sample within AMF families in the glasshouse and field study Table S3 ITS2 base pair ranges among AMF families in the MaarjAM ITS database Table S4 Identification of putatively AMF and non-AMF ITS sequences based on ≥ 90% similarity and coverage against INSDC in all four datasets Table S5 Observed and extrapolated taxonomic richness (VTs and OTUs for SSU and ITS, respectively) observed in the glasshouse and field study Methods S1 Detailed description of the field survey and glasshouse experiment, PCR amplification and associated bioinformatic and statistical analyses. Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.