Diversity of rumen microbiota using metagenome sequencing and methane yield in Indian sheep fed on straw and concentrate diet.

An in vivo study aiming to investigate the rumen methanogens community structure was conducted in Mandya sheep fed on straw and concentrate diet. The ruminal fluid samples were collected and processed for unravelling the rumen microbiota and methanogens diversity. Further, the daily enteric methane emission and methane yield was also quantified using the SF6 tracer technique. Results indicated that the Bacteroidetes (∼57%) and Firmicutes (25%) were two prominent affiliates of the bacterial community. Archaea represented about 2.5% of the ruminal microbiota. Methanobacteriales affiliated methanogens were the most prevalent in sheep rumen. The study inveterate that the ruminal archaea community in sheep is composed of 9 genera and 18 species. Methanobrevibacter represented the largest genus of the archaeome, while methylotrophs genera constituted only 13% of the community. Methanobrevibacter gottschalkii was the prominent methanogen, and Methaobrevibacter ruminantium distributed at a lower frequency (∼2.5%). Among Methanomassiliicoccales, Group 12 sp. ISO4-H5 constituted the most considerable fraction (∼11%). KEGG reference pathway for methane metabolism indicated the formation of methane through hydrogenotrophic and methylotrophic pathways, whereas the acetoclastic pathway was not functional in sheep. The enteric methane emission and methane yield was 19.7 g/d and 20.8 g/kg DMI, respectively. Various species of Methanobrevibacter were differently correlated, and the distribution of hydrogenotrophic methanogens mainly explained the variability in methane yield between the individual sheep. It can be inferred from the study that the hydrogenotrophic methanogens dominate the rumen archaeal community in sheep and methylotrophic/aceticlastic methanogens represent a minor fraction of the community. Further studies are warranted for establishing the metabolic association between the prevalent hydrogenotrophs and methylotrophs to identify the key reaction for reducing methane emission.

Introduction

The feed fermentation in the rumen is accomplished by millions of anaerobic microbes, which also assist in the removal of waste by-products (waste) such as carbon dioxide and methane. Bacteria are the most diversified and abundant microbes in the rumen and play a significant role in the degradation of fibre components of the diet. The anaerobic fermentation yield the production of microbial protein, volatile fatty acids, fermentation gases, and organic acids. The unique feature of the ruminal microbiota is their existence in the syntrophic fashion where the end product of one microbial consortium is being used by other. The by-products of anaerobic fermentation, such as CO2 & H2, are utilized by the ruminal archaea to form methane. Ruminal methanogenesis is affected by various factors such as feed ingredients, diet composition, feeding frequency, methanogens community structure, host species, and geographical region. Despite of the limited representation (3–5%) in the rumen microbiota (Henderson et al., 2015, Yanagita et al., 2000), archaea play a significant role in maintaining the H2 pressure within desirable limits (Gagen et al., 2015). Diverse substrate requirements and inter-linking with other ruminal microbes make the rumen methanogens community very complex. So far 155 species of methanogens isolated from various ecosystems found associated to the six orders and twenty-nine genera (Holmes and Smith, 2016). To date, 40 different strains of the methanogens have been isolated from all environments (Seedorf et al., 2015); wherein, only seven species have been isolated from the rumen ecosystem (Janssen and Kirs, 2008). The ruminal methanogens isolated till date represent only 1/10th of the archaeal community (Kim, 2012), and a larger fraction (∼90%) is yet to be identified. Cultivation based studies are usually unsuccessful in providing holistic insights to the rumen archaeal community. The methanogens community in sheep has been previously studied in Australia (Wright et al., 2006, Wright et al., 2004), Japan (Yanagita et al., 2000), China (Huang et al., 2016), and Latin America (Wright et al., 2008) using conventional techniques. With the development of cultivation-independent techniques, recent bioinformatics tools and precise taxonomic framework (Seedorf et al., 2015, Seedorf et al., 2014) had increased our understanding and identification of methanogens in the rumen. A successful methane mitigation strategy from livestock is not feasible until the whole ruminal methanogens community is carefully analyzed. Geographical locations may also have some influence on the methanogens community structure (Hook et al., 2010). The present study reports the community structure of ruminal archaea in native Mandya sheep. Our results will contribute to understanding the archaeal community in tropical sheep and developing effective strategies for manipulating rumen methanogens and methanogenesis. The findings will also help in comparing the rumen methanogens community between tropical and temperate world.

Materials and methods

Animals and feeding

The animal experiment was conducted at the Experimental Livestock Unit of the ICAR-National Institute of Animal Nutrition and Physiology. The animal handling was performed as per the regulations of the Institute Animal Ethics Committee (approval F. No. NIANP/IAEC/1/2019). Nine adult Mandya sheep of 30.61 ± 1.63 kg body weight (mean ± SD) were fed on a basal diet comprising finger millet straw and concentrate in a 55:45 ratio. The chemical composition of the diet and ingredients is given in Table 1. The sheep were dewormed in the beginning with fenbendazole (5 mg kg−1) and housed in tail to tail orientation in a ventilated shed. The animals were fed total mixed ration in the morning at 09.00 h; whereas the drinking water was freely accessible throughout day. The dried and ground samples (1 mm) of total mixed ration were analyzed for crude protein (N × 6.25) and ash content as per AOAC (AOAC, 2012). The acid detergent fibre (ADF) and neutral detergent fibre (NDF) in the diet were estimated as per Van Soest et al. (1991). The animals were fed on the basal diet for 45 days before collecting the rumen fluid samples to study the archaeal community composition.

Table 1

Chemical composition of feed ingredients and diet fed to sheep.

Particulars	OM	CP	NDF	ADF	TA
Finger millet straw	903	39.6	704	517	97
Concentrate	932	208	426	128	68
Total mixed ration	921	120	593	366	79
Concentrate ingredients
Maize	980	88.1	506	84	20
Soybean meal	930	470	251	126	70
Groundnut cake	939	444	201	135	61
Wheat bran	945	150	584	140	55

DM- dry matter; OM- organic matter; CP- crude protein; EE- ether extract; NDF- neutral detergent fibre; ADF- acid detergent fibre; TA- total ash; concentrate mixture was formulated with maize grain (32 parts), groundnut cake (12 parts), soybean meal (13 parts), wheat bran (40 parts), mineral mixture (2 parts) and salt (1 part).

Rumen fluid collection and DNA isolation

After the 45th day of feeding, rumen fluid samples were collected from the animal (30 ml) using a stomach tube as described previously (Malik et al., 2015b). The ruminal fluid samples were filtered through muslin cheese cloth and transported on ice to the laboratory. Each ruminal fluid sample was divided into two equal subsets. The first subset was processed to estimate volatile fatty acid (VFA), whereas another subset of ruminal fluid was used for the DNA isolation. In brief, the feed particles dissolved in the ruminal fluid were removed by short spinning at 1500 rpm for five minutes. The supernatant was collected into another tube and 2 ml of supernatant was centrifuged at 13000 rpm, 4 °C for 10 min. After centrifugation, the supernatant was discarded and 1 ml lysis buffer was added to the pellet for dissolving. The whole content was then transferred to a sterile bead beating tube containing 0.5 g pre-sterilized zirconia beads (Biospec USA). The genomic DNA was isolated by following the RBB + C method of Yu and Morrison (2004). Thereafter, the metagenomic DNA was checked on 0.8% agarose gel electrophoresis and quantified using Qubit 4.0 (Invitrogen).

Whole rumen metagenome analysis

The genomic DNA was cut into fragments using restriction enzymes, and amplicon libraries were prepared using Covaris M220 (Covaris lnc., Woburn, MA, USA) to generate 300 bp DNA fragments. Metagenome libraries were constructed using TrueSeq DNA PCR-Free Library Prep Kits (Illumina, San Diego, CA, USA) as per the manufacturer’s protocol. Sequencing data generated on HiSeq 2000 (Illumina) in the format of FASTQ files containing paired-end raw reads (2 × 150 bp) were processed online with Meta Genome Rapid Annotation using Subsystem Technology (MG-RAST, Meyer et al., 2008, Wilke et al., 2013). A dereplication step was performed using a simple k-mer approach to remove artificial duplicate reads. The RefSeq database hits were used for generating the phylum and genus level representation. The correlation matrix was computed and visualized in R (“R Core Team, 2021) using Corr package (Wei and Simko, 2021). The Subsystem database hits were used to identify different Subsystem categories and KO hits related to methanogenesis were visualized on KEGG Mapper color tool (Kanehisa et al., 2022).

Archaeal diversity

The metagenomic DNA were processed to prepare the 16S rRNA libraries. Amplicon libraries were prepared using the Nextra XT kit (Illumina Inc., San Diego, United States). Archaea specific primers Arch-344F (Wemheuer et al., 2012) and Arch-806R (Takai and Horikoshi, 2000) were used for preparing initial amplicons. The PCR amplification was performed as follows: initial denaturation at 95 °C for 3 min followed by 25 cycles of denaturation at 95 °C for 30 s, annealing at 55 °C for 30 s, extension at 72 °C for 30 s, and final extension at 72 °C for 5 min. A reaction without template was used as a negative control, whereas the Methanobrevibacter smithii DNA template was used as a positive control during the PCR amplification. After purification of the PCR products using AmPure XP beads Nextra XT index primers (i5 and i7) were added using Nextra XT index kit. All the indexed amplicon libraries were purified with AMPureXP beads (Beckman Coulter Life Sciences, United States) and analyzed individually on a 4,200 Tape Station (Agilent Technologies, United States). The libraries were multiplexed (10–20 pM of each) and sequenced on an Illumina MiSeq platform (Illumina, San Diego, CA, United States) using MiSeq reagent kit v3, and 2 × 300 bp paired-end reads were generated to obtain approximately 0.1 million sequences per library.

Bioinformatics analysis

Amplicon sequence reads generated using MiSeq were processed on CLC Genomics Workbench and CLC Microbial Genomics Module (V20; CLC, Qiagen, Aarhus, DK). Quality filtering was performed using a Phred score at a cutoff of < 20. The resultant filtered reads obtained from the samples were paired, and chimeric sequences were removed with a crossover cost of 3. Closed reference OTU pickup was performed using Rumen and Intestinal Methanogens Database (RIM-DB; (Seedorf et al., 2014)) at a similarity threshold of ≥ 97%. The representative sequences from each OTU cluster were taxonomically annotated up to the species level. The rarefaction curve was plotted using Vegan package 2.5.6 in R 3.5.1 (Oksanen et al., 2012). Rumen methanogens diversity was estimated using Shannon index by phyloseq package 1.26.1 (McMurdie and Holmes, 2013). The OTUs not present in any of the samples were pruned for computing the Shannon index. After normalization of OTU counts per 10,000 reads (Suppl. File 1), the abundance was studied at different taxonomic ranks and plotted using ggplot2 (Wickham, 2011). Core microbiome analysis was performed using the package microbiome V1.4.1 (Lahti and Shetty, 2012) in R with minimum prevalence in 50% of the samples. The microbiome analysis for the rumen archaeal prevalence in the samples was performed at different taxonomic levels.

To explore the co-abundance and co-exclusion among the rumen archaea in sheep at the genus level, network analysis was done using SParse InversE Covariance Estimation for Ecological Association Inference (SPIEC-EASI; (Kurtz et al., 2015)). Archaeal abundance data were processed and transformed into centered log-ratio. Neighborhood selection method was employed for graphical model inference with a minimum lambda ratio of 0.02 and iteration of 99 times. Thereafter, stability selection was performed using Stability Approach to Regularization Selection (StARS) as per the algorithm of Liu et al. (2010). General network attributes including network stability, degree distribution, betweenness centrality, and closeness centrality were retrieved from the SPIEC-EASI model. The undirected graph pertaining to the network connectivity was obtained from SPIEC-EASI and visualized in Gephi V0.9.2 (Bastian et al., 2009) using the Fruchterman-Reingold algorithm layout.

Methane emission

Daily enteric methane emission from the individual sheep was measured after 45 days of feeding using sulfur hexafluoride (SF6) technique (Berndt et al., 2014). The methane measurement and daily dry matter intake from the individual animal was recorded for seven consecutive days. The brass permeation tubes (34 mm long, 8.5 mm dia.), bored with 30 mm deep and 4.8 mm wide blind hole, were fitted with a Swagelok nut. The 7.0 mm diameter hole in the nut provided the permeation window for SF6 release. A Teflon septum (0.24 mm PTFE) and stainless steel frit (3/8″ OD, 2 µm pore size) were used for regulating the release rate from permeation tubes. Tubes were initially kept in liquid nitrogen for three days and thereafter charged with SF6 (750 ± 49.48 mg) while held in liquid nitrogen. The permeation tubes filled with SF6 were placed in an incubator at a temperature similar to the rumen. To calculated the daily release rate of SF6, the tubes were weighed daily at the same time continuously for 65 days. The calibrated tubes were inserted into the rumen seven days prior to the commencement of methane emission. An average release of 3.58 ± 0.13 mg/d (mean ± SE) SF6 was recorded from the permeation tubes. Following the guidelines of Williams et al. (2014), the PVC canisters were assembled with the nylon tube, capillary tube, and male–female Quick connectors. The methane and SF6 concentration in air was quantified using the background sample collected in the PVC canister hung on the ventilated iron mesh fixed. The gas samples in the canister were initially diluted with N2 and the successive sub-samples were drawn in an airtight syringe (Hamilton, 1 ml). After considering the local elevation and atmospheric pressure, the dilution factors were derived as per Lassey et al. (2014).where GS is the calculated concentration of methane (ppm) or SF6 (ppt) at average atmospheric pressure of 90 kPa at an elevation of 920 m. τf (kPa) was the final vacuum in canister after the addition of N2, τs (kPa) was the vacuum in the canister after the sample collection, τe was the vacuum in the evacuated canister before use, GA was the concentration of methane (ppm) or SF6 (ppt) in the sample presented to the GC.

Estimation of methane and SF6 was carried out using a gas chromatograph (GC 2010 plus, Shimadzu, Japan) fitted with flame ionization detector (FID) and electron capture detector (ECD), respectively. The gas chromatographic conditions described previously were used for the analysis of gases in diluted breath samples (Malik et al., 2021). The daily enteric methane emission was calculated using the formula of Moate et al. (2014).

RCH4 is daily CH4 emission (g/d); RSF6 represents the SF6 release rate; MWCH4 is the molecular weight of CH4; MSF6 is the molecular weight of SF6.

Daily enteric methane emission from the individual sheep was also compared uniformly by calculating the methane yield (MY, g/kg DMI). MY was calculated by considering the daily methane emission (g/d), and the dry matter intake (DMI) recorded over the measurement period.

VFA estimation

Total volatile fatty acid (TVFA) concentration in the ruminal fluid was estimated by Markham’s distillation apparatus, whereas individual VFA were estimated using the procedure of Filipek and Dvorak (2009). In brief, metaphosphoric acid (25%) was added to the supernatant of ruminal fluid in the ratio of 1:4 (v/v) and processed for the VFA estimation using a gas chromatograph (Agilent 7890B) by following the conditions previously described (Malik et al., 2021).

Principal component analysis

To reduce the dimensionality of large datasets, the Principal Component Analysis (PCA) was performed in GraphPad Prism version 9 by considering the distribution of individual methanogen species as variables and establishing the influence on methane emission (principal component). Total 18 variables (methanogens distribution) were taken, whereas daily methane emission was considered as the dependent variable.

Results

Rumen microbiome

Whole metagenome analysis revealed that the rumen microbiota in Indian Mandya sheep at the domain level was dominated by the bacteria, followed by archaea and the Eukaryota. The bacterial community at the phylum level (Fig. 1) was dominated by Bacteroidetes (∼57%), followed by Firmicutes (∼25%). The Firmicutes to Bacteroidetes ratio in sheep was 0.44, and both were negatively correlated. The third most abundant bacterial phylum was Proteobacteria, which constituted 6% of the total rumen microbiota, whereas the Fibrobacteres and Actinobacteria represented 1.75 and 2.29%. At the genus level, the top 10 most prominent genera (Suppl. Fig. 1) were Prevotella (35.1%), Bacteroides (14.3%), Clostridium (5.6%), Ruminococcus (3.0%), Eubacterium (2.4%), Parabacteroides (1.8%), Fibrobacter (1.7%), Butyrivibrio (1.7%), Methanobrevibacter (1.7%) and Bacillus (0.9%). Euryarchaeota were most prominent among archaea, constituting 2.49% of the microbiota. The correlation of Euryarchaeota with major bacterial phylum revealed that the Bacteroidetes were negatively correlated (R = −0.92), whereas Firmicutes were positively correlated (R = 0.92). Similarly, the Actinobacteria were also positively correlated (R = 0.91) with the Euryarchaeota (Fig. 2).

Fig. 1

Distribution of rumen microbiota in sheep at the phylum level.

Fig. 2

Correlogram showing the microbe-microbe correlations at the phylum level in sheep rumen.

The functional annotation (Suppl. Fig. 2) of the whole metagenome revealed that the genes related to the carbohydrate metabolism were most abundant (14.5%), followed by protein metabolism (10.6%) and amino acid and derivatives (9.9%). At the sub-system level 2, the monosaccharides-related genes were most abundant (27%). After that, central carbohydrate (23%) and oligosaccharides (18%) related genes were abundantly distributed. The genes related to one-carbon metabolism, including methanogenesis, were ranked fourth most abundant (∼8%; Fig. 3).

Fig. 3

Classification of rumen whole metagenome at subsystem level 1, carbohydrate metabolism and one carbon metabolism representing the abundance of genes in different pathways.

16S rRNA sequencing of ruminal fluid samples from nine sheep generated 3,390,112 raw reads with an average size of 376,679 reads per sample. After quality filtering and chimera removal, 1,836,553 reads were considered for further analysis (Suppl. File 1). Clustering of filtered reads yielded 142 OTUs at 97% similarity. Rarefaction analysis confirmed that the rumen methanogens diversity and species richness were adequately covered in this study (Suppl. Fig. 3). Alpha diversity assessed through Shannon index (Suppl. Fig. 3) revealed that the quantitative measurement of microbial diversity across the samples was similar with an average of 2.89 ± 0.21 (mean ± SD).

Rumen methanogens diversity characterized by 16S rRNA sequencing revealed that all reads obtained in this study belonged to the phylum Euryarchaeota. To overcome the limited resolution of taxonomic assignments beyond the genus level, a recently developed precise taxonomic framework, RIM-DB was used in the present study. Taxonomic assignments revealed that the rumen methanogens affiliated to Methanobacteria, Methanomicrobia, and Thermoplasmata classes were present in the sheep rumen (Fig. 4a). Methanobacteria associated methanogens were the most prevalent, which alone represented 85% of the total archaeal community. Methanogens affiliated to the Thermoplasmata were second most abundant in the sheep rumen and constituted about 13% of the total archaea. The remaining (∼1.5%) methanogens in Indian sheep were assigned to the Methanomicrobia class. In the present study, the order level taxonomic assignment confirmed the presence of Methanobacteriales, Methanomicrobiales, Methanosarcinales, and Methanomassiliicoccales methanogens in the sheep rumen. Among all, Methanobacteriales were the most prominent, while Methanomassiliicoccales were found to be the second most abundant. Although methanogens affiliated with Methanomicrobiales and Methanosarcinales were identified, aggregately, both represented only ∼2% of the archaea.

Fig. 4

Methanogens taxonomic assignments at class, order, genus, and species levels; a) denotes the archaea distribution at class and order level, the inner circle represents the class level abundance; while outer circle pertaining to the order level distribution of rumen methanogens; b) Top and bottom bar plots represent the genus and species level distribution of methanogens, respectively, each genus and species has been shown in a different color.

Further archaeal community analysis confirmed that the methanogens affiliated to nine genera were present in the Indian sheep (Fig. 4b). Among these, Methanobrevibacter represented the most considerable fraction (∼82%) of the rumen archaeome. Methylotrophic genera illustrated the second highest abundance (∼12.8%), and Group 12 methanogens were highly abundant among them. Other methanogens associated to Methanomicrobium and Methanosphaera genera constituted only 1.30 and 1.10% of the archaeome, respectively. In our study, 18 species of methanogens using RIM-DB taxonomic framework were identified in the sheep rumen. Methanobrevibacter gottschalkii was the most abundant species of rumen methanogens in sheep and constituted more than half (∼56%) of the total archaeal community (Fig. 4b). Methanobrevibacter millerae was found to be the second largest species of methanogens. Both hydrogenotrophic Methanobrevibacter wolinii and methylotrophic Group 12 sp. ISO4-H5 constituted an almost equal fraction (∼11%) of the rumen archaeal community in Indian sheep (Fig. 4b, Suppl. File 1); while, Methaobrevibacter ruminantium and Methanobacterium sp. were individually distributed at a lower frequency of ∼2.5%. Methanomicrobium mobile, another important species of hydrogenotrophic methanogens, had a lower abundance (<1.5%) in the sheep rumen. Similarly, Methanosphaera sp. ISO3-F5 also constituted approximately 1% of the archaeal community. The remaining 10 species of methanogens with the individual contribution of <1% (range 0.001–0.946%) aggregately constituted <2.5% of the community (Suppl. File 1). The functional annotation under one-carbon metabolism revealed that methanogenesis-related genes after the serine-glyoxylate cycle were the second most abundant in sheep (Fig. 3). Functional analysis also revealed that the genes for converting methylated compounds into methane contributed 2.7% to the total, and the remaining 97.3% were attributed to the remaining two pathways. KEGG reference pathway for methane metabolism (00680) indicated the formation of methane via hydrogenotrophic (M00567) and methylotrophic (M00356) pathways, whereas the acetoclastic pathway (M00357) was not detected in the present study (Fig. 5).

Fig. 5

KEGG metabolic pathways by functional microbial communities. The highlighted pathways in red color.

Core microbiome analysis

Core microbiome analysis revealed that Methanobrevibacter, Group 12, Methanobacterium, Methanosphaera and Methanomicrobium were prevalent genera with a minimum prevalence, and detection threshold of 50%, and 0.01, respectively (Fig. 6a). Similarly, at the species level, Methanobrevibacter gottschalkii clade, Group12 sp. ISO4-H5, Methanobrevibacter millerae, Methanobrevibacter wolinii, Methanobrevibacter ruminantium clade, Methanobacterium sp., Methanosphaera sp. ISO3-F5 and Methanomicrobium mobile constituted the core microbiome of rumen archaea in sheep (Fig. 6b). However, further analysis with a minimum prevalence of 70% reduced the number of genera representing the core microbiome to three (Methanobrevibacter, Group 12, and Methanobacterium). Similarly, the number of species was also reduced from eight to six at 70% prevalence (Suppl. Fig. 4).

Fig. 6

Core microbiome methanogens identified using microbiome package with a minimum 50% prevalence at a) genus, and (b) species level. The prevalence of methanogens has been shown with color/arrows (hydrogenotropic and methylotropic) represent the genes associated with methane metabolism in sheep.

Network analysis

Network wiring of OTUs belonging to the top genera of the methanogens in sheep has been presented in Fig. 7. Network stability obtained through SPIEC-EASI across the methanogens community was 0.049. Overall, a positive correlation of 71% was observed among the methanogens. The majority of intra-genus correlations were positive; while, OTUs belonging to different genera have shown both positive and negative correlations. The most dominated genus Methanobrevibacter had a strong positive correlation with Methanomicrobium and Methanobacterium. On the other hand, the correlation of Methanomassiliicoccales at the genus level (Group10 and Group 9) had a negative correlation with Methanobrevibacter. Methanobrevibacter was positively and negatively correlated with Group 12 and Methanosphaera, respectively.

Fig. 7

Methanogens associated network of sheep rumen. Nodes representing OTUs belonging to the top methanogens genera are color coded. Correlation between OTUs can be positive (green edges) and negative (red edges). The scale of correlation and thickness of the edge is proportional.

Methane emission and VFA production

The mean dry matter intake (g/d) in sheep was 946 ± 2.91 g/d, whereas the daily methane emission was 19.7 ± 1.11 g/d (Table 2). The DMI and daily methane emissions were weakly correlated in the present study (P = 0.433). Methane yield calculated using daily methane emission and dry matter intake from the sheep was within a range of 14.5 to 26.1 g/kg DMI. The acetate constituted the most considerable fraction (69.5%) of the total VFA followed by propionate (16.2%) and butyrate (10.9%). The acetate to propionate ratio in the present study was 4.3 ± 0.08.

Table 2

Ruminal methane yield and VFA production in sheep fed on straw and concentrate diet.

Attributes	Animal									Mean	SE
	S1	S2	S3	S4	S5	S6	S7	S8	S9
DMI (g/d)	937	938	952	947	958	945	956	933	950	946	2.91
Methane (g/d)	22.6	21.8	24.8	17.0	18.9	21.1	17.4	19.5	13.8	19.7	1.11
MY (g/kg DMI)	24.1	23.2	26.1	18.0	19.7	22.3	18.2	20.9	14.5	20.8	1.19
TVFA (mM)	62.6	70.0	77.8	67.5	65.3	78.0	85.2	70.3	60.2	70.8	2.71
VFA (mM)
Acetate	44.7	48.5	55.4	46.4	45.0	53.0	59.5	48.7	41.4	49.2	1.92
Propionate	11.3	10.6	12.4	11.6	10.9	12.7	13.1	11.4	9.31	11.5	0.39
Isobutyrate	0.52	0.72	0.64	0.86	0.81	0.90	0.67	0.81	0.82	0.748	0.04
Butyrate	5.03	8.69	7.96	6.87	6.85	9.34	10.41	7.44	7.09	7.743	0.53
Isovalerate	0.60	0.95	0.79	1.28	1.20	1.30	0.97	1.37	1.07	1.060	0.08
Valerate	0.42	0.51	0.54	0.57	0.52	0.69	0.57	0.58	0.45	0.538	0.03

DMI- dry matter intake; MY- methane yield; TVFA- total volatile fatty acids; VFA- volatile fatty acids; SE- standard error.

Methanogens impact on methane variability

The Methanobrevibacter wolinii, Methanobrevibacter smithii and Methanobrevibacter sp. ABM4 have a strong influence on PC1, whereas Methanosphera stadtmanae and Group 8 had more weightage on PC2 (Fig. 8). However, methanogens having a major influence on PC1 and PC2 were not correlated as Methanomicrobium mobile and Methanosphaera sp. ISO3-F5/Group 12 were in opposite orientation, therefore, considered negatively correlated in sheep archaeome. Similarly, Methanobrevibacter millerae and Methanobrevibacter gottschalkii/ Methanimicrococcus blatticola were also negatively correlated. Interestingly, the various species of Methanobrevibacter demonstrated a different correlation to methane emission. On the other hand, Methanobrevibacter smithii and Methanobrevibacter sp. ABM4, Methanobrevibacter gottschalkii and Methanimicrococcus blatticola, Group 10 and Methanomicrobium mobile; Methanosphaera sp. ISO3-F5 and Group 12 sp. ISO4-H5 had a close angle on PC loadings and therefore were positively correlated.

Fig. 8

PCA biplot depicting the influence of ruminal methanogens distribution on enteric methane emission. The length of the arrows in blue color is equivalent to the contribution of corresponding metadata to daily methane emission. The orientation on the graph and the degree of angles represent the intensity of correlation.

Discussion

Recent high-throughput sequencing technology enabled researchers to uncover the microbial capabilities and phenotypes yet to be described. A study in Vechur and crossbred cattle fed on a comparable diet with our study also reported the dominance of Bacteroidetes and Firmicutes (Sadan et al., 2020). Similarly, the Bacteroidetes and Firmicutes ratio in this study was also similar (0.45 Vs 0.44). The comparable results among the cattle and sheep on similar feeding regimes indicated that the host species do not impact the distribution of Bacteroidetes and Firmicutes in the rumen. From a previous study, Henderson et al. (2015) concluded that host is less influential on microbial community composition. Even the distribution of Bacteroidetes and Firmicutes in Holstein calves at a quite young age (1–8 weeks) and contrast diet (calf starter) were similar to the finding in the present study (Meale et al., 2016). The higher abundance of Bacteroidetes in this study can be attributed to the high fibre diet as Bacteroidetes mainly possess a strong capability to disintegrate the protein and polysaccharides (Huo et al., 2014, Pitta et al., 2016) and are net H2 utilizers (Stewart et al., 1997). Similarly, Firmicutes are efficient in the degradation of hemicellulose, cellulose and lignin, whereas Bacteroidetes are involved in the depolymerization of lignin (Gavande et al., 2021). Present study indicated an overall greater abundance of fibrolytic microbiota in the sheep rumen fed on finger millet straw and concentrate diet.

All the OTUs in this study were assigned to the phylum Euryarchaeota. The bacterial community has been reported taxonomically far rich than the ruminal archaea (Henderson et al., 2015), which is somewhat represented with similar universality and limited diversity. Most ruminal methanogens are associated with the phylum Euryarchaeota (Kim et al., 2011); however, Crenarchaeota and Nanoarchaeota affiliated methanogens are also reported in the rumen (Abecia et al., 2014, Shin et al., 2004, Xue et al., 2019). The most significant proportion of Methanobacteriales among the ruminal methanogens in this study was in good agreement with the previous studies (Huang and Li, 2018, Seedorf et al., 2015, Wright et al., 2004). Similarly, the Methanomassiliicoccales distribution was also in consonance with the global datasets (Huang and Li, 2018, Seedorf et al., 2015, Seshadri et al., 2018). The study indicated that the hydrogenotrophic methanogens dominate the archaeal community in sheep rumen, whereas aceticlastic methanogens had limited distribution frequency (∼0.10%). This was also supported by the KEGG reference pathway analysis revealing presence of hydrogenotrophic (M00567) and methylotrophic (M00356) pathways (Fig. 5). This is in agreement with the previous studies reporting the prevalence of hydrogenotrophic methanogens (Hungate, 1966, Hungate et al., 1970) and particularly of Methanobrevibacter (Henderson et al., 2015). A previous report also confirmed the methane metabolism (ko00680) as the most enriched pathway (Shi et al., 2014).

The methanogens community in Indian sheep appears to be more diversified than in Merino sheep (Wright et al., 2006, 2004), Finn-Dorset sheep (Snelling et al., 2014), and goats (Cheng et al., 2009). All the above studies reported the prevalence of Methanobacteriales and Thermoplasmatales methanogens, whereas the present study revealed the existence of Methanobacteriales, Methanomicrobiales, Methanosarcinales and Methanomassiliicoccales. The greater diversification in present study is also evidenced by the presence of methanogens from 9 genera and 20 species (Suppl. File 1).

These studies reported the limited distribution of Methanomicrobiales or Methanococcales, while the proportion of Methanomicrobiales in our study was within the global range of 5–15% (Henderson et al., 2015, Janssen and Kirs, 2008, Jeyanathan et al., 2011, Wright et al., 2006). On the other hand, few studies unexpectedly reported a very high distribution (52–100%) of Methanomicrobiales in the rumen (Chaudhary and Sirohi, 2009, Singh et al., 2012, Yanagita et al., 2000). Methanogens affiliated to nine genera with an overall dominance of the Methanobrevibacter (∼82%) were identified in this study. The dominance of Methanobrevibacter genus in archaeal community was in agreement with the various studies, reported the abundance in the range of 61–90% (Danielsson et al., 2017, Janssen and Kirs, 2008, Seedorf et al., 2015, Snelling et al., 2014). In contrast, few studies reported that the ruminal archaea community is dominated by the rumen cluster C- Thermoplasmatales (Sundset et al., 2009, Wright et al., 2008). Methanomassiliicoccales associated genera constituted 12.8% of the total archaea in sheep, where Group 12 being the dominant (Fig. 4b, Suppl. File 1). The methylotrophic methanogens distribution in Indian sheep was consistent with the previous studies (Jin et al., 2017, Seedorf et al., 2015) reported the abundance between 0.5 and 12%. The prevalence of hydrogenotrophic methanogens and Group 12 indicated that there could be a metabolic association between them. The positive correlation in network analysis also strengthens this hypothesis; however, this needs to be confirmed by experimentation.

Methanobrevibacter gottschalkii was the most prevalent methanogen in the sheep rumen at the species level (Fig. 6b), whereas Wright et al. (2006) from a study in Queensland sheep reported a limited distribution of Methanobrevibacter gottschalkii (∼10%). Other methanogens such as Methanobrevibacter millerae, Methanobrevibacter wolinii, Group 12 sp. ISO4-H5, and Methanobrevibacter ruminantium were also the constituents of core archaeome. Results on the dominance (33–81%) of Methanobrevibacter gottschalkii was consistent with the previous reports (Henderson et al., 2015, Janssen and Kirs, 2008). Huang and Li (2018) reported Methanobrevibacter millerae as one of the prominent methanogens (41.3%) in the Tibetan sheep, however, in our study it constituted ∼13% of the total rumen methanogens in Indian sheep. This deviation in the abundance of Methanobrevibacter millerae could be due to the difference in diet from our study (alpine meadow pasture Vs roughage based). Moreover, the geographical region (Huws et al., 2018, Seshadri et al., 2018), environmental conditions (Zhang et al., 2020), genotype, and developmental age of the animals may also influence the community composition. Feeds and diet composition are known to affect the rumen microbiota. However, the interaction of diet and ruminal methanogens is not fully understood. Diet is one of the major influencing factors determining the rumen microbial community composition (Kim et al., 2017). A remarkable difference in the rumen archaeal community was observed among the feeding of three diets based on pasture, oat hay and lucerne hay (Wright et al., 2007). However, no impact of the increasing proportion of concentrate on methanogens density was noticed previously (Hook et al., 2010, Lillis et al., 2011, Popova et al., 2011). The methanogens diversity in this study was comparable with a previous study by our group (Baruah, 2019), having similar dietary ingredients (finger millet straw and concentrate), but a different ratio (50:50 Vs 45:55). Minor differences were noted in the abundance of Methanobrevibacter gottschalkii clade, Methanobrevibacter ruminantium, and Group 12 sp. ISO4-H5 between two studies.

Among these, feed is one of the most important factors influencing the daily methane emission (Malik et al., 2015a). Though the mean MY (20.8 g/kg DMI) in this study was in agreement with the global datasets (Charmley et al., 2016, Hristov et al., 2013, Malik et al., 2021, Pinares-Patino et al., 2014); nevertheless, a variation in daily methane yield despite the same diet and intake was recorded between the sheep. Previous reports (Blaxter and Clapperton, 1965, Grainger et al., 2007) also confirmed a considerable variation in methane output among individuals fed on the same diet. This variation between the animals could be confounded by the variation in the methanogens community structure. This was also evidenced from the PCA biplot, where ∼30% variation in methane yield was primarily attributed to the hydrogenotrophic methanogens and about 19% to the methylotrophic methanogens.

Conclusions

Whole metagenome analysis revealed that the bacterial community dominated the rumen microbiota in sheep. Three dominant phyla were Bacteroidetes, Firmicutes, and Proteobacteria constituted 57, 25, and 6% of the microbiota. Functional analysis of the whole metagenome revealed the abundance of carbohydrate and protein metabolism genes. The greater abundance of methanogenesis genes characterized the one-carbon metabolism. KEGG reference pathway for methane metabolism indicated the formation of methane through hydrogenotrophic and methylotrophic pathways, whereas the acetoclastic pathway was not functional in sheep. It is inferred from the results that Methanobrevibacter is the most prominent genus in the Indian sheep, whereas Methanobrevibacter gottschalkii is the largest distributed species of methanogens.

Data accessibility

The archaeal metagenome sequencing reads from the experiment are accessible at the NCBI Sequence Read Achieve (SRA; https://www.ncbi.nlm.nih.gov/subs/sra) accession numbers SAMN15366540- SAMN15366548 under BioProject PRJNA641793. The OTUs abundance and taxonomical assignment data are available in the supplementary files.

Authors’ contributions

PKM, RB and HR conceived, designed, coordinated the study and helped in drafting the manuscript; ST, APK, and VS carried out molecular laboratory work, bioinformatic analysis, data visualization. All authors gave final approval for publication.

Funding

“This work was supported by the International Livestock Research Institute [window III of CGIAR, 2019] and Indian Council of Agricultural Research, New Delhi through a collaborative project on Methane Emission and its Mitigation”.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.