Diversity of microbes colonizing forages of varying lignocellulose properties in the sheep rumen

In situ rumen incubation and sample collection of the forages

All experimental procedures relevant to animals were approved by the Ethics Committee for Animal Experiments of the Animal Science Research Institute of Iran (Approval code: IR124-05-05-007-95004-95006). Rumen cannulation was performed according to the guidelines and regulations of American College of Veterinary Surgeons (ACVS) using a two-stage rumen cannulation technique as described previously (Martineau et al., 2015).

Six common lignocellulosic forages including camelthorn (Alhagi persarum, AP; both stem and leaves), common reed (Phragmites australis, CR; both stem and leaves), date palm (Phoenix dactylifera, DP; leaves), kochia (Kochia scoparia, KS; both stem and leaves), rice straw (Oryza sativa; cultivar Hashemi, RS; both stem and leaves), and salicornia (Salicornia persica, SC; both stem and leaves) were selected for in sacco rumen incubation and degradation analysis. Three adult male fistulated Shal (Iranian native) sheep of body weight between 62–68 kg were used for the study. The animals were housed in a barn and fed a diet of 70% wheat straw and 30% barley-based ad libitum twice daily at 08:00 and 16:30. The animals had free access to drinking water. Feeds were introduced into the rumen via the fistula; the material was first baked at 55 °C for 48 h in an air-circulating oven and was then ground finely enough to pass through a 2 mm sieve. An aliquot of the ground material (5 ±  0.05 g) was packaged within a heat-sealed nylon bag (5 × 10 cm; 50 µm pore size). The six feed samples were simultaneously incubated in duplicate per sheep and replicated in the three animals, which were placed within the rumen shortly after a morning meal; two bags per each feed were removed after 24 h, and a further two after each of 48 h, 72 h, and 96 h.

After removing the bags, in order to analyze fiber breakdown, one of the pair of bags was placed immediately in cold distilled water, severely rinsed for about 5 min, and then transferred to the laboratory where they were washed in a fully automatic washing machine with cold water for 15 min. The non-incubated nylon bags (0 h) were washed in the same manner using washing machine. The other bag rinsed in running distilled water, while rubbing gently between thumb and fingers until the water became clear, and finally squeezed by hands with sterile gloves to remove excess water. The bags were then snap-frozen in liquid nitrogen and stored at −80 °C until required for subsequent DNA extraction and the other one was processed for detergent fiber analysis.

Microbial dissociation and DNA purification

Microbial cells adhering to the digested material were recovered by suspending the digested feed particulates with a 1:2 (w/v) ratio in dissociation buffer containing 0.1% (v/v) Tween 80, 1% (v/v) methanol and 1% (v/v) tertiary butanol (pH 2), following Pope et al. (2011). The samples were vigorously vortexed every 1–3 min and this step was repeated at least five times. After that, the plant material was sedimented by imposing three rounds of low speed centrifugation, with the supernatant being transferred to a new tube after each spin. Finally, the microbial content was recovered by a rapid centrifugation (12,000 g for 5 min). Microbial DNA was extracted using the QIAamp^® DNA Stool Mini Kit (Qiagen, Hilden, Germany). The concentration of the DNA preparations was determined using a NanoDrop spectrophotometer (Thermo Scientific, Wilmington, DE, USA) and its integrity was checked by inspection of an electrophoretically separated (1% agarose gel) aliquot.

In sacco feed digestion analysis

The neutral detergent fiber (NDF), acid detergent fiber (ADF), and acid detergent lignin (ADL) contents of the digested feed samples were determined following Van Soest, Robertson & Lewis (1991). The measured loss in dry matter (DM), NDF, ADF, and ADL were used as indices of the forage’s degradation. The cellulose and hemicellulose contents of forages were calculated, by subtracting ADL from ADF and ADF from NDF, respectively. In sacco disappearances of the measured contents over sampling periods were calculated as the difference between the amount of those components in the substrates before and after incubation.

16S rRNA library preparation and sequencing

A two stage PCR amplification was employed to both enrich for 16S rRNA sequence and add Illumina sequencing adaptors and dual index bar-codes to the sequencing template. Initially, the universal primer pair Bact-0341F (5′-CCTACGGGNGGCWGCAG) / Bact-0785R (5′-GACTACHVGGGTATCTAATCC) was used to amplify a 460 bp fragment corresponding to V3 and V4 regions of the bacterial 16S rRNA sequence (Klindworth et al., 2013). Each 50 µL PCR contained 50 ng DNA (2 µL), 25 µL 2x PCR Master Mix (Thermo Scientific), 21 µL nuclease-free water, 1 µL forward primer (50 pM) and 1 µL reverse primer (50 pM). The PCR condition consisted of an initial denaturation at 94 ° C for 4 min followed by 25 cycles of 94 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s, and a final extension at 72 °C for 5 min. PCR products were purified and 2 µL of each reaction was used as a template for the second round of PCR during which the Illumina adaptors and barcode sequences were incorporated to the 5′-end of the amplified products. The second PCR was also performed in triplicate under the same running condition except the number of cycles were limited to 15. The amplicons were recovered from a 2% agarose gel and purified using the QIAquick Gel Extraction Kit (Qiagen), quantified fluorometrically, pooled in equimolar quantities, and each sheep individually paired-end sequenced (PE300) using the Illumina MiSeq System at Macrogen Inc. (Seoul, South Korea).

Sequence analysis

Trimmomatic v0.36 software was used to edit the raw paired-end sequence data, removing adaptor sequences and bases associated with a Phred score <3 (Bolger, Lohse & Usadel, 2014). A four base sliding window was used to scan the reads in order to delete when the average quality per base fell below 15. Paired-end reads were merged using FLASH software (Magoč & Salzberg, 2011) and reads which could not be assembled were discarded. Raw fastq files were de-multiplexed using the split_libraries_fastq.py script of the QIIME pipeline (version 1.9.1) (Caporaso et al., 2010b). Chimeric sequences were removed using the VSEARCH with the most recent version of RDP database (release 11) as a template (Rognes et al., 2016). Sequences were then clustered into operational taxonomic units (OTUs) using UCLUST software (Edgar, 2010), applying a 97% similarity cut-off. The most abundant sequences belonging to each OTU cluster were considered as representative and were aligned against the Greengenes core set (DeSantis et al., 2006) using the PyNAST aligner (Caporaso et al., 2010a). The Ribosomal Database Project naïve Bayesian classifier was used to obtain taxonomic information, applying a minimum support threshold of 80% (DeSantis et al., 2006; Wang et al., 2007).

The OTU table was filtered for low abundant OTUs using filter_otus_from_otu_table.py script with –min_count_fraction option was set to 0.00001 (discarding OTUs represented by <0.001% of the sequences). The OTU table was rarefied to the sequencing depth 3460 reads which corresponded to the number of reads in the sample with minimum number of reads as the sub-sampling depth, then rarefaction curves of the clustered OTUs were drawn using the QIIME pipeline (Caporaso et al., 2010b). Rarefaction curves through alpha diversity indexes were also plotted using high sequencing depth (50000). Alpha diversity indices including Shannon, Simpson, Good’s_coverage, and Chao1 were calculated using core_diversity_analyses.py script in the QIIME pipeline. Beta diversity indices including weighted and unweighted Unifrac phylogenetic distance matrices were constructed with rarefied OTU table as input and visualized through principal coordinate analysis (PCoA) plots made by R/ggplot2 (Wickham, 2016). In the present study, due to uneven sequencing depths across samples, we needed to take low rarefaction value to keep more samples and avoid losing any treatment. The extent to which that could affect the diversity results was assessed through rerunning the alpha and beta diversity analysis using higher sequencing depth (50000).

Statistical analysis

Statistically significant differences in physicochemical data including DM, NDF, ADF, ADL, cellulose, and hemicellulose contents were analyzed by one-way ANOVA and the general linear model (GLM) procedure of the SAS software v9.3 (SAS Institute Inc., Cary, NC, USA). The main and interaction effects of forage and incubation time were measured based on chemical compositions of six experimental forages. The orthogonal contrast was performed to assess the distributional trend of ruminal microbiome (genus level) adhered to the six lignocellulosic forages during their rumen incubation. In situ degradability kinetics for the chemical composition of the experimental forages was evaluated by the exponential model (Orskov & McDonald, 1979).

Differences in taxa abundances between forages and sampling intervals were estimated using analysis of composition of microbes (ANCOM) based on relative abundances of OTUs summarized at various taxonomic levels (Mandal et al., 2015). Means were compared by Duncan post-hoc test in PAST v4.1 given Bonferroni P-value cutoff <0.05 (Hammer, Harper & Ryan, 2001). Error correction was made based on the number of groupwise comparisons performed at each taxonomic level. Permutational multivariate analysis of variance (PERMANOVA) was performed using the adonis function of R-package vegan v2.5-5, to test for significant differences between community compositions of forage attached microbial communities. In addition, permutation multivariate analysis of group dispersions (PERMDISP), using the betadisper function of vegan, was used to test for the homogeneity of dispersions (variances). The Spearman correlation analysis was performed using the corr.test function of the psych package v1.8.12 and visualized using R/ggcorrplot package. The p-values were corrected using Bonferroni method based on the total number of correlations calculated for each variable separately. For all tests, p-values less than 0.05 were considered statistically significant.

Sequence deposition

All raw amplicon sequences were submitted to NCBI short read archive under BioProjectID: PRJNA642227.

Physicochemical characteristics of the incubated feeds

The selected forages were analyzed with respect to their contents for NDF, ADF, ADL, cellulose, and hemicellulose before and after rumen incubation (Table 1, Table S1). They showed a significantly high diversity (p < 0.05) concerning their NDF contents within which CR showed the highest (4.39 ± 0.02) but AP and RS showed the lowest (3.93 ± 0.04 and 3.89 ± 0.02, respectively) mean NDF contents (Table 1). The ADF contents of forages varied with AP (3.07 ± 0.02) had the highest and RS (2.33 ± 0.02) but SC (2.37 ± 0.03) had the lowest amounts. With respect to the ADL contents, the differences between forages were also statistically significant (p < 0.05). AP (1.25 ± 0.02) showed the highest but RS (0.20 ± 0.006) the lowest amounts of initial ADL while the difference between CR and KS was not significant in this respect. The contents of forages for cellulose and hemicellulose also varied, in which CR and RS showed the highest cellulose while SC and CR showed the highest hemicellulose contents, whereas AP displayed the lowest amounts of both cellulose and hemicellulose among other forages.

In sacco rumen biomass degradation

The six incubated forages were monitored for changes in the contents of DM, NDF, ADF, ADL, cellulose, and hemicellulose during their rumen incubation (Table S2 and Fig. S1). DM degradation was the fastest in KS but slowest in RS (33.8% vs 20.17%, p < 0.05) during the first 24 h of rumen incubation. The amount of DM degradation after 48 h of rumen incubation was significantly less than the first 24 h; RS (26.3%) had the lowest degradation but DP (41.3%) had the highest degradation in the second 24 h. With respect to NDF fermentation, RS showed the lowest (16.73%) while KS the highest amount of NDF loss (31.00%) after 24 h of rumen incubation. 48 h after rumen incubation, the highest rate of NDF degradation occurred in DP while the lowest amount was observed in RS. This pattern continued until the end of 96 h of incubation in the rumen. In the first 24 h, besides RS, in which the lowest rate (9.85%) of ADF degradation occurred, all other forages showed a constant rate of degradation (average of 21.2%) and their difference was not significant. 48 h later, a significant difference was observed between DP and RS. After that, up to 96 h, the degradation rate was slight and the differences between the forages were not significant. Regarding ADL degradation, forages were divided into two groups containing AP, DP and KS (average of 25.6%) compared with SC, CR, and RS (average of 7.11%), which had the highest and lowest amounts of average for ADL degradation, respectively. These homogenous subsets showed steadily pattern of degradation over ruminal incubation. In general, by extending the rumen incubation beyond 24 h, the degradation of NDF, ADF, and ADL was significantly declined (Table S2 and Fig. S1). It was apparent that the initial differences in fiber-related physicochemical properties significantly affected the rumen digestion of the six forages. Moreover, forage and incubation time effects significantly affected all fiber-related parameters, however, an interaction effect between forage and incubation time influenced the rumen degradability of ADL and hemicellulose (Table S3).

Degradation kinetics of chemical compositions

In situ DM, NDF, and ADF degradation kinetic parameters of forages, calculated by an exponential model, are presented in Table S4. There was large variability in soluble fraction (a), slow degraded fraction (b), and degradation rate (C). DM degradability of the forages revealed a low range of the effective degradability (ED) over all experimental forages, on an average of 29.7%, which was the highest in KS (37.50%) following DP (32.70%) but the lowest amount of 16% for RS (p < 0.05). The potential degradability also showed a similar tendency with effective degradability, which was over 60% for KS but under 32% for RS. This pattern was in agreement with DM disappearance illustrated in Fig. S1. In general, there was a high range of soluble fraction for NDF (17.51 ± 8.3) and ADF (14.7 ± 7.1) degradation data. Effective degradation (ED2) of NDF ranged from 19.90% to 32.30%, and that of ADF ranged from 15.30% to 28.10%, however, there was no true difference between the forages for ADF and NDF degradability.

16S rRNA sequencing data analysis

MiSeq-based sequencing of 16S rRNA gene amplicons resulted in 5,306,118 sequence pairs (average 73,696 pairs of sequences/sample). After quality trimming, sequences were merged using FLASH which resulted in the assembly of 3,967,162 sequences (ranged from 2,729 to 1,34,454 on average of 55,099 sequences/sample) with a mean length 456 bp. A further quality filtering during demultiplexing step discarded 715,635 low quality and short sequences. A total of 839,990 sequences were identified as chimera and discarded. The remaining 2,295,869 non-chimeric sequences were clustered at the 97% sequence identity. After filtering for low abundant features 2,861 clean OTUs representing 2,163,318 sequences were remained for downstream diversity analyses. To test whether our sequencing effort provided sufficient coverage to recover the species diversity of fiber-attached microbiota a rarefaction analysis was performed. The rarefaction curves for most samples approached to plateau but did not reach to the saturation, indicating that our sequencing effort was not enough to fully explore the diversity of forage attached microbiota (Figs. S2 and S3).

Microbial diversity analysis

Alpha diversity analysis showed no statistically significant difference among the forages but significant variations were observed among sampling intervals (Fig. 1). Chao1 index differed between 72 and 96 h samples, while both Shannon and Simpson indices were significant between 24 and 96 h time points (p < 0.05). Nevertheless, diversity indices were not affected by the incubation length within different forages. All samples showed a high Good’s coverage (>0.92, Fig. 1) at the sampled time points, indicating that our sequencing effort allowed to detect >90% of the microbial diversity associated to the forages.

Beta diversity analysis through weighted and unweighted Unifrac divergence indices also revealed no significant difference between forages regardless of their incubation time (Fig. 2). However, when the microbial communities of forage-attached microbes were compared across the incubation length, significant differences were observed in weighted Unifrac (PERMANOVA P < 0.001) and unweighted Unifrac (p < 0.001) indices (Fig. 2). Testing for homogeneity of group dispersions also revealed no significant difference in dispersions across forages and sampling times (Fig. 2B, PERMDISP p > 0.6 and p > 0.8, respectively).

In the present study, due to uneven sequencing depths across samples, we needed to take low rarefaction value to keep more samples and avoid losing any treatment. The extent to which that could affect the diversity results was assessed through rerunning the diversity analysis using no-rarefied samples. The outcomes appeared to be identical, so that, significant differences were just observed between the incubation times (Figs. S4 and S5)

Community composition of forage-attached microbes

The microbial community associated to forages were assigned to 16 bacterial phyla and one archaeal phylum (Fig. 3). Bacteroidetes (CR and DP) and Firmicutes (CR and AP) showed significant differential abundances between forages, over 96 h ruminal incubation (p < 0.05). Bacteroidetes comprised the most abundant phylum accounting for 42.1% of all sequences. Firmicutes were the second most abundant phylum accounting for 30.3% of sequences followed by Proteobacteria (18.5%), Fibrobacteres (3.4%), Spirochaetes (3%), Actinobacteria (0.84%), Synergistetes (0.18%), Verrucomicrobia (0.17%), Tenericutes (0.15%), and TM7 (0.14%). Comparing the prevalence of major bacterial phyla in forage-attached microbes regardless of the incubation time revealed some differences between forages. Bacteroidetes were particularly overrepresented in CR–(48.7%) but were underrepresented in DP-associated microbial community (34.9%). Firmicutes were less represented in CR-attached microbiota (22.6%) while they were equally distributed in other forages (on average 31.9%). Proteobacteria were less represented in the community of KS but overrepresented in that of DP-attached microbes (14.9% and 22.6%, respectively, vs average of 18.3% in other feeds). An interesting finding was the higher prevalence of fiber-degrading bacteria of Fibrobacteres phylum in the microbial community associated with CR (5.4%) and RS (5.5%), two forages with the highest initial cellulose contents, while they accounted for less than 3.3% of reads in other forages.

At the family level, the fiber-attached sheep rumen microbiota was affiliated to 61 bacterial and 2 archaeal families (Fig. S6). Prevotellaceae was the most predominant bacterial family representing about 27.3% of the attached bacterial community followed by Succinivibrionaceae (18%), Lachnospiraceae (13.8%), unclassified Bacteroidales (8.9%), Veillonellaceae (6.2%), Ruminococcaceae (4.1%), Unclassified Clostridiales (4%), and Fibrobacteraceae (3.4%). Bacteroidaceae and Lachnospiraceae were the two families differed between the forages (ANCOM p < 0.05; Fig. S7). Bacteroidaceae occurred abundantly in CR and RS-associated microbiota. Members of Lachnospiraceae were less represented in CR- and RS- (p < 0.001), while more highly represented in AP-associated microbiota (p < 0.02).

At the genus level, forage-attached microbes were distributed in 106 bacterial genera (Fig. 4), of which, six including Butyrivibrio, Lachnospira, Pseudobutyrivibrio, BF311, unclassified Lachnospiraceae and Shuttleworthia displayed differential abundances among the forages (ANCOM p < 0.05, Fig. 4). All members of Lachnospiraceae family showed least significant abundance in CR compared to the other forages (p < 0.05). However, unclassified BF311 showed high abundance in CR and RS and they also differed significantly between CR and AP (p < 0.05).

The Random Forest machine learning analysis enabled to identify the most predictive taxa discriminating forages from each other during the course of rumen incubation (Fig. 5). During the first 24 h of rumen incubation, Treponema was the most predictive taxa for which relative abundances were differed between forages. At this time point, the relative abundances of Oscillospira and Lachnospira were also helpful for discriminating microbial communities associated to the forages but at a lower degree of accuracy. At the 48 h time point, Treponema was also the most predictive taxa followed by Pseudobutyrivibrio, Anaerovibrio, and Lachnospira which their relative abundances could be predictive of microbial communities associated to different forages. At the 72 h time point, Butyrivibrio was the only taxa that could be used to separate forages, they showed a significantly low abundance in the forages with the highest initial cellulose contents (e.g., CR and RS). By extending the incubation time to 96 h the prevalence of unassigned sequences was the most predictive one followed by Butyrivibrio which was also informative at the previous time point (e.g., 48 h).

Changes in forage-associated microbes during their rumen incubation

Variations in abundance of genera were also surveyed over the incubation length. The abundance of 14 genera changed with incubation (Fig. S8). Among which, there were several highly abundant rumen bacterial genera which are known to play key roles in plant lignocellulose degradation, including Fibrobacter, Ruminococcus, unclassified Lachnospiraceae, unclassified S24-7, unclassified Ruminococcaceae, unclassified Veillonellaceae, Lachnospira, Coprococcus, Dialister, Roseburia, Selenomonas, Desulfovibrio, Anaerovibrio, and RFN20. The abundance of unclassified Lachnospiraceae gradually increased with incubation length particularly in AP (the forage with the highest initial ADL concentration). The unclassified Ruminicoccaceae were increased with incubation and peaked at 72 h. Members of unclassified Veillonellaceae were more highly represented at 24 h in AP (2.1%) and RS (1.7%) and at 96 h in CR and RS the two forages with the highest initial cellulose concentrations.

Based on orthogonal contrast results (Table S5), Fibrobacters showed a significant linear increase in abundance in RS, KS, DP, and SC during incubation (Fig. 6 and Table S5). The prevalence of Ruminococcus and Coprococcus cubically changed over nRS rumen incubation, however, Coprococcus showed the same trend during DP incubation in rumen (p < 0.05). Abundance of Selenomonas showed a linearly (AP) and quadratically (CR, DP, and SC) decreasing trends during ruminal incubation (p < 0.05). Moreover, Butyrivibrio and Roseburia indicated linearly decreasing trend in SC (p < 0.01) and RS (p < 0.05) during ruminal incubation, however, Roseburia revealed quadratic variation in abundance during SC incubation (p < 0.05) in the rumen. Un_Veillonellaceae and Treponema abundances changed quadratically during RS and SC incubation (p < 0.01). Un_S24-7 genus linearly increased over SC incubation in rumen, however, its abundance quadratically changed in RS and KS during ruminal incubation.

The connection between fiber-attached microbiota and lignocellulose degradation

A spearman correlation analysis was performed to infer the association between digestibility indices and the community composition of fiber-attached microbiota (Fig. 7). This analysis revealed that cellulose content of forages negatively correlated with the abundance of unclassified S24-7 (r =  − 0.37, p = 0.001), unclassified Paraprevotellaceae (r =  − 0.47, p = 0.001), unclassified Christensenellaceae (r =  − 0.32, p = 0.002), and unclassified Lachnospiraceae (r =  − 0.53, p = 0.002) but positively correlated with the prevalence of Bifidobacterium (r = 0.37, p = 0.001), Prevotella (r = 0.38, p = 0.001), unclassified Veillonellaceae (r = 0.38, p = 0.001), and Selenomonas (r = 0.36, p = 0.002). The hemicellulose content of the forages was also inversely correlated with the population of unclassified Lachnospiraceae (r =  − 0.34, p = 0.001), Lachnospira (r =  − 0.40, p = 0.001), and Pseudobutyrivibrio (r =  − 0.33, p = 0.001). Members of unclassified Veillonellaceae (r = 0.37, p = 0.002; r = 0.46, p = 0.001) and Selenomonas (r = 0.39, p = 0.001; r = 0.42, p = 0.001) positively correlated with NDF and ADF concentrations but Treponema did inversely (r =  − 0.34, p = 0.001; r = 0.46, p = 0.001). Moreover, unclassified Lachnospiraceae (r =  − 0.45, p = 0.001) and Pseudobutyrivibrio (r =  − 0.39, p = 0.001) correlated negatively with NDF content. The ADL content of forages, however, indicated no statistically significant association.