The success of microbiome studies (composition, structure, diversity, and function) is primarily ascribable to the development of bioinformatics tools embedded in creative algorithms specially tailored to overcome the technical challenges posed by the analysis of massively paralleled, high-throughput sequencing data (Simon and Daniel, 2011; Siegwald et al., 2017). These bioinformatics tools make use of several techniques (e.g., read mapping, k-mer alignment, and composition analysis) (Piro et al., 2017) and can be categorized into two distinct groups: (1) programs that use all available genome sequences (Lindgreen et al., 2016), also called assignment-first approaches (Siegwald et al., 2017) (e.g., CLARK – Ounit et al., 2015; GOTTCHA – Freitas et al., 2015; KRAKEN – Wood and Salzberg, 2014; MG-RAST – Meyer et al., 2008), and (2) programs that target a set of marker genes (Lindgreen et al., 2016), also known as clustering-first approaches (Siegwald et al., 2017) (e.g., QIIME – Caporaso et al., 2010; MOTHUR – Schloss et al., 2009; MetaPhlAn – Segata et al., 2012; mOTU – Sunagawa et al., 2013). In the assignment-first tools, all reads are assigned to the lowest taxonomy unit (lower common ancestor-LCA) within a reference database based on their annotations, while in the clustering-first approaches the reads are grouped into Operational Taxonomic Units (OTUs) using different OTU picking strategies (closed or open reference) to assign reads to a taxonomic group based on their sequence similarities (Siegwald et al., 2017).

However, most of the above studies are focused on demonstrating how single analytical steps (e.g., sequence pre-processing, OTU clustering or taxonomic assignment) generated by the existing tools impact the microbial classification in real or simulated datasets derived from the Human Microbiome Project (Siegwald et al., 2017). Comparison of methodologies to comprehensively classify the rumen microbiome is lacking which may be in part due to its complexity, as the rumen microbial community consists of bacteria, archaea, protozoa and fungi (Russell and Rychlik, 2001). A recent study by Li F. et al. (2016) developed a Mothur (Schloss et al., 2009) based pipeline to assess active rumen microbiota from data generated from total RNA sequencing. Later, the same researchers applied this pipeline to investigate linkages between the active rumen microbiome (structure and function) and feed efficiency in beef cattle using metatranscriptomics (Li and Guan, 2017). Using the developed mothur-based pipeline for taxonomic assignment, the authors identified that the active microbial taxa differed in the rumen of cattle with differing feed efficiency and suggested that the active rumen microbiome is one of the biological factors that may contribute to variations in feed efficiency in beef cattle (Li and Guan, 2017). There were two steps employed in taxonomic classification by Li F. et al. (2016): bacterial sequences belonging to V1–V3 regions were extracted from the aligned Greengenes database, and archaeal sequences belonging to the V6–V8 regions were aligned with a rumen-specific archaeal 16S rRNA gene database (Janssen and Kirs, 2008). Despite the efficacy of this pipeline, it still remains a challenge for researchers to determine which approach (assignment- or clustering-first methods) of taxonomic classification delivers the most realistic representation of rumen microbial ecology.

In the current study, we propose a comparative analysis of the outcomes of Kraken (Wood and Salzberg, 2014) and the pipeline of Li F. et al. (2016) with a focus on the biological interpretation of the rumen microbial classification from the perspective of two conceptually different software packages. Unlike the pipeline developed by Li F. et al. (2016), Kraken algorithms can make multiple comparisons of single or assembled k-mers against any hypervariable region, providing useful information regarding a particular species detected in a region of the 16S rRNA gene that is different from the targeted internal conserved region initially sequenced (Wood and Salzberg, 2014; Valenzuela-González et al., 2016). Although Kraken algorithms have been originally designed to assign taxonomic identity to short DNA reads (Wood and Salzberg, 2014), studies have shown that Kraken is also useful to provide taxonomic classification for long (up to 1352.1 ± 153.72 bp) metagenomic DNA sequences (Valenzuela-González et al., 2016). Therefore, our objectives were (i) to compare and contrast the pipeline of Li F. et al. (2016) and Kraken to assess the taxonomic profiles of rumen bacteria and archaea and (ii) to investigate the impact of the comparative analysis of both analytical approaches on the biological interpretation of the rumen microbial classification obtained from cattle exhibiting different feed efficiencies.

In this study, the comparison of taxonomic outcomes of two pipelines, Mothur (developed by Li F. et al., 2016) and Kraken (developed by Wood and Salzberg, 2014, and adapted to the conditions of this study), was performed to determine which is a better approach in rumen microbial classification when total RNA-seq data were used. The advent of high-throughput sequencing has greatly advanced our knowledge of the ecology and functional capacity of rumen microbes and their role in converting low-quality and unusable feedstuffs into energy sources for host productivity (McCann et al., 2017). As a result, an assiduous effort has been made to unveil the linkage between the rumen microbiota and phenotypic traits of interest such as feed efficiency (Li and Guan, 2017), enzyme discovery (Qi et al., 2011) and methane emissions (Kittelmann et al., 2014; Shi et al., 2014; Kamke et al., 2016). Metagenomic studies have shown that the host may regulate the microbiota and its metabolic activity in relation to feed efficiency (FCR) through host-microbiome cross talk genes such as TSTA3 (GDP-L-fucose synthetase) and Fucl (L-fucose isomerase), suggesting that the relative abundance of these genes could be used as a predictor for host feed efficiency (Roehe et al., 2016). Although the number of rumen metagenomics and metatranscriptomics studies has grown enormously over the last couple of years (McCann et al., 2017), the functional outcomes and biological interpretation of omics data strongly depend on the computational methods used (Simon and Daniel, 2011; Siegwald et al., 2017). In this study, both Mothur and Kraken pipelines showed the rumen of the bulls to be dominated by Prevotella, Treponema, Ruminoccocus, Fibrobacter, and Butyrivibrio, which are considered as part of a “core bacterial microbiome” (Henderson et al., 2016). In addition to the mutual “core microbiome” shared by the two pipelines at the genus level, Kraken detected a relatively high abundance of (1) Prevotella ruminicola (Supplementary Table 2), which is involved in the ruminal digestion of hemicellulose and pectin (Marounek and Duskova, 1999); (2) Fibrobacter succinogenes (Supplementary Table 2), a gram-negative, fiber degrader species (Suen et al., 2011); and (3) non-motile species within the Ruminoccocus genus (Supplementary Table 2) that share different niches (La Reau et al., 2016): R. bicirculans, which selectively utilizes hemicelluloses but not cellulose or arabinoxylan (Wegmann et al., 2014), and R. albus, which is capable of digesting cellulose and xylan (Christopherson et al., 2014).

Interestingly, both methods identified about 1% of Cyanobacteria (Supplementary Tables 1, 2), corroborating the findings of previous studies that have reported low abundances of these oxygenic phototrophic bacteria in the rumen of dairy (Scharen et al., 2017) and beef cattle (Li and Guan, 2017), and of camels (Gharechahi et al., 2015). Cyanobacteria are aerobic bacteria that can perform carbohydrate fermentation in a deficient N₂ concentration (heterocystous) or in a combination of N₂ deficiency and anoxic conditions (non-heterocystous) (Nandi and Sengupta, 1998). Although the ruminal environment is widely considered to be anaerobic, significant concentrations of O₂ (60 and 100 nmol/min per mL) can be detected in the rumen fluid (Newbold et al., 1996), indicating that the presence of Cyanobacteria in the rumen may be related to O₂ scavenging and sugar fermentation performed under restrict aerobic conditions. It is important to mention that although Cyanobacteria has been widely detected in aqueous and soil environments (Williams et al., 2004; Cruz-Martinez et al., 2009), the identification of this phylum in the mammals’ gut has raised critical questions on what roles these organisms may play in aphotic and anaerobic habitats (Soo et al., 2014) like the rumen. Recent researches have reported that gut Cyanobacteria are highly conserved but their 16S rRNA gene phylogenetic tree differed from the photosynthetic Cyanobacteria, which led to the designation of a new candidate class called Melainabacteria (Soo et al., 2014) whose members are capable of fermenting a range of sugars (e.g., glucose, fructose, sorbitol) into acetate and butyrate in the gut (Di Rienzi et al., 2013). Neither Kraken nor Mothur identified Melainabacteria in the samples, demonstrating that further studies are needed to disentangling its role in the rumen.

However, the two methods (Kraken and Mothur) generated microbial classification at different taxonomic levels for rumen bacteria. To completely understand the function of the rumen microbiota, it is essential to identify organisms at the species level since different species, within the same genus, can have varied functions and niches. The Mothur based method was useful to identify a diverse bacterial microbiota from the RNA-seq datasets, but it was not able to classify any of the bacterial sequences further than the genus level (Tables 1, 2). Microbial classification up to the species level is a major challenge for clustering-first approaches based on targeted regional 16S rRNA when short (up to 250 bp) or even longer reads generated from total RNA-seq are used to identify environmental microbes (Xiang et al., 2017). Most existing tools (for bacteria and archaea) lack solid probabilistic-based criteria to evaluate the accuracy of taxonomic assignments to determine the best-matched database hits to distinguish multiple species from the targeted sequence region of the 16S rRNA gene (Xiang et al., 2017). To identify bacteria at the species level, sequencing of full length of 16S rRNA is desired and thus future studies need to increase the sequence length to enhance the resolution for microbial identification. For the Kraken based approach, the reference database was built based on all known microbial genomes and as a result it generated a higher resolution (to the species level) of the rumen microbiota, enabling the program to annotate each microbial sequence to the LCAs (Wood and Salzberg, 2014). In this process, k-mer paths formed by Kraken assign a specific weight to each node (equal to the number of sequences associated with the node’s taxon) while increasing the sensitivity of the species classification even if regions (for example, V3–V5) of the 16S rRNA gene were analyzed (Wood and Salzberg, 2014; Valenzuela-González et al., 2016). Consequently, the generation of chimeric trees using short or long input sequences is improbable with Kraken as unlike other programs (such as Ribosomal Database Project classifier and Mothur), it leaves out specific sequences if there is insufficient evidence for classification and they are designated as unclassified (Valenzuela-González et al., 2016). Therefore, inputting short or long environmental sequences (containing most of the 16S hypervariable regions or mRNA sequences) into Kraken may generate a more representative profile of complex microbiomes (Valenzuela-González et al., 2016) like the rumen. However, the lack of reference genomes for rumen microorganisms also limits Kraken. For example, the classification of Xenorhabdus (Table 2) and Xenorhabdus doucetiae [data not shown; relative abundance (%): H-FCR, 0.1 ± 0.10 found in 6 samples; L-FCR, 0%], a motile, gram-negative soil bacterium usually described as being part of entomopathogenic nematode/bacterium symbiotic complex (Furgani et al., 2008) has not been previously reported in amplicon based sequencing (Li F. et al., 2016) or metagenomic/metatranscriptome sequencing (Li and Guan, 2017) of rumen contents. The classification of this bacterial species may indicate that Kraken did not properly identify the microbe since the reference genome information was built mostly from all microbial genomes annotated in the NCBI database. However, these organisms may have been actually detected in the rumen since cattle can consume soil, raising the possibility that their detection was transitory.

It is noteworthy that Methanobrevibacter (family Methanobacteriaceae) was identified in both databases (Supplementary Tables 1, 2, and 4). This genus has been reported to be the most abundant archaeal population in the rumen based on DNA datasets (Kittelmann et al., 2013; Henderson et al., 2016), but it had a lower abundance than Methanomassiliicoccaceae at the RNA level in this study. This result is consistent with the research conducted by Li F. et al. (2016), who reported a predominance of Methanomassiliicoccaceae over Methanobrevibacter in RNA-based datasets when compared to DNA Amplicon-seq outcomes, suggesting that Methanomassiliicoccaceae may be more active in the rumen than Methanobacteriaceae. However, further studies are needed to determine whether the differences in abundance between those two archaeal populations have a methodological influence or are controlled by diet, host animal or management strategies. Unlike bacterial classification, Kraken and Mothur generated contrasting results on archaea identification (Table 3), which reflects the divergent taxonomic profiles at the species level. For example, certain archaeal genomes, such as Methanobrevibacter wolinii and Methanobrevibacter woesei, were only found in the rumen-specific archaea database, as the Kraken standard database lacked these complete genomes. However, Kraken was able to detect Candidatus Methanoplasma termitum and Candidatus Methanomethylophilus alvus, which were not identified by Mothur pipeline. Li Y. et al. (2016) isolated the archaeon ISO4-H5 (member of the order Methanomassiliicoccales) from the sheep rumen and discovered that this archaeal taxon exhibited genome size (1.9 Mb) and GC content (54%) similar to Candidatus Methanoplasma termitum (enriched from the termite gut) and Candidatus Methanomethylophilus alvus (enriched from human feces). These two species encode pathways required for hydrogen-dependent methylotrophic methanogenesis by reduction of methyl substrates, without the ability to oxidize methyl substrates to carbon dioxide (Li F. et al., 2016). Thus, it is possible that these microbes reside in the rumen. Future analysis with archaeon ISO4-H5 sequences included in the databases of both pipelines as well as its isolation, culture and characterization may provide further evidence of this possibility.

To further verify how these two methods affected data interpretation, the rumen microbiota of H-FCR and L-FCR bulls were compared based on the taxonomic outcomes generated by the two software packages. Both computational pipelines revealed differences in microbial abundance between H- and L-FCR groups at all taxonomic ranks, with Mothur exclusively identifying a higher abundance of poorly characterized bacterial phylotypes (e.g., R4-41B) in L-FCR bulls (Table 2). It has been reported that the abundance of R4-41B was negatively correlated with production traits over the first 12 weeks postpartum in dairy cows (Lima et al., 2015), suggesting that it may have undesirable impacts on the function of the rumen microbiome of L-FCR cattle. Although Kraken identified a relatively higher abundance of Xenorhabdus in H-FCR bulls (Table 2), this result could be erroneous with further validation needed as described above. However, researchers have enumerated and identified a high number (15.7 × 10⁴ Most Probable Number/g) of chlortetracycline resistant Enterobacteriacea in cattle feces that largely consisted of Xenorhabdus doucetiae (Watanabe et al., 2016). Since antimicrobial agents (e. g., chlortetracycline) are typically administered subtherapeutically to beef cattle (Inglis et al., 2005), our results suggest that H-FCR animals may be more susceptible to harbor chlortetracycline resistant bacteria than L-FCR animals in the event of a therapeutic administration of this antibiotic. Further investigations aiming to evaluate the effects of antimicrobial agents (e. g., chlortetracycline) on the development of antimicrobial resistance in Xenorhabdus recovered from less efficient cattle (H-FCR) are warranted. Kraken also detected a higher (P = 0.09) abundance of Methanococcaceae [relative abundance (%): H-FCR, 13.6 ± 8.96; L-FCR, 4.1 ± 4.85] in the rumen of H-FCR bulls, indicating that Methanococcaceae may play a potential role in the linkages between methanogenesis and reduced feed efficiency in cattle. Although RNA-targeted DNA probes and genomic DNA sequencing have revealed a significant population of this archaeal family residing in the rumen (Janssen and Kirs, 2008) and exhibiting a positive correlation with increased forage content in the diet (Pitta et al., 2016), members of this methanogenic archaea family still need to be cultured from the rumen to test our findings.

Finally, our study demonstrated that both pipelines (Mothur and Kraken) were effective in detecting a lower bacterial diversity in efficient (L-FCR) cattle (Table 4), corroborating the recent findings by Li and Guan (2017) and (Shabat et al., 2016) that the rumen microbiota of efficient cattle is less complex and more specialized in harvesting energy from the diet through simpler metabolic networks (e.g., acrylate pathway) than inefficient cattle. However, only Kraken identified a significantly lower diversity in the archaeal community in L-FCR bulls, but this result should be carefully interpreted as many archaea phylotypes classified by Kraken are environmental organisms that have not yet been described in the rumen. For example, the methane-producing archaeon Methanothermococcus okinawensis (the third-most abundant archaea taxon classified by Kraken, Supplementary Table 2) was first isolated from a deep-sea hydrothermal vent system (Takai et al., 2002), Picrophilus torridus and Acidilobus saccharovorans (the fourth and fifth-most abundant archaea taxa detected by Kraken, Supplementary Table 2) were isolated from a dry solfataric field (Fütterer et al., 2004) and a terrestrial acidic hot spring (Mardanov et al., 2010), respectively. Thus, it is worth mentioning that, in spite of the Kraken’s promising results, this pipeline is severely limited when studying a microbiome that is not well described in its standard database (like the rumen), indicating that Mothur (using a specific archaea database described by Li F. et al., 2016) could be more suited for identifying archaeal taxonomic profiles.

The current study is the first to compare the molecular-phylogenetic outcomes of Mothur and Kraken using transcriptomic sequence data (∼140 bp in length) of rumen samples. The Kraken pipeline has been adapted to include reference genomes for rumen specific organisms, which has led to the identification of rumen bacteria at species level and more bacterial phylotypes. However, the results of the archaeal classification as well as some of the bacterial species identified by Kraken should be carefully interpreted as many detected phylotypes have not yet been described in the rumen, highlighting the importance of strengthening the Kraken database through the inclusion of more genomes annotated by single cell sequencing of rumen cultures/isolates to enable a more accurate classification. As to the future directions, we plan to include new sequenced genomes (410 draft bacterial and archaeal genomes) by Hungate1000 project (JGI database) into Kraken standard database and the recently developed Rumen and Intestinal Methanogen Database (Seedorf et al., 2014) for further analysis, with the goal of improving the accuracy of the results. We also propose the configuration of a joint pipeline using both Kraken and Mothur simultaneously to improve the resolution of taxonomic profiling of the rumen microbiome. This joint pipeline will produce a final rumen microbial profile obtained from the combination of multiple results generated from different bioinformatics tools as outlined by Piro et al. (2017), who published a computational method called MetaMeta that executes and integrates results from six metagenomic analysis tools (CLARK – Ounit et al., 2015; DUDes – Piro et al., 2016; GOTTCHA – Freitas et al., 2015; KRAKEN – Wood and Salzberg, 2014; KAIJU – Menzel et al., 2016; and mOTUs – Sunagawa et al., 2013). If the rumen microbiome datasets are strengthened to the same level as the human databases, the joint pipeline will generate more sensitive and reliable results than those of the best single profile (generated separately by each tool) (Piro et al., 2017). We also believe that a joint pipeline supported by a collection of tools could be useful to control sources of variation present in any metagenomics/metatranscriptomic analysis (e.g., analytical pipelines, related databases and software parameters), which will ultimately lead to standardized results and more reliable biological interpretations. In addition, although Kraken has improved the taxonomic assessment at species level, the high number of unclassified sequences (65%) suggests a need for identifying rumen microbes with a more resolved taxonomic classification. Regardless of the approach we undertake, the only way for improvement is through a continued strengthening of the databases by including additional information of whole genome sequencing of rumen isolates as well as single cell sequencing of unculturable rumen microbes, as the ability to culture rumen microorganisms is still limited.

AN and LG conceived and designed the experiment. AN and BG executed Kraken, and FL executed the pipeline based on Mothur. AN analyzed the data and performed the statistical analysis. AN, FL, and BG executed the experiment and wrote the manuscript. LG and TM contributed to the experiment and revised the manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.