Investigating Microbial Eukaryotic Diversity from a Global Census: Insights from a Comparison of Pyrotag and Full-Length Sequences of 18S rRNA Genes

Ten samples were collected from five different marine environments at various depths (5 to 2,500 m) between August 2000 and November 2005 as part of a Global Protistan Survey (Table 1). Water samples were collected using Niskin bottles mounted on a rosette equipped with conductivity, temperature, and depth (CTD) sensors.

Water samples of 2 to 20 liters were prefiltered through 200-μm and 80-μm mesh to remove most metazoa and were then filtered (<10 mm Hg) onto 45-mm GF/F (Whatman) filters. The filters were placed in 2 ml of 2× lysis buffer (100 mM Tris [pH 8], 40 mM EDTA [pH 8], 100 mM NaCl, 1% SDS) and were stored frozen (−20°C) on shipboard or flash frozen in liquid nitrogen for later DNA extraction. DNA was extracted according to the procedures detailed by Countway et al. (15).

Pyrosequencing of the V9 hypervariable region of the 18S rRNA gene was performed according to the procedures described by Amaral-Zettler et al. (29). The primers used to sequence the V9 region were 1380F (5′-CCCTGCCHTTTGTACACAC-3′), 1389F (5′-TTGTACACACCGCCC-3′), and 1510R (5′-CCTTCYGCAGGTTCACCTAC-3′). The emulsion PCR was performed using standard Roche protocols, and sequencing was carried out on the Genome Sequencer FLX system (Roche, Basel, Switzerland) using the GS LR70 long-read sequencing kit (Roche).

Sequence reads that did not have (i) an exact match to the proximal primer (1380F or 1389F) and (ii) the presence of the distal primer (1510R) were removed, as were reads that had one or more ambiguous bases (N′s). Further quality filtering was done using the free software package mothur, version 1.30.0 (http://www.mothur.org/), by following the standard operating procedures for 454 sequences (30). These procedures included removing pyrotags with an average quality score of <35 using a 50-bp sliding window and those with >8 homopolymers present. Preclustering with a 1-bp mismatch allowance was performed to minimize the effects of sequencing or amplification errors on OTU clustering (31), and chimeras were removed using the UCHIME algorithm (32). Sequences are available at the Visualization and Analysis of Microbial Population Structure (VAMPS) website, hosted by the Marine Biological Laboratory (VAMPS.mbl.edu), and in the NCBI Sequence Read Archive under SRA number SRP001225.

Full-length 18S rRNA gene sequences were amplified by PCR using the universal eukaryote primers Euk-A (5′-AACCTGGTTGATCCTGCCAGT-3′) and Euk-B (5′-GATCCTTCTGCAGGTTCACCTAC-3′). The PCR and cloning procedures are based on those of Countway et al. (33) with slight modifications, including the use of AmpliTaq Gold for PCR amplification, and amplicons were cloned using a TOPO-TA kit (Invitrogen). Sequencing of the full-length 18S rRNA gene was conducted by the Joint Genome Institute (JGI, Walnut Creek, CA) using two vector primers (T3 and T7) and an internal primer (Euk-570F).

Sequences that did not have both the proximal (Euk-A) and distal (Euk-B) primers, or that could not be assembled due to failure of the internal sequencing primer, were removed from the data set. Sequences were then passed through the pintail algorithm to remove chimeras (34).

The mothur software package (version 1.30.0) was also used for calling OTUs for both the pyrotag and full-length sequence data sets (average neighbor method) to ensure that the OTU-calling procedures for the two sequence data sets were comparable. An aligned SILVA eukaryotic full-length 18S rRNA gene reference database (version 102; aligned and provided by the mothur software package) was used to aid in the alignment of full-length and pyrotag sequences using mothur prior to OTU calling. All full-length sequences were pooled, and OTU calling was performed at a sequence similarity of 97% based on the observation that OTUs called using full-length sequences at 97% by mothur yielded a number of OTUs most comparable to that called by the Microbial Eukaryote Species Assignment (MESA) program at 95% sequence similarity (20; S. K. Hu and D. A. Caron, unpublished data). Calling OTUs at 95% sequence similarity in MESA was designed to create OTUs of full-length or nearly full length 18S rRNA gene sequences at approximately species-level distinctions using a data set of well-curated and morphologically well described protistan species and their sequences (20).

The pyrotag and full-length sequence data sets were first compared by examining the number of OTUs formed by the two data sets. Pyrotag sequences were randomly subsampled to yield the same number of sequences in each sample as the full-length data set (i.e., to standardize to the full-length sequence data set). This subsampled pyrotag data set was used to call OTUs at different sequence similarities (95 to 100%) using mothur, and the resulting numbers of OTUs were compared to the numbers of OTUs called for the full-length sequences using mothur at 97% sequence similarity.

The analysis described above allowed the comparison of OTU formation between the pyrotag and full-length sequence data sets using the same number of sequences, in order to avoid the potential bias inherent in using a much larger pyrotag sequence data set. One of the strengths of NGS, however, is the substantially greater number of sequences that can be generated from a sample. Subsequent comparisons between the pyrotag and full-length sequence data sets, therefore, included all pyrotag sequences. Richness for the two data sets (all full-length or all pyrotag sequences) was examined through the construction of rarefaction curves for both data sets. Rarefaction curves of the observed numbers of OTUs were generated for OTUs called at 97% sequence similarity for both data sets and additionally at 98% and 99% similarity for the pyrotag data set. Rarefaction curves were constructed with and without the inclusion of singletons in the two data sets.

The best BLAST result for a representative sequence from each OTU using the SILVA database (version 111) was used to provide taxonomic information for pyrotag and full-length sequence OTUs. The representative sequence from each OTU was selected by mothur and was the sequence with the smallest distance to all other sequences within the OTU. If more than one sequence matched this condition, then the sequence with the smallest average distance to other sequences in the OTU was selected.

All pyrotag sequences from the 20 pyrotag OTUs that had the most pyrotag sequences were mapped to the full-length sequences in order to investigate how an OTU generated with pyrotag sequences related to OTUs generated with full-length sequences. The pyrotag sequences were mapped to full-length sequences using the Burrows-Wheeler Aligner algorithm with no mismatch allowance (35), and the number of full-length OTUs matched perfectly by the pyrotag sequences in a single pyrotag OTU was tallied. In the event that a pyrotag sequence mapped perfectly to multiple full-length sequences (i.e., the full-length sequences had identical V9 regions), one of these matching full-length sequences was randomly selected as the perfect match to the pyrotag sequence (i.e., each pyrotag sequence matched with only one full-length sequence for subsequent analysis).

The results of nonmetric multidimensional scaling (nMDS) analysis using the pyrotag and full-length sequence data sets were compared in order to investigate potential differences in community similarities inferred from the two types of sequences. The number of sequences in each data set was standardized separately for the two data sets (450 for the full-length data set and 4,736 for the pyrotag data set) by randomly subsampling the sequences in each sample down to the number of sequences in the sample with the least number of sequences for each data set. This standardization procedure was performed to avoid the potential bias of having different numbers of sequences in each sample. In addition to the analyses using the full-length and pyrotag sequence data sets that included singletons, a third analysis was performed on the pyrotag data set with the singletons removed (yielding a data set of 4,472 pyrotags per sample). All subsampling of sequences was performed using mothur. The standardized data sets were square-root transformed to down-weight highly represented OTUs. The PRIMER software package (version 6) was used to calculate the Bray-Curtis community composition similarity values that were based on the OTU composition of each sample. A matrix of pairwise Bray-Curtis similarity values was constructed for each data set (i.e., full-length sequences, pyrotag sequences, and pyrotag sequences without singletons) and was then used to perform CLUSTER (group average mode with the SIMPROF test for significance) and nMDS analyses. The results of the SIMPROF significance test from the CLUSTER analysis were overlaid on the nMDS plots, which are visual representations of the similarities among communities.

The full-length 18S rRNA gene sequences analyzed in this study have been submitted to GenBank under accession numbers KJ757035 to KJ759741, KJ760393 to KJ762454, and KJ762829 to KJ764638.

We compared the number of OTUs formed from pyrotag sequences called at different sequence similarities (95 to 100%) to the number of OTUs formed from full-length sequences called at 97% sequence similarity in order to investigate which sequence similarity used to call OTUs for the pyrotag data set resulted in a number of OTUs most comparable to that generated by the full-length sequence data set. The pyrotag data set for this analysis was standardized by randomly subsampling the number of pyrotags in each sample to match the number of full-length sequences in each sample (Table 1) in order to avoid potential biases from using a much larger pyrotag data set. OTU calling for all full-length sequences (6,579 sequences) at 97% similarity in mothur resulted in a total of 1,400 OTUs, while the same number of pyrotag sequences called at various sequence similarities, ranging from 95 to 100% sequence similarity, resulted in 1,113 to 1,695 OTUs (i.e., ≈80 to 120% of the number of OTUs called for the full-length sequences) (Fig. 1). OTU calling of the pyrotag sequences at 98% sequence similarity yielded a total of 1,419 OTUs, which was most similar to the number of full-length OTUs called at 97%. Singleton OTUs (i.e., OTUs composed of only one sequence in the data set) constituted a large portion (>50%) of the OTUs in all treatments of both data sets. Calling OTUs using full-length sequences at 97% sequence similarity generated the highest proportion of singletons (72% of the total number of OTUs) compared to the number of singletons generated from the standardized pyrotag sequence data called at different levels of similarity (range, 54 to 64% of the total number of OTUs). The pyrotag singletons in this analysis included some “secondary” singletons that were generated as a consequence of the subsampling process. The use of the standardized data set yielded an average of 59% singleton OTUs over the range of sequence similarities used to call OTUs (Fig. 1), whereas calling of OTUs using the entire pyrotag data set yielded an average of 51% singleton OTUs.

Rarefaction curves were constructed using all sequences available for the pyrotag (190,696 sequences) and full-length (6,579 sequences) sequence data sets in order to examine sequence coverage (Fig. 2). Full-length OTUs were called at 97% similarity, while pyrotag OTUs were called at 97%, 98%, and 99% sequence similarities. Calling pyrotag sequences at 98% sequence similarity resulted in a number of OTUs that was most comparable to that from the full-length sequence data set in the previous standardized comparison (see the preceding section), and OTU definitions of 97% and 99% were included in order to investigate how differences of 1% in sequence similarity would affect the analysis results. The entire pyrotag data set generated a large number of OTUs (at 97% sequence similarity, 6,414 OTUs; at 98%, 7,775 OTUs; at 99%, 11,234 OTUs), while the full-length sequence data set generated 1,400 OTUs, as noted above. The rarefaction curve for the full-length sequence data set was similar to the pyrotag rarefaction curves for the portion of the curves where the data sets (i.e., number of sequences sampled) overlapped (Fig. 2B). The curve for the full-length sequences appeared to have a slightly different trajectory than those for the pyrotag sequences (Fig. 2B), but the curve for the full-length sequences was always bracketed by the curves generated by the pyrotag OTU data sets called at 97%, 98%, and 99% sequence similarities.

The removal of singleton OTUs altered the form of the rarefaction curves dramatically. The effects of removing singletons were investigated because there have been suggestions that some singleton OTUs may be artifactual (31, 36, 37) (for more detail, see Discussion). The number of OTUs observed decreased substantially for both data sets with the removal of singletons, as expected (Fig. 2C and D). The number of OTUs observed for the pyrotag data set was reduced to 40 to 52% of the total number of OTUs when singletons were removed (from 11,234 to 4,530 OTUs at 99% sequence similarity; from 7,775 to 3,828 OTUs at 98%; from 6,414 to 3,349 OTUs at 97%). The number of OTUs observed in the full-length sequence data set without singletons was reduced to 28% of the total number of OTUs (from 1,400 to 387 OTUs). The rarefaction curve for the full-length data set without singletons had a much more shallow slope and a lower maximum than the pyrotag rarefaction curves for the same number of sequences sampled (Fig. 2D).

The higher-level taxonomic compositions (supergroup level) of the entire pyrotag and full-length sequence data sets were tallied in order to (i) compare the taxonomic compositions obtained by calling pyrotag OTUs at different levels of similarity, (ii) compare the taxonomic compositions of the pyrotag and full-length sequence data sets, and (iii) examine the effect of the removal of singletons on the taxonomic information inferred using both data sets. This analysis focused on pyrotag OTUs called at 97%, 98%, and 99% sequence similarities for the reasons given in the preceding section.

Calling OTUs using the pyrotag data set at three different levels of sequence similarity had very little effect on the higher-level taxonomic composition of the sequences (Fig. 3). Similarly, the presence or removal of singleton OTUs did not substantially influence the taxonomic composition derived from either the full-length or the pyrotag sequence data set. However, the full-length sequence data set yielded a larger proportion of sequences identified as rhizarians than the pyrotag sequence data set (14 to 15% of full-length sequences and ≈3% of pyrotag sequences) and smaller proportions of sequences identified as stramenopiles (1 to 3% of full-length sequences and 12 to 13% of pyrotag sequences) and as cryptophytes and haptophytes (≈0% of full-length sequences and ≈3% of pyrotag sequences). Other groups (e.g., alveolates) constituted similar proportions of the two data sets.

All pyrotag sequences from a single pyrotag OTU were mapped to full-length sequences at 100% sequence similarity in order to investigate how individual pyrotag OTUs related to full-length OTUs. Full-length OTUs for the analysis were called at 97% sequence similarity. This analysis focused on 20 pyrotag OTUs that had the most pyrotag sequences. Pyrotag OTUs were called at 95%, 97%, and 99% sequence similarities (Table 2).

Pyrotag sequences from nearly all of the 20 pyrotag OTUs with the most sequences mapped to full-length sequences from multiple full-length OTUs. The results were similar among the three sequence similarities examined for the pyrotags, but lower sequence similarity tended to hit more full-length OTUs. In only three cases did the pyrotags from a well-populated pyrotag OTU match a single full-length OTU (the 17th, 18th, and 19th most populated OTUs when pyrotag sequences were called at 99%, 97%, and 95% sequence similarities, respectively). The pyrotag OTUs with the most sequences formed at 95%, 97%, or 99% sequence similarity consisted of high numbers of pyrotags (34,787 to 49,709), and pyrotags from these three OTUs mapped perfectly to full-length sequences that were distributed among 127, 109, or 75 full-length OTUs, respectively (Table 2). While sequences from a single pyrotag OTU aligned perfectly to sequences from multiple full-length OTUs, the full-length OTUs were generally from the same higher-level (e.g., class or order) taxonomic group (e.g., Dinoflagellata) (Table 3).

nMDS plots were constructed using the pyrotag and full-length sequence data sets to (i) compare community similarities among samples using pyrotag and full-length sequences (Fig. 4A and B) and (ii) examine the effect of the removal of singletons on the clustering of communities for the pyrotag data set (Fig. 4B and C). OTUs for the pyrotag data set were called at 98% sequence similarity, while OTUs for full-length sequences were called at 97%, for the reasons provided above. Each data set was standardized to the number of sequences in the sample with the lowest number of sequences (for the full-length data set, 450 sequences/sample; for the pyrotag data set, 4,736 sequences/sample; for the pyrotag data set without singletons, 4,472 sequences/sample). The Bray-Curtis similarity values estimated for the full-length data set were generally lower (average, 5%) than those estimated for the pyrotag data sets (average with singletons, 14%; average without singletons, 13%), but the overall patterns and the clustering of some communities in the nMDS plots were similar for the three data sets (Fig. 4). For example, communities from the Ross Sea (RS) consistently clustered together, and communities from deep water at the East Pacific Rise (EPR) (sampling depths, 1,500 m and 2,500 m) consistently clustered together as well in both the pyrotag and full-length sequence data sets. One difference between the pyrotag and full-length sequence data sets concerned the community from the Arctic Ocean (AO) collected at 35 m. The two Arctic samples (35 and 500 m) were not significantly different (P, >0.05 by the SIMPROF test) in the full-length data set (Fig. 4A) but were significantly different in the pyrotag data sets. There were also some differences in the clustering of communities collected from the Gulf Stream (GS), EPR, and Eastern North Pacific (ENP) among the three data sets as indicated by the results of the SIMPROF significance test (i.e., subclusters on the nMDS plots), although these samples occupied the same positions on the plots in all three analyses.

The removal of singletons from the pyrotag data set had little effect on the resulting nMDS plot except that the community collected from the ENP at 5 m was statistically similar to the community collected from the GS at 15 m in the nMDS plot constructed with the pyrotag data set after the removal of singletons (Fig. 4C). On the other hand, the community collected from the ENP at 5 m was not statistically different from the community collected from the EPR at 20 m in the nMDS plot constructed with the entire pyrotag data set (Fig. 4B).

Pyrotag OTUs that were formed at a higher sequence similarity (98%) yielded a similar number (∼1,400 OTUs) as full-length OTUs called at 97% sequence similarity (Fig. 1). Our results indicate that the utility of the V9 hypervariable region of the 18S rRNA gene to adequately distinguish eukaryotic taxa may be less than the taxonomic resolution afforded by the entire gene. Dunthorn et al. (38) concluded that the genetic distances among species of ciliates that were estimated using the full-length 18S rRNA gene were more similar to the genetic distances estimated with the V4 region than to those estimated with the V9 region, similarly suggesting that the V9 region of the 18S rRNA gene may not be the ideal choice for distinguishing between ciliate species. Youssef et al. (24) also reported that some hypervariable regions of the 16S rRNA gene underestimated bacterial richness in comparison to full-length sequences called at the same level of sequence similarity.

The choice of a sequence similarity for calling OTUs heavily influenced the overall species richness observed for the entire pyrotag data set (Fig. 2A). A 2% difference in sequence similarity for OTU calling resulted in an almost 2-fold difference in the number of OTUs formed. The numbers of OTUs formed from the pyrotag data set called at 97% and 99% sequence similarities were 6,414 and 11,234, respectively, when all sequences were included for OTU formation (Fig. 2A). The removal of singletons greatly reduced the overall number of OTUs observed, as well as the difference in the number of OTUs observed between the two sequence similarities (3,349 OTUs for 97% sequence similarity and 4,530 OTUs for 99%). Stoeck et al. (23) also observed that a small difference in the sequence similarity used for calling OTUs resulted in a large difference in the number of OTUs formed. Caution is therefore needed in assessing the ecological information contained in sequence data sets, because the sequence similarity level used for calling OTUs can significantly affect the species richness. There is currently no consensus on a single sequence similarity for distinguishing between species, even for a closely related and well-defined taxonomic group such as the ciliates (39), let alone for all the different groups of protists.

Singletons contributed significantly to the number of OTUs formed in all treatments of the data set (Fig. 1) and heavily influenced the observed species richness (Fig. 2). Singletons would also heavily influence the estimation of total species richness using richness estimators (40–42). Questions regarding the potentially artifactual nature of pyrotag singletons have been raised previously, and a few studies employing collections of known taxa as artificial communities have demonstrated that some singletons were sequences that were not a part of the original community (31, 36, 37). Dickie (41) noted that the number of sequencing artifacts for pyrosequencing increased as the number of pyrotags increased. Therefore, some studies have removed singletons to avoid overestimating species richness, despite the risk of removing real taxa present at very low abundances (see, e.g., references 22, 43, 44, and 45). While the removal of singletons resulted in asymptotes of the rarefaction curves in this study (Fig. 2), the removal of singletons did not always lead to asymptotes (22, 44).

Full-length sequences are less likely to contain a significant percentage of sequencing artifacts that would lead to the formation of artifactual singleton OTUs than are very short sequences (e.g., sequences with ∼100 bp in this study). Nonetheless, singleton OTUs still contributed significantly to the full-length sequence data set (72%) (Fig. 1). This result would argue that a substantial “rare biosphere” existed within the data set generated from these globally distributed samples, regardless of the nature of the pyrotag singletons. We also compared rarefaction curves generated from the pyrotag data set with singletons removed to the curve generated from the full-length sequence data set with singletons included (data not shown). The curve for the full-length sequence data set had a steeper slope and indicated a higher observed richness than the pyrotag curves for a given number of sequences sampled. Taken together, these results imply that the pyrotag data set, with singletons removed, would result in a lower number of species observed (assuming that most full-length singletons represented valid OTUs). We speculate that the true observed richness of the samples (as constrained by the current methodologies for filtering, DNA extraction, and sequencing) most likely lies between the full-length curve with singletons and the pyrotag curves without singletons.

Singletons had a large effect on observed OTU richness, but they had a negligible effect on the overall taxonomic composition for both data sets (Fig. 3). This result implies either that the singletons were real and had proportionally the same taxonomic composition as the entire data set or that they were artifactual sequences generated at proportions similar to those of the populated OTUs from which they were generated. Regardless of whether singletons are artifactual sequences or sequences from rare organisms, they are not biased toward specific or unique groups.

OTUs obtained using the pyrotag data set at different sequence similarities also did not lead to substantive changes in higher-level taxonomic composition (Fig. 3). Therefore, the coalescing of pyrotag OTUs as the imposed sequence similarity was relaxed resulted in OTUs generally converging to the same high-level taxonomic groups. On the other hand, the taxonomic composition differed between the full-length and pyrotag data sets (Fig. 3). The full-length sequence data set yielded a higher proportion of rhizarians and smaller proportions of stramenopiles, cryptophytes, and haptophytes. This difference in taxonomic composition was not due to differences between the samples or DNA extraction procedures, because both the full-length data set and the pyrotag data set were sequenced from the same DNA extracts. These differences are more likely due to differences in the primers used for the two types of sequencing approaches (46). Other studies comparing full-length and pyrotag sequences have also demonstrated differences in the relative proportions of certain groups within the community, with certain groups or species detectable only by either Sanger sequencing or pyrosequencing (26, 27). In the present study, the full-length sequence data set yielded a very small proportion of cryptophytes and haptophytes (0.38% from the data set with singletons; 0.25% from the data set without singletons), possibly due, in part, to biases against haptophytes in the primers used for the full-length sequences (11, 47).

We also compared the taxonomic compositions of the two data sets at a higher taxonomic resolution (i.e., genus level) for sequences identified as ciliates (data not shown). The results were similar to those at the supergroup level, in which the presence/absence of singletons and the sequence similarity threshold used to call pyrotag OTUs had little effect on the genus composition. Certain genera were detected only by one of the two sequencing methods, but unsurprisingly, the pyrotag data set detected a much larger number of genera (43 to 58 genera for pyrotag sequences called at 97 to 99% sequence similarity, with and without the presence of singletons) than the full-length sequence data set (25 genera for full-length sequences called at 97% sequence similarity with singletons; 8 genera without singletons). Despite the discrepancies in the numbers of genera identified by the two data sets, both data sets had the most sequences identified as members of the genus Strombidium. Similarly, Santoferrara et al. (28) compared the investigation of ciliate diversity using Sanger sequencing, pyrosequencing, and microscopy and found that the dominant taxa obtained by the three methods were generally comparable.

Pyrotags from a single pyrotag OTU that had many sequences generally mapped at 100% sequence similarity with full-length sequences associated with multiple full-length OTUs (Table 2). This result could be explained by an inability of the short pyrotags to resolve many OTUs or by the possibility that the OTUs formed from full-length sequences were overestimating diversity. We speculate that the former situation is the more likely explanation. The majority of full-length sequences with identical V9 regions were found to be placed in the same full-length OTU (called at 97% sequence similarity) in the present study, but there are rare cases in which full-length sequences with exactly the same V9 region can be <90% similar across the entire gene (data not shown). Similarly, Sogin et al. (48) randomly compared the pairwise distances of 5 × 10⁶ full-length 16S rRNA (bacterial) gene sequences and found that >90% of full-length sequences with identical V6 hypervariable regions were at least 95% similar across the entire gene. This suggested that ∼10% of full-length sequences with identical V6 regions can be >95% different across the gene and that a single hypervariable region might not contain sufficient information to differentiate many OTUs. If this speculation is true, our results indicate that estimates of community richness obtained from pyrotag-based OTUs should be viewed with caution.

Despite this caveat, full-length sequences from multiple OTUs that were matched perfectly by pyrotags from a single pyrotag OTU were generally from the same higher taxonomic group (e.g., Dinoflagellata [Table 3]). This was the case despite the fact that pyrotags from a single pyrotag OTU matched as many as 127 full-length OTUs (Table 2). Thus, investigations of the higher taxonomic composition (Fig. 3 presents an example) should be less affected by the level of resolution afforded by the pyrotags. Huse et al. (49) performed an in silico analysis that mapped V3 and V6 hypervariable regions to the full-length 16S rRNA gene and obtained similar results in that 99% of the short sequence tags matched full-length sequences from the same class, order, or phylum level.

The number of full-length OTUs matched by pyrotags in a single pyrotag OTU decreased gradually as the number of pyrotags in each pyrotag OTU decreased (Table 2). This finding may simply reflect the fact that less-populated pyrotag OTUs had fewer sequences to match to full-length OTUs. Our analysis yielded only 3 pyrotag OTUs (one for each sequence similarity threshold) among the 20 most populated pyrotag OTUs that contained pyrotags aligning perfectly to only a single full-length OTU each. There were also a few pyrotag OTUs (one for each sequence similarity threshold) that did not match to any full-length sequences at 100% sequence similarity. The latter situation is not surprising given the limited number of full-length sequences in the data set and our requirement for perfect sequence similarity between a pyrotag and a full-length sequence.

Despite differences in the absolute numbers of OTUs obtained from the pyrotag and full-length sequence data sets, community analysis using Bray-Curtis similarity values and nMDS plots indicated that the conclusions derived from the two data sets were generally similar (Fig. 4). Liu et al. (50) demonstrated that short sequences (100 bp) were able to provide similar resolution for the comparison of communities as full-length 16S rRNA genes using UniFrac (a analysis of the distances between two environments based on the evolutionary history inferred from phylogenetic trees) in an in silico study of short sequences extracted from near-full-length 16S rRNA gene sequences. In addition, the removal of singletons in our data set had little effect on the results of the nMDS plots (Fig. 4), indicating that this beta diversity analysis may be robust toward the different sequencing approaches (Sanger versus NGS), as well as different processing of sequences (with or without singletons). Egge et al. (51) also showed that different quality filtering procedures for pyrotag sequences of the V4 hypervariable region of the 18S rRNA genes of haptophytes had little effect on the clustering of communities on nMDS plots.

Next-generation sequencing approaches offer tremendous potential for conducting studies of microbial ecology, but only if sequence data can be translated into generally accepted, taxonomically meaningful data. To date, the accuracy of these approaches (i.e., their ability to provide accurate taxonomic characterization of a microbial assemblage) has not been adequately investigated (50, 51). Community analysis using pyrotags of approximately 100 bp of the V9 region of the 18S rRNA gene in the present study was shown to yield overall higher-level taxonomic information similar to that yielded by an analysis of full-length 18S rRNA gene sequences from the same samples and also provided similar patterns of community similarity among samples (e.g., using the Bray-Curtis similarity index) (Fig. 4). However, OTUs called using short pyrotag sequences generally matched to multiple OTUs called using full-length sequences, confusing the meaning of OTUs generated using these two approaches. We conclude that the V9 region of the 18S rRNA gene is insufficient for distinguishing between low-level taxonomic groups (i.e., species) and may be more appropriate for identifying higher taxonomic levels, such as phyla or classes. Our results indicate that at present, caution is warranted when one is attempting to use NGS sequence data sets of the V9 hypervariable region alone to characterize community richness across the whole domain of microbial eukaryotes. Microbial ecologists are moving toward the goal of adopting a molecular taxonomy for protists that incorporates molecular data into holistic assessments of microbial diversity and community composition (7), but vetting these approaches and the resulting conclusions continues to be essential.