Methane yield phenotypes linked to differential gene expression in the sheep rumen microbiome.

Ruminant livestock represent the single largest anthropogenic source of the potent greenhouse gas methane, which is generated by methanogenic archaea residing in ruminant digestive tracts. While differences between individual animals of the same breed in the amount of methane produced have been observed, the basis for this variation remains to be elucidated. To explore the mechanistic basis of this methane production, we measured methane yields from 22 sheep, which revealed that methane yields are a reproducible, quantitative trait. Deep metagenomic and metatranscriptomic sequencing demonstrated a similar abundance of methanogens and methanogenesis pathway genes in high and low methane emitters. However, transcription of methanogenesis pathway genes was substantially increased in sheep with high methane yields. These results identify a discrete set of rumen methanogens whose methanogenesis pathway transcription profiles correlate with methane yields and provide new targets for CH4 mitigation at the levels of microbiota composition and transcriptional regulation.

Methane (CH4) accounts for 14% of total global greenhouse gas emissions and is the second largest contributor to global warming (Intergovernmental Panel on Climate Change 2007). Almost a third (28%) of anthropogenic CH4 emissions are due to enteric fermentation in livestock (Yusufa et al. 2012), an impact predicted to rise further due to an increased worldwide demand for meat, milk, and other animal products. The dominant source of CH4 emissions from livestock is from ruminants (Naqv 2011), where CH4 is formed as a byproduct of feed fermentation in the fore-stomach (rumen) by CH4-producing archaea, known as methanogens (Boone et al. 1993). Methanogens use a limited range of substrates, including CO2/H2, formate, acetate, and methyl compounds (Zinder 1993; Hook et al. 2010). Only a few rumen methanogens have been cultivated or characterized in detail, and their respective contributions to CH4 production under in vivo conditions in livestock remain poorly defined (Buddle et al. 2011). Measurements of ruminant CH4 emissions are mainly from animal trials in which the effects of particular diets or inhibitors of CH4 formation were assessed (Machmüller et al. 2003; Lila et al. 2005; Nkrumah et al. 2006; Denman et al. 2007; Hegarty et al. 2007; Martínez-Fernández et al. 2014). However, a program investigating natural variation in CH4 emissions from sheep is underway in New Zealand, and measurements made using both tracer gas techniques and open-circuit respiration chambers suggest there is repeatable and heritable variation between individual animals in CH4 yield (CH4 produced per unit of feed consumption) (Lassey et al. 1997; Pinares-Patiño et al. 2011a,b, 2013). In the present study, CH4 production in a cohort of New Zealand sheep was examined under controlled experimental conditions, followed by deep metagenomic and metatranscriptomic sequencing of their rumen contents to examine the microbial contribution to this variation. By comparing the microbiota from low- and high-CH4-yielding animals, we identify and characterize specific archaeal clades and transcriptional characteristics that appear to explain CH4 yield differences in sheep.

Results

Measurement of variation and reproducibility of CH4 yield in sheep

Twenty-two age-matched crossbred rams (Fig. 1A) fed on a pelleted lucerne (alfalfa) diet had CH4 yields measured using open-circuit respiration chambers (Fig. 1B) at two time points separated by 2 wk. We found a high concordance of CH4 yields from the same sheep at different time points (P = 0.85 for differences) but substantial variation between sheep (P = 0.0001, one-way ANOVA test) (Fig. 1C), findings consistent with previous observations (Grainger et al. 2007; Dengel et al. 2011; Pinares-Patiño et al. 2011b). The lowest- and highest-yielding sheep groups differed by 4.41 g of CH4 produced per kilogram of dry matter intake (38.5%, P = 0.0002) (Fig. 1C).

Figure 1.

The measurement of CH4 yields in sheep. (A) New Zealand sheep used for this study. (B) CH4 yields from the sheep in grams of CH4/kg dry matter intake (DMI) were measured using open-circuit respiration chambers (http://www.globalresearchalliance.org). (C) CH4 yield measurements from 22 sheep (each with two time points) sorted by mean values. Four high (red) and four low (blue) emitters are selected for further study. P-value indicates the statistical significance of the differences in CH4 yield between the two selected groups.

Comparison of relative abundance of different microbial populations in low and high CH4 yield sheep

We performed deep metagenome and metatranscriptome sequencing (Supplemental Fig. 1) on rumen content samples from four rams with the highest mean CH4 yields, four rams with the lowest mean CH4 yields, and two rams with intermediate CH4 yields at two separate times (20 samples total). We generated ∼50 Gb of unamplified metagenome whole-genome shotgun (WGS) sequencing data from each rumen sample, totaling 1020 Gb (Supplemental Table 1). To compare the microbial community structures of low- and high-CH4-yielding sheep, we aligned the metagenome WGS reads to bacterial and archaeal 16S and eukaryotic 18S ribosomal RNA (rRNA) gene sequences from the comprehensive SILVA database (Pruesse et al. 2007). None of the microbial domains showed a significant change in relative abundance between low and high CH4 yield sheep (Fig. 2A). We also specifically quantified the overall proportion of methanogens in low and high CH4 yield sheep, which also did not reveal any significant differences between low and high CH4 emitters (Fig. 2B), and validated this by quantitative polymerase chain reaction (qPCR) using universal primers targeting the 16S rRNA genes of methanogens (Supplemental Fig. 2). Comparison of the methanogen community structure at the class level using the SILVA and Greengenes databases (see Methods) (DeSantis et al. 2006; Pruesse et al. 2007) displayed a similar community composition and structure between the high and low CH4 yield sheep, with Methanobacteria the dominant class, followed by Thermoplasmata and Methanomicrobia (Fig. 2C; Supplemental Fig. 3A,B). Detailed analysis of the community structure at the genus and subgenus level within Methanobacteria showed elevated Methanosphaera spp. in the low CH4 yield sheep and higher relative abundances of organisms belonging to the Methanobrevibacter gottschalkii clade in the high CH4 yield sheep (Supplemental Fig. 4A), which was independently confirmed by pyrotag sequencing, which provided similar results (Supplemental Fig. 4B).

Figure 2.

Comparison of relative abundance of different microbial populations in low and high CH4 yield sheep. (A) Relative abundance of microbial domains in low and high CH4 yield sheep. (B) Relative abundance of methanogenic and nonmethanogenic archaea in low and high CH4 yield sheep. (C) Relative abundance of classes of CH4-producing Euryarchaeota in low and high CH4 yield sheep. (NS) No statistical difference in Wilcoxon rank-sum test in each subgroup.

Metagenomics analysis and methanogenesis genes abundance comparison

To examine possible differences in the presence of functionally relevant methanogenesis genes, three major methanogenesis pathways were examined: the predominant CO2/H2 pathway (Fig. 3A) of hydrogenotrophic methanogens, as well as the alternative aceticlastic and methylotrophic pathways, which have genes in common with other nonmethanogenic pathways (Supplemental Fig. 5A; Morgavi et al. 2010; Glass and Orphan 2012). While a total of 297 KEGG genes show significant enrichment in low or high CH4 yield sheep by multiple statistical analyses (Supplemental Fig. 6A), these included only one gene (fdhA, K05299) that is part of the CH4 production pathway (ko00680) (Supplemental Fig. 6B). Given that different KEGG genes may have the same biochemical functionalities (Kanehisa and Goto 2000), we then grouped the KEGG genes based on the biochemical reactions they catalyze and repeated the statistical analyses. Except for acetyl-CoA synthetase (Enzyme Commission [EC]: 6.2.1.1), an enzyme involved in many metabolic pathways other than methanogenesis, none of the gene encoding enzymes involved in methanogenesis were significantly enriched in high CH4 yield sheep (Fig. 3B; Supplemental Fig. 5B). In summary, these results suggest that although there are some shifts in subpopulations of methanogens, the high CH4 yield levels in sheep are unlikely to be due to an increased relative abundance of methanogens and methanogenesis pathway genes.

Figure 3.

Comparisons of gene and transcript abundance for enzymes involved in methanogenesis between high and low CH4 yield sheep. (A) Diagram of CO2/H2 methanogenesis pathway shows enzymes involved in each biochemical reaction. (B,C) Gene (B) and transcript (C) abundance for each enzyme. (D) Transcriptions per gene for each enzyme. (RPM) Reads per million; (NS) no statistical significance in Wilcoxon rank-sum test; (*) P < 0.05; (**) P < 0.01. Error bars, SE.

Metatranscriptomics analysis and methanogenesis pathway genes expression comparison

To explore the possibility that changes in methanogenesis-related gene expression might be responsible for the differences in CH4 yields, we generated 6.6 Gbp of metatranscriptome sequences on average from each mRNA-enriched sample (Supplemental Fig. 7; Supplemental Table 1). Using the same method as above to quantify gene abundance, transcripts from 349 KEGG genes were significantly enriched in high or low CH4 yield sheep (Supplemental Fig. 8A). Notably, three of the top 10 KEGG genes with increased transcripts in high CH4 yield sheep code for enzymes in the methanogenesis pathway (Supplemental Table 2). These three genes belong to an operon encoding the subunits of methyl coenzyme M reductase (mcr, EC: 2.8.4.1), an enzyme that catalyzes the final and rate-limiting step during CH4 biogenesis (Cedervall et al. 2010). Furthermore, the CH4 metabolism pathway (ko00680) was the most significantly enriched pathway among ∼300 known KEGG pathways (Supplemental Fig. 8B).

The observation of an increased number of mcr transcripts in high CH4 yield sheep prompted an examination of the transcript abundance for each of the enzymes involved in methanogenesis (Fig. 3A). In contrast to gene abundance (Fig. 3B), transcript abundance of all of the genes encoding enzymes involved in the CO2/H2 pathway was significantly increased in high CH4 yield sheep (Fig. 3C), whereas the transcript abundance for the two alternative methanogenesis pathways (aceticlastic and methylotrophic) was not altered, other than the final rate-limiting step that is shared with the CO2/H2 pathway (Supplemental Fig. 5C). Specifically, genes encoding MCR showed the largest difference in transcript levels (Fig. 3C), suggesting that methanogenesis gene expression was elevated in high-CH4-yielding sheep. Importantly, increased transcript, but not gene abundance, was observed for every step in the pathway (Fig. 3B,C). To rule out that significant differences in transcript levels arise from subtle, subsignificant differences in gene abundance, we normalized transcripts by gene counts. For seven of the 10 steps in the CO2/H2 methanogenesis pathway, we observed significant, twofold to fourfold increases in high-CH4-yielding sheep (Fig. 3D). These findings indicate that genes involved in the CO2/H2 methanogenesis pathway were expressed at significantly higher levels in high-CH4-yielding sheep. Furthermore, metagenome and metatranscriptome data from rumen samples of two sheep intermediate in their CH4 yields had methanogenesis transcripts/genes intermediate between those from low and high animals, further supporting this relationship (Supplemental Table 3).

Assembly and validation of mcr/mrt operons

To identify and further characterize the methanogen species responsible for the increased methanogenesis transcript levels, we focused our analysis on mcr operons. The mcr genes are categorized into two related operons (mcr and mrt) that encode isoenzyme complexes (McrABG and MrtABG), which are thought to be differentially regulated by H2 concentration (Reeve et al. 1997). All methanogen genomes encode either the mcr or the mrt operon, and some encode both. We combined the metagenomic and metatranscriptomic sequence data from all sheep and assembled 35 distinct mcr/mrt operons, of which 28 (80%) were considered to be full-length based on the subunits they contained (Supplemental Table 4; Cedervall et al. 2010). The authenticity of these assembled operons was confirmed by PCR amplification of the predicted sequences directly from the metagenomic DNA samples (Supplemental Fig. 9) and by sequencing the products (Supplemental Table 5), which showed that nearly all of the mcr/mrt operons predicted on the basis of short-read assemblies extracted from the metagenomic and metatranscriptomic data represent authentic operons present in the sheep rumen methanogens.

Phylogenetic analysis of mcr/mrt genes in sheep rumen

Among the subunits of mcr/mrt operons, the genes encoding the alpha subunits (mcrA/mrtA) have been established as reliable phylogenetic markers for methanogens (Hallam et al. 2003; Paul et al. 2012). Phylogenetic analysis of the assembled mcr/mrt operons based on their alpha subunits and 146 previously identified full-length McrA/MrtA subunits from the protein database of NCBI indicated that methanogens in the rumen of sheep in this study clustered into three groups (Sheep Rumen MCR Groups I–III, SRMR1–3) (Fig. 4A). Two groups (SRMR2 and SRMR3) represent substantial expansions within the Methanobacteria of mrtA and mcrA genes, respectively. The SRMR2 group includes new mrtA genes within both Methanobrevibacter spp. (encompassing the Methanobrevibacter ruminantium [Leahy et al. 2010] and M. gottschalkii clades [Hansen et al. 2011]) and Methanosphaera spp. (Fricke et al. 2006), whereas SRMR3 contains new mcrA genes within Methanobrevibacter spp. (Fig. 4A; Supplemental Fig. 10). In contrast, SRMR1 represents a cluster of mrtA genes whose sequences show considerable divergence from all other known full-length genes, suggesting they belong to a group of poorly characterized rumen methanogens (Fig. 4A). The existence of a seventh order of methanogenic archaea, called Methanoplasmatales, has previously been proposed, but evidence of rumen representatives has been limited to 16S rRNA gene sequences and partially sequenced mcrA/mrtA genes (Paul et al. 2012). More recent evidence suggests that these organisms are methylotrophic, using methylamines, and possibly methanol, as substrates (Poulsen et al. 2013). Our analysis, based on full-length mcrA/mrtA genes, provides direct support that SRMR1 represents a phylogenetically divergent cluster of rumen methanogens and enables insights into the biology of these species through metatranscriptomic analyses.

Figure 4.

Phylogenetic analysis of methanogens in sheep rumen. (A) A phylogenetic tree constructed based on full-length methyl coenzyme M reductase alpha subunit (McrA/MrtA) protein sequences. Known McrA/MrtA proteins from NCBI are shown in black; new ones from this study, in color. (B) Genes and transcripts for three groups of identified sheep rumen methanogens. (RPM) Reads per million; (NS) no statistical significance in Wilcoxon rank-sum test; (*) P < 0.05; (**) P < 0.01. Error bars, SE. (C) Relative contribution of each group of sheep rumen methanogens to the overall abundance (RPM) of genes and transcripts in low and high CH4 yield sheep. The sizes of each pie indicate the abundance of genes/transcripts.

To understand the contribution of the methanogen species in the three SRMR clades to high CH4 yield, we aligned the metagenome and metatranscriptome WGS reads to the alpha subunit genes of 35 mcr/mrt operons to quantify the gene and transcript abundance, respectively. The methanogens in these three groups combined are likely to represent the majority of methanogens in the sheep rumen, since >90% of the total mcrA/mrtA metagenomic reads and >97% of the total mcrA/mrtA metatranscriptomic reads derived from methanogens can be mapped to one of the three groups (Fig. 4C). To obtain insight into the general function of each of the three SRMR groups, we considered all mcrA/mrtA genes within each group cumulatively and compared the total counts from the low and high CH4 yield animals. Consistent with our analyses based on 16S rRNA genes and KEGG enzymes, none of the groups showed significant differences in mcrA/mrtA gene abundance (Fig. 4B). However, there were significantly more transcripts from SRMR1 and SRMR3 in high CH4 yield sheep (Fig. 4B). Furthermore, metatranscriptome data from the two sheep with intermediate CH4 yields had SRMR1 and SRMR3 transcripts at levels between the low and high animals, providing an additional link between expression of these SRMR groups and CH4 yields (Supplemental Fig. 11). Overall, transcripts from SRMR1 and SRMR3 contributed to >90% of the total mcrA/mrtA transcripts in both low and high CH4 yield sheep and were 2.84- and 2.85-fold more abundant in high CH4 yield sheep for SRMR1 and SRMR3, respectively (Fig. 4C). SRMR1 and SRMR3 accounted for 34.6% and 63.1%, respectively, of the overall transcript increase in high CH4 yield sheep. In contrast, transcript levels from SRMR2 were very low, were only slightly increased in low CH4 yield sheep, and did not contribute to the overall transcript increase in high CH4 yield animals (Fig. 4B). To validate these findings, we used qPCR and gene-specific primers (Supplemental Table 6) to quantify five randomly selected assembled mcrA/mrtA genes from DNA and RNA from the four high CH4 yield, four low CH4 yield, and two intermediate CH4 yield sheep used for metagenome and metatranscriptome sequencing, as well as from two sheep with intermediate CH4 yields. Gene and transcript abundances (Supplemental Fig. 12C,D) were consistent with the results from metagenomic and metatranscriptomic analysis (Supplemental Fig. 12A,B), which support transcriptional up-regulation as the primary microbial mechanism contributing to higher CH4 yield levels among sheep.

Discussion

CH4 emitted from sheep is formed by methanogenic archaea in the rumen as an end product of microbial degradation of forage material. It is therefore likely that the ruminal microbiome contributes to the host CH4 yield phenotype. The exact mechanism causing the high and low CH4 yield phenotypes observed in sheep is still unclear, and our understanding of the microbial contribution to differences in CH4 yields among sheep has been limited by the low throughput of previous cellular and molecular manipulations (Warnecke et al. 2007; Pope et al. 2010; Hess et al. 2011). Also, previous studies have suggested that microbial-derived phenotypes, including CH4 production levels, are primarily determined by microbial abundance profiles (Kao-Kniffin et al. 2011; Fox 2012). In contrast, the deep metagenomic and metatranscriptomic sequencing of rumen content in the present study revealed that increases in CH4 output are primarily associated with increases in the expression of methanogenesis pathway genes.

A possible mechanism explaining CH4 yield differences between animals is based on the amount of time that feed particles are retained in the rumen (Benchaar et al. 2001), with longer particle retention times leading to higher CH4 yields. Particle retention time in ruminants is known to be a heritable trait (Orskov et al. 1988; Smuts et al. 1995) and may explain at least some of the CH4 yield variation observed in sheep (Pinares-Patiño et al. 2003). Recently, CH4 yield in sheep in Australia has been directly correlated with the retention time of feed particles and liquid and with the total amount of feed particles and rumen volume (Goopy et al. 2013), further supporting this view.

Differential particle retention time may explain our findings of altered expression of methanogenesis pathway genes in sheep via a substrate-mediated effect. Differences in the passage rate of particles through the rumen is predicted to affect ruminal H2 levels according to a model based on microbial growth kinetics and fermentation thermodynamics (Janssen 2010). In this model, an increased particle passage rate is associated with higher rumen H2 concentrations, a thermodynamic negative feedback of H2 that results in less H2 formation by the fermentative microbes and, hence, less CH4 formation. Conversely, slower particle passage results in lower H2 concentrations, enhanced H2 formation during fermentation, and more CH4. These hypotheses are consistent with the finding that low CH4 yield sheep have fewer H2-producing bacteria and high CH4 yield sheep have more H2-producing bacteria in their rumens (Kittelmann et al., unpubl.).

Under ruminal conditions of slower particle passage rate and lower H2 concentrations, there will be a higher turnover rate of a smaller H2 pool through the methanogenesis pathway to account for the elevated CH4 formed. The lower ruminal H2 concentration means that methanogens have to increase expression of methanogenesis genes to maintain the H2 turnover rate. This is because enzyme concentrations as well as substrate concentrations can limit the flux through a pathway, and increasing enzyme expression partially overcomes the limitation of lower substrate concentrations (Morgan et al. 1997; Enoki et al. 2011; Walker et al. 2012; Browne and Cadillo-Quiroz 2013) Conversely, a high particle passage rate and high H2 conditions would require a lower level of expression of methanogenesis pathway genes to permit the same flux.

While there have been few studies on characterizing rumen microbial populations associated with natural variation in ruminant CH4 yields (Kittelmann et al. 2013), there have been numerous investigations on feedlot cattle selected for efficiency of feed conversion (also known as residual feed intake [RFI]), for which some CH4 yields' data are also available. Low-RFI animals are considered to be feed efficient and have lower CH4 yields compared with high-RFI, or feed-inefficient, animals (Nkrumah et al. 2006; Hegarty et al. 2007). Comparisons of ruminal microbiomes between low- and high-RFI animals using a variety of methods have shown differences in bacterial and archaeal community profiles correlated with RFI, although these associations are often influenced by the energy content of the diet (Guan et al. 2008; Zhou et al. 2010; Carberry et al. 2012; Hernandez-Sanabria et al. 2012). Methanogen-related differences observed in these studies included a specific high-RFI–related PCR-DGGE band associated with Methanobrevibacter smithii PS (Zhou et al. 2010), an elevated abundance of Methanosphaera stadtmanae, and Methanobrevibacter sp. strain AbM4-like sequences in high-RFI animals (Zhou et al. 2009), and a higher abundance of M. smithii genotypes in high-RFI animals (Carberry et al. 2014). Where measured, total methanogen densities in the rumen contents did not differ between the feed efficiency groups, indicating that the composition of the methanogenic community was the important difference. These observations are generally consistent with our findings of no changes in total methanogen numbers and an increase in the relative abundance of the M. gottschalkii group within the Methanobacteria. However, the elevated levels of Methanosphaera spp. in high-RFI cattle relative to low-RFI animals (Zhou et al. 2009) differ from our observation of elevated Methanosphaera in the low-CH4-yielding sheep. This may be due to the large difference between the diets fed (high grain feedlot diet for cattle [Zhou et al. 2009] vs. pelleted lucerne diet for sheep [this study]) or to innate differences between ruminant species (cattle vs. sheep).

The main findings of this study indicate that there are strong correlations between the expression levels of the hydrogenotrophic methanogenesis pathways in rumen methanogens and CH4 yields in sheep, in the absence of significant changes in methanogen community structure or relative abundance. This indicates a response of methanogenesis functions of the resident methanogens to the supply of their main substrate, H2. We predict that these gene expression changes are indirectly controlled by particle retention time or digesta passage rate in sheep. This is an avenue for future investigation within New Zealand’s sheep CH4 screening program, with the long-term goal of selecting animals with lower CH4 yields without compromising their productivity or reproductive ability. Furthermore, the identification of specific groups of methanogens that encode up-regulated methanogenesis genes correlated with high CH4 yield in sheep confirms current gene targets under investigation and provides new microbial and pathway targets for CH4 mitigation technologies in ruminants.

Methods

Sheep CH4 yield measurements and rumen contents sampling

Based on previous CH4 yield data and Central Progeny Testing breeding values (Pinares-Patiño et al. 2013), 11 high-CH4-yielding and 11 low-CH4-yielding rams were selected from the Woodlands Research Station progeny flock. Remeasurement of CH4 yields from these rams was conducted twice in respiration chambers (Fig. 1B) at the New Zealand Ruminant Methane Measurement Center, AgResearch Grasslands, Palmerston North, New Zealand, after adaptation to pelleted lucerne diet (for composition, see Supplemental Table 7) for 2 wk. Rumen contents were collected from all 22 sheep on two occasions (June 13 and June 28, 2011) immediately after CH4 measurements were made, by stomach intubation 4 h after morning feeding. The pH of the rumen contents was measured, and the samples were immediately snap-frozen as pellets in liquid N2 and stored at −85°C for DNA and RNA extraction.

DNA and RNA extraction

Based on CH4 yields, frozen rumen samples were selected from four high, four low, and two intermediate sheep (each at two time points, 20 samples total) and used for DNA extraction. For metagenomic and pyrotag sequencing, DNA was extracted using a “Repeated Bead Beating and Column (RBB+C) purification” method (Yu and Morrison 2004). For metagenomic sequencing of large paired-end insert libraries, high-molecular-weight DNA was extracted using a nonmechanical lysis DNA extraction method (Rosewarne et al. 2011). The extracted DNA was checked for concentration and molecular weight (Supplemental Material).

RNA was extracted from the same batch of rumen contents for transcript analysis using a hot lysis-acid phenol extraction method (Supplemental Material).

SSU rRNA gene sequencing and analysis

Amplicons of archaeal ssrRNA genes were generated according to the method previously described (Supplemental Material; Kittelmann et al. 2013) using Qiagen Taq PCR MasterMix (Qiagen) and were sequenced using 454 pyrosequencing (Roche). Sequence reads were quality filtered, assigned to samples by their nucleotide barcodes, size and quality filtered (≥220 bp, with <3% low-quality bases), and clustered by pyroclust via the PyroTagger pipeline (Kunin and Hugenholtz 2010).

A representative sequence from each cluster was BLAST-searched against the Greengenes database (DeSantis et al. 2006) and assigned a taxonomic identity or an operational taxonomic unit (OTU). Singletons were removed; diversity metrics were calculated; and the taxonomic and phylogenetic assignments were used to compare the types of microbial communities present. The abundance of Euryarchaeota classes was compared using the Wilcoxon rank-sum test. For the detailed analysis within the class Methanobacteria, archaeal 16S rRNA gene sequencing reads were clustered into OTUs at a 97% sequence similarity by uclust using the QIIME pipeline (Caporaso et al. 2010). A representative sequence from each cluster was BLAST-searched against an in-house rumen archaea reference database (Janssen and Kirs 2008) and assigned a taxonomic identity. OTUs were summarized at clade level, and clades were tested for statistical significance between high- and low-CH4 animals using the Wilcoxon rank-sum test using R software (R Development Core Team 2014).

Metagenomic library construction and sequencing

For each short insert metagenomic library, 2 μg of DNA extracted from each rumen contents sample was used as a template for a WGS sequencing unamplified library (tight insert size of 250 bp for high-throughput sequencing from both ends by 2 × 150 bp using an Illumina HiSeq 2000 instrument). Two additional libraries with insert sizes of 8 kb for combined samples with low and high CH4 yield traits, respectively, were constructed and sequenced to facilitate genome assembly. Four Illumina HiSeq 2000 runs were conducted (∼1 Tbp of sequence), and raw sequence data (plus initial titration runs) were passed through a JGI-developed filtering program to filter out known Illumina sequencing and library preparation artifacts and host contaminations. For the 8-kb libraries, duplicated read pairs derived from polymerase chain reaction (PCR) amplification during library preparation were identified and consolidated into a single consensus read pair. Artifact-filtered sequence data from 250-bp tight insert libraries were joined by FLASH (Magoc and Salzberg 2011) and then combined, screened, and trimmed according to k-mer analysis. The jointed reads were used for digital counts-based comparative analyses and gene assembly.

Ribosomal RNA depletion and cDNA library generation for metatranscriptomic analysis

For cDNA library construction, ∼2.0 μg of total RNA per rumen sample (four high, four low, and two intermediate CH4 yield sheep at two time points) was enriched for mRNA, using the Ribo-Zero rRNA removal kit (Meta-Bacteria, Epicenter Biotechnologies), and mRNA-enriched RNA and total RNA from one untreated sample were fragmented using mRNA fragmentation reagents (Ambion). Double-stranded cDNA (ds cDNA) was synthesized using SuperScript II reverse transcriptase (Invitrogen) using random hexamers (MBI Fermentas). The cDNA sequencing libraries were generated and amplified using the Illumina TruSeq genomic sample prep kit (Illumina) following the manufacturer’s instructions. The amplified libraries were purified and size-selected, and the pooled library was sequenced using the Illumina HiSeq 2000 platform. Full details for the RNA and cDNA manipulations can be found in the Supplemental Material.

Annotation of metagenome and metatranscriptome WGS reads

Artifact-filtered metagenome and metatranscriptome WGS reads were annotated by comparison with the KEGG database (release 58.1, June 1, 2011) (Kanehisa and Goto 2000) using USEARCH 6.0 (Edgar 2010).

The abundance of ssrRNA genes in metagenome data was quantified by aligning the jointed WGS reads to the SILVA (Pruesse et al. 2007), Greengenes (DeSantis et al. 2006), and RDP databases and via BLAST searches of an AgResearch in-house rumen archaea reference database (Supplemental Figs. 3C,D; Supplemental Material).

Statistical analyses

CH4 yield measurements from sheep were analyzed using one-way ANOVA, while individual KEGG genes or gene-encoding ECs differing in gene or transcript abundance were identified using the Wilcoxon rank-sum test in R, with the P-value computed using 10,000 permutation tests. All other statistical comparison tests used the Wilcoxon rank-sum test without permutation. The P-values for multiple tests were corrected using the Benjamini-Hochberg approach (Benjamini and Hochberg 1995).

KEGG pathway enrichment analyses used the Fisher’s exact test based on the differentially enriched genes in gene or transcript levels and known KEGG pathways (n = 293). Hierarchical clustering analyses and heatmaps visualization were performed using the R package.

The reconstruction of mcr/mrt-containing operons of methanogens from metagenome sequence data, the validation of the assemblies, and the phylogenetic analysis of mcrA/mrtA genes are described in the Supplemental Material.

Quantification of gene and transcript levels of the 35 mcrA/mrtA genes

The metagenome and metatranscriptome WGS reads from the eight low CH4 yield, eight high CH4 yield, and four intermediate CH4 yield sheep rumen samples were aligned to mcrA/mrtA genes of the assembled 35 mcr/mrt operons using BWA-based JGI in-house developed gene counting software with a cutoff of 97% identity. The abundance of genes and transcripts were normalized to reads per million (RPM) for further analysis.

Quantitative PCR

To validate methanogen abundance within rumen samples, 16S rRNA gene copies were enumerated by qPCR using methanogen-specific primers and a LightCycler 480 SYBR green I master kit (Roche Applied Science) with real-time amplification on a Rotor-Gene 6000 real-time rotary analyzer (Corbett Life Science) as described previously (Supplemental Material; Jeyanathan et al. 2011).

The qPCR and reverse-transcription qPCR (RT-qPCR) used to verify the abundance and expression of five mcrA/mrtA genes in rumen samples are described in the Supplemental Material.

Data access

Raw sequence data reported in this study have been submitted to the National Center for Biotechnology Information Sequence Read Archive (SRA; http://www.ncbi.nlm.nih.gov/sra/) under accession no. SRA075938. The sequenced samples include run identifications SRR873595–SRR873602 (low CH4 yield metagenome samples); SRR1206671, SRR873604–SRR873610 (high CH4 yield metagenome samples); SRR1138235, SRR873603, SRR1138232, SRR1138234 (intermediate CH4 yield metagenome samples); SRR873450–SRR873457 (low CH4 yield metatranscriptome samples); SRR1206249, SRR873459–SRR873465 (high CH4 yield metatranscriptome samples); and SRR1138694, SRR1138697, SRR1138702, SRR873458 (intermediate CH4 yield metatranscriptome samples). The Whole-Genome Shotgun project has been deposited at DDBJ/EMBL/GenBank under accession no. AUXO00000000. The version described in this article is version AUXO010000000. This transcriptome Shotgun Assembly project has been deposited at DDBJ/EMBL/GenBank under accession nos. GALQ00000000 and GALT00000000 for low and high CH4 yield sheep, respectively. The version described in this article is the first version, GALQ01000000 and GALT01000000. The alpha subunits and assembled operons encoding the methyl coenzyme M reductase have been deposited in GenBank under accession nos. KF214817–KF214824 and KF312304–KF312338.