Between 5 and 30 November 2011, 34 composite soil samples were collected from woody plant-dominated ecosystems in Papua New Guinea (PNG) following a standard protocol (Tedersoo et al. 2014). Briefly, two soil cores (5 cm diam. to 5 cm depth) were collected at a distance of 1-1.5 m from 20 randomly selected tree trunks in each of the 2500-m² plots. These 40 soil cores from each plot were pooled and subjected to DNA extraction, analysis of soil chemistry, and assignment of climatic variables as described in Tedersoo et al. (2014). For in-depth molecular analysis, we selected 26 samples that represented various ecosystems in PNG and that had relatively well-preserved and inhibitor-free DNA (Suppl. material 1).

To address barcode and primer biases, we selected seven barcodes in SSU (variable domains V4 and V5), ITS (ITS1 and ITS2), and LSU (variable domains D1, D2, and D3) of the nuclear rDNA (Figure 1). These short barcodes represent the most common marker choices in cloning and metabarcoding studies of fungi either alone or combined within a longer stretch of rDNA. To amplify these barcodes, we selected the most commonly used primers for each respective marker, targeting primer pairs that produce amplicons of 250-400 nucleotides on average. For the ITS1 and ITS2 barcodes, we used alternative forward primers to further quantify potential biases arising from primer choice. Due to a lack of suitable primers between the variable D1 and D2 domains of LSU, we designed the LF402 primer and its reverse complement LF402F that are relatively specific to fungi compared to other eukaryotes (Appendix 1). For the V4 domain of SSU, we generated a reverse primer Euk742R that covers a wide range of eukaryotes including the most common fungal groups. In addition to previous modifications to the ITS3 and ITS4 primers (ITS3tagmix and ITS4ngs, respectively; Tedersoo et al. 2014), we also optimized the commonly used SSU515F, SSU1196R, ITS1, ITS1F, and LR0R primers to improve their taxonomic coverage and to balance their melting temperature compared to the respective reverse primers. All primers are detailed in Table 2 and their relative positions are indicated in Figure 1 and on the UNITE homepage (https://unite.ut.ee/primers.php). We further evaluated the suitability of these primers through in silico-based assays on SSU and LSU data sets of SILVA release 119 (Quast et al. 2013), the SSU, ITS and LSU data sets of Bellemain et al. (2010), and the ITS reference data set of UNITE (Abarenkov et al. 2010). In addition to manual comparison of primer and template sequences, we used EcoPCR (www.grenoble.prabi.fr/trac/ecoPCR) for SSU and LSU. The mismatched primer variants were further BLASTn-searched against the International Nucleotide Sequence Database consortium (INSDc) to determine the natural frequency of inferred mismatches and to distinguish these from rare sequencing errors. The primer mismatches with fungal templates are outlined in Appendix 1.

Figure 1. 

Map of ribosomal DNA indicating variable regions as well as primers used and/or discussed in this study. Primers pairs used for HTS are highlighted.

A reverse or forward primer for each barcode was supplemented with one of the sixteen 10-base identifier tags (Table 1). We used a mock community sample with known composition of 24 species (Suppl. material 2) and a soil sample from Alaska as positive controls. We also added negative controls and spiked two soil samples (G2655 and G2658) with 2% genomic DNA from Archaeorhizomyces finlayi (isolate Ny10) and A. borealis (isolate 600) to screen accumulation of their sequences in primer pairs that have no or two mismatches (ITS4ngs and LR0Rngs) to the Archaeorhizomycetes. The PCR cocktail consisted of 0.6 µl DNA extract, 0.5 µl each of the primers (20 pmol), 5 µl 5xHOT FIREPol Blend Master Mix (Solis Biodyne, Tartu, Estonia), and 13.4 µl double-distilled water. PCR was carried out in four replicates using the following thermocycling conditions: an initial 15 min at 95 °C, followed by 30 cycles at 95 °C for 30 s, 55 °C for 30 s, 72 °C for 1 min, and a final cycle of 10 min at 72 °C. Four replicate PCR products were pooled and their relative quantity was estimated by running 5 µl amplicon DNA on 1% agarose gel for 15 min. DNA samples yielding no visible band were re-amplified using 35 cycles in an effort to obtain sufficient PCR product, whereas samples with a very strong band were re-amplified with only 25 cycles. Amplicons were subjected to quantity normalization with a SequalPrep Normalization Plate Kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s instructions. Normalized amplicons were subjected to ligation of Illumina adaptors using two variants of the TruSeq DNA PCR-free HT Sample Prep kit (Illumina Inc., San Diego, CA, USA). All samples were sequenced in a single Illumina MiSeq 2×300 paired-end run. HTS sequences are available in Short Read archive under accession SRP055957.

Selected soil samples (Suppl. material 1) were also subjected to shotgun sequencing of the full metagenome. DNA was fragmented to 300–400 bases on average with the Covaris S2 instrument (Covaris Ltd, Woburn, MA, USA). Samples were prepared for sequencing with the TruSeq Nano DNA HT Library Prep Kit (Illumina Inc.) according to the instructions of the manufacturer. These samples were sequenced using the Illumina HiSeq 2×250 paired-end protocol covering 80% of a full plate.

We amplified DNA from two soil samples (G2655 and G2658 that were spiked with A. finlayi and A. borealis DNA) using the ITSOF and LR7 primers (Hopple and Vilgalys 1994; Tedersoo et al. 2008) and Pfu DNA Polymerase (Fermentas, Kaunas, Lithuania) to produce 2000–2200 base-pair (bp) amplicons for improved classification. The PCR mix included 2.5 µl 10× Pfu buffer, 0.4 µl Pfu polymerase, 0.25 µl dNTP mix (20 mM each), 0.7 µl each primer (20 µmol), 2 µl template DNA, and 18.45 µl ddH₂O. Thermocycling conditions included an initial denaturation at 95 °C for 3 minutes, 30 cycles of denaturation for 30 seconds, annealing at 53 °C for 30 seconds, extension at 72 °C for 5 minutes, and a final extension for 10 minutes. The amplicons were cloned using StrataClone Blunt PCR Cloning kit (Agilent Technologies, La Jolla, CA, USA) following the manufacturer’s instructions. Positive inserts were re-amplified using primers ITSOF and LR5 (Hopple and Vilgalys 1994). Amplicons were sequenced bidirectionally using primers ITS5 (White et al. 1990) and LR5. Unique sequences are available in the UNITE database (accessions UDB014801-UDB014812).

Paired-end sequencing (2×300 bp) in the Illumina MiSeq sequencer resulted in 12,771,565 reads. LSU and SSU amplicons were paired, quality filtered, and demultiplexed using the LOTUS pipeline (Hildebrand et al. 2014) with the following options: overlap >9 bases, average quality >28, no DNA ambiguities allowed, no differences to primer+identifier allowed, homopolymer length <9, and minimum sequence length >199 bases. Primers and tags were removed and reads were trimmed to 200 bases covering the variable domain, and clustered with cd-hit 4.6.1 (Fu et al. 2012) at 97% sequence similarity into operational taxonomic units (OTUs). Chimeric OTUs were detected using the UCHIME de novo and reference-based algorithms as implemented in USEARCH 7 (Edgar et al. 2011). The SILVA 119 data set was used as reference corpus.

The ITS reads were quality filtered using MOTHUR 1.33.3 (Schloss et al. 2009) with the following options: average quality over 15 nucleotide window >30, no ambiguities allowed, and homopolymer length unrestricted. The remaining reads were paired using the PANDASEQ assembler (Masella et al. 2012) with minimum overlap of 25 bases and demultiplexed using mothur (no differences to primer+identifier allowed, minimum sequence length >200 bases). ITS1 and ITS2 sequences were processed in ITSx 1.0.9 (Bengtsson-Palme et al. 2013) to remove flanking 18S, 5.8S and 28S rRNA gene fragments. Potential chimeras were removed as described above, using the UNITE-UCHIME release as reference (Nilsson et al. in press). The few reads containing primer sequences were excluded using Unix commands. The quality-filtered ITS1 and ITS2 sequences were clustered as described above.

Following exclusion of singletons from all HTS data sets (cf. Brown et al. 2015), representative sequences were selected for each OTU based on the greatest similarity to the consensus sequence. These sequences were subjected to taxonomic assignments using the NAÏVE BAYESIAN CLASSIFIER (NBC; Wang et al. 2007) and BLASTn (options: word size=7, gap penalty=-1, gap extension penalty=-2 and match score=1) searches retrieving 10 best BLASTn matches. For both methods, The SILVA 119 eukaryote data set and the INSDc fungal data set served as references for SSU and LSU. The UNITE 7.0beta data set (https://unite.ut.ee/repository.php) and the INSDc eukaryote data set were used as a reference for BLASTn- and NBC-based taxonomic assignment of ITS, respectively.

Because of differential taxonomic resolution among the barcodes, we used the taxonomic assignments of both NBC and BLASTn searches to complement each other as both methods alone provided no assignment for ca 40% of OTUs due to poor representation of fungal data in SILVA, obvious misidentifications in INSDc, and great abundance of taxonomically unassigned sequence data (resulting in poor resolution using NBC). To optimize classification, we therefore combined and verified results of different methods and determined approximate E-value and sequence similarity thresholds for robust identification to phylum or class level for each barcode. For the SSU barcodes, we determined that sequences with a BLASTn E-value <e^-80 and >90% similarity to verified fungal sequences could be reliably assigned to the fungal kingdom. Sequence similarities at or greater than 95% served to provide class-level assignments. Taxa with short introns in the SSU V5 barcode (primers pair SSU817F-SSU1196Rngs) such as some members of Pezizomycotina and Agaricomycetes exhibited more variation, and these groups were therefore manually evaluated for class-level affinity. For ITS1, ITS2, and LSU D2 (primers LF402Fmix-TW13), we determined that BLASTn E-values <e^-50 and sequence similarity >75% over >70% sequence length allowed robust assignment to the fungal kingdom. For individual classes, we relied on an 80% sequence similarity threshold, except the early diverging lineages, Archaeorhizomycetes, and Cantharellales, where we used 75% sequence similarity. For the LSU D1 barcode (LR0Rngs-LF402 primers), BLASTn E-values <e^-80 and sequence similarity above 80% were sufficient to consider sequences to be fungal, whereas similarity >85% was indicative of class-level affiliation. For the D3 barcode of LSU (LR3R-LR5 primer pair), BLASTn E-value <e^-100 and sequence similarity >85% to fungi was strongly suggestive of fungal origin, and sequence similarity >90% allowed placement to classes. Many Leotiomycetes and Eurotiomycetes possessed introns close to the LR5 primer site, which allowed identification of these groups at >80% similarity level. For all barcodes, phylogenetic placement of BLASTn matches and NBC were highly concordant, except several instances in Glomeromycota, Microbotryomycetes, and Tremellomycetes. For these groups, we relied on the 10 best BLASTn matches.

We followed the taxonomy of INSDc, except raising several early diverging lineages to phylum rank (cf. Tedersoo et al. 2014). The OTUs of ITS1 and ITS2 barcodes were assigned to EcM lineages following Tedersoo and Smith (2013). For several common lineages such as /hebeloma-alnicola, /cortinarius, /cenococcum and /meliniomyces, SSU and LSU barcodes provided insufficient discrimination from non-EcM relatives such that it was impossible to reliably assign ecological categories to SSU and LSU data.

For the shotgun metagenome data, samples were demultiplexed, and LSU and SSU regions of all organisms were extracted using SORTMERNA (Kopylova et al. 2012). Ribosomal DNA sequences were classified against the SILVA 119 reference database (SSU, ITS, and LSU) and UNITE 7.0beta (ITS region). We also generated sets of 10 best BLASTn matches against the INSDc fungal data set. We manually checked conflicting assignments as well as those with <95% similarity to fungal sequences to improve phylum- and class-level classification.

To understand potential amplification biases related to sequence length in the ITS1 and ITS2 barcodes, we downloaded all ITS sequences of the 16 most common fungal classes (based on our amplicon data) from UNITE 7.0beta data set. Ribosomal RNA genes flanking the ITS1 and ITS2 barcodes were trimmed using ITSx 1.0.9. Average and median values and standard deviations were calculated for each group.

For OTU-based statistical analyses, we removed all non-fungal sequences and rarefied all amplicon samples to a depth of 8609 sequences using MOTHUR. This depth represents the median number of sequences of the ITS1 (ITS1ngs-ITS2 primers) barcode that was the second lowest among all markers (Table 2). The LSU D2 barcode comprised several times less sequences, and this marker was therefore removed from comparisons of taxonomic richness and statistical performance of barcodes, but kept for comparisons of taxonomic coverage. Soil sample G2653 was poorly covered by sequence data for three barcodes (<1000 sequences) and was therefore omitted from all analyses. Samples G2666 (SSU V5 barcode), G2677 (SSU V4 barcode), and G2678 (LSU D3 barcode) were variously poorly represented in single data sets. Their values were treated as missing data in the statistical analyses.

OTU accumulation curves and their 95% confidence intervals were computed for ITS1 and ITS2 barcodes using ESTIMATES 9.1.0 (Colwell 2013). To further address whether ITS1Fngs and ITS1ngs or ITS3tagmix and gITS7 primers exhibit additive performance (i.e., recover more OTUs when combined), we pooled the primer pair data by samples and rarefied these combined samples to 8609 sequences, followed by inference of OTU accumulation curves as described above.

Differences in OTU richness among samples and barcode-primer combinations were evaluated based on two-way main-effect ANOVAs supplemented with Unequal n HSD tests for multi-level comparisons. To address the relative performance of the eight barcode-primer pair combinations in recovering the role of spatial, edaphic, floristic, and climatic predictors on fungal community composition, we performed multivariate permutational ANOVAs as implemented in the ADONIS routine of the vegan package of R (R Core Development Team 2013). Geographical coordinates were translated into Principal Coordinates of Neighbouring Matrices (PCNM) vectors with soil element concentrations being logarithm-transformed prior to analyses. All four categories of variables were included in separate matrices for these analyses following Põlme et al. (2014). We used both Bray-Curtis dissimilarity and Hellinger distance in analyses of community composition. We also analyzed the relative importance of primer choice for the ITS1 and ITS2 barcodes in ADONIS. Because data sets differed substantially in matrix connectance, we tested whether this matrix property affects the adjusted determination coefficients (R²_adj) that are widely used to interpret model performance quantitatively through the explained variance. Therefore, we re-analyzed the ITS1Fngs-ITS2 data set at 12 levels of connectance by selectively removing species based on their frequency from one (4%) to 12 (48%) and calculated the correlation between the logarithm of connectance and R²_adj as revealed from ADONIS analyses.

To further test changes in phylogenetic community structure among samples and barcode-primer combinations, we assigned the OTUs to fungal classes. Sequence number-based proportions of classes were log-ratio transformed in relation to the proportion of non-Agaricomycetes and unassigned fungi to ensure statistical independence of these groups. The proportion of EcM and arbuscular mycorrhizal (AM) to non-mycorrhizal sequences was similarly log-ratio transformed but only analyzed for the four ITS data sets. Using the Hellinger distance and Bray-Curtis distance, the effects of barcode-primer pair combination were analyzed using ADONIS. Global Nonmetric Multidimensional Scaling (GNMDS) graphs were drawn in parallel using the same options. 95% confidence ellipses were calculated using ORDISURF routine in the Vegan package of R.

Potential analytical biases in the recovery of fungal classes by ribosomal DNA region (SSU, ITS, and LSU) and analysis method (metagenomics and amplicon) were addressed based on the average values of 14 shared samples using two-way main-effect ANOVAs neglecting interactions. rDNA region-based biases in amplicon and metagenomics data sets were further tested using two-way ANOVAs including the samples and regions as fixed factors.

To understand possible amplification biases related to sequence length in the ITS1 and ITS2 barcodes, we downloaded all ITS sequences of the 16 most common fungal classes from UNITE 7.0beta. Ribosomal RNA genes flanking the ITS1 and ITS2 barcodes were trimmed using ITSx 1.0.9. Average and median values and standard deviations were calculated for each group for illustrative purpose.

Although combinations of samples and primer pairs were normalized separately, sequences assigned to each primer pair were differentially represented in the raw and final data sets (Table 2). Much of the data loss from a total of 11,691,106 demultiplexed sequences was attributable to unpaired reads and stacked reads, where certain reads consisted of two repeats of itself. Some of these incomplete repeats could not be successfully filtered out by bioinformatics tools, but these sequences were easily recognized due to their abnormal length (>480 bases) and were removed manually. Remarkably, one such artefactual OTUs formed the second most abundant taxon in the ITS2 data set. The unequal occurrence of barcode-marker combinations was not related to the average marker length as the least represented LSU D2 (LF402Fmix-TW13) barcode was the shortest among all barcodes. Therefore, we recommend not to multiplex different barcodes in the same libraries unless their roughly equal representation is known.

Barcodes generated by the universal primer pairs SSU515Fngs-Euk742R (SSU V4) and LR3R-LR5 (LSU D3) exhibited the distinctly lowest proportion of fungal sequences (67–68%), suggesting that fungi account for roughly two thirds of eukaryote ribosomal DNA in the studied soils on average. The proportion of fungal sequences was the greatest for the primer pair LR0Rngs-LF402, reaching 99.9% of all sequences. Classifications based on both BLASTn searches and NBC individually assigned >90% of the reads to fungi, indicating that this primer pair could indeed be the most fungus-specific of those tested. Of fungal classes, Agaricomycetes was the dominant group in all data sets.

The barcode-primer combinations exhibited five-fold differences in the number of fungal OTUs recovered in total and on the basis of samples rarefied to 8906 sequences (Figure 2; Suppl. material 3–6). As expected, the most conserved barcode - SSU V5 (SSU817F-SSU1196Rngs) - revealed the lowest number of OTUs, followed by LSU D3 (LR3R-LR5) and SSU V4 (S515Fngs-Euk742R) barcodes. Of rDNA genes, the LSU D1 barcode recovered greater OTU richness compared to other SSU and LSU barcodes. Overall, however, the ITS1 and ITS2 barcodes recovered the greatest number of OTUs. Both primer pairs targeting ITS2 sequences produced significantly more OTUs than those used to amplify ITS1, although the latter is considered somewhat less conserved (Nilsson et al. 2009; Wang et al. 2015). For the ITS2 barcode, the forward primer ITS3tagmix produced significantly more OTUs per sample compared with the gITS7 primer (especially Leotiomycetes and Lecanoromycetes), but their richness estimates were strongly correlated (n=25; R=0.96; P<0.001; Table 3). Similarly, of the two primers targeting the ITS1 barcode, the ITS1Fngs forward primer produced significantly more OTUs compared to the ITS1ngs primer (especially Sordariomycetes), but again these estimates were strongly correlated (n=25; R=0.83; P<0.001). OTU accumulation curves of separate and mixed ITS1 data sets indicated that combining the results of both primer pairs enabled recovery of 6.6% more OTUs, but the fraction of gITS7 data did not add to the ITS3tagmix data (Figure 3). Overall, primer pairs for the ITS1 barcode were the least correlated with other barcodes in predicting OTU richness and the median value of richness. The LSU D1 barcode and ITS2 barcode (ITS3tagmix-ITS4ngs) had the strongest average correlation with other combinations.

Figure 2. 

Sample-based OTU richness as recovered by different barcode-primer pair combinations. Error bars denote standard error; different letters indicate statistically different groups.

Figure 3. 

Rarefied OTU accumulation curves for samples based on the (a) ITS1 and (b) ITS2 barcodes and their 95% confidence intervals.

Across all samples, barcode-primer pair correlations were strongly correlated in recovering the relative abundance of fungal classes (Table 4). These correlations were slightly stronger for primer pairs targeting the same barcode. The ITS3tagmix-ITS4ngs and gITS7-ITS4ngs primer pairs exhibited the strongest correlations to the median values across all barcodes. On average, the most conserved LSU D3 (LR3R-LR5) and SSU V5 (SSU817F-SSU1196Rngs) barcodes had the distinctly lowest correlations to other barcodes and the median values of all barcodes taken together.

Ecological analyses using all barcodes consistently revealed that vegetation structure was the strongest predictor of fungal communities (Table 5). Contrary to our predictions, the most conservative barcodes SSU V5 (SSU817F-SSU1196Rngs) and LSU D3 (LR3R-LR5) indicated that floristic variables explain the greatest amount of variation in the fungal community among all barcodes. Conversely, floristic variables accounted for the lowest amount of variation among all barcodes in all four ITS data sets. We suspected that lower matrix connectance could result in relatively lower detection power in the more sensitive and larger data sets. Successive removal of infrequent species indicated that the logarithm of connectance was extremely strongly correlated with the coefficient of determination up to connectance values of 0.45 (Bray-Curtis distance: n=8; R=0.980; P<0.001; Hellinger distance: n=8; R=0.995; P<0.001). Beyond this point, the adjusted determination coefficients reached a plateau or began to decline (Figure 4).

Figure 4. 

Relationship between connectance and adjusted coefficient of determination (R²_adj) for floristic variables across different barcode-primer pair combinations based on (a) Bray-Curtis distance and (b) Hellinger distance. Pointed line indicates correlation in the ITS1Fngs-ITS2 data set (filled circles), covering eight connectance classes (C<0.45). Open circles, other ITS1 and ITS2 primer pairs; triangles, SSU barcodes; rectangles, LSU barcodes.

In the ITS1 barcode data set, the choice of primers (ITS1Fngs vs. ITS1ngs) explained 2.0% and 1.8% of the community variation based on Hellinger distance (partial analysis: F_1,33=3.12; P<0.001) and Bray-Curtis dissimilarity (F_1,33=3.12; P<0.001), respectively. This amount of variation was marginal compared to the effects of vegetation (16.0–20.0%) and edaphic parameters (5.5–6.4%) but comparable to climatic effects (2.8–3.8%). In the ITS2 barcode data set, the primer pair (ITS3tagmix vs gITS7) accounted for <1% variation in fungal community composition based on both distance measures (P<0.001); this was negligible compared to the effects of vegetation (19.6–22.2%), spatial vectors (5.5–6.9%), soil variables (5.6–5.8%), and climatic predictors (3.3–3.4%). The class-level analysis of all nine barcode-primer pair combinations revealed that primer pair was the strongest overall predictor, explaining 38.1% of the community variation (partial analysis: F_8,182=38.3; P<0.001; Figure 5), overshadowing the environmental effects (altogether 22.7%).

Figure 5. 

Global Nonmetric Multidimensional Scaling (NMDS) graph demonstrating the relative placement of samples (lower case letters, encoded in Suppl. material 1) in the ordination space. 95% confidence ellipses are indicated for each barcode-primer pair combination. For two-dimensional solution, stress=0.191 (R²=0.875).

Of the 290,779,313 high-quality metagenome sequences from the PNG soil samples, 1,309,342 (0.45%) were assigned to ribosomal DNA of prokaryotic and eukaryotic organisms. Bacterial sequences and eukaryote sequences unassigned to any kingdom dominated the rDNA subset. Only 16,833 (1.29%) of these sequences were determined to represent fungal nuclear rDNA. Across all samples and regions, Agaricomycetes (41.4%), Eurotiomycetes (10.4%), and Dothideomycetes (5.5%) dominated, whereas 13.2% sequences could not be assigned to any fungal class (Suppl. material 7). The fungal class Mixiomycetes was only recovered based on metagenome analysis, but many other uncommon groups in the amplicon data set were absent from the metagenome data set, presumably due to the comparatively low number of fungal metagenome sequences. The relative abundance of fungal classes in the amplicon-free soil metagenome was significantly correlated with that of all barcodes; the two most conserved barcode-primer pair combinations exhibited the weakest correlations (Table 4).

Within the soil metagenome, there were substantial differences in the recovery of fungal classes based on SSU, ITS, and LSU (Figure 6). Both in the metagenome and amplicon data set, 11 out of the 16 most common fungal classes exhibited significant differences in recovery by the rDNA markers. In spite of PCR and primer biases, the taxonomic composition of fungi was more similar based on the same marker using either metagenomics- or amplicon-based HTS. Chytridiomycota and Pezizomycetes were significantly less common in the metagenome compared to the amplicon-based HTS data, but Leotiomycetes and the unassigned category were more abundant in the metagenomics data set (Figure 6).

Figure 6. 

Relative abundance of fungal classes in the amplicon and metagenomics data sets divided into SSU, ITS, and LSU subsets averaged over different barcodes (amplicon data) and 14 shared samples. Asterisks in the margins indicate significant differences in recovery of classes among SSU, ITS, and LSU of metagenomics (right) and amplicon (left) data sets. Asterisks in the center indicate significant differences between the metagenomics and amplicon-bases approaches.

Metagenomic analysis of the mock community enabled us to estimate class-level biases in identification of fungi. While other groups were roughly evenly represented, Archaeorhizomyces spp. accounted for 17.0%, 5.1%, and 10.1% of metagenomics sequences in the SSU, ITS, and LSU data subsets, respectively, which is consistent with trends in the metagenomics and amplicon data sets of soil samples. Across all three markers, 6.2% and 2.4% of sequences respectively belonged to the unknown category or were affiliated with classes that were not included in the mock community (Suppl. material 7). The proportion of unknowns and incorrectly identified sequences were 4.3 and 2.5 times greater for the LSU compared with ITS data subset. The proportion of incorrect identifications is probably much greater, because our mock community represents sequences from the most common classes and misidentifications among these classes cannot be inferred. Therefore, use of positive controls representing a mock community and single isolates is warranted for evaluation of identification biases in metagenomics data sets.

In-depth analysis of primer mismatches to fungal templates revealed potential systematic biases inherent to different primer pairs (Appendix 1). The fungus-specific ITS1F primer family has central mismatches to multiple fungal groups as well as 3’ terminal mismatches to Chytridiomycota, Saccharomycetes, and several taxa within Dothideomycetes. The “universal” ITS1 primer has two central mismatches to most Sordariomycetes. These two primers complement each other very well in terms of taxon coverage and avoidance of introns, except for the fact that both miss most Microsporidia. The ITS2, ITS3, and gITS7 primers exhibit single mismatches to multiple fungal groups especially within the early diverging lineages. The use of a mixture of multiple modified ITS3 variants accounts for the vast majority of these mismatches (Appendix 1; Tedersoo et al. 2014). Due to extreme variation in the 5.8S rDNA, many lineages in the Tulasnellaceae family are missed by all these primers. The same thing goes for the Microsporidia, whose unusual configuration of the ribosomal genes makes this a very challenging group to target in molecular ecology efforts. The ITS4 and LR0R primer families possess two mismatches to Archaeorhizomycetes and one mismatch to most Chaetothyriales, the latter of which was accounted for by a degenerate position in the ITS4ngs and LR0Rngs primers. Both primers are ill-suited for amplification of a large proportion of Tulasnellaceae and Microsporidia due to multiple central mismatches (Taylor and McCormick 2008; Oja et al. 2015). The LF402 and LF402F primers discriminate among most non-fungal groups and cover most of the fungal kingdom, except for a number of taxa within Cantharellus, Tulasnellaceae, Ustilaginaceae, Candida, and some early diverging lineages. These can be accounted for by adding degenerate positions or constructing primer mixes (Appendix 1). The LR3R, TW13, and LR5 primers are truly universal, with perfect match to nearly all eukaryotes. Some Saccharomycetes taxa exhibit a single central mismatch to all these primers. The SSU primer SSU515F is nearly universal across the three domains of life. However, Microsporidia and Archaeorhizomycetes exhibit a single central mismatch to this primer. This mismatch can be accounted for by a degenerate position in the original variant of this primer (Turner et al. 1999). The original primer contains an intron site common to many lineages in Saccharomycetes and Pezizomycotina. Shortening the SSU515Fngs primer served to remove the intron position from the priming site in this study. The Euk742R and SSU817F primers possess multiple mismatches to several fungal groups, especially Cantharellus, Tulasnellaceae, and Saccharomycetes but also to certain other groups as well. Therefore, these primers are not recommended for fungal metabarcoding. The SSU1196R primer has multiple (including terminal) mismatches to Cantharellus, Tulasnellaceae, Saccharomycetes, some groups of Pezizomycetes, and Glomeromycota. Shortening the primer as in SSU1196Rngs ameliorated the mismatches only in part. Of other commonly used primers, the ITS5 primer is universal to most eukaryotes with a single central mismatch to some fungal and non-fungal groups. The LR5F primer has a perfect match to nearly all fungal groups, animals, and stramenopiles, but at the same time it exhibits a terminal mismatch to plants. Of basidiomycete-specific primers, the ITS4-B has multiple mismatches to nearly half of the major orders, but it effectively excludes other fungal phyla and non-fungal groups. The LB-W primer covers nearly all Basidiomycota (except Cantharellus and Ustilaginomycetes) and early diverging lineages, but it discriminates strongly against Ascomycota.

Fungal taxa differed roughly three-fold in the length of the ITS1 and ITS2 barcodes (Figure 7). The Archaeorhizomycetes exhibited both very short ITS1 and ITS2 barcodes, with the full ITS region averaging 370 bp. On average, Mortierellomycetes, Glomeromycetes, Tremellomycetes, and Wallemiomycetes had a disproportionately short ITS1 barcode compared with the ITS2 barcode. Except for Wallemiomycetes, the abundance of these groups was lower based on the ITS1 barcode in spite of smaller amplicon size. However, critical mismatches to the ITS1Fngs and ITS2 primers (Appendix 1) may have caused relatively lower recovery of Glomeromycota and Tremellomycetes based on the ITS1.

Figure 7. 

Differences in sequence length in the ITS1 and ITS2 barcodes of 16 most abundant fungal classes as revealed based on amplicon libraries in this study. Columns, asterisks, and error bars represent mean and median values and standard deviation, respectively. Numbers inside bars indicate the number of sequences analyzed (n). Taxa are ordered by average length of the ITS1 region.

Our analyses of seven barcodes indicate that markers differ substantially in their ability to recover OTUs at the 97% sequence similarity threshold, a threshold value that is almost universally used in HTS studies. Consistent with the lower species-level discrimination power of SSU and LSU compared with ITS (Schoch et al. 2012), both SSU and LSU barcodes recovered fewer OTUs compared to ITS1 and especially ITS2 barcodes. Because the 97% sequence similarity threshold is considered too low for separation of closely related species for the ITS region (Kõljalg et al. 2013), the true taxonomic richness is probably underestimated with all barcodes used here, but relatively more so with SSU and LSU. Clustering of the mock community sequences indicated that the relatively more species from different fungal orders are inadvertently lumped together based on sequences of conservative domains of SSU and LSU (Suppl. material 2). Although OTU richness estimates of different barcode-primer pair combinations exhibited strong or moderate linear correlations, we caution against using SSU and LSU for metabarcoding. These genes may nonetheless be suitable for class or phylum-level analyses and inference of community phylogenies when functional assignment of taxa is not of concern.

In spite of the low taxonomic resolution of SSU and LSU, these barcodes were relatively more efficient in recovering trends in community composition in terms of greater proportion of variance explained. There are two alternative and perhaps additive explanations to this observation. First, phylogenetic niche conservatism among fungi may reinforce this pattern (Cavender-Bares et al. 2009) as taxonomic richness of fungal classes is driven by different environmental variables such as soil pH and climate (Tedersoo et al. 2014). Second, lumping of closely related rare species increases the connectance of community matrices, i.e. results in greater matrix fill. Lower connectance associated with relatively more taxon-rich data sets adversely affects the performance of many community association metrics (Olesen and Jordano 2002). We demonstrated that greater connectance may improve statistical properties of multivariate analyses by reduction of noise. Importantly, qualitative aspects of these analyses were not altered as no other variable beside vegetation became significant with removal of infrequent OTUs. These results suggest that possibilities for optimal performance of multivariate analyses of large and sparse data sets warrant further exploration. It is important to assess the effects of connectance and removal of rare OTUs separately, because these two effects were not disentangled in our study.

Although all barcode-primer pair combinations revealed that floristic variables account for the strongest effects in fungal community composition, there were statistically significant primer biases that were not reported in previous studies (Arfi et al. 2012; Daghino et al. 2012; Blaalid et al. 2013). At the level of OTUs, the choice of primer pairs for the ITS1 barcode explained 2% of the community composition, which is less than any environmental predictor group. Tedersoo et al. (2010) found that selection of either ITS5 or ITS1F forward primers explains 15% of the community variation in ectomycorrhizal roots. Differences in performance can be explained by primer-template mismatches of the ITS1Fngs, ITS5, and ITS1ngs forward primers (Appendix 1) and by the common occurrence of an intron between ITS1Fngs/ITS5 and ITS1ngs primer sites in many ascomycete groups (Bhattacharya et al. 2001). Lower taxonomic resolution or mismatches to primers in the 3´ end of the original ITS2 primer compared with optimized gITS7 and ITS3tagmix primers may have contributed to lower taxonomic richness of the ITS1 compared with ITS2 barcode. Taxa falling into Basidiomycota (except Agaricomycetes), Glomeromycota, Rozellomycota, Dothideomycetes, and Geoglossomycetes were disproportionately more OTU-rich in the ITS2 data set (Tedersoo et al. 2015b).

Currently, the Illumina HiSeq technology enables generation of 4×10⁸ DNA sequences per run. In our soil samples, fungal nuclear rDNA represented ca. 0.005% of all DNA molecules, resulting in an average sequencing depth of approx. 1000 sequences per sample (n=16). In spite of relatively even representation of all metagenomics sequences across samples, the coefficient of variation for the number of fungal sequences was roughly three times higher for metagenome samples compared with amplicons. Such uneven sequencing depth further complicates downstream statistical analyses.

PCR and primer biases in the amplicon data sets are well addressed, but there are also certain biases inherent to metagenomics approaches that are related to base composition and replication of DNA fragments (Gomez-Alvarez et al. 2009). Because the metagenomics sequences exhibit very low overlap across the rDNA, it is impossible to assign these sequences to OTUs and recover taxonomic richness (cf. Bengtsson et al. 2012).

Except for some minor but statistically significant taxonomic biases, the metagenomics data set covering SSU, ITS, and LSU provided highly comparable results to that of all barcodes taken together but especially to the results from ITS1 and ITS2. The metagenomics analyses confirmed that the Agaricomycetes, Eurotiomycetes, and Sordariomycetes are the key players in tropical soils of PNG. Furthermore, the ITS data set of both amplicons and metagenome exhibited a similar proportion of mycorrhizal fungi (Tedersoo et al. 2015b). We caution that some care should be taken in interpreting our comparison of metagenomes and amplicons, because both analyses were performed in a single run, such that there may be differences between the two Illumina models, MiSeq 2×300 and HiSeq 2×250, and their chemistry.

Given the low and uneven recovery of fungal rDNA sequences and difficulties in correct taxonomic assignment (see below), metagenomics with the sole purpose of metabarcoding is clearly a waste of financial and computational resources. Enrichment of targeted molecules such as mitochondrial DNA may improve the problems associated with insufficient sequencing depth (Zhou et al. 2013). Nonetheless, vast amounts of functional data produced in the metagenomics analyses promise a new paradigm in community ecology (Zhou et al. 2015).

Our in silico analysis of primer mismatches extends the results of Bellemain et al. (2010) and Op De Beeck (2014) by covering many additional primers and pointing to particular mismatching positions and taxa. This has enabled us to fine-tune primers for targeting the ITS region in fungi and other soil eukaryotes for cross-kingdom studies (Tedersoo et al. 2015a). We also demonstrate that there is no universally good solution for covering all fungi, because groups with rapidly evolving nuclear rDNA such as Microsporidia, certain groups within Saccharomycetes, Tulasnellaceae, and Cantharellus commonly escape detection by these primers for all of SSU, ITS, and LSU. Similarly, excessively long ITS barcodes and introns within barcodes may hamper amplification and sequencing of additional groups (Lindahl et al. 2013; Tedersoo et al. 2015b). Besides primer bias, unequal representation of amplicons may result from extreme GC to AT ratio (Aird et al. 2011). Although there was no evidence for length bias in our study, Ihrmark et al. (2012) demonstrated that barcode length has a strong inverse correlation with recovered abundance in a mock community based on 454 pyrosequencing technology. In particular, Archaeorhizomycetes stand out as a group with exceptionally short ITS1 and ITS2 barcodes. Because of the very strong size selection in certain cloning methods and PCR (Kanagawa 2003; Carninci et al. 2011), the unprecedented dominance of Archaeorhizomycetes in molecular studies based on cloning of long DNA fragments (Schadt et al. 2003; Castro et al. 2010) may represent an extreme example of a cloning bias.

The technical biases of sample preparation steps are poorly understood and these may be largely specific to platforms, models, and chemistry. Nonetheless, depending on the base composition and size of DNA fragments, unequal competition for adaptors may occur (Aird et al. 2011). Besides potential competitive disadvantage in ligation, the relatively longer reads in the 300–600 bp frame are less likely to be recovered because of reduced read quality in the last quarter of the sequence and problems with pairing low-quality sequences.

In addition to PCR and primer biases inherent to amplicon-based analyses, a bias related to incompleteness of reference database and uncertainty of taxonomic assignment is common to both amplicon and metagenomics data sets. A part of this so-called “identification bias” results from differential quality and abundance of reference data that may affect the probability of identification of particular taxa (Kõljalg et al. 2013). For example, reference sequences covering the full LSU are available for most early diverging lineages and the most common groups of Ascomycota and Basidiomycota, but not for “minor” groups such as Orbiliomycetes, Geoglossomycetes, and Archaeorhizomycetes. Cryptomycota and Microsporidia were poorly represented in the ITS and LSU reference data sets, because taxonomic studies in these groups are largely focused on the SSU. Similarly, many other early diverging lineages such Zoopagomycota, Neocallimastigomycota, and Blastocladiomycota were underrepresented in the ITS reference data sets. Although the early diverging lineages represented <0.5% relative sequence abundance in all data sets, many of the OTUs not assigned to phyla probably represent one of these lineages. Another aspect of the identification bias was evident when comparing taxonomic assignments of reference databases and results from BLASTn vs. Naïve Bayesian Classifier. Combining the databases and methods enabled us to identify a vast majority of OTUs and metagenome sequences to class level. The identification bias is very difficult to quantify or explicitly test among different barcode data sets and among studies, because of confounding primer effects, PCR conditions as well as subjectively established similarity thresholds and manual re-evaluation based on previous experience and knowledge about sequence variability.

Differential representation of taxa in the reference data sets certainly plays a much greater role in metagenomics data sets targeting all genes (Korsakovsky Pond et al. 2009) or only ribosomal DNA, because the rDNA sequences used for identification may fall into any part of the SSU, ITS, or LSU, or into intergenic spacers (IGS) that are usually not considered. Since the SSU, ITS, and LSU all include variable as well as highly conserved regions, it is very difficult to provide correct taxonomic assignments to short fragments (100-450 bases) at lower taxonomic levels. For several fungal groups, full-length sequences of LSU are lacking, which results in partial ability to identify these. Typically, fungal orders can be distinguished based on specific similarity thresholds and NBC classifications for ITS and variable domains of SSU and LSU. However, orders, classes, and even phyla can be lumped together at the same levels of similarity if the metagenomic sequence happens to cover a highly conserved region. Furthermore, our metagenomics rDNA data comprised a remarkable proportion of SSU and LSU sequences that clearly belonged to the ophistokonts, but could not be reliably linked to animal and fungal kingdoms or minor lineages related to these groups. Therefore, the uncertainty and potential mistakes in taxonomic assignment are much more relevant to the metagenomics approach compared to DNA barcodes. To ameliorate identification biases in metagenomics data, reference databases should be updated with full-length rDNA sequences originating from genomics initiatives such as the 1000 Fungal Genomes Project (http://1000.fungalgenomes.org) and covering all classes of fungi and closely related organisms.