The animal experimental protocol of this research was reviewed, discussed and approved by the institutional ethics committee of the Institute of Genome Research Institutional Review Board (IGR IRB), Vietnam (Approval number: No. 03/QD-NCHG).

In previous study, The goats Co (five animals) and Bach Thao (five animals) living on natural hay in high rocky mountains at private goat farms at Ninh Binh, Thanh Hoa in Vietnam were selected and used for harvesting bacteria in the goats’ rumen [18].

The 39,579 ORFs from the metagenomic data that could be linked to NCBI taxonomy [18] were used for taxonomic assessment of the bacterial community structure in Vietnamese goats’ rumen. The ORFs related to lignocellulose degradation including 821 ORFs for CEs and PLs, 816 ORFs for cellulases, and 2,252 ORFs for hemicellulases [18] were used for assessment of cellulolytic bacteria as described in the following section.

In previous study [18], a classified group of 39,579 ORFs was obtained from comparison of 164,644 ORFs from the metagenomic data of bacteria in goats’ rumen against the NCBI non-redundant protein (NR) database using the BLASTX algorithm with an E-value of less than 10⁻⁵. In this study, the ORFs were subjected to the Metagenome Analyzer program (MEGAN) (version 4.6) [19] for taxonomic assignment into NCBI taxonomy using lowest common ancestor (LCA) algorithm. The LCA-based algorithm assigns genes to taxa in the way that the taxonomical level of the assigned taxon reflects the level of conservation of the gene. So sequences placed at the deeper nodes of the phylogeny represent more conserved genes relative to the genes on the leaves the tree that are more specific to that species. Thus all ORFs (represented in gene code) in the classified group were assigned in the detail in every taxon.

From 39,579 ORFs classified in MEGAN, we separated the taxon level: kingdom, phylum, class, order, family, genus, species in separated files. After that, we counted for abundance of each taxonomy level and then summarized for the most abundances in every levels. We also used the ratio of Firmicutes/Bacteroidetes to assess the bacterial structure in our goats’ rumen in comparison with the ratio in other bacterial consortia in different animal rumens in previous studies.

From 821 ORFs encoding CEs and PLs involved in lignocellulose pretreatment, 816 ORFs encoding cellulases, 2,252 ORFs encoding for hemicellulases had been annotated from 164,644 ORFs [18], we filtered all these ORFs accompanied with the corresponding gene codes and functional annotation to a new file. Based on the gene codes, we used formula Vlookup in exel file to assign the genes with the corresponding taxonomic levels from taxonomic profile of total 39,579 ORFs to generate a file with gene code, functional and taxonomical information combined. The every taxonomic level related to the lignocellulolytic genes were separated in to each file for counting abundance bacteria associated with GHs for cellulases, hemicellulases and CEs, PLs for lignocellulose degradation. We also used the ratio of Firmicutes/Bacteroidetes for assessment of bacterial role in lignocellulose degradation.

This project was deposited in the DNA Databank of Japan with the accession ID PSUB006562.

In this study, from ~9 Gb of rumen bacterial metagenomics data, 39,579 ORFs (accounted for 24.03%) were classified. The majority of the organisms belonged to bacteria (99.8%). The 39,501 bacterial ORFs were affiliated to 28 phyla, 41 classes, 95 orders, 181 families, 571 genera and 1,634 species (Table 1, Supplementary Table S1). The most dominant phylum was Bacteroidetes (63.5%) followed by Firmicutes (22.6%), Proteobacteria (7.5%), and Synergistetes (3.1%) (Table 1). Typical ratio of Firmicutes versus Bacteroidetes in the the bacterial community was 0.36:1 (Figure 1A). Further analysis showed that, within Bacteroidetes phylum, the most abundant order was Bacteroidales (60.8%). The most represented order within the Firmicutes phylum was Clostridiales (16.4%) followed by Selenomonadales (2.8%) (Figure 1). Synergistales—a unique member of Synergistetes phylum represented about 3.1% of the taxonomic classified ORFs and Aeromonadales—a member of Proteobacteria, were also abundant in this data. With regard to the taxonomy at genera level, Prevotella (35.3%) and Bacteroides (16%) were most abundant which belong to Bacteroidales. Analysis of the 13 most abundant species resulted in eight species belonging to Prevotella, one species belonging to Bacteroides, two species assigned to Clostridiales and only one species of Synergistales could be identified (Figure 1B).

In the previous publication, we indicated that lignocellulolytic genes were highly diverse in Vietnamese goats’ rumen [18] if compared to another study on goat rumen microbial diversity [4]. We supposed that the diversity is an adaptation of the rumen bacteria to goat’s diet. As such, we affiliated all the ORFs related to lignocellulose degradation in the goats’ rumen. As a result, among 816 ORFs genes for cellulases, 2,252 ORFs for hemicellulases and 821 ORFs encoding enzymes involved in lignocelluloses pretreatment, that were annotated by different databases [18], 221 genes for cellulases, 544 genes for hemicellulases and 226 genes for lignocellulose pretreatment were affiliated by NCBI taxonomic classification. That means about 24% to 27% lignocellulolytic genes were successfully assigned to taxon and a large portion of the genes could not be affiliated to taxonomy. In general, these gene groups were mostly found to be associated with the Bacteroidetes phylum (854 ORFs of the total 991 ORFs; accounting for 86.2%). The second most abundant phylum was Firmicutes (94 ORFs; accounting for 9.5% of ligocellulolytic genes) (Figure 2A). The ratio of Firmicutes/Bacteroidetes is therefore quite low (0.11:1) and lower than the ratio in bacterial structure (0.36:1). Looking at order level, the most lignocellulolytic associated ORFs were affiliated to Bacteroidales, and Clostridiales (Figure 2, Supplementary Table S2). Analysis in genus level, 398 ORFs were affiliated to Provotella, 300 ORFs were belonged to Bacteroides, 39 ORFs were predicted to be derived from Paraprevotella, 39 ORFs were originated from Ruminococcus. Specially Aerophaga was predicted to harbour ORFs for only hemicellulase (7 ORFs), while Algoriphagus was predicted to carry ORFs for only CEs and PLs.

When we analysed specific classes of lignocellulolytic genes, we observed divergent patterns of taxonomic distribution, depending on enzyme function. For pretreatment enzymes (CEs and PLs), the most abundant phylum is again Bacteroidetes (201 ORFs), followed by Firmicutes (14 ORFs) (Figure 2A). In other words, Bacteroidetes hold 88.9% genes for lignocellulose pretreatment, Firmicutes kept 6.2%. Thus the ratio between Firmicutes versus Bacteroidetes is very low (0.07:1). However, genes encoding for esterase enzymes showed a different taxonomic pattern. The 56 ORFs encoding CE1 were affiliated to many bacteria orders such as Bacteroidetes, Firmicutes, Planctomycetes, Spirochaetes, while CE1 domains collocated with a carbohydrate binding domain (CBM) CBM48 (10 ORFs) were only assigned to Bacteroidetes. This pattern was similar for CE6 and PL1, CE6 (21 ORFs), which were predicted to be produced in Bacteroidetes and Verrucomicrobia, while PL1 (45 ORFs) could be linked to Bacteroidetes and Firmicutes. However, the bi-fuctional enzymes CE6-GH95 and PL1 carrying CBM6 domain or a fibronectin like 3 domain (FN3), which help for increase the affinity between enzyme and corresponding substrate, were again assigned to Bacteroidetes only. Moreover, several enzymes such as CE12 (25 ORFs), CE7, CE7-CH26, CE8, CE8-Ig (immunoglobulin like domain), PL9, PL10 were predicted to be produced only in Bacteroidetes, while CE13 was only generated in Firmicutes (Figure 2B). At species level, Prevotella buccae was the most abundant species among lignocellulose pretreatment enzymes.

The most cellulase ORFs were affiliated to Bacteroidetes (153 ORFs; 69.2%), followed by Firmicutes (53 ORFs, ~24.0%), and then Proteobacteria (7 ORFs) (Figure 2A). Thus ratio of Firmicutes/Bacteroidetes is again low (0.35:1) but higher than the ratio in CEs and PLs producers (0.07:1). We found that Bacteroidetes and Firmicutes were predicted to produce many enzymes such as GH3, GH3-FN3, GH5, GH9, GH9-Ig. More specifically, GH3 is predicted to have beta-glucosidase activity, while GH5 and GH9 exhibit endoglucanase activity. Within cellulases, we also identified genes uniquely assigned to Firmicutes. For instance, GHs (GH6, GH5, GH9) accompanied with CBMs, such as CBM2, CBM3, CBM4, CBM63 were only identified among Firmicutes, while no GHs collocated with CBMs in this groups could be assigned to Bacteroidetes. Furthermore, GH1, GH44, GH5-FN3, GH74 were also uniquely associated with Firmicutes (Figure 2C). Further analysis of species harbouring cellulolytic genes showed that, Bacteroides uniformis, Ruminococcus bicirculans, Eubacterium siraeum were the most dominant to bear putative cellulases (Figure 3, Supplementary Table S2).

The hemicellulase group is the most diverse compared to cellulases and lignocellulose pretreatment genes, with 20 different GHs activity domains (Figure 2D). Within this group, the most abundant phylum was again Bacteroidetes with 500 ORFs accounting for 91.9%, followed by Firmicutes (27 ORFs), and then by Proteobacteia (7 ORFs) (Figure 2A). Therefore, the ratio of Firmicutes versus Bacteroidetes was again low being 0.05:1. Here, GH2 was estimated to be produced from the most diverse phyla such as Actinobacteria, Bacteroidetes, Firmicutes, Proteobacteia, Thermotogae, followed by GH10, GH15, GH53. Firmicutes and Bacteroidetes were also associated with genes encoding six GH enzymes (GH2, GH28, GH36, GH6, GH51, GH53). However, many other GHs including bi-functional enzymes containing also a CE domain or a CBM domain, and were predicted to be produced solely associated with Bacteroidetes, while GH113, GH26-CBM35 were only linked only to Firmicutes (Figure 2D). At the species level, Prevotella bergensis, P. bryantii, P. copri, Paraprevotella clara, B. graminisolvens, P. albensis, P. paludivivens, B. cellulosilyticus. P. buccae, B. eggerthii, P. ruminicola were the most abundances for producing hemicellulases.

Potential degraders of the cellulosic biomass usually have both beta-glucosidase and endoglucanase [10], while the potential opportunists possess only beta-glucosidase to hydrolyze cellobiose or short polysaccharides to glucose for their own consumption. In accordance with this, the most potential opportunist bacteria in Vietnamese native goats’ rumen were proposed to be Bacteroides fluxus, Lachnospiraceae bacterium JC7, Bifidobacterium breve, Bacteroides caccae, Prevotella nceiensis, Hungatella hathewayi, Pseudomos syringae and Treponema sp. JC4 (Supplementary Table S3) that produce beta-glucosidase only.

Many species were classified as potential lignocellulose degraders in this study, of which 56 species harboured at least 2 different catalytic domains. For cellulose degradation, R. bicirculans and R. flavefaciens were associated with six different catalytic enzymes. R. bicirculans harbored GH3-FN3, GH44, GH5, GH5-CBM2, GH5-FN3, GH74, while R. flavefaciens was predicted to contain CBM63, GH16-CBM4, GH5, GH9, GH94, GH9-CBM3 (Supplementary Table S3).

For effective hemicellulose decomposition, B. graminisolvens possesses nine different catalytic domains, P. copri and P. paludivivens had eight different catalytic domains (Supplementary Table S3). In the bacterial group for lignocellulose pretreatment, Prevotella sp. MSX73 was predicted to the most potential degrader (Supplementary Table S3) to harbour six different catalytic domains. It is interesting that so many species harboured genes for so many functions of lignocelluose degradation. For instance, Prevottella albensis was predicted to bear both cellulases and hemicellulases when this bacteria was predicted to possess both cellulase domains GH3-FN3, GH9-Ig and hemicellulase domains GH10, GH2, GH30, GH43-CBM13, GH43-CBM6, GH97. Seven species contained both GHs for hemicellulases and CEs and PLs for lignocelluose pretreatment. More specifically, B. cellulosilyticus and P. bryantii had 10 genes for different catalytic domains (Supplementary Table S3).

Remarkably, seven species (Draconibacterium orientale, P. buccae, P. copri, P. oryzae, P. paludivivens, R. bicirculans, R. flavefaciens) contained all domains conducting all steps for lignocellulose degradation: pretreatment, cellulose and hemicellulose saccharification. Most domains (14 domains) were identified in one species P. paludivivens, while four other species harboured 13 catalytic domains (Supplementary Table S3).

To get an insight of the bacterial community and their role in lignocellulose digestion in Vietnamese native goats’ rumen, the enriched bacterial metagenomic DNA was sequenced by Illumina platform. The result from NCBI taxonomic classification showed significant amount of unclassified genes that accounted for 50.25%. This suggests that the major part of bacteria in the goats’ rumen was novel or the threshold set for classification was too conservative.

Among 39,579 ORFs affiliated to NCBI taxonomic classification, we found Bacteroidetes and Firmicutes were the most abundant phyla, accounted for 63.5% and 22.6% respectively (Figure 1). In general, the phylogenetic structure of our sample rumen appeared to be similar to the structures described in already published studies of animal rumens. Also, the ratio between Firmicutes and Bacteroidetes was comparable. In our study, ratio Firmicutes/Bacteroidetes was 0.36:1, that is the same ratio as in the rumen of Svalbard reindeer [7], sheep [20], Jinnan cattle [13] and relevant to the ratio 0.3:1 found in goat rumen at 110 days old after adaptation to plant materials [15] (Table 2). However this ratio is strongly dependent on diet throughout the life cycle of goats, and the low ratio between Firmicutes and Bacteroidetes only establishes later in life when the diet consists of mainly plant materials. This was elegantly shown in a previous study [15], they analyzed changes of bacterial community in goats’ rumen after weaning at 60-day old and getting adaptation with plant materials and showed that the ratio Firmicutes/ Bacteroidetes changed from 2.1:1 at 80 days old to 1.7:1 at 100 days old and to 0.3:1 at 110 days old [15]. This demonstrates that for adaptation to plant diet the bacterial structure in rumen typifies a greater abundance of Bacteroidetes if compared to Firmicutes. So, the ratio of Firmicutes/Bacteroidetes changes from higher to lower level. Thus in this study, the low ratio of Firmicutes/Bacteroidetes (0.36:1) indicates that bacterial structure in Vietnamese goats’ rumen adapted well with the diverse of lignocellulose materials in high mountain pastures and dry crop residues. That also explains the diverse lignocellulolytic genes observed in our previous publication [18] compared with Korean goats’ rumen [4]. And this diversity is apparently associated with the increase in Bacteroidetes.

In contrast, the ratio of Firmicutes/Bacteroidetes in our sample is quite different from the observation in rumen of Indian goats (1:1) adapted to Pasture grass, green bajra or rumen of hay-fed cow (0.8–1.7:1) [20], giraffe (1.6:1) [21], as well as rumen from jak (1.7:1) [13] (Table 2). This confirms that bacterial structure in rumen depends strongly on diet.

While, Firmicutes, Bacteroidetes are the two most abundant groups, Proteobacteria and Synergistetes were ranked third and fourth in abundance in our goat rumen. Synergistetes was observed to be abundant in goat rumen in a previous study but not abundant in other rumens [15]. Our study confirms this abundance, and supports the hypothesis that Synergistetes may be unique in goat rumen, potentially reflecting the influence of the host genotype on the bacterial structure.

Regarding the taxonomy at genera level, the two most abundant genera in our sample were Prevotella (35.3%) then Bacteroides (16%) belonging both to Bacteroidales. Prevotella was also the most abundant genus in cattle rumens [22,23] and other goat rumen studies [15,22]. However in that particular case, the Prevotella was reduced during adaptation to other plant matters [15] and went down to lower than 40%.

In previous study, the low ratio Firmicutes/Bacteroidetes was defined to be a parameter for efficient design of a technical system to degrade lignocelluloses on an industrial scale [17], as well as for assessing the optimality of natural cellulolytic systems by herbivorous animals. However, very limited studies showed evidence for this. Here, we elucidated the association between the ratio of Firmicutes versus Bacteroidetes in a functional context to provide scientific evidence for the use of Firmicutes/Bacteroidetes as an indicator of lignocellulose breakdown capacity.

This study indeed confirms that the ratio of Firmicutes/Bacteroidetes was very low (0.11:1) in the bacterial community responsible for lignocellulose breakdown. In other words, the number of lignocellulases associated with Bacteroidetes is ~10 times when compared to Firmicutes or Bacteroidetes plays an important role in lignocellulose degradation, especially in the Vietnamese rumen. In contrast to this, biogas fermenters [17] show a ratio of 2.6:1, while hay-fed cow rumens [20] show a ratio of ~1.4:1. Overall, the ratios Firmicutes/Bacteroidetes in these studies were higher (Table 2).

The lowest ratios of Firmicutes/Bacteroidetes were observed in communities predicted to carry ORFs for hemicellulases and CEs, PLs (0.05:1 and 0.07:1 respectively) (Table 2). Again, the ratios differ between earlier studies done on other biomass degrading environments (biogas fermenter [17] and hay-fed cow rumens [20]) (Table 2). This implicates that the number of CEs and PLs for lignocellulotic pretreament, hemicellulases, are mainly found among Bacteroidetes (Figure 2B). Additionally, Bacteroidetes also harboured very high diverse GHs for hemicellulases (14 families only produced in Bacteroidetes) accompanied with CBMs (Figure 2D). This is the first time CBMs allocated with catalytic domains were so abundantly associated with Bacteroidetes.

Although Bacteroidetes is the most dominant group associated with cellulases, Firmicutes harboured some very specific cellulases GHs domains, accompanied with many CBMs including CBM2, CBM3, CBM4, and CBM63. That is contrary to the CBMs situated in hemicellulases catalytic domains. CBM63 domain found in this study also contained expansin domain that disrupts hydrogen bonding between cellulose microfibrils and enhance the accessibility of cellulases to cellulose substrates [24]. The CBM63 is usually seen in various bacteria [24] but this is the first time observed in R. flavefaciens. CBM63 was highly expressed in Escherichia coli and showed activity to increase activity of cellulase (CAS No. 9012548, Sigma-Aldrich, Toluca de Lerdo, Mexico) [18].

We also investigated whether typical species were associated with lignocellulase activity in the goat rumen, but the species distribution turned out to deviate from a previous study investigating goat rumen microbial diversity [4].

In the present study we identified B. uniformis, R. bicirculans, E. siraeum to be the most dominant species bearing cellulases. In contrast, Prevotella ruminicola, Butyrivibrio proteoclasticus were the two most abundant species associated with cellulases activity in the Korean goat rumen study [4]. However, the dominant species producing cellulases were not abundant species in overall bacterial structure in Vietnamese goats’ rumen.

In the group producing hemicellulases, we found a correlation in abundant species when compared to the overall microbiome in the Vietnamese goats’ rumen. As such, the most six abundant species producing hemicellulases also were in the list of 13 most abundant bacteria in the goat rumen community. Clostridium is known to inhabit the gastro-intestinal tract of ruminants and produce a range of lignocellulolytic enzymes in a multi-enzyme complex designated as cellulosome. However in this study, we did not observe the presence of Clostridium sp, and confirms the absence of this species in another study of goat rumen [15].

In this study, we found many potential lignocellulolytic degraders harbouring many different catalytic domains (up to 14 different catalytic domains) that function for all lignocellulose pretreatment and hemicellulose, cellulose saccharification. These degraders were affiliated in both Bacteroidetes and Firmicutes phyla. Remarkably, R. flavefaciens belonged to Firmicutes is the only gut bacterium so far shown to produce a cellulosome-type enzyme complex [25] was also found in the present study.

Recent research has uncovered that bacteria belonging to Bacteroidetes contain very versatile polysaccharide utilization loci (PULs), that play an important role in cellulose degradation in rumen, gut and faecal samples of herbivores [26]. A metagenomic data investigating the bacteria producing lignocellulolytic enzymes in cow rumen showed several PULs that contain up to 10 domains at one locus. All of these PULs were affiliated with Bacteroidetes [27]. Thus we hypothesize that the species harbouring many catalytic domains in this study may produce many enzymes in the same way as described for PULs previously. However, in this study we also emphasize the significance of lignocellulolytic degrading catalytic domains in Firmicutes as well. This novel insight will be further elucidated in future studies.