Ssu Rrna Analysis Essay

Previous SectionNext Section


Microbial phylogeny, identification, and evolution studies were revolutionized by the introduction of small-subunit (SSU) rRNA analysis 25 years ago (1), and with the advent of PCR and high-throughput sequencing, community structure studies now are commonplace (2–5). The growing sizes of SSU rRNA gene databases provide a rich ecological and phylogenetic context for SSU rRNA gene-based community structure surveys (6, 7). However, the accuracy of PCR-based amplicon approaches is reduced by primer bias and chimeras (8, 9).

Unlike gene-targeted amplicon sequencing, shotgun sequencing takes samples from the entire community by sequencing randomly sheared fragments of DNA (10, 11). Hence, while amplicon sequencing can provide far deeper coverage of SSU rRNA genes with the same amount of sequencing, shotgun sequencing may provide a more accurate characterization of microbial diversity, including functional diversity (12). In particular, shotgun sequencing may provide an improved means to detect divergent sequences not recovered by standard SSU rRNA gene primers, such as those of Verrucomicrobia, as well as eukaryotic members of the community (8, 12–14). Note that both approaches remain prone to sequencing error and bias from environmental DNA extraction (9).

The challenges for using shotgun DNA for rRNA analyses are in efficiently searching for these fragments in large sequence data sets and the subsequent analysis of the matching short reads. Several methods have been developed for SSU rRNA retrieval in large data sets (15–18), but speed improvements still are needed to match the growth in data size; moreover, none of them provide further community analysis using the identified rRNA gene sequences. In addition, traditional community analysis tools (6, 19, 20) are largely designed to handle sequences that are amplified by PCR primers. There are two primary types of community analyses: reference based (supervised) and operational taxonomic unit (OTU) based (unsupervised). The reference-based method assigns SSU rRNA gene sequences to bins based on the taxonomy of their closest reference sequences, while OTU-based methods assign overlapping gene sequences to bins based on de novo clustering with a specified similarity cutoff (e.g., 97%). The reference-based method can be applied easily to shotgun data once SSU rRNA gene fragments are retrieved (21) and several tools are available for this (22–26), but the OTU-based approach still remains challenging with shotgun data because reads are from randomly sheared fragments.

The main goal of this study is to enable unsupervised OTU-based analysis of large shotgun metagenomic data sets from soil. We improved speed and memory efficiency with a hidden Markov model (HMM)-based method, which already has been shown to be fast and accurate for SSU rRNA searches (16–18), using a well-curated and up-to-date training reference sequence collection from SILVA (7). Our unsupervised clustering method first was tested on a synthetic community with shotgun data of 100-bp reads. We next applied the method to soil data sets, where we assembled longer reads from the overlapping paired-end Illumina HiSeq reads and mapped those to 150-bp small hypervariable regions of SSU rRNA genes for de novo clustering and further diversity analysis. We retrieved and analyzed the large-subunit (LSU) ribosomal gene for confirmatory analysis. Finally, we went beyond traditional primer evaluation (in silico database search) by evaluating primer biases using the paired shotgun and amplicon data produced from the same DNA extract (27, 28).

Previous SectionNext Section


Soil samples, DNA extraction, and sequencing.Two sets of soil samples were used. The first sample, which was used to develop the method, was a bulk (non-root-influenced) soil sample (SB1) taken in 2009 from between rows of switchgrass. The method then was applied to the second sample set taken in 2012, which consisted of seven replicate rhizosphere samples from both corn (C) and Miscanthus (M) plots. All samples were from the Great Lakes Bioenergy Research Center (GLBRC) Cropping System Comparison Site at the Kellogg Biological Station in southwest Michigan ( The rhizosphere samples were closely associated with the roots (<1 mm).

DNA extraction and SSU rRNA gene amplification methods were described previously (29). The SSU rRNA gene amplicons from the first sample were sequenced by the Joint Genome Institute (JGI) in their standard work flow, which used 454 GS FLX and Titanium platforms and a primer set (926F, AAACTYAAAKGAATTGACGG; 1392R, ACGGGCGGTGTGTRC) that targeted the V6-V8 variable region of bacteria, archaea, and eukaryotes. The second set also was sequenced at the JGI but at a later time, so the Illumina MiSeq platform and a primer set (515F, GTGCCAGCMGCCGCGGTAA; 806R, GGACTACHVGGGTWTCTAAT) that targeted the V4 variable region were used. Shotgun sequencing also was done by JGI using Illumina GAII platforms for the first set and HiSeq 2500 with 250-bp insert libraries and two 150-bp reads for the second set. We had about 8 Gb of data for the first set and about 300 Gb of data each for corn and Miscanthus for the second set.

Data preprocessing.Data preprocessing is necessary for both shotgun and amplicon data due to sequencing errors. However, it is not included as part of the core pipeline, because users have their own preferences. We trimmed trailing bases with quality score 2, called the read segment quality control indicator (encoded by ASCII 66 “B” in Illumina GAII or ASCII 35 “#” in Illumina HiSeq shotgun data), and discarded reads shorter than 30 bp and with “N.” The reads then were quality trimmed with fastq-mcf (version 1.04.662) ( with “-l 50 -q 30 -w 4 -x 10 -max-ns 0 -X.” The paired-end reads overlapping by more than 10 bp were assembled into one long read by FLASH (version 1.2.7) (30) with “-m 10 -M 120 -x 0.20 -r 140 -f 250 -s 25.” Roche 454 pyrotag amplicon data were processed using the RDP Pipeline PIPELINE INITIAL PROCESS and CHIMERA CHECK tools (6). Reads were sequenced from the reverse primer end to the forward primer end. Since the targeted region is about 467 bp (926F/1392R), most reads were not long enough to reach the forward primer. Thus, only the reverse primer product was used for quality trimming. The minimum length was set to 400 bp, and defaults were used for other parameters.

Building SSU and LSU rRNA gene models.We quality trimmed SILVA (7) SSU and LSU Ref NR database (version 115) sequences by discarding all sequences with ambiguous DNA bases and converting U to T. We then clustered them at a 97% similarity cutoff using (default with UCLUST) and from QIIME (version 1.8.0) (20). We chose the longest representative in each OTU to be further clustered at an 80% similarity cutoff. We collected the longest sequence in each OTU, resulting in 4,027 representative sequences for the SSU rRNA gene and 1,295 for the LSU rRNA gene to obtain a phylogenetically diverse set of reference genes. We further grouped these sequences into two groups, one combining Bacteria and Archaea and the other containing only Eukaryota (see Discussion). Each group was used to make two HMMs (hidden Markov models), one with sequences from a previous step and the other with reverse complement, using hmmbuild in HMMER version 3.1 (31). Finally, the HMM files were concatenated into a single file for each gene. This step is not part of the pipeline, and the resulting HMMs were included in the database of this pipeline.

Identification of rRNA gene fragments from metagenomic data.The analysis framework is shown in Fig. S1 in the supplemental material, as are the reasons for our choices of pipeline components. We searched Illumina shotgun metagenomic data with hmmsearch in HMMER version 3.1 (31) using the LSU and SSU HMMs. For testing the sensitivity of newly built models, we analyzed our tool with meta-rna (16) and metaxa (18). We used an E value of 10 for hmmsearch with the newly built HMMs. Meta-rna ( (the package was not assigned a version number; the most recent version update, 21 October 2011, was used) was run with flags “-k euk,bac,arc -e 0.00001” (16), and metaxa (metaxa_x) version 2.0.2 was run with flags “–allow_single_domain 1e-5,0 -N 1 -E 1e-5.” A bulk soil (SB1), a rhizosphere soil (M1), and a synthetic community sample (12) were used as test data. We aligned the HMMER hits from the E value cutoff of 10 using the multiple-sequence aligner align.seqs in mothur (version 1.33.3) (32). For SSU rRNA gene fragments, 18,491 full-length SSU rRNA gene sequences (14,956 from Bacteria, 2,297 from Archaea, and 1,238 from Eukaryota) from the SILVA database (release 102) (7) were downloaded from the mothur website ( and used as the template with flags “threshold = 0.5” and “flip = t.” For LSU rRNA, multiple-sequence alignments (MSA) of representative sequences of the SILVA LSU Ref NR database, clustered at a 97% similarity cutoff, were used as the template with the same flags as those for SSU. Based on alignment information provided in the align.seqs output report file, those shotgun reads with more than 50% mapped to a reference gene were designated SSU rRNA or LSU rRNA gene fragments.

Testing the effect of target region size and variable region on clustering.We used shotgun data of a synthetic community comprised of 64 species, which were sequenced by the paired-end 100-bp method on an Illumina HiSeq 2000 (12). For testing the effect of target region size on clustering, we picked V4 with starting position 577 in Escherichia coli. Sizes from 50 to 180 bp with a 10-bp increment were chosen. The minimum read length was set to the target region size minus 5 bp if the region size was less than 100 bp, and it was set to 95 bp when the region size was greater than 100 bp. We used the pre.cluster command in mothur (19) with 1 edit distance to collapse reads with errors and their original reads. De novo clustering then was achieved by RDP McClust with an algorithm for unweighted pair group method using average linkages (UPGMA) and a minimum read overlapping length of 25 bp (33). We chose McClust as the clustering tool due to its speed and memory efficiency (33, 34). We chose V2, V3, V4, V5, V6, and V8, starting at positions 127, 427, 577, 787, 987, and 1227 in E. coli, respectively (35), to test the hypervariable region effect on clustering results. Target region sizes of 80 bp and 120 bp and a distance cutoff of 5% were chosen. Further, the analyses described above also were applied to 16S rRNA genes from the 64 species comprising the community to obtain the true OTU numbers.

Community structure comparison based on OTUs from de novo clustering.For the clustering analysis of shotgun and amplicon data, 150 bases corresponding to the V8 region (E. coli positions 1227 to 1377) were aligned. Reads shorter than 100 bp were removed from the alignment and the remainder clustered using McClust with a minimum overlap of 25 bp and the UPGMA method (6). The clustering result was converted to the mothur format, and community structure comparisons were done using make.shared with a label of 0.05, dist.shared(calc = thetayc), and the pcoa command. E. coli positions 127 to 277, 577 to 727, and 997 to 1147 were chosen for V2, V4, and V6, respectively, for comparisons of different regions (35). Procrustes analysis as implemented in QIIME (20) was used to transform V2, V6, and V8 PCoA (principal-coordinate analysis) results and to minimize the distances between corresponding points in V4. The bulk soil sample (SB1) was sequenced in six lanes from one Illumina plate using DNA from the same extraction. We used these as technical replicates for testing the reproducibility of de novo OTU-based analysis on shotgun data. Since sequencing depth is critical for reproducibility testing, we pooled these into two samples of three lanes each, with the first three lanes as SB1_123 and the remaining three lanes as SB1_456.

Comparison of OTU-based microbial community structures inferred from shotgun and amplicon SSU rRNA gene sequences.The abundance of each OTU in shotgun data and amplicon data (V6-V8 for SB1 and V4 for M1) from the same DNA extraction were compared to check the consistency between the two sequencing approaches. Pearson's correlation coefficient and linear regressions were used to evaluate the correlation between the two types of data and between technical replicates. All of the abundances of each OTU were increased by a pseudocount of one to allow display on a log scale (avoiding zeros).

Comparison of taxonomy-based microbial community structures inferred from shotgun and amplicon SSU rRNA gene sequences.The SSU rRNA fragments from shotgun data and amplicon data were classified using RDP Classifier (21). The reference SSU rRNA genes from RDP and SILVA are provided on the mothur website and were used as training sets, with a bootstrap confidence cutoff for classification of 50%. Representative sequences of SILVA LSU Ref NR clustered at a 97% similarity cutoff were used as a training set with taxonomy information built from the sequence file for the LSU rRNA gene. The bacterial taxonomy profiles from shotgun data and amplicon data were compared at the phylum level.

Copy number correction.We used the SSU rRNA gene copy number database in CopyRighter (36), which provides the copy number for each taxon in the Greengenes database (37). In the taxonomic summary, the abundance of each taxon is weighted by the inverse of its SSU rRNA gene copy number. Similarly, in OTU-based analysis, the abundance of each OTU is weighted by the inverse of SSU rRNA gene copy number of its consensus taxon. Unclassified sequences are weighted by the inverse of the average copy number of all taxa in the data set.

Accession numbers.The amplicon data for C1 to C7 and M1 to M7 have been deposited in the JGI genome portal under project identifier (ID) 1025756 with library ID M2094 and M2113, respectively, and SB1 was deposited in the NCBI SRA under accession number SRX902929. The shotgun data for the same three data sets were deposited in the JGI portal (C1 to C7 are under project ID 1023764, 1023767, 1023770, 1023773, 1023776, 1023779, and 1023782, and M1 to M7 are under project ID 1023785, 1023788, 1023791, 1023794, 1023797, 1023800, and 1018623; SB1 is under project ID 402775).

Previous SectionNext Section


We developed an optimized pipeline that readily analyzes large data sets (see Fig. S1 in the supplemental material). The pipeline has two major steps: SSU rRNA gene fragment search and unsupervised OTU analysis. HMMER-based methods search with HMMs and scale with increasing sizes of SSU rRNA gene databases (6), so we chose them for the first search step. We used meta-rna (16) but could not run it on large data sets due to its poor memory management. Therefore, we simplified and optimized the approach used by meta-rna. Since the search step is still the computational bottleneck (see Fig. S1), our interest here was to make an improvement on search speed and memory efficiency while retaining accuracy. Our implementation is about 4 times faster and 100 times more memory efficient than meta-rna and is 10 times faster and 15 times more memory efficient than metaxa (18) (Table 1). The speed improvement is realized from two modifications. (i) We reduce the number of HMMs to search with by merging Bacteria and Archaea models. SSU rRNA genes are highly conserved; thus, the merged model still has high sensitivity. Even more so, we can increase the sensitivity by using a more relaxed E value cutoff, since false positives are tolerable in this initial search step. (ii) We use reverse-complement HMMs rather than reverse complementing the reads, because the latter scales poorly with large data sets. The newly built HMMs for Bacteria and Archaea together and HMMs for Eukaryota cover most hits found by separate HMMs of the three domains (Bacteria, Archaea, and Eukaryota) in meta-rna in all three test data sets, the two soil data sets in this study, and one synthetic community (Table 1). Our method identified 15,566 (0.03%) of 44,787,632 quality-trimmed and paired-end merged sequences as SSU rRNA gene fragments in the bulk soil sample (SB1) and 112,402 (0.04%) of 274,060,925 reads in the first Miscanthus replicate (M1).

Unsupervised OTU analysis with shotgun data is not available in any current pipeline, so we developed a method for OTU clustering around a small region where all reads overlap. To show the validity of our unsupervised method, we did tests on effects of target region sizes and different variable regions with shotgun data from a synthetic community. We found all region sizes from 50 bp to 160 bp in V4 had an OTU number that approached the species number at a distance cutoff of 4% or 5% (see Fig. S2 in the supplemental material) when testing target region size effect on OTU number. We also did a similar test with only full-length SSU rRNA genes from 64 species to make sure the OTU number is close to the species number and confirmed that the OTU number was close to the species number (at a range of 50 to 60) when a cutoff of 0 to 0.06 was chosen (see Fig. S3). When the target region size was larger than 170 bp, the clustering tool (McClust) (33) did not cluster because the percentage of nonoverlapping reads exceeded its threshold. On the basis of the results described above, we chose 80 bp or 120 bp as the target region size and a 5% distance cutoff for testing different hypervariable regions. The number of OTUs created in all variable regions is close to the real number of species in the synthetic community except when using V3.

Reproducibility between technical replicates is important and is a basic feature of a sequencing method (39, 40). We evaluated it by comparing the correlation of OTU abundance between technical replicates from the bulk soil sample and found high correlation between them (Pearson's correlation coefficient of 0.997). The consistency was better for the more abundant OTUs (see Fig. S4A in the supplemental material). Log transformation could reduce the effect of high-abundance OTUs on the overall correlation. Even with log-transformed abundance, the two replicates have a Pearson's correlation coefficient of 0.91 and linear regression (R2) of 0.88 when OTUs with less than 25 total counts were discarded (see Fig. S4B). The choice of a cutoff of 25 was chosen to compare data to those from another study on amplicon reproducibility (41) and as mentioned in Discussion.

We also compared OTU-based microbial community structures inferred from shotgun and amplicon SSU rRNA gene sequences. OTU abundances in shotgun and amplicon data, however, do not correlate as well (Pearson's correlation coefficient of 0.87 for bulk soil sample SB1 [see Fig. S4C in the supplemental material] and 0.58 for Miscanthus rhizosphere sample M1 [see Fig. S4D]), showing that the amplicon and shotgun methods do not provide the same information. The classification of OTUs with total abundance higher than 10 and with a ratio between two data types higher than 5-fold shows bias against Verrucomicrobia in the bulk soil sample (SB1_PT) amplified by the V6-V8 primer, while it shows a favorable bias for Verrucomicrobia in the rhizosphere sample (M1_PT) amplified by the V4 primer. The classification showed a bias against Actinobacteria in M1_PT amplified by the V4 primer (see Fig. S5). These results were consistent with the taxonomy-based comparison of the two data types (see below) that suggested primer bias in amplicon data.

We applied ordination analysis to OTU tables from an unsupervised analysis of corn and Miscanthus rhizosphere samples. OTUs from shotgun and amplicon data both showed separation of rhizosphere communities of corn and Miscanthus (Fig. 1, horizontal dimension) as well as a significant difference between the two data types (Fig. 1, vertical dimension), confirming the difference between shotgun and amplicon data. Significant separation (P < 0.001 by analysis of molecular variance [AMOVA] test in mothur) of corn and Miscanthus samples also was observed when V2, V4, V6, and V8 shotgun data were used for clustering (Fig. 2), but the sample groupings were the same for all variable regions. Figures 1 and 2 showed that the dispersion among the seven corn replicates was much higher than that for the Miscanthus replicates. Miscanthus samples had higher alpha diversity than corn samples, as shown for each of the V2, V4, V6, and V8 regions, although there were variations among these regions (Fig. 3).


Principal-coordinate analysis (PCoA) of amplicon- and shotgun-derived data from seven field replicates. There are significant differences between amplicon (filled)- and shotgun (unfilled)-derived data (along the y axis) and between corn (circle) and Miscanthus (square) rhizosphere samples (along the x axis) (P < 0.001 by AMOVA [analysis of molecular variance]). PCoA was applied to the OTU table that resulted from clustering with shotgun data and amplicon data using 150 bp of the V4 region. Labels with suffix “_PT” indicate amplicon data, and others indicate shotgun data.


Comparison of ordination analysis with different variable regions. PCoA of OTUs from different SSU rRNA variable regions (V2, V4, V6, and V8) was applied on corn (“_C”) and Miscanthus (“_M”) rhizosphere samples. Different colors indicate the seven replicates. OTU tables from the clustering of shotgun data using 150 bp of V2, V4, V6, and V8 regions were used for PCoA, and Procrustes analysis in QIIME was used to transform the PCoA results from different regions and plot them in the same figure.


Alpha diversity comparisons between corn and Miscanthus samples using V2, V4, V6, and V8 regions. All variable regions showed Miscanthus samples are more diverse than corn samples, even though there was variation of diversity among rRNA gene regions. Alpha diversity was calculated with inverse Simpson (Invsimpson) using OTUs resulting from clustering with different variable regions.

We compared the taxonomy-based microbial community structures inferred from shotgun data with those from amplicon SSU rRNA gene sequences (12,163 amplicons for SB1 and 60,148 amplicons for M1), and we confirmed known primer biases and revealed a new bias. Before comparing two data types, we looked at the taxonomy profile of shotgun data using different variable regions. Shotgun data mapped to different variable regions show similar taxonomies at the bacterial phylum level (Pearson's correlation coefficient of >0.96), except that V6 has more unclassified sequences (see Fig. S6 in the supplemental material). Since different variable regions may provide different levels of taxonomic precision for certain groups (42), taxonomy information from all regions may better represent the taxonomy profile. Thus, we used all SSU rRNA gene fragments for taxonomy comparison with amplicon data. Both shotgun and amplicon data show Actinobacteria, Proteobacteria, and Acidobacteria as the three most abundant phyla, as is expected for soil (43). Since shotgun data are more accurate at estimating community structure than other methods, we accepted the shotgun data as the reference (9, 12). The 926F/1392R (V6-V8) primer set is biased against Verrucomicrobia (0.3% in amplicon data versus 5.8% in shotgun data by RDP database) in bulk soil sample SB1 (Fig. 4), and the 515F/806R (V4) primer set is biased against Actinobacteria (11.6% in amplicon data versus 26.6% in shotgun data) and in favor of Verrucomicrobia (5.9% in amplicon data versus 3.2% in shotgun data) in rhizosphere sample M1 (Fig. 4).


Taxonomy profiles of bacterial phyla of shotgun fragments from SSU and LSU rRNA genes and of amplicon reads. Classifications were done using both RDP and SILVA reference databases. The suffix “_PT” indicates amplicon data. The lower graph depicts bulk soil sample (SB1); the amplicon data used the V6-V8 primer set (515F/806R) and showed fewer Verrucomicrobia detected using both databases. The upper graph depicts the rhizosphere sample (M1), and its amplicon data used the V4 primer set (515F/806R) and showed fewer Actinobacteria and more Verrucomicrobia using both databases.

To take advantage of the fact that shotgun data are untargeted, we retrieved and classified the LSU rRNA genes, which are cotranscribed with SSU rRNA genes. Their taxonomy profile was similar to that of SSU rRNA genes (Pearson's correlation coefficient of 0.87 for SB1 and 0.91 for M1), except that more reads (19.6%) remain unclassified (Fig. 4). This is expected because of the much lower number of reference LSU rRNA genes in the SILVA database. The two genes show consistent community profiles at the bacterial phylum level, and they also confirm the known primer bias against Verrucomicrobia in the 926F/1392R (V6-V8) primer set and that against Actinobacteria in the 515F/806R (V4) primer set (Fig. 4). Further, both the LSU and SSU HMMs show the ability to identify members of the Eukaryota and give about the same taxonomy profile at the domain level (see Fig. S7 in the supplemental material). It is worth noting that both LSU and SSU shotgun data show Archaea to be twice as numerous (2% versus 1%) in the bulk soil as in the rhizosphere and that Eukaryota were much more numerous (6% versus 1%) in the Miscanthus rhizosphere soil than in the bulk soil (see Fig. S7). A higher fungal percentage (2.54% versus 0.39%), fungus/bacterium ratio (0.027 versus 0.004), and arbuscular mycorrhizal fungus (AMF) percentage in fungi (0.18% versus 0) were found in the rhizosphere sample (M1) than in the bulk soil sample (SB1) (see Table S1 in the supplemental material). Copy correction was applied to two soil samples (SB1 and M1). Both samples showed that Firmicutes and Bacteroidetes had the highest fold change after copy number correction (see Fig. S8). Despite copy number correction, the clustering of our soil samples did not change (see Fig. S9) compared to that without copy number correction (Fig. 1), probably because of the low proportion of taxa with large rrn number corrections.


Comparison of search results from SSUsearch, metaxa, and meta-rnaa

Previous SectionNext Section


We present, characterize, and validate an efficient method for retrieving and analyzing SSU rRNA gene fragments from shotgun metagenomic sequences. The pipeline enables unsupervised diversity analysis with copy number correction on multiple variable regions, has the scalability to handle large soil metagenomes, is expandable to other phylogenetic marker genes, and is publicly available on GitHub.

We apply a two-step approach for retrieving SSU rRNA gene fragments, a loose HMM filtering step followed by a more stringent step that screens by identity to the best-match reference. The first step leverages HMMER (31); thus, it should have better scalability than existing shotgun analysis pipelines that use BLAST-like tools such as MG-RAST (44). MG-RAST annotates shotgun reads by BLAT search against rRNA databases, and the taxonomy of reads is inferred from the best hit or least common ancestor of several top hits (24, 45, 46). BLAT or BLAST-like tools are not scalable for large data sets and must be run in parallel on large computer clusters, because BLAST-like tools typically do pairwise comparisons of reads against large and growing rRNA databases, such as RDP, SILVA, and Greengenes (6, 7, 37), while HMM-based methods compare reads to only a fixed number of models (commonly one for each domain) and are more scalable (44). Moreover, these pipelines lack unsupervised community analysis. An HMM-based search has been used before, as it is fast and sensitive for rrn retrieval (16–18, 22), and current existing implementations, such as meta-rna, RNASelector, and metaxa, all are wrappers around HMMER (31). We chose meta-rna and metaxa for comparison, because RNASelector can run only in a graphic interface that is not suitable for large data sets.

The second step evaluates hmmsearch results based on identity to their best-match reference and also prepares the alignment of SSU RNA gene fragments for clustering. Since there is no clear sequence identity threshold for SSU rRNA genes, the choice of identity cutoff is arbitrary (the default is 50%). This is also a common practice for amplicon analysis platforms, where reads with low identity to reference sequences are discarded prior to clustering (19, 20). For consistency in comparison, sequences in our amplicon data sets with less than 50% identity to reference sequences also are discarded. An alignment of the SSU rRNA gene fragments is essential for the later unsupervised analyses. Compared to methods that use only the 16S rRNA gene in E. coli as an alignment template (47), our method takes advantage of the rich phylogenetic diversity of SSU rRNA genes provided by the SILVA database. Increasing the number of reference sequences can improve the quality of alignment but also linearly increases the memory required (32, 48).

Our unsupervised analysis with shotgun data is a novel and important part of the pipeline. Our tests show that regions as small as 50 bp can be applied to clustering. Thus, short reads around 50 bp can be applied to this method as long as there are sufficient numbers of reads aligned to the target region (sequencing depth is a limiting factor; see below). This is consistent with pilot studies from 454 amplicon sequencing (5, 49). When sequencing depth is limited, there is flexibility to control the number of reads to include for clustering by adjusting the target region size and read length cutoff within certain limits (see Fig. S2B in the supplemental material). The caveat of using very short or very large target regions is that the overlapping portion of reads will decrease; thus, there will be a decrease in the accuracy of clustering and the impact of sequencing error will increase. For example, an error in a 50-bp read can cause 2% distance; accordingly, we need to set a larger distance cutoff for clustering. In addition, we can obtain longer sequences from overlapping paired ends, as shown in our soil data (see Fig. S10). Thus, reads from Illumina shotgun data (ranging from 75 to 250 bp) can be used for unsupervised analysis. Note that the flexibility on the choice of variable region for analysis is another advantage of shotgun data (see Fig. S2C).

We also found good reproducibility of OTU abundance between technical replicates, which is critical for the validity of our method (see Fig. S4A in the supplemental material). Generally, OTU-based analysis provides higher resolution than taxonomy-based diversity analysis for community comparison, largely due to the databases lacking reference sequences from uncultured microbes (50). The high correlation of OTU abundance in two technical replicates shows the reproducibility of the analysis of shotgun data, which is comparable to the reproducibility of amplicon data shown in another study in terms of Pearson's correlation coefficient and linear regression (R2) (41). Further, comparison of OTU abundances in shotgun data and amplicon data sequenced from the same DNA extraction also show that many OTUs have inconsistent abundances between the two types of data (see Fig. S4C and D), which agrees with the differences seen in the taxonomy-based comparison (Fig. 4).

Community comparison by ordination methods such as PCoA and NMDS (nonmetric multidimentional scaling) is one of the most common analyses in microbial ecology. To the best of our knowledge, the methods used in two previous studies (47, 51) are the only existing tools that are designed to deal with the clustering of SSU rRNA gene fragments from Illumina shotgun data. The method used in the first study could result in poor alignment by using only the 16S rRNA gene in E. coli as the alignment template. The method used in the second study (PhylOTU) was applied to larger shotgun sequences from Sanger sequencing. It determines the OTU clustering of SSU rRNA gene fragments aligned over the whole gene length, which can be problematic, because fragments aligned to different regions do not overlap and the clustering results are not reliable, even though the reference sequences included in the clustering process can improve the results. Since our tests show that a hypervariable region as small as 50 bp can be used for unsupervised analyses (see Fig. S2A and C in the supplemental material), we made sure all of the sequences overlapped by picking one small region (150 bp), and all sequences included in clustering have lengths longer than 100 bp in our clustering method. In addition, longer reads obtained by assembling overlapping paired-end reads (see Fig. S10) can make use of more overlap among reads and are more suitable for clustering. As read lengths increase with the improvement of sequencing technology, larger regions can be chosen and the clustering results will be even more reliable. Also, shotgun data provide the flexibility to choose any variable region (Fig. 2 and 3; also see Fig. S2C and S6), and the consistency of results from different variable regions provides more confidence in the biological conclusions as well as the method itself.

Primer bias is a major limitation of amplicon methods, and it was apparent in our comparisons of community profiles from amplicon versus shotgun data (Fig. 4). Commonly, it is difficult to tell if a bias is caused by primer or DNA extraction. The paired amplicon and shotgun data from the same DNA extract provide us a new opportunity to evaluate this issue, since the difference is only from the sequencing step. We used two main SSU rRNA gene databases, RDP and SILVA, to make sure the taxonomy distribution was not biased by the choice of reference databases, and we further confirmed by taxonomy the distribution of the LSU rRNA gene. The bias against Verrucomicrobia with V6-V8 primers is consistent with other studies showing that the abundance of Verrucomicrobia in soil samples is underestimated due to primer bias (8). Meanwhile, the bias toward Verrucomicrobia with V4 primers agrees with studies showing that the V4 primer set has better coverage of Verrucomicrobia (8). Furthermore, the V4 primer set shows bias against Actinobacteria. The V4 primer set has been reported to cover 92.4% of Actinobacteria in reference databases, the lowest level of coverage among nine common bacterial phyla (8). Bias against Actinobacteria also has been reported in a study on a synthetic community where no primer mismatch was found with members from Actinobacteria (12), and another study on environmental samples using Sanger sequencing (24F/1492R) (52) suggested that the in silico evaluation of primers is not sufficient and that factors other than primer mismatch are causing the bias, for example, competition between primers and melting temperature (53). Thus, primer bias detection by comparing paired amplicon and shotgun data is superior to methods that only search primers in the reference databases.

Another advantage of our method over amplicon approaches is that we can identify the SSU rRNA gene from Bacteria, Archaea, and Eukaryota. Fungi are of special interest in microbial ecology due their critical roles in ecosystems (54, 55). The fungus/bacterium ratio is an important indicator of C/N ratio and soil health (56, 57). Our shotgun metagenome (DNA based) shows a fungus/bacterium ratio of 0.004 in bulk soil (SB1) and 0.027 in rhizosphere soil (M1) (see Table S1 in the supplemental material), while studies using phospholipid fatty acid analysis (PLFA) at the same sampling site (KBS) typically show ratios of 1 to 1.3 (58). The difference can be explained by the higher biomass of fungi relative to that of DNA, since some fungal hyphae may not be filled with nuclei. In addition, we also found higher percentages of AMF of rhizosphere soil (M1) than of bulk soil (SB1), which is consistent with their symbiotic relationship with grass roots.

We also show that copy number correction can be achieved in our pipeline. Gene copy number is another source of bias that limits one's ability to accurately profile microbial communities. There are up to 15 SSU rRNA gene copies in some bacteria and up to 5 in archaea (59). This pipeline utilizes the SSU rRNA copy database in CopyRighter (36). As expected, Firmicutes in soil samples have the highest fold change (see Fig. S8 in the supplemental material). However, due to their low proportion in these soils, the impact of copy number correction on the overall community profile was minor (see Fig. S9). Copy number correction, however, is still an unresolved issue, because SSU rRNA gene copy number data for most species and/or OTUs are lacking and copy number can be incorrect even for species with complete genome sequences because of the misassembly of these repeated regions.

Sequencing depth is another important factor in considering this method for diversity analysis. The percentage of SSU rRNA gene fragments in shotgun data varies depending on the SSU rRNA gene copy number and genome size of each member. In our bulk soil sample (SB1) and the Miscanthus rhizosphere soil sample (M1), we classified about 0.03% and 0.04%, respectively, of the total shotgun data as SSU rRNA. In an ideal situation we want to obtain enough SSU rRNA gene fragments to see saturation of the rarefaction curve in OTU-based analysis, which is difficult for soil samples because of their high diversity and the presence of sequencing error. However, studies have shown that near-saturation sequencing of SSU rRNA amplicons is not necessary for beta-diversity analysis (4, 60). Thus, the empirical fold coverage of 3,000, based on the whole length (about 1,500 bp) of the SSU rRNA gene, is suggested for surface soil samples, which require 11.2 Gb [i.e., 1,500 bp × 3,000 bp ÷ 0.04%] of shotgun data, assuming the SSU rRNA gene comprises about 0.04% of total data.

In this study, the LSU rRNA gene was used mainly as confirmation for SSU rRNA gene-based diversity analysis (Fig. 4; also see Fig. S7 in the supplemental material). However, the LSU rRNA gene offers additional stretches of variable and characteristic sequence regions due to its longer sequence length and yields better phylogenetic resolution (61). For this reason and because there are more available references for fungi, the LSU rRNA gene is used more commonly for fungal community studies (62–65). Currently, the use of the LSU rRNA gene is limited by reference sequences and available universal primer sets (61), but its increased resolution should not be overlooked, since a too-limited resolution of the SSU rRNA gene is a barrier in many ecological studies (54). In the future, other single-copy genes with phylogenetic references and finer resolution, such as rplB, gyrB, and recA, also could be recovered from metagenomic sequences and used for community structure analysis by a similar pipeline (66).

Conclusions.We developed a fast and efficient pipeline that enables unsupervised diversity analysis with Illumina shotgun data. The pipeline has the scalability to analyze large data sets (5 central processing unit hours for 38 Gb data, with 4.8 Gb peak memory) and can be run on most desktops with more than 5 GB of memory. Since SSU rRNA-based community analysis is an important method in microbial ecology, this method can save projects with existing shotgun sequence data from the additional cost of SSU rRNA amplicon sequencing. Moreover, shotgun sequencing is not as affected by primer bias and chimeras as amplicon sequencing; thus, it can improve the measurement of microbial community structure. As read length and sequencing depth increase, longer and more SSU rRNA gene fragments can be recovered. Thus, clustering and diversity analysis by this pipeline will become even more reliable.


This work was funded in part by the U.S. Department of Energy (DOE) Great Lakes Bioenergy Research Center (DOE Office of Science BER DE-FC02-07ER64494) and by DOE Office of Science grants BER DE-FG02-99ER62848 and DE-SC0004601.

Previous SectionNext Section


    • Received 26 August 2015.
    • Accepted 13 October 2015.
    • Accepted manuscript posted online 16 October 2015.
  • Address correspondence to James M. Tiedje, tiedjej{at}
  • Citation Guo J, Cole JR, Zhang Q, Brown CT, Tiedje JM. 2016. Microbial community analysis with ribosomal gene fragments from shotgun metagenomes. Appl Environ Microbiol 82:157–166. doi:10.1128/AEM.02772-15.

  • Supplemental material for this article may be found at


  1. 1.↵
  2. 2.↵
  3. 3.↵
  4. 4.↵
  5. 5.↵
  6. 6.↵

Citation: Taib N, Mangot J-F, Domaizon I, Bronner G, Debroas D (2013) Phylogenetic Affiliation of SSU rRNA Genes Generated by Massively Parallel Sequencing: New Insights into the Freshwater Protist Diversity. PLoS ONE 8(3): e58950.

Editor: Stefan J. Green, University of Illinois at Chicago, United States of America

Received: July 13, 2012; Accepted: February 11, 2013; Published: March 14, 2013

Copyright: © 2013 Taib et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Funding: Funding came from Le sonseil régional d'auvergne "". The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.


The development of molecular ecology was prompted by indisputable evidence that, for most environments on Earth, the majority of existing organisms had not yet been cultured. This evidence came from the analysis of sequences recovered directly from environmental samples. Vast new lineages of microbial life were uncovered by this approach, changing our picture of the microbial world and yielding a phylogenetic description of community membership [1], [2]. More precisely, the sequencing of the small sub-unit (SSU) rRNA genes highlighted new monophyletic groups or clades in the environment, such as SAR11 [3] or MGI [4] among the Bacteria and Archaea respectively. Similarly, several new lineages of protists have been discovered in oceanic systems during the last decade [5]. Recent studies conducted in lakes have also highlighted numerous phylogenetic groups, especially putative parasites (Fungi and Perkinsozoa), and this finding is modifying our view of the microbial loop and therefore, the functioning of aquatic ecosystems [6], [7].

Recent advances in next-generation sequencing (NGS) technologies are spurring progress in determining the microbial diversity of various ecosystems by highlighting, for example, the rare biosphere and the activity of these low abundance organisms [8], [9]. Currently, the pyrosequencing of amplified SSU rRNA gene variable regions is mainly used to determine bacterial and archaeal diversity and structure in various ecosystems, such as soil [10], ocean [11] or gut microbiota [12]. The recent results obtained regarding the composition and structure of the microeukaryote communities using high-throughput amplicon sequencing performed with the Roche 454 pyrosequencing platform in freshwater systems [13], [14] have fuelled the current debates on the biogeography of these microorganisms and on the role of the rare biosphere. The taxonomic assignment of such data is often inferred from supervised classification with the Ribosomal Database Project Classifier (RDP) [15], sequence similarity with BLAST [16]–[18] or both [19], [20]. Pairwise identity scores via BLAST remain the most commonly used tool for large eukaryotic datasets [14], [21]–[26]. However, as claimed by Bik et al. [26], assigning accurate taxonomy to eukaryotic operational taxonomic units (OTUs) is more difficult than the approaches used for Bacteria; the relative paucity of sequences in public eukaryotic databases results in many sequences without significant top BLAST matches [26]. Furthermore, the best BLAST match assigns a single organism as the most likely phylogenetic neighbor, without specifying the level of relatedness (class, order or phylum) of the compared sequences [27].

Phylogenetic methods assess relatedness among various groups of sequences by inserting unknown OTU sequences within a known phylogeny. On the one hand, these methods allow query sequences to be affiliated with their relatives. Tree-based assignment is, therefore in theory, a more robust approach [28] and current FLX Titanium longer reads now make it possible to extract phylogenetic information with a high degree of reliability [29]. On the other hand, phylogenetic analyses allow for the description of clades, which may lead to new insights into the structure and functioning of ecosystems, as previously mentioned. Moreover, these phylogenetic analyses are not limited to the taxonomic assignment of an individual sequence as implemented in bioinformatic pipelines dedicated to NGS and used in microbial ecology studies (mainly on 16S rRNA gene amplicons): phylogenies can also be used to compare environments (beta-diversity) using methods based on tree topology and/or branch length such as the popular tool UNIFRAC (unique fraction metric) [30]. Although more robust, these methods are less frequently used than BLAST or probabilistic classifiers, as they require more computing resources (Table S1). Though large computational capacity is now more accessible (e.g., QIIME [20] can be implemented on a cloud), massively parallel sequencing projects that seek to elucidate the phylogenetic structure of microbial populations are still faced with the attendant computational challenges of classifying the sequences obtained.

In this work, we introduce a tree-based treatment designed for analyzing massively parallel sequencing outputs that automatically affiliates sequences from SSU rRNA gene amplicons and builds phylogenetic trees composed of very large numbers of sequences. As short-read sequence data (e.g., 100 base sequences generated by the Illumina sequencing platform) provide limited phylogenetic resolution [29], our work is focused on the treatment of moderately long (∼ 450 bp obtained for example with Titanium platforms) to near full-length sequences. Designed for the analysis of any microorganism (protists, Bacteria and Archaea), the value of this treatment is highlighted here on the protist diversity as the pipelines dedicated to the study of eukaryotic pyrotags are still scarce. Indeed, 16S rRNA gene reads were widely investigated in previous studies [31]–[33] to assess bacterial diversity, which enhanced the development of specific 16S rRNA gene analytical tools. However, 18S rRNA gene surveys and tools allowing for the accurate and rapid taxonomic affiliation of protists from NGS data are needed because the number of studies dealing with protists diversity is currently increasing (e.g., [14], [34]). We first tested the accuracy and speed of phylogenetic affiliation on large fragments of well-annotated 18S rRNA gene sequences (>1,200 bp) and on short sequences that simulate pyrosequencing outputs. Secondly, the different methods of taxonomic assignment (i.e., tree-based, similarity and probabilistic approaches) were compared with each other, in a first attempt to determine the best method for affiliating protists in the context of massively parallel sequencing of amplicons. Thirdly, the accuracy of phylogenetic affiliation was compared on amplicons covering different variable regions (V1 to V9), and finally, a dataset of original pyrosequencing data obtained from lacustrine small protists was analyzed by the tree-based approach that was developed.


Evaluation of performance on reference sequences

In the analysis of near full-length reference sequences of 18S rRNA gene, taxonomic groups were found in similar proportions to those initially present in the samples. Our phylogenetic affiliation method, referred to as PANAM (Phylogenetic Analysis of Next-generation AMplicons), was more accurate using LCA (lowest common ancestor) assignment for the different taxonomic ranks, ranging from 99.1% to 90.8% versus 98.6% to 86.7% for PANAM using the NN (nearest neighbor) method (Figure 1.A). For comparison, when refining affiliations from kingdom to genus, the accuracy of the standard phylogenetic affiliation using ClustalW [35] and PHYML [36] as implemented in STAP, ranged between 96.1% and 74.6%. At the finest phylogenetic level studied (i.e., genus), BLAST and RDP allowed for the affiliation of 62.3% and 68.4% of reference sequences. Thus, our phylogenetic affiliation method outperformed the other methods on near full-length sequences. However, as environmental sequences are generally quite divergent from referenced ones and their affiliation needs to be checked manually, sequences belonging to freshwater clades [6], [7] were also processed by our phylogenetic affiliation method to evaluate how it behaved on these datasets. The phylogenetic analysis of these environmental sequences (Sanger, >1,200 bp) enabled us to retrieve the affiliations obtained by other authors together with the delineation of freshwater clades corresponding to Cercozoa clade [6] and Perkinsea clades 1 and 2 [7] (Figure S1).

Figure 1. Accuracy of the phylogenetic affiliation of PANAM compared to different approaches and on different regions

. 1.A. Accuracy of the phylogenetic affiliation of PANAM-LCA, PANAM-NN, STAP, BLAST and RDP Classifier. 1,000 near-full-length sequences were randomly picked from the reference database and removed from it for the simulations. For PANAM, simulations were repeated 5 times and the standard variation is less than 0.03. 1.B. Accuracy of the phylogenetic affiliation in relation with the variable region targeted. The specificity was tested with PANAM-LCA and a sequence length equal to 400 bp.

Different 18S rRNA gene regions were targeted by simulating amplicons with lengths of 200 and 400 bp starting from a conserved region given by the following forward primers: NSF4, NSF370, NSF573 NSF963, NSF1179 and NSF1419 (Table S2). Because the V8–V9 region is often missing in public databases, the results obtained from this region were based only on 300 sequences included in the reference database. The affiliation results at the genus rank differed according to length, variability within the studied region and method used for taxonomic affiliation (Table 1). For the six regions tested, the accuracy increased with amplicon length for both affiliation methods implemented in PANAM, LCA and NN. Considering the affiliation methods, LCA specificity was higher than that of NN for fragments of 200 bp only for the V1 and V8 regions, and LCA specificity was always better for fragments of 400 bp. The comparisons with the other affiliation tools implemented in pipelines dedicated to pyrosequencing results showed that at 200 bp, BLAST outperformed RDP, STAP and PANAM, with the exception of the V8 region, for which PANAM (LCA) gave the highest result (68.7%). In contrast, for 400 bp amplicons, the most accurate affiliations were obtained with PANAM, with the exception of the V5–V6 amplicon. In this last region, we observed a decrease in the accuracy of the affiliation, coupled with a sharper decline for the phylogeny-based affiliations. The specificity therefore varied between 64.2% (V5–V6) and 79.2% (V8–V9) at the genus level.

In addition to the accuracy of assignment, this phylogenetic affiliation method was developed to optimize processing time for large datasets. Thus with a 2 GHz Intel(R) Xeon(R) and 24 GB RAM and with a single 32-bit CPU, PANAM can process the phylogenetic analysis of 1000 eukaryotic OTUs of 400 bp in approximately 20 minutes, regardless of the affiliation method. The run time increased with the number of OTUs, regardless of the length. For example, for 400 bp, the run time ranged from 24 minutes for 5000 OTUs to 6 days and 14 hours for 1 M eukaryotic OTUs. For near full-length sequences, PANAM was able to process 1 M sequences in 16 days (Figure S2).

Reliability of the phylogenetic affiliation in relation to the region targeted and the taxa of interest

The reliability of affiliations was compared for 400 bp reads spanning the 18 S rRNA gene for four taxonomic groups: Alveolata, Stramenopiles, Fungi and Viridiplantae at the genus level (Figure 1.B). Generally, the fragment affiliation depended on the taxonomic group and the region considered. According to previous results, the regions from V5 to V6 gave, on average, the weakest accuracy. Another general trend observed in this analysis was a poor taxonomic restitution for sequences belonging to Viridiplantae compared to other groups, between 52.4% and 59.1% regardless of the region targeted. The best specificity values for Stramenopiles, Alveolata and Fungi were obtained in different regions: V1–V2 (89.5%), V4 (81%), and V8–V9 (82.2%) respectively. The taxonomic affiliation for these three groups from the V8–V9 region was relatively similar, from 78.5% to 82.2%.

Tree-based analysis of pyrosequencing data from small lacustrine protists

In silico simulations have shown that primers NSF573 and NSR1147, used to target the V4 region of the 18S rRNA gene captured the greatest diversity (data not shown) and that the region amplified by these primers is suitable for taxonomic affiliation (Table 1). The reads were clustered at 95% similarity, and 6% of the OTUs (4% of reads) defined from this pyrosequencing run matched with Metazoa sequences and were not processed further. The diversity and richness indexes obtained for each environment are shown in Table 2. The lowest and highest richness indexes (Chao1) were found on Anterne Lake and Villerest Lake respectively, whereas the normalized indexes (based on 3759 sequences) showed that Bourget Lake harboured the largest number of species (Table 2). This normalization also had an effect on the richness estimates in Godivelle Lake and Geneva Lake.

In the lakes studied, regarding level 2 and 3 from EMBL taxonomy (displayed in Table S3, a PANAM table output, including number of sequences, OTUs and diversity indexes), the major phylogenetic groups were Fungi, Alveolata and Stramenopiles representing 73.2% of OTUs and 78.6% of sequences (Figure 2). These mean values mask some disparities between lakes. Thus, Anterne Lake harboured mainly reads affiliated to Fungi (99.4% of total), whereas the main phylum in Geneva Lake was Alveolata (Figure 2; Table S3). Sequences belonging to the phylum Cryptophyta were the most abundant in Pavin Lake and Sep Lake.The results highlighted the presence of freshwater clades delineated in previous studies [6], [37] such as Cryptophyta_2 to Cryptophyta_4, Rhizophydium or Cryptomycota (previously known as LKM11) among Fungi (Table S3). Sequences derived from Fungi, which were very abundant in sequence libraries from Anterne Lake and Aydat Lake, belonged to this last Cryptomycota clade (Table S3, Figure S3). These data demonstrate the presence of Chlorophyta and Haptophyta in all of the lakes studied, with the exception of Anterne Lake, which is characterised by an over-representation of Fungi and an absence of Haptophyta. This tree-based approach allows for the study of beta-diversity from phylogenies. The UNIFRAC metric showed that Bourget, Aydat and Anterne Lakes differed from other ecosystems regardless of the phylogenetic level (total Eukaryotes, Stramenopiles and Fungi) at which the analysis was performed (Figure 3).

Figure 3. Principal coordinate analysis computed using a Unifrac distance metric from the phylogenies of the Stramenopiles,

Fungi and the total eukaryotes. This analysis permit to differentiate environments according to their taxonomic composition. For example, Lake Godivelle seems to be different from the other lakes for the Stramenopiles, while it is similar for all eukaryotes.

In a comparison of the OTUs found in this study to those present in previous studies on the small protists, only 4.8% were previously detected in lakes. If only the dominant OTUs (>1% of reads) are taken into account, then the proportion of OTUs similar to specific lacustrine sequences increased to 19.7%. Moreover, new light is shed on putative clades of small protists. Specifically, these clades include the chlorophycean group of Mamiellophyceae, represented in Figure 4; Foraminifera (Rhizaria); Dictyochophyceae (Stramenopiles); and Euglenida (Euglenozoa). These clades were supported by high bootstrap values (> 0.8), included 23, 14, 17 and 23 OTUs respectively, and were found in at least three of the eight lakes. The novel clade within the Euglenozoa was composed only of OTUs present at less than 1% of reads.

Figure 4. Main putative clades detected among Mamiellophyceae (Chlorophyceae) based on 18S SSU reads (425 bp ± 114).

The OTUs affiliated to Chlorophyceae were generated at 95% similarity. A profile alignment was processed using HMMalign and the phylogeny was built by FASTTREE2 with 100 bootstraps. The distribution of the OTUs among different lakes shows a main presence of clade 1 in Lake Pavin while clade 2 is mainly present in Lake Godivelle.


As the interplay between evolution and ecology receives more attention in ecosystem studies [38], there is greater interest in phylogenetic approaches for deciphering the mechanisms that govern the diversity and functioning of communities and ecosystems. However, the phylogenetic methods that are typically applied to Sanger-sequenced SSU rRNA are computationally expensive and cannot be readily used to handle NGS datasets; therefore, pyrosequencing reads are mainly analyzed by other approaches. The method described in this study is a response to the challenge of analyzing hundreds of thousands of SSU rRNA genes in a phylogenetic framework, inferring taxonomies from sister sequences and describing clades. This method has been implemented and tested for microorganisms with an emphasis on protists, which are not well served by bioinformatics tools dedicated to NGS data, although the early focus on bacterial and archaeal diversity has recently broadened to include eukaryotic microorganisms [39], [40]; thus, the database provided in PANAM includes reference sequences from protists, Bacteria and Archaea and can be used for taxonomic assignment of all microorganisms.

Accuracy of affiliation methods for protist sequences

Our taxonomic affiliations were compared with BLAST, a tool commonly used for the identification of microorganisms especially microeukaryotes (e.g., [22]); RDP, which is currently used to classify bacterial and archaeal SSU rRNA sequences and fungal LSU rRNA sequences; and STAP implemented in WATERS [41]. This method, based on ClustalW alignments and PHYML phylogenies, is a standard method for taxonomic affiliations based on phylogenetic analyses. The RDP Classifier [42] is often considered to be restricted to bacterial and archaeal taxa [26] and therefore, is not used for eukaryotic classification of SSU rRNA genes after amplicon pyrosequencing. We used this tool for the first time for taxonomic affiliation of 18S rRNA gene amplicons generated with high-throughput pyrotag sequencing. The affiliation of simulated amplicons were obtained by the RDP Classifier trained on the near full-length sequences of the reference database used in PANAM. Surprisingly, trimming the reference database to the primer region did not result in an improvement of classification for 18S rRNA gene sequences (data not shown), in contrast to the results of Werner et al. [43] on 16S rRNA gene sequences. As noted by these authors, a naïve Bayesian classification depends on the training set size. The weak performance on the truncated sequences could thus be explained by the limited number of 18S rRNA gene sequences in public databases compared with 16S rRNA gene sequences, particularly for the V9 region (see the discussion below).

The comparison of the tree-based method proposed with these tools in the context of taxonomic affiliation of 18S rRNA gene amplicons shows that regardless of the method that is used, taxonomic reliability depends on the sequence length and amplicon location on the SSU rRNA gene sequence. These results, which to our knowledge have not been examined for 18S rRNA gene sequences, are consistent with observations of 16S rRNA gene sequences from Bacteria and Archaea [44].

Our results mostly illustrate the impact of sequence length on phylogenetic methods, which appears to be the main limitation of this approach. According to Liu et al. [31], it is possible to use short fragments from the 16S rRNA gene to draw the same conclusions as with full-length sequences. However, by comparing different affiliation methods, they also noted that the short reads generated by pyrosequencing (i.e., 200 bp) were likely to be problematic for inferring phylogeny due to their small number of bases; similarity and probabilistic methods are therefore the most accurate. However, our analysis, similar to the one proposed by Jeraldo et al. [29] for 16S rRNA gene sequences, demonstrates that with the current average length achieved by the pyrosequencers (Titanium generation; > 400 bp), phylogenetic methods are reliable and offer an advantage over other methods such as RDP. From 400 bp amplicons, the phylogenetic affiliation method implemented in PANAM outperforms the classical tools dedicated to NGS analysis at the genus level with the exception of amplicons sequences covering the V5–V6 region of the SSU rRNA gene. Phylogenetic methods are generally considered superior to other approaches for taxonomic affiliation [45] as they assess relatedness between a set of sequences. They are also considered to be difficult to automate as i) their reliability greatly depends on the quality of the alignments, which need to be validated by experts in the field, and ii) they use intensive, time-consuming methods for tree building.

In this study, we use the curated alignments sequences provided by SILVA, which is, at least for eukaryotic sequences, the only up-to-date curated database. All high-quality and near full-length aligned sequences suitable for in-depth phylogenetic analysis were selected. However, the guide-tree for eukaryotes provided by SILVA, in contrast to the other domains, represents only an approximate phylogeny. Tree-based approaches can implement other tools based on the tree-insertion methods like pplacer [46] as proposed by Bik et al. [28]. Similarly to STAP, this tool analyzes one sequence at a time. Thus, clades may be, at best, approximated from a frozen backbone tree, while the addition of distant taxa, as can be expected from environmental sequences, may require a re-evaluation of the phylogenetic tree [46]. In terms of processing time, we demonstrated that the tree-based method described here can process 1 M sequences in a reasonable (about three hours) time scale. For comparison, while pplacer processes 10,000 sequences in ∼0.5 hour, PANAM can process 30,000 sequences in the same amount of time with the same computational resources. However, while a pyrosequencing run can produce up to 1.2 M reads, the raw sequences first go through a quality control stage that eliminates poor quality reads and replicates. Additionally, in diversity studies, the raw sequences are first cleaned (i.e., quality trimmed) and clustered, and phylogenetic analyses are applied to the representatives of each OTU and not to all of the raw reads from a run. Consequently, in current studies of diversity, the effective number of sequences to be affiliated is on the order of tens of thousands, which can be processed in a few hours on a personal computer.

Accuracy of the protist affiliation in relation to the region targeted

The primers used for the taxonomic assignment of Bacteria traditionally span the regions V3, V6 and V9 of the SSU rRNA gene [12], [47]. However, some studies [32], [48] suggest that the V6 region is not optimal for taxonomic affiliation as it overestimates richness and the number of OTUs at different cut-offs [49]. In the microeukaryotic field, the regions V2–V3 [13], V3 [14], [34], V4 [22], [23], [39] and V9 [21], [22], [24], [25], [39] were investigated with limited in silico analysis. Behnke et al. [39] partially addressed this concern because they compared the V4 and V9 regions for analyzing sequencing errors; V4 amplicons are likely more prone to an increased frequency of Roche 454 pyrosequencing homopolymer errors relative to the V9 region [22]. However, the inclusion of at least some part of the variable regions of the SSU RNA gene is necessary for the methods to retrieve sufficient signal for taxonomic affiliation. Liu et al. [32] stressed that tree-based methods are more sensitive to the 16S rRNA gene region targeted than are similarity-based methods because of different rates of evolution among regions [44], and/or the difference of homopolymer incidence and length between the regions [48]. The same conclusions can be drawn from our results from 18S rRNA gene amplicon sequences, because the accuracy of the phylogenetic affiliation for the region V5–V6 dropped for both phylogenetic methods used in this study (STAP and PANAM). Interestingly, the accuracy of the taxonomic affiliation of the main phyla varied with the region analyzed, but regardless of the variable region analyzed, simulated amplicons from Viridiplantae were always difficult to affiliate reliably at the genus level. Thus, the bias observed between variable regions [22] could be due to primers that may not anneal uniformly to all groups, but also to the bioinformatic process used for the taxonomic identification. In summary, with the exception of Viridiplantae, the V8–V9 region appears to be a good candidate for the study of protist diversity because the reliability of the taxonomic affiliation did not differ according to the phyla considered (i.e., Stramenopiles, Fungi, Alveolata). However, sequence databases such as GenBank contain many fewer sequences that include the V9 region than other variable regions.

New insights into the small protist composition of the lacustrine ecosystem

In this analysis, our goal was not to explain the spatial pattern of the protist community composition (PCC) but to characterize the structure of these communities (richness, diversity and composition) by high-throughput SSU rRNA gene amplicon sequencing and sequence affiliation utilizing a tree-based method. We focused on the optimization of processing environmental data and on the description of the general picture of protists diversity obtained for these lakes.

For an in-depth analysis of this PCC from lacustrine ecosystems, we introduced environmental sequences and taxonomies in the reference database to delineate specific clades as defined in previous publications (e.g., [6], [37], [50]). The introduction of “environmental reference” sequences reflecting the taxonomies of protists originating from specific environments can enhance the affiliation of poorly represented environmental sequences. Phylogenetic methods provide a clear edge in describing under-studied and complex communities. However, as with other methods, the precision of sequence mapping falls off when experimental sequences lie distant from reference SSU rRNA gene sequences [51]. This observation is particularly true for environmental sequences, for which the availability of close relatives and well-annotated sequences in reference databases is limited, as is the case for the V9 region. If the referenced trees do not include known relatives branching close to experimental reads, divergent lineages form long-branch taxa with no close reference sequences at relatively deep internal nodes. This phenomenon results in a less precise taxonomic affiliation of these sequences; however, clades of interest could still be drawn, as very similar sequences (i.e., sequences with low pairwise distance) are very well preserved among tree searches from de novo phylogenies [29].

Most eukaryotic species are defined on morphological differences, however, as the majority of existing microorganisms on Earth have not yet been cultured, their phenotypic traits can hardly be described. Thus, environmental microbial species are delineated according to a sequence similarity cut-off based on comparisons of SSU rRNA gene sequences to demarcate operational taxonomic units [52]. Although they do not technically represent species, OTUs composed of multiple sequences can be used to describe novel species, using the provisional designation of “Candidatus”, when the SSU rRNA gene sequences are sufficiently different from those of recognized species [53]. In this study, after dataset cleaning and sorting, the reads left for the affiliation were clustered at a 95% identity threshold as proposed by Caron et al. [54] to delineate eukaryotic taxa. These authors defined this similarity threshold after studying the distribution of intra- and inter-specific variations of the 18S rRNA gene in protistan communities. However, as they pointed, this cut-off is a conservative estimator of species richness, and may mask considerable physiological diversity in some OTUs. In other studies, taxon clustering is performed at sequence similarity from 90% to 100% [23]. As the error rate of many NGS platforms in any case is ∼1% it is recommanded to cluster at a lower threshold than 99%. Some authors chose a similarity of 97% because this value is commonly used to define OTUs in Bacteria (e.g., [22]). However, this value has been defined for delineating a species from the full-length 16S rRNA gene. Thus, from in silico analysis of 16S rRNA genes, Kim et al. [33] showed that the clustering threshold must be chosen according to the variable region amplified and the domain studied (i.e., Archaea or Bacteria). A less conservative cut-off could overestimate the richness and diversity because in some phyla, such as diatoms, the level of intragenomic polymorphism in the SSU rRNA gene can reach 2% [55]. Finally, in a previous study, Mangot et al. [56] defined a threshold of 95 % by adding an internal standard (a clonal sequence derived from a copy of the 18S rRNA gene in Blastocystis subtype 4 genome) before amplifying and sequencing the DNA samples. Indeed, all the amplicons derived from this sequence clustered in one OTU at this cut-off.

Our tree-based treatment applied to NGS sequences demonstrated that few OTUs have been previously described by the traditional cloning-sequencing (CS) method. As these OTUs represent taxa present in relatively low abundance in many environments, little information is available about them. These novel OTUs were contained in a broad range of higher level taxa, including i) well-established clades such as Cryptomycota, ii) in phyla rarely detected by cultivation-independent sequencing (e.g., Ichthyosporea) and iii) in novel clades previously undescribed in lacustrine ecosystems, such as Foraminifera.

Thus, according to this study, the OTUs representing the most abundant sequences were found among Fungi, Alveolata, Stramenopiles, Cryptophyta and Rhizaria. More precisely, the phylogenetic affiliation allows to delineate three of the four previously defined freshwater Cryptophyta clades [6]. Within the Fungi, numerous OTUs were associated with Cryptomycota [57] or Chytridiomycota, which include both parasitic and saprotrophic organisms [58]. The presence of Chlorophyta and Haptophyta was confirmed in most of the lake environments sampled in this study. By the CS method used for describing PCC, Chlorophyta and Haptophyta were often absent [59], [60] or found at a very low proportion [6], [37], whereas these phyla represented a significant proportion of PCC when counting methods such as FISH were used [61]. Such a bias has also been highlighted in marine environments since epifluorescence microscopy reveals a dominance of phototrophic or mixotrophic cells over heterotrophic cells [62]. Another example of phyla rarely described yet detected here is the Ichthyosporea phylum, which was found only in hyper-eutrophic conditions [63]. Finally, some clades supported by high bootstrap values in our phylogenies, e.g., Mamiellales or Foraminifera, seem original because they have not been detected by CS with 'universal' eukaryotic primers. To our knowledge this is the first time that a clade closely associated to Mamiellales, as defined by Marin and Melkonian [64], has been detected in lakes. Present but scarce in our pyrosequencing data, these microalgae constitute the dominant photosynthetic group among the picoplankton 18S rRNA gene sequences in marine surveys (∼ 1/3 of the sequences), especially in coastal waters, and have been shown to account for 45% of the picoeukaryotic community, as targeted by TSA-FISH in these waters [65], [66]. The freshwater counterpart of this group, the Monomastigales, is rarely recovered from environmental samples and likely requires new molecular approaches that will specifically target photosynthetic organisms in the environment [64]. Freshwater Foraminifera, a group of granuloreticulosan protists largely neglected until now have already been detected by using specific primers in one study of freshwater ecosystems [67]. Thus, a NGS sequencing analysis with a moderate depth (∼ 10,000 cleaned read per sample for Eukaryota) allows for the detection of the main phylogenetic phyla but also rarely detected phyla or phyla only detected by specific primers which act similar to massively parallel sequencing by focusing on one clade. Among the biases commonly assigned to CS, other than the variability in the cell lysis efficiency, the rRNA gene copy number, which range from 1 to 12,000 [68] is certainly the most important and may result in an over-representation of heterotrophic organisms notably of the alveolate taxa [34]. However, even if these differences in copy number distort the interpretation in number of reads and OTUs for both the CS and NGS methods, the massively parallel sequencing can at least increase detection of rare lineages or organisms with low gene copy numbers thanks to the increased depth of sequencing. We can hypothesize that this copy number could be more homogeneous at a specific lower taxonomic level (for example Alveolata), and the various indexes were therefore computed for each phylum instead of considering the whole protistan community (Table S3).


These results show that phylogenetic methods provide a clear edge in describing under-studied and complex communities, allowing the taxonomic affiliation of experimental sequences within an evolutionary framework; the study of relatedness among both environmental and reference sequences; and the evaluation of proximity of experimental sequences (“binning”). Thus, the tree-based method presented in this work, applied to the whole spectrum of microorganisms diversity (i.e., Eukaryota, Bacteria and Archaea), makes it possible to seek typical clades, allowing for the discovery of new putative lineages that are rarely or never recovered by classical sequencing approaches and the investigation of specific features within ecosystems considering sampling depths and periods. This feature cannot be inferred with a similarity search, a naïve Bayesian classification (RDP) or tree-based methods that process one sequence at a time.

Materials and Methods

The data originating from simulations and pyrosequencing were processed by a pipeline, referred to as PANAM (Phylogenetic Analysis of Next-generation AMplicons) that is based on publicly available programs. In addition to the phylogenetic analysis, this pipeline allows for the complete analysis of a full pyrosequencing run, including raw data processing, sequence clustering into OTUs and generating phylogenies for the taxonomic affiliation. The description of the procedure is detailed in the following sections (“Processing of raw pyrosequencing reads and OTU picking”; “Phylogenetic affiliation”; “Richness and diversity indexes”). It is written in Perl and can be run on Linux. The package comprises a reference sequence database, a taxonomy file and reference profile alignments and can be obtained from

Processing of raw pyrosequencing reads and OTU picking

The pyrosequencing reads can be cleaned according to different methods commonly used in the field of molecular microbial ecology. Pyrosequencing errors can therefore be reduced by removing the primers (e.g., [69]), defining a minimal score and length of the reads (e.g., [14]) or removing reads with unidentified bases (Ns).

Short sequences and sequences with low-quality scores are removed using PANGEA scripts [16] and only sequences with a primer match percentage above a defined threshold are selected using Fuznuc [70]. Alternatively, other quality filtering methods can be implemented; the platform does not depend upon the filtering approach described above. When several samples are analyzed, the checked sequences are split into different files depending on their bar code or tag. Then, generated files are clustered using USEARCH [71] at a user-defined threshold, and representative sequences from OTUs are selected for the phylogenetic assignment.

Phylogenetic affiliation

For the phylogenetic affiliation, a dedicated database of reference sequences, verified taxonomy and alignments was built using sequences extracted from the SSURef 108 database of the SILVA project [72]. For this purpose, all the sequences (16S and 18S rRNA genes) with more than 1,200 bp, quality score > 75%, and a pintail value > 50 were extracted. The sequence quality score defined by SILVA is a combination of the percentages of ambiguities, homopolymers longer than 4 bases and possible vector contaminations, and the pintail value corresponds to the probability that the rRNA sequence is chimeric. The complete database, after filtering according to the criteria above, contains 164,353 sequences (Archaea: 11,092; Bacteria: 131,428; and Eukaryota: 21,833) together with their taxonomy. To speed up the phylogenetic processing, the 3 domains were split into 37 phyletic groups of unicellular organisms corresponding to the first monophyletic clade after domains, as annotated in the guide-tree of SILVA (ARB format), and clustered at 97% identity.

Each profile corresponds to the first rank beneath that of domain. As the taxonomy of Bacteria and Archaea follow standardized taxonomic paths, the monophyletic profiles of these two domains correspond to phylum, the first level occurring after the domain. For Eukaryota domain, the taxonomy does not necessarily fit this organization, and the position of the taxon in the taxonomic hierarchy does not imply rank as it is the case with Bacteria and Archaea. Therefore, for the eukaryotic profiles, we opted for the rank position (the first one after the eukaryotic domain) and the monophyly, regardless to the taxonomic level.

For each of the 37 phyletic groups, an outgroup containing one sequence from each other group belonging to the same domain plus 2 external sequences were added to the alignment to root the phyletic tree to be produced and to specify the relatedness of early diverging sequences from the root of the group. To broaden the targeted diversity, the user can add specific environmental sequences to the database and the profiles.

Using this dedicated database, the phylogenetic affiliation is carried out following the different stages described in the Figure 5.

Figure 5. Flow chart describing the phylogenetic affiliation.

A primary classification, sorts and splits reads into groups according to the taxonomy of their best USEARCH hit (1). Next, a file containing aligned reads and sequences from the corresponding group is generated by processing a profile alignment by HMMER. This file is used by FASTTREE to build a phylogenetic tree (2), which is then parsed to assign a taxonomy to each read and to report putative clades (3).

1- First, OTUs are compared against the reference database described above with USEARCH [71]. As this first step does not intend to provide an exact affiliation, but rather to give a first approximation to perform a rapid and accurate phylogenetic analysis, the query sequences are sorted according to the taxonomy of their best hits, whatever their similarity score. Several files are generated, each containing the reads and their 5 best hits, assigned to one of the 37 specific phyletic groups.

2- After reads have been assigned to phyletic groups, they are aligned to the reference sequences of the corresponding profile alignment for that group using hmmalign from the HMMER package [73]. Synthetic files, which include the reference sequences and the aligned experimental reads, are generated.

3- Using FASTTREE [74], a bootstrapped phylogenetic tree (100 iterations) is built for each phyletic profile, including OTUs associated with their 5 best hits and the reference sequences. The trees are then parsed to generate files containing the taxonomy of the inserted sequences and files reporting the clades that could be identified from reads forming monophyletic groups. Two methods for taxonomy assessment are implemented: lowest common ancestor (LCA) and nearest neighbor (NN). In this last method, for each query sequence, all the nodes containing the sequence are scanned from the most recent to the deepest. The closest neighbor is defined as the first referenced sequence starting from the lowest node. The query sequence will acquire the complete taxonomy of its nearest neighbor. For LCA [32] each node holds only the common taxonomy between all of its descendants and thus may be incomplete. Each query sequence will inherit the taxonomy of its lowest node. The final taxonomy assignment is based on the phylogeny. The relatedness between all sequences (both experimental and referenced) are re-evaluated, and the similarity based assignments proposed on stage 1 are therefore revised to provide a more phylogeny-driven affiliation. Regarding the clades, their definition differs according to authors (e.g., [75], [76]), although in general, a new clade is declared when the cluster contains environmental sequences from at least 3 different sources and is supported by bootstrap values generally higher than 70%. The files generated describe monophyletic clusters with all the information required for experts in the field to define a putative environmental clade: a bootstrap value, a list of all the experimental sequences affiliated to it and the nearest reference neighbour together with its taxonomy. The implementation of PANAM (files generated) is extensively described in the documentation associated with the pipeline.

Richness and diversity indexes

After the cleaning step, richness (Chao1 and ACE), diversity (Shannon) indexes, and coverage are computed for each sample [77]. Subsequently, sequence library sizes are equalized to avoid biases associated with different sampling depths (e.g., [78]). Briefly, the same number of sequences (i.e., the number of sequences in the smallest sample) are randomly sampled from each library, and diversity indexes are calculated for these equalized datasets. After phylogenetic affiliation, Chao1 and the Shannon diversity indexes are computed for levels 2 and 3 from the EMBL classification (e.g., Stramenopiles and Bacillarophyta ).

Analysis of sequencing data obtained from simulations

PANAM was first tested on near full-length sequences with known taxonomy using 5 sets of 1000 sequences randomly picked from the reference database and removed from it for evaluations to be re-affiliated. The reliability of PANAM taxonomic affiliations was evaluated for specificity defined as the proportion of ranks correctly affiliated among the detected ones. A pyrosequencing simulation was also performed with pseudo-reads being generated by clipping the 5 × 1000 full-length sequences datasets from 6 universal forward primers for Eukaryotes [79] (Table S2). Clipped sequences were extended 200 and 400 bp from the forward primer positions defined on the Saccharomyces cerevisiae sequence (V01335), thus covering regions with different variability along the 18S rRNA gene. As emphasized, this pipeline allows taxonomic affiliations within an evolutionary context: its performance was thus primarily compared with that of STAP (Small Subunit rRNA Taxonomy and Alignment Pipeline) [51], the phylogenetic affiliation method used in WATERS (Workflow for the Alignment, Taxonomy, and Ecology of Ribosomal Sequences) [41], but was also compared with non-phylogenetic methods, including BLAST and the RDP Classifier implemented in MOTHUR [19] trained on the near full-length and trimmed sequences of the reference database.

The computational load of the phylogenetic analyses using PANAM was also tested with increasingly large datasets to evaluate processing time on a personal computer and to detect any scaling issues.

Analysis of sequencing data obtained from environmental studies

The PANAM tree-based method was run on environmental sequences, namely i) a set of environmental sequences originating from published studies on the diversity of protists and belonging to described environmental lacustrine clades of Perkinsozoa and Cercozoa [6], [7] and ii) from an environmental survey of the lacustrine protist diversity performed in eight freshwater ecosystems.

For this purpose, eight lakes or reservoirs, described in Table 2 (Lakes Anterne, Aydat, Bourget, Godivelle, Geneva, and Pavin, and Reservoirs Sep and Villerest), were sampled once during their thermal stratification (from May to August according to the lake). Water samples from the epilimnion (1 to 5 m) were collected with a Van Dorn bottle at a permanent station (the deepest zone of the lake). Water samples (from 100 to 120 ml) were successively filtered through 5 µm-pore-size and 0.2 µm-pore-size polycarbonate filters (Millipore), and the membranes were stored at-80°C until nucleic acid extraction. All samples were extracted following the protocol described previously by Lefranc et al. [37].

The V4-V5 variable region of eukaryotic 18S rDNA was amplified with primers Ek-NSF573 and Ek-NSR1147 (Table S2). To discriminate each sample, a 5 bp multiplex tag was coupled with the Roche 454 pyrosequencing adaptor A. The amplification mix (30 µl) contained 30 ng of genomic DNA, 200 µM of deoxynucleoside triphosphate (Bioline, London, UK), 2 mM MgCl2 (Bioline), 10 pmol of each primer, 1.5 U of Taq DNA polymerase (Bioline) and the PCR buffer. The cycling conditions were an initial denaturation at 94°C for 10 min followed by 30 cycles of 94°C for 1 min, 57°C for 1 min, 72°C for 1 min and 30 s and a final 10-min extension at 72°C. Finally, the amplicons of all of the samples were pooled at equimolar concentrations and pyrosequenced using a Roche 454 GS-FLX system (Titanium Chemistry) by GATC (Konstanz, Germany). The reads, alignments and trees have been deposited in Dryad ( The reads used in this study were selected from a full run, separated into bins according to the tags, analyzed by PANAM, using trimming criteria of quality score > 22 and sequence length > 200 bases and clustering into OTUs with a 95% similarity threshold. UNIFRAC metrics [30] and a principal coordinate analysis were used to compare the small protist community between the lakes based on phylogenetic information obtained by PANAM using the packages Picante and ade4 implemented in the R software [80].

To broaden the covered diversity, more specifically regarding the environmental and pyrosequencing datasets processed in this study, and to build phylogenies with more similar sequences for the studied environment, 173 sequences from eukaryotic clades specific to lacustrine ecosystems, defined in previous works (e.g., [6], [37]), were introduced in the eukaryotic reference database and the corresponding groups.

Author Contributions

Conceived and designed the experiments: NT GB DD. Performed the experiments: NT JFM ID GB DD. Analyzed the data: NT GB JFM ID DD. Contributed reagents/materials/analysis tools: NT ID GB DD. Wrote the paper: ID GB DD.

Categories: 1

0 Replies to “Ssu Rrna Analysis Essay”

Leave a comment

L'indirizzo email non verrà pubblicato. I campi obbligatori sono contrassegnati *