Improved taxonomic assignment of rumen bacterial 16S rRNA sequences using a revised SILVA taxonomic framework

Rumen samples and microbial community data used in analyses

Sequencing data generated using 454 GS FLX Titanium chemistry from the 684 samples in the GRC study that generated usable sequence data for bacteria (Henderson et al., 2015), which covered a wide range of ruminant species on different diets, were used in this study. The samples were all processed for sequencing of bacterial 16S rRNA sequences as previously described (Henderson et al., 2015; Kittelmann et al., 2013). Raw sequencing reads (BioProjects PRJNA272135, PRJNA272136, and PRJNA273417 in NCBI’s Sequence Read Archive; Leinonen, Sugawara & Shumway (2011)) were processed in QIIME (Caporaso et al., 2010) using standard parameters unless indicated otherwise. The primers used generated amplicons containing the V1-V3 regions of the 16S rRNA gene (Henderson et al., 2015). During the split library process, sequences that were at least 400 bp long were retained (input arguments -l 400, -L 1000 -r–z truncate_remove), covering at least the V1 and V2 regions of the 16S rRNA gene. The resulting 4,557,252 sequencing reads were concatenated and grouped into 774,769 operational taxonomic units (OTUs) using UCLUST with a 97% similarity definition criterion. Taxonomic identities were assigned to the repset sequences in QIIME using a BLAST (Altschul et al., 1990) search against either the SILVA database version 119 (Quast et al., 2013), or version 119Rum (this study, see below), Greengenes (version 13_8; McDonald et al., 2012) or RDP training set versions 14 and 16 (obtained from http://www.mothur.org/wiki/RDP_reference_files). Data were summarized at the genus level. Additionally, taxonomic identities were assigned against RDP’s bacterial 16S rRNA gene dataset using Classifier release 11.4 (Cole et al., 2014; Wang et al., 2007) and a confidence cut-off of 80% to summarize data. The number of sequence reads assigned to each OTU was used to assign abundance to each repset sequence.

Improving the resolution of the SILVA taxonomic framework

Operational taxonomic units that contained at least 50 sequences at the 97% similarity cut-off were selected for further analysis (i.e., 8,985 or 1.16% of all OTUs, representing 2,833,335 or 62.2% of all GRC sequences). These representative OTU sequences were aligned with SINA version 1.2.11 (Pruesse, Peplies & Glöckner, 2012) using the SILVA SSURef database version 119 as a reference alignment. This contained 16S rRNA (gene) sequences ≥1,200 nt. The aligned representative OTU sequences were imported into ARB version 6.0.2 (Ludwig et al., 2004), together with information on their abundance and prevalence in rumen samples, and they were added to the guide tree using the inbuilt Escherichia coli filter for positions 1,044–11,892 (corresponding to E. coli 16S rRNA gene positions 28–514). The SILVA databases and associated guide trees are rigorously curated, and sequence quality inclusion criteria, guide tree construction, and maintenance are described in detail elsewhere (Pruesse et al., 2007; Glöckner et al., 2017). The 16S rRNA gene sequences of isolates included for genome sequencing in the Hungate1000 project (Seshadri et al., 2018) were also aligned and added to the guide tree as described above. We assigned taxonomic designations down to genus-equivalent level, where possible, to all predefined monophyletic bacterial groups in the guide tree that were found to contain rumen bacterial sequences. We used the following similarities between sequences in a monophyletic radiation to define a genus-level (>94.5%) or family-level (>86.5%) group (Yarza et al., 2014). Clusters that did not also contain sequences of described species were named to the lowest taxonomic level that could be assigned and a strain, clone, or uncultured genus-level group (UCG) identifier included in the name as appropriate (e.g., Bacteria; Bacteroidetes; Bacteroidia; Bacteroidales; Prevotellaceae; Prevotellaceae NK3B31 group or Bacteria; Firmicutes; Clostridia; Clostridiales; Ruminococcaceae; Ruminococcaceae UCG_001). Taxonomic designations were independently reviewed by three of the authors, and once agreed upon by all were incorporated into the SILVA 119Rum database used for testing (see below), and then incorporated into SILVA SSURef database release 123 (http://www.arb-silva.de/documentation/release-123/).

Assessing the taxonomic assignments generated by different taxonomic frameworks

The 50 most abundant and 50 most prevalent bacterial OTUs were used as a test subset to compare the differences in taxonomic assignments made using different frameworks. Abundance was defined as the average proportion of an OTU in all 684 samples. Prevalence was defined as the proportion of the 684 samples that the OTU occurred in. This subset contained 77 unique OTUs because some were both abundant and prevalent. Sequence alignments of these OTUs and the closest type strain and other cultured relatives based on sequence similarity were manually curated and used to generate similarity matrices. Sequence similarity cut-offs recommended by Yarza et al. (2014) were used to identify likely identities of OTU sequences at genus (94.5%), family (86.5%), order (82.0%), class (78.5%), and phylum (75.0%) levels, as these provide a unified basis for the classification and nomenclature of uncultured bacteria that is compatible with the taxonomy of cultured bacteria (Yarza et al., 2014).

Revising the nomenclature of genus-level clusters for rumen bacteria

We first identified the radiations of abundant rumen bacteria using the SILVA 119 framework. To do this, we placed the repset sequences from the 8,985 largest OTUs, that is, those with the greatest number of sequence reads assigned to them, from the GRC dataset into the SILVA 119 tree. These fell into 241 pre-defined clusters in the guide tree provided with the SILVA 119 database. Later (see below), when the entire GRC dataset was reanalyzed after these clusters had been given unique taxonomic identifiers (again, see below), 97.8% of all sequences in the GRC dataset were assigned to these 241 clusters, showing that these potentially covered a large part of rumen bacterial diversity (Table S1).

Next, we examined these clusters based on only the sequences in the SILVA 119 database, that is, excluding the short reads from the GRC dataset. These clusters therefore contained 16S rRNA (gene) sequences ≥1,200 nt. Of these 241 clusters, 99 had no unique or had inadequate taxonomic identifiers in the SILVA 119 framework. For example, there were 12 distinct radiations named “Lachnospiraceae incertae sedis,” 16 named “Ruminococcaceae uncultured,” and three named “Coprococcus” (Table S1). This means that sequences matched to them would be assigned to an undefined or unclassified group at higher than the genus level, or they would be given a genus assignment that would be incorrect because the named reference sequence was >5.5% different to that of the type species of the genus. Here, we applied the genus-level sequence to sequence similarity cut-off of 94.5% identity proposed by Yarza et al. (2014). Each of these 99 clusters was given a unique nomenclature based on named species in the cluster, current trivial names used in the literature, or the names of isolates or other sequences in the cluster. Based on the level of sequence similarity and the level of separation from other named groups, these were provisionally given genus- or family-level status. Examples include the definition of one of many radiations designated as “Lachnospiraceae uncultured” as Lachnospiraceae group NK4B4, named after an isolate that falls into this genus-level group, and the distinction of multiple radiations previously all designated as “Erysipelotrichaceae uncultured” as Erysipelotrichaceae UCG001 to Erysipelotrichaceae UCG010, where UCG indicates a genus-level cluster without recognized cultured isolates (Table S1).

The names of an additional 11 clusters (containing 13.9% of GRC sequences) were modified to clarify their taxonomic positions, but the original names were unique and so these changes had no material impacts on the classification of sequences assigned to them. No changes were made to the nomenclature of the remaining 131 clusters, containing 20.1% of GRC sequences. Where a group contained named isolates that were clearly not in the same genus-level group as the type species of the genus, the genus name was retained but written in square brackets to indicate this. An example is the [Eubacterium] ruminantium group (Table S1), which is in the family Lachnospiraceae, and not close to the type species of the genus Eubacterium, E. limosum, which is in the family Eubacteriaceae (Ludwig, Schleifer & Whitman, 2009).

All of these new taxon designations were linked to sequences available in the SILVA SSURef database release 123 (SILVA 123 framework, http://www.arb-silva.de/documentation/release-123/). A summary of these designations is given in Table S1.

Comparison of taxonomic assignments generated using different taxonomic frameworks

To allow us to assess the impact of the improved taxonomic resolution, we produced a temporary version of SILVA119 that included these new designations, which we called SILVA 119Rum. We compared the taxonomic assignments made using this new framework with those made using four other frameworks, namely Greengenes 13_8, RDP release 11.4 (using RDP Classifier), RDP training sets 14 and 16, and the parent, SILVA 119. For our comparisons, we used all 774,769 OTUs generated from the 4,557,252 partial bacterial 16S rRNA gene sequences that had been generated from 684 rumen samples in the GRC dataset, rather than just the 8,985 largest OTUs that we used to identify the bacterial groups that needed refinement. The repset sequences from this larger set of OTUs were assigned to between 723 and 1,441 uniquely taxonomic strings using the different frameworks (Fig. 1A), although not all were resolved to the genus level (Fig. 1B). The proportion of sequences in the GRC dataset that were assigned to a taxon at any rank that was uniquely labelled within the preceding rank was also calculated (Fig. 1B). For example, genus-level groups labelled uncultured or incertae sedis without further qualifiers were not considered unique, since it is not possible to tell how many genus-level groups with the same identifier there might be in a family. There were between 557 and 1,355 identifiable genus-level taxa (Fig. 1A) using the different frameworks.

Overall, assignments made using RDP release 11.4 were the most conservative, returning the smallest number of taxa (Fig. 1A) and assigning the lowest proportion of sequences to a taxon at any taxonomic rank (Fig. 1B). The Greengenes taxonomic framework resulted in the second most conservative assignment of taxonomic identities to sequences, followed by the SILVA 119 and SILVA 119Rum taxonomic frameworks. Greengenes and RDP release 11.4 resulted in fewer <50% and <30% of sequences being assigned to a named genus, respectively. The SILVA 119 taxonomic frameworks allowed 55.4% of sequence data to be assigned to 751 uniquely-named genera. The improved SILVA 119Rum framework resulted in 87.1% of all sequences in the GRC dataset being assigned to one of 871 uniquely identifiable genus-level groups.

The RDP training sets (RDP versions 14 and 16) resulted in >95% of sequences being assigned at the genus level when they were used as BLAST databases (Table S2). The RDP training sets predominantly contain sequences of characterized organisms whose taxonomic strings are nearly all resolved to the genus level. It contains few reference sequences from unnamed or taxonomically poorly-resolved groups. For this reason, these frameworks are not suitable for use with the QIIME “parallel_assign_taxonomy_blast.py” script. This script always assigns a match to a sequence, and by default that match will almost always return a genus name because the RDP training sets contain few sequences without a valid genus name. Sparse databases like the RDP training sets will result in an assignment to a match that may not be a close one. This results in an over-assignment of sequences to named genera to which the sequences do not belong; however, if these databases were used with Classifier as intended, it is likely these assignments would have been given low confidence scores (Wang et al., 2007) that would not result in a genus level assignment, as we observed (Fig. 1A).

Comparison of apparent microbial community composition using different frameworks

Apparent rumen microbial community compositions were broadly comparable at phylum level, regardless of the taxonomic framework used (Fig. 2A). However, at the genus level the apparent make-up within the dominant phyla Bacteroidetes (Fig. 2B) and Firmicutes (Fig. 2C) differed considerably. This was due to the different naming conventions, different nomenclature, and finer taxonomic resolution of groups in some of the databases. This means that data analyzed using different frameworks cannot be directly compared, as can be seen in Figs. 2B and 2C, where the same dataset gave very different apparent community structures using the different taxonomic frameworks. The impact of these differences is illustrated by the taxonomic identities assigned to individual OTUs using the different taxonomic frameworks. As examples, we show the taxonomic assignments of the 77 most abundant and prevalent OTUs in the GRC dataset made using the different frameworks (selected data are shown in Table 1; full data in Table S2). In some cases, the differences were in nomenclature, but in others it was due to the lack of resolution near the genus-level. Other examples of the taxonomic equivalence of groups with different names in the different frameworks are shown in Fig. 2.

Increased taxonomic resolution with SILVA 119Rum

More OTUs, and hence more sequences otherwise assigned only to the family level (such as Lachnospiraceae, Ruminococcaceae, and Succinivibrionaceae), were assigned to multiple new genus-level groups within these families using the SILVA 119Rum framework (Table 2). Other frameworks only assigned these OTUs to order or family levels, or to genera that are not monophyletic or contain sequences with sequences that have <5.5% identity. These finer scale subdivisions of families and genera are reflected in the greater number of genus-level taxa found in the GRC dataset classified using SILVA 119Rum (Figs. 1A, 2B and 2C).

The great diversity of rumen Lachnospiraceae and Ruminococcaceae that is not yet formally described and named has been reported previously (Creevey et al., 2014). Sequences from members of Lachnospiraceae made up >17% of the GRC dataset. Refinement of the taxonomy increased the genus level separation, with three quarters of those sequences assigned to more exact genus-level taxa (Table 3; Table S1). This included separation of undefined taxa named Lachnospiraceae incertae sedis and uncultured Lachnospiraceae. Named genera like Butyrivibrio, Blautia, and Coprococcus also appeared to contain sequences originating from multiple genus-level groups. Similarly, Ruminococcaceae contained >13% of all GRC sequences, and over half of these were classified in SILVA 119 as Ruminococcaceae incertae sedis or as uncultured Ruminococcaceae. These were separated into 20 new genus-level groups, one of which, Ruminococcaceae NK4A214 group, contained nearly 2.8% of all GRC sequences (Table S1). This group is named after isolate NK4A214, the first recognized isolate that belongs in this apparently new genus (Kenters et al., 2011).

We divided some genera, such as Prevotella and Ruminococcus, into multiple genus-level groups based on their degree of sequence divergence. Sequences assigned to Ruminococcus were split between two different genus level groups, Ruminococcus 1, characterized by R. flavefaciens and R. albus and containing 67% of all Ruminococcus sequences in the GRC dataset, and Ruminococcus 2, which contains R. bromii and 33% of all Ruminococcus sequences in the GRC dataset (Fig. 2C; Table S1). The potential separation of R. bromii from the species that fall into Ruminococcus 1 was suggested by an early 16S rRNA gene sequence analysis of this genus (Rainey & Janssen, 1995).

The genus Prevotella contained 22.0% of all sequences in the GRC dataset when originally analyzed (Henderson et al., 2015). Finer scale subdivision of this taxon into genus-level groups still resulted in 18.3% of all sequences falling into a genus we designated Prevotella 1, which contained the type species of the genus, P. melaninogenica, and also P. ruminicola. Four additional genus-level groups of sequences that previously fell into Prevotella were identified and given unique designations: Prevotella 2, Prevotella 6, Prevotella 7, and Prevotella 9. Six further genus-level clusters, all previously named “uncultured,” were given unique designations, and so increased the resolution within the family Prevotellaceae (Fig. 2B; Table S1).

Other groupings defined in the SILVA 119Rum framework allowed assignment of sequences to taxa that were not included in the other databases. The RC9 gut group was absent from the Greengenes and RDP release 11.4 frameworks, meaning this group, making up on average 5.7% of sequences, was only identified when the SILVA 119 or SILVA 119Rum databases were used. Otherwise they were assigned to the order Bacteroidales. Members of the family Christensenellaceae were prominent in the GRC dataset, making up 10.7% of all sequences when reanalyzed using the SILVA 119Rum framework. This group was poorly resolved until a cluster containing isolate R-7 was defined and given a unique identifier (Christensenellaceae R-7 group) in SILVA 119Rum (Table 2). This genus-level group contained 10.6% of all GRC sequences and 99.5% of all of the sequences that were assigned to the family Christensenellaceae. In some cases, sequences that were previously assigned to poorly-defined groups were added to taxa that had names, because the reference sequences they matched to were given the same family-level designation. Examples include the family-level Bacteroidales BS11 gut group and S24-7 group.