|Home | About | Journals | Submit | Contact Us | Français|
Microbes commonly exist in milieus of varying complexity and diversity. Although cultivation-based techniques have been unable to accurately capture the true diversity within microbial communities, these deficiencies have been overcome by applying molecular approaches that target the universally conserved 16S ribosomal RNA gene. The recent application of 454 pyrosequencing to simultaneously sequence thousands of 16S rDNA sequences (pyrotags) has revolutionized the characterization of complex microbial communities. To date, studies based on 454 pyrotags have dominated the field, but sequencing platforms that generate many more sequence reads at much lower costs have been developed. Here, we use the Illumina sequencing platform to design a strategy for 16S amplicon analysis (iTags), and assess its generality, practicality and potential complications. We fabricated and sequenced paired-end libraries of amplified hyper-variable 16S rDNA fragments from sets of samples that varied in their contents, ranging from a single bacterium to highly complex communities. We adopted an approach that allowed us to evaluate several potential sources of errors, including sequencing artifacts, amplification biases, non-corresponding paired-end reads and mistakes in taxonomic classification. By considering each source of error, we delineate ways to make biologically relevant and robust conclusions from the millions of sequencing reads that can be readily generated by this technology.
Microbes in nature typically exist in complex communities, which can be highly variable in both the composition and the abundance of their constituent species. Initial surveys of microbial community diversity applied cultivation-based methods, but these were often biased and of limited value owing to the inability to cultivate many (more likely, the majority of) naturally occurring species. Such methods were supplanted by molecular approaches that did not rely on cultivation but were based rather on the PCR amplification, cloning and Sanger sequencing of universally conserved molecules, usually the 16S ribosomal RNA gene. Current sequence databases contain over a million full-length 16S rRNA sequences spanning a broad phylogenetic spectrum than can serve as a benchmark for assessing the bacterial taxa (also referred to as ‘phylotypes' or ‘ribotypes') present in environments worldwide (Cole et al., 2009).
Although the traditional 16S-cloning-and-sequencing approach has the potential to provide an in-depth view of the richness and evenness of bacterial species within a community, its application is somewhat laborious and costly, with the result that most applications have assayed on the order of only 100 sequences per sample, which may not be sufficient to fully characterize all but the simplest communities. This has led to the development of alternative techniques for assessing rRNA variation (for example, ARISA, DGGE, tRFLPs) (Muyzer et al., 1993; Liu et al., 1997; Fisher and Triplett, 1999); however, the most revolutionary innovation has been the application of 454 pyrosequencing to survey hundreds of thousands of 16S rRNA gene amplicons (termed ‘pyrotags') in a single sequencing run (Sogin et al., 2006). Since its initial introduction, pyrotag analysis has been widely implemented (for example, Dethlefsen et al., 2008; Bowers et al., 2009; Cardenas et al., 2010; McLellan et al., 2010) and further optimized in order to obtain thousands of sequences from multiple samples in parallel (Parameswaran et al., 2007; Hamady et al., 2008; Costello et al., 2009; Engelbrektson et al., 2010).
Advancements in sequencing technologies, by offering greater numbers of sequencing reads at much lower costs, might eventually supplant pyrotag analysis for the characterization of microbial communities. Recent interest has focused on the applicability of other sequencing methodologies, most notably Solexa/Illumina (Lazarevic et al., 2009; Claesson et al., 2010; Gloor et al., 2010; Caporaso et al., 2011; Zhou et al., 2011), which is currently less than 1/100 the cost per read than 454 pyrosequencing. Although realizing much shorter read lengths, the Illumina technology can be tailored to yield sequences of increased lengths, as can be obtained by merging the paired-end reads generated from the same amplicon (Gloor et al., 2010; Rodrigue et al., 2010; Zhou et al., 2011). By integrating sample-identifying barcodes into the amplification primers, the Illumina platform, like 454 pyrosequencing, is amenable to a high level of multiplexing, which further increases its utility for examining large and complex sets of samples (Gloor et al., 2010).
Here, we use the Illumina sequencing platform for 16S rRNA amplicon analysis (called ‘iTags') and assess its generality, practicality and potential complications. We first fabricate and sequence a paired-end library of hyper-variable 16S rDNA fragments amplified from samples that varied in contents from a single bacterium, an artificial community and a highly complex natural community. The inclusion of identifying barcodes and multiple primer pairs, combined with the selected amplicon lengths, allowed us to evaluate the level of sequencing error, amplification biases and association between paired-end reads, among other variables. Despite recovering over 30 million reads, our imposition of stringent quality filters resulted in up to an 85% reduction in the number of potentially informative reads. Nevertheless, tens to hundreds of thousands of informative reads were returned for each of the multiplexed samples. These results uncovered both the major advantages and obstacles in the application of these methods, and provide insights into the ways that such impediments can be ameliorated prior to the application of these procedures.
During the development of this project, the Illumina GAIIx sequencer could produce reads of up to 100 nucleotides (nt) in length from both ends of single DNA fragments or amplicons (that is, paired-end reads). We exploited this technology to assemble reads spanning a >100-nt region of the 16S rDNA. To accomplish this, we identified candidate universal primer pairs spanning the 16S rRNA variable loops (Neefs et al., 1991) that would maximize sequence information and be short enough to allow adequate read overlap.
The phylogenetic distribution of each primer sequence was analyzed using the RDP Probe Match (Cole et al., 2009). The goal was to obtain broadly distributed primers that amplify regions of approximately 100–150bp in length such that the paired-end reads from an amplicon would show some overlap. Of primer sets meeting these criteria, two flanking the variable loop-6 (V6) were chosen (XXX denotes sample-specific barcodes): 967F 5′-XXX CAACGCGAAGAACCTTACC-3′ and 1046R 5′-XXX CGACAGCCATGCANCACCT-3′ (Short 1 V6 (S1V6); 98bp in Escherichia coli) (Sogin et al., 2006), and 917F 5′-XXX GAATTGACGGGGRCCCGC-3′ and 1061R 5′-XXX CACGRCACGAGCTGACGAC-3′ (Long V6 (LV6); 163bp in E. coli) (Keijser et al., 2008). We further modified the S1V6 primer pair to increase both the extent of overlap between paired-reads and the number of potentially informative sites as follows: 970F 5′-XXX CGCGAAGAACCTTACC-3′ and 1050R 5′-XXX ACGACAGCCATGCANC-3′ (Short 2 V6 (S2V6); 96bp in E. coli).
Sample-identifying barcode sequences were included at the 5′ ends on both the forward and the reverse primers. The barcodes were designed (1) to be short, 3 or 4bp, to minimize loss of sequence information; (2) to have two or more differences from one another (such that a single sequencing error could not convert one into another); and (3) to have the 3′ nucleotide of the barcode be a low-frequency match of nucleotide upstream from the 5′ end of the priming site. The barcodes for the primer pairs were as follows: S1V6–ATG, CAG, TCG, ACT, CGT, GCA, GAT, and TGA; S2V6–AGC, ATT, GTC, TCC, CCAT, CGGC, GCGT and TTAT; and LV6–GTC, AAC, CCA, TAA, AGGA, GCGC, TTGC and ATAC.
We tested 16S primers on sets of samples representing a range of microbial diversity. First, we used two well-defined samples containing either DNA extracted from a single sequenced strain of E. coli K-12 (MG1655) (Blattner et al., 1997) or from a mixture of genomic DNAs of known concentrations from 19 cultivated strains of bacteria (Table 3). Next, we analyzed a sample of unknown microbial diversity: a fecal sample from a laboratory mouse. The QIAamp DNA stool mini kit (Qiagen, Valencia, CA, USA) was used following the manufacturer's protocol to extract total DNA from the feces of an individual lab-reared mouse, Mus domesticus strain WSB.
PCR amplifications were performed in triplicate 30-μl reactions using 30–50ng of template DNA and one of two DNA polymerases. Samples amplified with the barcoded S1V6 primer pairs were performed using 1.2U of Taq DNA polymerase (5 PRIME) and final concentrations of 1X polymerization buffer, 1m dNTPs and 0.4μ of each primer. The reaction mixtures were subjected to an initial denaturation cycle of 95°C for 2min, followed by 25 cycles at 95°C for 30s, 55°C for 45s and 72°C for 1min, and a final extension step at 72°C for 10min. The samples amplified with the barcoded S2V6 or LV6 primer pairs used the high-fidelity Advantage 2 polymerase mix and the Advantage 2 SA PCR buffer (Clontech, Mountain View, CA, USA) following the manufacturer's instructions. Amplifications with the S2V6 and LV6 primer pairs were performed in a laminar-flow hood using screw cap tubes to minimize contamination. The thermal cycling conditions were identical to those detailed above except that the annealing temperature was 58°C.
Triplicate reactions for each sample were pooled, visualized on 1% agarose gels, purified with the MinElute PCR purification kit (Qiagen) and eluted in 11μl of EB buffer. The concentrations of the amplification products were quantified on a Bioanalyzer using the DNA 1000 LabChip kit (Agilent Technologies, Santa Clara, CA, USA). The barcoded products of the S1V6 and S2V6 primer sets were combined in equimolar amounts, re-purified with the MinElute kit, eluted in 30μl of EB buffer and submitted for 2 × 75-nt paired-end sequencing at the Yale Center for Genomic Analysis. Similarly, the LV6 products were mixed in equimolar amounts, but these longer amplicons were subjected to 2 × 100-nt paired-end Illumina sequencing.
Initial base-calling using Bustard (v 1.7) and quality filtering of paired-end reads were performed by the Yale Center for Genomic Analysis. Upon receipt, these reads were evaluated for quality (mean quality scores; presence of sites denoted as ‘N'), and identifiable primers and barcode sequences. Read pairs with perfectly matching primer and barcodes were extracted and binned by barcode. The resulting Illumina fastq files, containing quality values equivalent to Phred scores, were converted to Sanger fastq files using Maq (‘maq ill2sanger', v 0.7.1) (Li et al., 2008). Individual read pairs were converted to fasta and qual files, and then assembled using Phrap (v 1.080812) using the following parameters: ‘-vector_bound 0 -trim_start 0 -forcelevel 3 -preassemble -bandwidth 10 -repeat_stringency 0.98 -ace -minmatch 5 -maxmatch 20 -minscore 10'.
To investigate the effect of quality scores on estimates of diversity, the assembled reads for each barcoded sample were further filtered by consensus quality scores. Barcodes were pruned from sequences and the data sets were aligned in the RDP pyrosequence aligner (Cole et al., 2009), which implements the Infernal algorithm (Nawrocki and Eddy, 2007). The resulting alignments were edited to remove primer sequences and clustered to generate operational taxonomic units (OTUs) using the RDP clustering algorithm. The taxonomic assignments of OTUs were derived from the RDP Classifier using a 50% bootstrap cutoff, as recommended for reads less than 100bp. Custom PERL scripts were written to generate a pipeline to analyze, organize and format these data. The PERL scripts and raw sequence data presented and analyzed here are available at www.yale.edu/ochman/data.
Over 30 million paired-reads were collected from each of the two multiplexed Illumina sequencing runs. Initial filtering removed 85% of the reads from run1 and 54% from run2 owing to errors stemming from several sources (Table 1), which is comparable to other studies that report filtering 40–70% of their initial reads (Gloor et al., 2010; Caporaso et al., 2011). We secured a total of 4.5 and 16.0 million paired-end reads for the two runs in spite of the high number of reads that were discarded because they did not meet our thresholds for accuracy; this, in turn yielded tens of thousands to millions of raw read-pairs per multiplexed sample (run1 =380253±176319; run2 =3936966±2299439) (Table 2). Because the runs were founded with equimolar concentrations of amplicons, some of the variation in the number of reads returned per sample is most likely attributable to amplification biases that occurred during the synthesis of the Illumina sequencing library. Future experiments might avoid this bias by implementing an amplification-free library approach (Kozarewa et al., 2009).
We used a strict set of quality filters to minimize the impact of erroneous reads on measures of microbial diversity (Kunin et al., 2010). Errors in the barcode/primer sequence of either of the paired-end reads necessitates the removal of both reads in the pair; and in the majority of read-pairs that were culled, either only one of the reads perfectly matched the primer sequence or the paired reads lacked identical barcodes. The primer and barcodes are not informative for sequence classification and together constitute over 20% of the sequencing read, but the inclusion of both sequences provide a useful filter for removing erroneous reads and estimating sequencing error rates.
Error frequencies were estimated for the barcode/primer regions of all reads with high sequence quality (q30) for the first 25nt and with single-nucleotide differences from the expected sequence. This included reads initially filtered out for not containing a perfect primer or barcode sequence (see above). Across the 14 barcode/primer combinations analyzed from run1, the error rates per nucleotide varied over three-fold (1.6 × 10−4 to 5.6 × 10−4), and rates varied only two-fold (4.1 × 10−4 to 8.4 × 10−4) among the four barcode/primer combinations from run2 (Supplementary Table 1). For two reasons, we ascribe these variants to errors generated during sequencing rather than to artifacts that arose during PCR amplification: (1) Oligonucleotide synthesis is typically very accurate, such that errors in primer regions are more likely generated during the sequencing step; (2) we performed independent experiments using two different Taq polymerases each having unique error rates, but obtained approximately the same rate (analysis of covariance; F1,15=1.07, P=0.317) and spectrum of mutations with each. The fact that error rates were independent of the source of Taq polymerase further suggests that errors are introduced during DNA sequencing.
There were large numbers of paired-end reads with mismatched barcodes, but owing to our use of barcodes that require at least two mutations to convert one barcode into another, not many of these mismatched read-pairs were attributable to sequencing errors. In fact, many read-pairs contained recognizable barcodes that differed between the corresponding reads, which could result from (1) over-clustering during the sequencing run; (2) the production of chimeric amplicons during the amplification step of the Illumina library or (3) primer contamination, the likelihood of each of which is discussed below.
Over-clustering can result from low nucleotide diversity in the first four sequenced nucleotides and high cluster density on the flow cell. Although the first four bases of our barcodes and primers did not contain equal frequencies of each of the four nucleotides, over-clustering is an unlikely cause of the mismatched barcodes detected. An analysis of the run showed that read signal intensities and the percentage of clusters passing filter were both above normal expectations (intensity >200 for first cycle, >80% clusters passing filter). Alternatively, chimeric amplicons are a common by-product of 16S rDNA amplifications (Qiu et al., 2001); however, our protocols used a low number of amplification cycles (10–12) and the high-fidelity Taq polymerase for library construction, both of which reduce the frequency of chimeric amplicons. This leaves primer contamination during the initial PCR amplifications as the principal contributor of mismatched barcodes. Although we implemented precautionary measures intended to minimize contamination (using screw cap tubes; removing primers and templates individually; setting up the reactions in a laminar flow hood), contamination occurred in these samples. With current interest in multiplexing scores (or even hundreds) of samples, we let these results serve to forewarn investigators that highly stringent procedures are necessary to prevent cross-contamination.
We first evaluated the efficacy and accuracy of the iTag approach by sequencing a single template, E. coli K-12 (MG1655). The E. coli genome contains seven rRNA operons, six of which have identical alleles at V6 and the seventh encoding a V6 allele with nine nucleotide differences (Figure 1a). At both 100% and 97% identity thresholds only two OTUs are expected.
The V6 loop of E. coli was amplified using three primer pairs (S1V6, S2V6 and LV6), and we performed a technical replicate of S2V6 (Table 2). Nearly all reads assembled with well-supported overlaps of their distal 3′ regions (S1V6 and S2V6, ≈50bp; LV6, ≈30bp). Alignment and clustering of the assembled iTag sequences showed that each of the E. coli samples contained much more than the two predicted phylotypes (n=4090–159224; Table 2).
Previous work has shown that quality trimming of pyrotags is essential for accurate prediction of community composition (Kunin et al., 2010); therefore, we filtered the assembled reads according to consensus quality scores estimated by Phrap. We compared two quality filters, requiring that all sites had error probabilities of either 10−2 or 10−3 (corresponding to Phred scores of q20 or q30). These measures reduced the number of unique phylotypes (aka 100% OTUs), but failed to remove all of the spurious OTUs (Table 2).
Among these confounding OTUs, we found that 3–7% of high-quality reads (q30) differed by only 1nt from the expected E. coli alleles (Figures 1b–e). Variants occurred at every nucleotide position regardless of the primer pair or the source of Taq polymerase. Given the numbers of variants, it is not likely that each represents a new mutation that occurred during growth in culture. Therefore, these phylotypes are attributable to errors that occurred during sequencing.
The very large number of sequence variants made them impossible to cluster into the two expected OTUs, even when applying very relaxed thresholds (for example, 90% OTUs). Most of the remaining OTUs (Figure 1) represent phylotypes that were detected in other samples and are best explained by primer contamination (discussed above). In addition, some of the contaminants in the LV6 samples were because of nonspecific amplification of other chromosomal loci; for example, we recovered a 123-bp fragment of uridylyltransferase (glnD) from E. coli and a 78-bp fragment of an ATP transport permease from Streptococcus sp.
To rid samples of artifactual phylotypes generated by sequencing errors and PCR contamination, we used a sequence abundance threshold in a manner similar to that of Caporaso et al. (2011). We based this threshold on the representation of a unique phylotype (100% OTU) given the number of q30 reads for that barcode. Through stepwise increments in the abundance threshold, we followed the degree to which spurious reads were eliminated (Figures 1 and and2).2). Eventually, by applying a 1.0% abundance threshold (that is, removing all OTUs occurring at frequencies less than 1.0%) we were able to recapitulate the two phylotypes expected to occur in the E. coli genome.
Having addressed the sequence variation detected within a single species by the iTag approach, we analyzed the data generated for a synthetic community of 19 strains of bacteria (Table 3). This sample was amplified using each of the three primer sets, which should yield 19 100% OTUs, and 16 (S1V6, S2V6) or 15 (LV6) 97% OTUs. A similar pattern showing an excess of erroneous and contaminating high-quality OTUs emerged even after stringent quality filtering (Table 3). Although the application of a 1% abundance threshold removed artifactual reads, it also removed several of the expected OTUs that occurred below this threshold (Figure 2). The initial template DNA concentrations of 12 of the 19 strains represented less than 1% of the total DNA mass, so this result was anticipated (Table 3); however, we also found that the actual read abundances for the 19 100% OTUs departed from the expected values (Figure 3). The likely cause of the deviation between the relative amount of input DNA and the relative number of reads could be primer specificity. For example, the Streptococcus cristatus OTU was similarly abundant in the CAG–S1V6 (1.7%=2754/162125) and ATT–S2V6 (1.8%=10207/556184) data sets, but nearly absent from the AAC–LV6 data set (0.0022%=9/407785).
Short 16S tags can be difficult to assign taxonomically and alone provide limited phylogenetic resolution (Liu et al., 2007). We found that the RDP Classifier accurately identified all of the 19 expected OTUs to taxonomic Class, but its success rate progressively decreased for lower taxonomic divisions, with only about 60–90% assigned correctly. Alternative BLAST-based classification approaches can increase accuracy provided related full-length 16S rRNA sequences are available (Huse et al., 2008). Of course, neither of these approaches are able to distinguish among species that have identical V6 loop sequences (for example, E. coli, Escherichia fergusonii and Salmonella enterica; Figure 3).
By way of analyzing a complex microbial community whose composition and constituents are uncharacterized, we assessed the V6 sequence diversity within a fecal sample from the house mouse M. domesticus strain WSB/EiJ. This sample was amplified with GCA–S1V6 and TCC–S2V6, generating, respectively, 384480 and 223708 binned reads. As before, even with rigorous quality filtering, a large number of low-frequency 100% OTUs remained (Table 2). We applied similar abundance thresholds (10−5 to 10−2), which resulted in a dramatic reduction in phylotypes (Figure 4). The utilization of the same threshold as that applied to the E. coli samples (that is, 10−2) leaves only 12 OTUs, which represent 67% (GCA–S1V6) and 74% (TCC–S2V6) of the high-quality reads. Thus, it appears that application of a 1.0% threshold is too stringent for actual biological samples and probably removes many tags that are likely to represent the actual microbial constituents of these communities. The two amplified samples show an extensive overlap in 100% OTUs (Sørensen distance, 0.010; Jaccard distance, 0.020), which results in a highly similar taxonomic distribution regardless of the abundance threshold applied (Figure 4).
Massively parallel sequencing platforms supply an efficient and sensitive means of characterizing microbial communities, but the sequencing depths offer by these methods are accompanied by artifacts that can affect the interpretation of community diversity. We have identified and quantified several sources of errors inherent to 16S amplicon sequencing (iTags) on the Illumina platform (contamination, amplification and sequencing errors), thereby enhancing the utility of these data for assessing community diversity. Most recent efforts to interrogate microbial communities have relied upon 16S rRNA gene sequences generated by 454 pyrosequencing (for example, Dethlefsen et al., 2008; Bowers et al., 2009; Costello et al., 2009; Cardenas et al., 2010; McLellan et al., 2010); however, Illumina sequencing technologies currently produce much larger numbers of much shorter reads for a fraction of the cost and will likely supplant the use of 454 pyrotags.
Our experimental iTag design used barcoded primers flanking the V6 segment of 16S rRNA, a region commonly amplified in pyrotag experiments (Sogin et al., 2006; Keijser et al., 2008; Galand et al., 2009). Individual samples were amplified, mixed and then used as templates to construct and sequence two standard Illumina paired-end libraries. Although Illumina sequence quality decays along the length of reads, our motivation for assembling paired-end reads was to increase the quality and confidence of the overlapping region (Gloor et al., 2010; Rodrigue et al., 2010; Zhou et al., 2011). Therefore the amplified V6 region was restricted to 100–160bp to ensure an adequate overlap of the forward and reverse paired-end reads. We subsequently filtered out low-quality tags, leaving only assembled reads that contained barcode and primer sequences at each end, which bounded a phylogenetically informative amplified region. This design allowed the simultaneous estimation of sequencing error rates and microbial diversity.
Over half of >30 million paired-end reads from each of Illumina runs were subsequently discarded owing to either (1) sequencing errors in one or both of the primer regions or (2) paired-end reads that contained mismatched barcodes (Table 1). Other Illumina-based 16S rRNA studies have encountered similarly high sequencing error rates, making such extreme read filtering necessary, but these studies did not report a high incidence of mismatched barcodes (Gloor et al., 2010; Caporaso et al., 2011). After ruling out over-clustering and sequence chimeras, primer contamination during the initial sample amplifications was viewed as the most likely cause of the mismatched barcodes (see above). For example, >60% of the contaminants detected in the E. coli libraries were present in one or more of the other samples, and frequently these contaminants represented the more abundant phylotypes in the other samples (data not shown). The sequencing depth afforded by the Illumina platform greatly increases the likelihood of detecting contaminants; however, we stress the utility of barcoding and sequencing both ends of multiplexed amplicons in order to detect this source of errors. In the absence of paired barcodes, removing reads that lack a significant match to a reference database can be implemented. Although this latter method is of limited utility when examining novel or poorly characterized communities.
To establish baseline levels of sensitivity of the iTag approach, we performed several controls, including the amplicons generated from a single bacterial genome and from a defined mixture of DNA from 19 strains. In both cases, we found an excess of high-quality, 100% phylotypes regardless of the primer pair, Taq polymerase or sequencing run. However, sequencing reads bearing a single mutation relative to the two expected E. coli phylotypes, or to the 19 phylotypes from the mixture of strains, represented 2.2–8.2% of the reads. Based on these reads we estimate a sequencing error rate of 1.3 × 10−3 to 4.1 × 10−4 per nucleotide (Supplementary Table 2), which is slightly greater that the rate estimated from errors confined to the primer regions (run1 =3.8 × 10−4; run2 =6.4 × 10−4). This difference is in line with expectations as, in the Illumina platform, error rates are known to increase with distance from the priming site (Rodrigue et al., 2010; Zhou et al., 2011). These errors impact taxonomic identification using the RDP Classifier, particularly reducing the ability to assign Order, Family and Genus identities to the shorter S1V6 and S2V6 tags.
To account for the large number of spurious reads in these samples, we applied strict abundance thresholds to our samples. We came the closest to recapitulating our expected diversity from these samples when we applied an abundance threshold at which all phylotypes at frequencies <1% were removed. This cutoff is based on the total number of high-quality reads for a given sample and more-or-less corresponds to the 0.01% threshold applied by Caporaso et al. (2011), whose threshold is based on the total number of read-pairs recovered for their entire experiment (10000 out of 87507177 reads). Application of such stringent thresholds is not particularly satisfying, but will remain a necessary measure until sequencing error rates are improved.
Mammalian guts house a diverse array of resident microbes and are numerically dominated by the Firmicutes and Bacteriodetes (Ley et al., 2008). A previous pyrosequencing effort identified 1000 bacterial phylotypes (97% OTUs) in healthy Mus musculus individuals using the same S1V6 primers (Antonopoulos et al., 2009); however, our sample from a lab-reared, wild-derived inbred strain of M. domesticus WSB yielded only 12 phylotypes (97% OTUs) after application of the equivalent quality and abundance filters derived from the contrived samples (q30 and 10−2). Reducing the stringency of the abundance threshold increased the number of recovered iTag phylotypes, but did not alter the taxonomic distribution of the iTags (Figure 4). (Lowering the abundance threshold from 10−2 to 10−4 yielded about 20% more iTags containing only a single-nucleotide difference relative to the 12 high abundance OTUs.) The source of the difference in the number of phylotypes detected in the Mus fecal samples is unknown, and direct comparisons between iTags and pyrotags can only be assessed from the corresponding analyses of 16S diversity within E. coli (Kunin et al., 2010; this study). Therefore, we suggest the amplification and sequencing of a control sample (for example, E. coli) during each multiplex experiment to empirically determine appropriate abundance thresholds and to account for differences in sequencing error rates between runs.
Interest in leveraging the Illumina sequencing platform to examine 16S rRNA gene diversity has led to the recent development of several strategies and observations about its implementation (Claesson et al., 2010; Gloor et al., 2010; Caporaso et al., 2011; Zhou et al., 2011). Until recently, the Illumina sequence read length was limited to ≈100nt. As a result, paired-end sequencing of single, short variable loops (V6) has been assembled (Gloor et al., 2010; Zhou et al., 2011; this study), and longer variable region(s) (255–465bp) were artificially merged or analyzed separately (Claesson et al., 2010; Caporaso et al., 2011). Assembly of overlapping, paired-end reads results in phylotypes of considerably better quality because it identifies and rectifies problems that arise from the deterioration of read quality as a function of read length (Claesson et al., 2010; Gloor et al., 2010; Zhou et al., 2011). But because PCR and sequencing errors persist, application of abundance thresholds are still necessary to avoid inflated estimates of species richness (Gloor et al., 2010; Zhou et al., 2011). Moreover, conservative estimates of error rates using contrived DNA samples vary six-fold, from 2.4 × 10−3 to 4.1 × 10−4 per nucleotide (Zhou et al., 2011; this study).
The iTag approach holds enormous promise as means to investigate microbial diversity, but in its current applications there are three major obstacles: (1) contamination, (2) the utility and classification of short read lengths, and (3) sequencing error rates. Contamination can likely be mitigated by changes in laboratory procedures, including use of liquid-handling robots and sterile work environments. Although not considered here, contamination derived from reagents or consumables is also possible and would require alternative strategies for its detection and elimination. Despite the prevalence of studies using the V6 region and evidence that short reads perform adequately for community analyses (Liu et al., 2007), short V6 tags appear to systematically overestimate species richness (Youssef et al., 2009). Current Illumina read lengths (≈150nt) will allow the recovery and assembly of larger (V4) or more (V6+V7) 16S variable regions that better reflect the microbial diversity obtained when analyzing the entire 16S rRNA molecule. Additionally, longer read lengths will facilitate the use of additional barcodes that will allow multiplexing hundreds of samples in a single run. Although the observed error rates are of little consequence to most genomic studies, as errors are resolved through sequencing depth, amplicon tagging approaches treat each unique sequence read as a novel OTU. Until there are improvements to the Illumina sequencing chemistry and analysis programs, erroneous reads will need to be filtered through the application of abundance thresholds (Gloor et al., 2010; Caporaso et al., 2011). The cost and capacity of the Illumina sequencing make it the preeminent platform for assessing microbial community diversity, and we have called attention to the errors that need to be resolved in its application and in the interpretation of results.
We thank Michael Nachman, Nancy Moran and Zakee Sabree for contributing the template material used for library construction, and Zakee Sabree for coining the term ‘iTag'. Funding for this project was provided by the National Institutes of Health to HO.
Supplementary Information accompanies the paper on The ISME Journal website (http://www.nature.com/ismej)