Received: 2 May 2019 | Revised: 12 July 2019 | Accepted: 13 July 2019 DOI: 10.1111/jbg.12427 ORIGINAL ARTICLE Identification of rumen microbial biomarkers linked to methane emission in Holstein dairy cows Yuliaxis Ramayo‐Caldas1,2 | Laura Zingaretti3 | Milka Popova4 | Jordi Estellé1 | Aurelien Bernard 4 | Nicolas Pons 5 | Pau Bellot 3 | Núria Mach 1 | Andrea Rau 1 | Hugo Roume 5 | Miguel Perez‐Enciso 3 | Philippe Faverdin 6 | Nadège Edouard 6 | Dusko Ehrlich 5 | Diego P. Morgavi 4 | Gilles Renand 1 1 UMR 1313 GABI, INRA, AgroParisTech, Université Paris‐Saclay, Jouy‐en‐Josas, Abstract France Mitigation of greenhouse gas emissions is relevant for reducing the environmental 2 Animal Breeding and Genetics impact of ruminant production. In this study, the rumen microbiome from Holstein Program, IRTA Torre Marimon, Caldes de cows was characterized through a combination of 16S rRNA gene and shotgun Montbui, Spain 3 metagenomic sequencing. Methane production (CH4) and dry matter intake (DMI) Department of Animal Genetics, CRAG, UAB, Bellaterra, Spain were individually measured over 4–6 weeks to calculate the CH4 yield (CH4y = CH4/ 4 VetAgro Sup, UMR 1213 DMI) per cow. We implemented a combination of clustering, multivariate and mixed Herbivores, INRA, Université Clermont model analyses to identify a set of operational taxonomic unit (OTU) jointly associ- Auvergne, Saint‐Genès‐Champanelle, France ated with CH4y and the structure of ruminal microbial communities. Three ruminotype 5 INRA METAGENOPOLIS Unit, Jouy‐en‐ clusters (R1, R2 and R3) were identified, and R2 was associated with higher CH4y. Josas, France The taxonomic composition on R2 had lower abundance of Succinivibrionaceae and 6 UMR 1348 PEGASE, INRA, Agrocampus‐ Methanosphaera, and higher abundance of Ruminococcaceae, Christensenellaceae Ouest, Saint‐Gilles, France and Lachnospiraceae. Metagenomic data confirmed the lower abundance of Correspondence Succinivibrionaceae and Methanosphaera in R2 and identified genera (Fibrobacter Yuliaxis Ramayo‐Caldas and Gilles and unclassified Bacteroidales) not highlighted by metataxonomic analysis. In addi- Renand, UMR 1313 GABI, INRA, AgroParisTech, Université Paris‐Saclay, tion, the functional metagenomic analysis revealed that samples classified in cluster Jouy‐en‐Josas, France. R2 were overrepresented by genes coding for KEGG modules associated with meth- Emails: yuliaxis.ramayo@irta.cat (Y. R.‐C.) anogenesis, including a significant relative abundance of the methyl‐coenzyme M and gilles.renand@inra.fr (G. R.) reductase enzyme. Based on the cluster assignment, we applied a sparse partial least‐ Funding information squares discriminant analysis at the taxonomic and functional levels. In addition, APIS‐GENE we implemented a sPLS regression model using the phenotypic variation of CH4y. By combining these two approaches, we identified 86 discriminant bacterial OTUs, notably including families linked to CH4 emission such as Succinivibrionaceae, Ruminococcaceae, Christensenellaceae, Lachnospiraceae and Rikenellaceae. These selected OTUs explained 24% of the CH4y phenotypic variance, whereas the host ge- nome contribution was ~14%. In summary, we identified rumen microbial biomark- ers associated with the methane production of dairy cows; these biomarkers could be This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2019 The Authors. Journal of Animal Breeding and Genetics published by Blackwell Verlag GmbH J Anim Breed Genet. 2020;137:49–59.  wileyonlinelibrary.com/journal/jbg | 49 50 |    RAMAYO‐CALDAS et al. used for targeted methane‐reduction selection programmes in the dairy cattle industry provided they are heritable. KEYWORDS metagenomics, metataxonomics, methane emission, microbial biomarker 1 | IN TRO D U C T ION 2 | M ATERIAL S AND M ETHOD S Cattle have the remarkable ability to digest and transform 2.1 | Phenotype and host genotype details non‐edible plant cell wall components into high‐quality pro- The experiment was carried out at the INRA experimental teins for human consumption. The digestion of feeds entails, farm in Méjusseaume (Le Rheu, France). Management of however, the production of the greenhouse gas CH4. The experimental animals followed the guidelines for animal possibility to mitigate CH4 emissions while simultaneously research of the French Ministry of Agriculture and other improving feed efficiency is highly relevant for the sustain- applicable guidelines and regulations for animal experi- ability of cattle production systems. Notwithstanding, the mentation in the European Union (European Commission, mechanisms to achieve it are not yet fully understood (Flay 2010). Approval number for ethical evaluation was et al., 2019). The contribution of the gastrointestinal micro- APAFIS:3122‐2015112718172611. Sixty‐five loose‐housed biota to feed digestion and enteric CH4 production is well lactating Holstein cows were used in this experiment. They established in ruminants (Delgado et al., 2019; Difford et were allocated to three pens equipped with individual troughs al., 2018; Huws et al., 2018; Ross, Moate, Marett, Cocks, and automatic gates detected by radio‐frequency identifica- & Hayes, 2013; Tapio, Snelling, Strozzi, & Wallace, 2017). tion tags attached on the cow ears. Each pen was of simi- However, enteric CH4 production is a complex trait de- lar size (n = 21–23). Parity was equilibrated between pens, termined not only by the rumen microbiome (Huws et al., with two‐thirds of the cows that were in their first lactation 2018; Ross et al., 2013; Tapio et al., 2017), but also by host (n = 42). Cows received the same total mixed ration (TMR) genetics (Difford et al., 2018; Ross et al., 2013), and envi- throughout the experiment consisting of maize silage (65%), ronmental factors (Gerber et al., 2013; Martin, Morgavi, soybean cake (18%) and energy concentrate (8%), composed & Doreau, 2009; McAllister, Cheng, Okine, & Mathison, of corn, wheat, barley, dehydrated beet pulp, dehydrated al- 1996). Exploring the abundance and composition of mi- falfa (8%) and minerals (Table S1). The TMR was offered ad crobial communities in the gastrointestinal tract of cattle in libitum and individually weighed every morning. Each morn- relation to the host genome is of great interest for quantify- ing, refusals from the previous day were weighed. Samples ing animal variability in feed digestibility and enteric CH4 of forages, concentrates and refusals were analysed for dry emission (Huws et al., 2018; Leahy et al., 2013; Li et al., matter content. Dry matter intake was calculated daily as the 2019). Several studies in ruminants have explored the re- difference between offered and refused dry matter weights. lationship of the rumen microbiota composition with feed Methane production was measured with two GreenFeed efficiency and CH4 emission (Delgado et al., 2019; Jami, emission monitoring (GEM) systems (C‐Lock Inc.). This sys- White, & Mizrahi, 2014; Jewell, McCormick, Odt, Weimer, tem automatically measures CH4 when animals visit a con- & Suen, 2015; Li et al., 2019; Myer, Smith, Wells, Kuehn, & centrate feeder equipped with a head hood and an extractor Freetly, 2015; Wallace et al., 2015). Pioneering rumen‐engi- fan for the capture of breath and eructation gases. The animals neering studies have suggested that microbial communities are attracted and kept attracted to the feeder with pellets that are highly resilient and host‐specific (Cole, 1991; Weimer, are distributed in small quantities. In the current experiment, 2015). Although these latter properties make it difficult to the cows were allowed to visit the GEM a maximum of four manipulate the ruminal microbial community, they also en- times per day, with each visit separated by at least 6 hr. At able the analysis of the covariation of these ecosystems with each visit, cows received six drops of concentrate separated host performance and jointly selecting both host genome and by 30 s for at least 3 min. At each visit, CH4 production rate is microbiome variants. calculated, combining the gas concentration (measured with The main goal of this study was to implement an in- a non‐dispersive infrared analyser) to the airflow in the pipe tegrative approach using 16S rRNA and shotgun metage- (measured with a flow meter). An algorithm developed and nomic sequencing data to identify microbial biomarkers applied by C‐Lock Inc. calculates the CH4 rate at each visit if linked to CH4y emission. In addition, we also explored the head of the animal is correctly positioned (controlled by the role of host genetics on the determinism of this a laser beam) for at least 2 min. With two GEM systems and phenotype. RAMAYO‐CALDAS et al.    | 51 three batches, the following experimental design was applied: during a first period (6 weeks in January–February 2017), 2.3 | Bioinformatics and statistical analysis cows in batches B1 (n = 21) and B2 (n = 21) were measured Sequences corresponding to the 16S rRNA gene data were and rumen liquid sampled. During a second period (the next analysed on an in‐house Galaxy‐based graphical user inter- 4 weeks in March 2017), cows in batch B1 were measured face for IM TORNADO (Jeraldo et al., 2014) and mothur again and sampled together with cows in batch B3 (n = 23). (Schloss et al., 2009) for bacteria and archaea, respectively. The daily dry matter intake (TMR and concentrate in the The workflow included a quality control step to remove se- GEM system) and the per‐visit CH4 emission rate measures quences with Phred scores of <33 and trimmed sequences were averaged over the testing periods to estimate the indi- based on expected amplicon length, as well as merge vidual dry matter intake (DMI), methane emission rate (CH4) paired reads, remove chimera and select OTUs (97% iden- and methane yield (CH4y = CH4/DMI). The traits were ad- tity). Finally, after removing doubleton OTUs, only those justed beforehand for the contemporary group mean in a sim- OTUs representing more than 0.001% of the total were re- ple linear model including the batch × period effect: tained. Taxonomic classification was based on the SILVA v123 database (Quast et al., 2012) for bacteria and RIMDB y = Xb + e (Seedorf, Kittelmann, Henderson, & Janssen, 2014) for ar- chaea. As previously mentioned, animals in the study were where y is the trait value vector, b the fixed effect vector of divided into groups due to experimental constraints (see the contemporary group (batch × period) effect and e the re- Phenotype2.1 section). Therefore, to estimate the stability sidual phenotype. of ruminal bacterial communities, we used the 16S rRNA Rumen fluid (~400 ml) was sampled via oesophageal gene data of 21 cows in batch B1 that were sampled twice tubing in the morning before feeding (last week of February in the two successive periods. In a first step, we estimated for period 1 and last week of March for period 2). The av- and contrasted diversity metrics among time points (T1 erage days in milk (DIM) in period 1 were 148.5 (SD 12.5) vs. T2) with vegan r package (Jari Oksanen et al., 2018). and 182.8 (SD 18.6) in period 2. Samples were filtered Alpha‐diversity was evaluated with the Shannon index through a polyester monofilament fabric (250 μm mesh ap- (Shannon, 1984), and beta‐diversity was assessed using the erture), and 2 ml of the filtrate was centrifuged at 20,000 g, Whittaker index (Whittaker, 1972). Subsequently, stability 20 min, 4°C. The supernatant was discarded, and the pellet between sampling points was estimated using the RV coef- was snap‐frozen in liquid nitrogen and stored at −80°C. ficient on the two OTU normalized abundance tables. The The Illumina BovineSNP50 v.2 BeadChip (Illumina Inc.) RV coefficient was calculated between times points as the was used to genotype the 65 Holstein cows. Quality control total co‐inertia (sum of eigenvalues of the product of two was performed to exclude single nucleotide polymorphisms cross‐product matrices) divided by the square root of the (SNPs) with minor allele frequencies (MAF) <5%, rates of product of the squared total inertia. A zero RV score indi- missing genotypes above 10%, as well as those that did not cates no similarity, whereas the RV score approaches 1 for conform to Hardy–Weinberg expectations (threshold set at a increasing co‐structure between two data sets. p‐value of .001). We also excluded SNPs that did not map The whole‐metagenome data were processed as fol- to the bovine reference genome (ARS‐UCD1.2 assembly) or lows. Gene abundance profiling was performed using the that were located on the X‐chromosome. 16.6 million gene integrated reference catalogue of the rumen microbiome (Junhua et al., 2019). First, low‐quality 2.2 | Rumen microbial DNA extraction, and host contaminant reads (from Bos taurus genome ARS‐ PCR amplification and sequencing UCD1.2) were removed using AlienTrimmer (Schubert, DNA from rumen liquid fraction was extracted with an estab- Lindgreen, & Orlando, 2016) and Bowtie2 (Langmead & lished protocol (Yu & Morrison, 2004). Extracted DNA was Salzberg, 2012), respectively. For a more detailed infor- sent to the University of Illinois Keck Center for Fluidigm mation which includes the creation of gene abundance and sample preparation and Illumina sequencing. Primers target- KEGG orthologous (KO) tables, as well as the assembly of ing the V3–V5 region (F357 and R926) were used to amplify metagenomic species (MGS) clusters, see the supplemen- a region of 570 base pairs of the bacterial 16S rRNA gene. tary material (Supporting Information). Archaea‐specific primers (349F and 806R) were used to am- plify a 457‐base‐pair 16S rRNA gene fragment. The amplicons were sequenced on one MiSeq flow cell for 251 cycles. The 2.3.1 | Structure of the ruminal ecosystem whole‐metagenome shotgun sequence of 30 samples collected To infer the structure of ruminal bacterial communities, rumi- from rumen of Holstein cows with low and high CH4y emission notype cluster detection was done using the genera abundance and distributed across the three batches were generated using a (101 genera based on the bacterial data) in each sample, as quantitative metagenomic pipeline (Supporting Information). previously described for human gut enterotypes (Arumugam 52 |    RAMAYO‐CALDAS et al. TABLE 1 Mean and standard deviation of production and reliability corresponding to the sPLS‐DA model was as- methane emission traits sessed as function of the prediction maximum distance Trait Unit Mean SD CV (%) between overall misclassification error rate and balanced error rate (BER) after fivefold cross‐validation repeated 500 Live weight kg 634 49 8 times. In addition, sPLS was also performed in regression Dry matter intake kg/day 21.2 2.2 10 model to identify OTUs associated with CH4y phenotypic (DMI) variation. Subsequently, we implemented a conservative Milk production kg/day 31.1 4.8 15 approach and retained only those OTUs found in common Milk efficiency kg/kg DMI 1.47 0.16 11 between the two approaches for downstream analysis. The GreenFeed visits n/day 2.50 0.53 21 sPLS‐DA approach using ruminotype‐like sample classifi- Visit duration s 224 16 7 cation was also implemented on the metagenomic data at Methane emission g/day 506 56 11 MGS and functional (KEGG) levels. rate Methane yield g/kg DMI 24.1 3.1 13 2.3.3 | Mixed model To estimate the proportion of CH4 and the CH4y pheno- et al., 2011). Briefly, sample clusters were detected using the typic variance explained by the host (heritability) as well as probability distribution distance metric related to Jensen– by the bacterial community (microbiability), the following Shannon divergence and partitioning around medoids. The Bayesian mixed model was implemented using the bglr r optimal number of clusters was determined following the package (Pérez & de los Campos, 2014): Calinski–Harabasz (CH) index (Caliński & Harabasz, 1974), and the statistical consistency of the corresponding partition y = 1𝝁 + g + b + 𝜺 was evaluated using the Silhouette coefficient (Rousseeuw, where y is the ( phenotype vector (i.e., CH4y), 1µ is the inter- 1987). Furthermore, sample stability within each cluster was ) cept, g ∼ N 0,G𝝈 g and b ∼ N 0,B𝝈 2b , G is the genomic 2 ( ) estimated through a cluster‐wise Jaccard bootstrap analysis (100 repetitions; Hennig, 2007). The association between relationship matrix (GRM) based on 38,872 autosomal SNPs, predicted ruminotype clusters with CH4y was obtained using and 𝝈 2g represents the additive genetic variance. B represents a least‐squares analysis as implemented in the lsmeans r the microbial relationship matrix, calculated based on the package. To identify genera with significantly different Bray–Curtis dissimilarities distance matrix, and ε is the re- abundance among the predicted ruminotype cluster groups, sidual term. A second model was employed using the method a differential abundance (DA) analysis was performed using proposed by Ross et al. (2013) which build the microbial re- the zero‐inflated Gaussian mixture model implemented in the lationship matrix based on the variance–covariance matrix metagenomeseq r package (Paulson, Stine, Bravo, & Pop, from the log‐transformed and standardized OTU table. In 2013), using a threshold for adjusted p‐values of 5%. In ad- both cases, the models were run using a Gibbs sampler with dition, a presence–absence (PA) test was also performed to 30,000 iterations and a burn‐in of 2,000 rounds; we used identify genera that were unique for each of the identified standard flat priors for the intercept. clusters. 3 | RESULTS AND DISCUSSIO N 2.3.2 | Multivariate analysis To identify a set of OTUs jointly associated with CH4y 3.1 | Phenotype description phenotypic variation and the structure of ruminal bacte- The descriptive statistics of production and CH4 traits are rial communities, a combination of multivariate analyses presented in Table 1. Methane emissions were calculated was performed using sparse partial least‐squares (sPLS) as as the average value of all visits throughout the study. With implemented in the mixomics r package (Rohart, Gautier, 2.50 visits per day, emission rates per each cow were cal- Singh, & Lê Cao, 2017). sPLS is a statistical approach em- culated for an average of 105.1 and 69.7 visits for periods ployed to identify a small subset of variables that maximize 1 and 2, respectively. This large number of visits to GEM the covariance between two different data sets (for instance, system ensured a precise measurement of the animal phe- a table of centred log ratio‐transformed OTUs and rumino- notype, as it largely exceeds the recommended 20–30 mini- type clusters or CH4y values; Lê Cao, González, & Déjean, mum visits (Barchia et al., 2017; Manafiazar, Zimmerman, 2009). In a first step, sPLS discriminant analysis (sPLS‐ & Basarab, 2016). The variability among cows for methane DA) was applied based on sample classification accord- production (g/day; CV 11%) was higher than previously ing to ruminotype cluster assignment. The classification reported (Ricci, Rooke, Nevison, & Waterhouse, 2013). RAMAYO‐CALDAS et al.    | 53 From period 1 to period 2, the repeatability of the methane Cows that clustered within R2 emitted more CH4y (p < .05) emission traits was 0.82 and 0.73 for CH4 and CH4y, re- than those clustered in R1 or R3, and no significant differ- spectively. These values were lower than the repeatability ence in CH4y emission was observed between the R1 and R3 of production traits (0.95–0.97 for DMI, milk production clusters (Figure 1c). and milk efficiency). Methane production rate was moder- We evaluated the taxonomic composition of rumino- ately correlated with DMI (Pearson's correlation r = .44) type clusters through a combination of presence–absence and milk production (r = .28) and was poorly correlated (PA) and differential abundance (DA) analysis. According with body weight (r = .16). to the PA test, the only genera showing significant differ- ences between clusters was Succinivibrionaceae_UCG‐001, which was absent in most cows classified as R2 (high CH4y 3.2 | Link between ruminotype emission) compared with R1 and R3. Only four cows out structure of ruminal bacterial communities and of the 16 in R2 had Succinivibrionaceae_UCG‐001, al- CH4y emissions beit with a lower abundance compared to the other clusters The 16S rRNA gene sequences from 65 rumen liquid frac- (Table 2). The DA analysis confirmed the PA test results, as tions were analysed. After quality control, bacterial 1,198 Succinivibrionaceae_UCG‐001 showed significantly differ- OTUs and 1,764 archaeal OTUs were identified (Supporting ent abundance among clusters (Table S2). We note that larger Information). The relative abundance of bacterial genera in differences in Succinivibrionaceae_UCG‐001 abundance, as each sample was used for cluster detection as described for well as a greater number of DA genera, were observed be- human gut enterotypes (Arumugam et al., 2011). This method tween samples classified as R2 and R3 (Table S2). However, first performs a sample stratification, followed by the identi- no significant difference in Succinivibrionaceae_UCG‐001 fication of the optimal number of clusters and the statistical relative abundance was observed between R1 and R3 (Table consistency of the predicted partition. Cows clustered into S2). These observations are consistent with previous reports three ruminotype clusters (R1, R2 and R3; Figure 1a,b); 30 describing members of the Succinivibrionaceae family as- cows were assigned to R1, 16 to R2 and 19 to R3. In concord- sociated with low CH4 (Danielsson et al., 2017; Kittelmann ance with our results, three ruminotype clusters have previ- et al., 2014; Wallace et al., 2015). Furthermore, this is in ously been reported in sheep (Kittelmann et al., 2014), while agreement with the higher CH4y emissions observed in in dairy cattle (based on a principal coordinate analysis at the R2. Members of families associated with high CH4 emis- OTU level), two clusters were reported by Danielsson et al. sion, such as Ruminococcaceae, Christensenellaceae and (2017). Therefore, we tested sample stability by comparing Lachnospiraceae (Kittelmann et al., 2014; Tapio et al., 2017), two (k = 2) or three (k = 3) putative clusters. Cluster‐wise exhibited DA patterns with higher abundant in R2 compared Jaccard bootstrap analysis revealed moderate sample stabil- to R1 and R3 (Table S2). In a similar way, we explored the ity, showing a more stable solution for three (K3: 0.63, 0.64 archaeal taxonomic composition of the ruminotype clusters. and 0.63) clusters compared with two (K2: 0.61 and 0.4). According to a DA analysis, significant differences were only In concordance with the aforementioned studies in sheep detected at the OTU level. Similar to the aforementioned pat- (Kittelmann et al., 2014) and dairy cattle (Danielsson et al., terns at the bacterial level, the largest number of DA OTUs 2017), we observed a significant association between rumi- was observed between samples classified as R2 and R3 (Table notype cluster assignments and CH4y emission (Figure 1c). S3). Of noted, OTU members of Methanosphaera showed (a) (b) F I G U R E 1 Structure of the ruminal bacterial community of lactating Holstein (c) cows displaying natural differences in methane emission. (a) Sample distribution. (b) Calinski–Harabasz (CH) indexes of potential numbers of clusters. (c) Pairwise comparisons of CH4y emission between ruminotype‐like clusters R1 (red), R2 (blue) and R3 (green) [Colour figure can be viewed at wileyonlinelibrary.com] 54 |    RAMAYO‐CALDAS et al. TABLE 2 Results from the presence– Adj p absence test between ruminotype‐like values clusters Comparison Genus Odds ratio p Values (BH) R1_R2 Succinivibrionaceae_UCG‐001 17.72 5.27E−05 .005 R2_R3 Succinivibrionaceae_UCG‐001 0.022 4.22E−05 .003 F I G U R E 2 Analytical framework employed to identify microbial biomarkers of methane emission in lactating Holstein cows [Colour figure can be viewed at wileyonlinelibrary.com] consistently lower abundance in the R2 cluster (higher CH4 cannot be neglected. Therefore, to conservatively focus on emitters) and were only detected when R2 was compared with the primary OTU markers associated with CH4y emissions, R1 or R3. Similar patterns were observed for OTUs classi- we combined the cluster analysis results with those of a pre- fied as Methanobrevibacter ruminantium, which showed dictive model. We applied a two‐pronged strategy using mul- lower abundance in R2 compared with R3. Meanwhile, no tivariate analysis, including a supervised sPLS‐DA (based on clear patterns were observed for OTU members of the genera sample cluster classification) and a sPLS regression model Methanobrevibacter (gottschalkii clade); some OTUs were that considers the joint covariation of OTU relative abun- most abundant in R2, while others showed higher abundance dances and CH4y emissions (Figure 2). For the sPLS‐DA, the in R1 or R3 (Table S3). first and second principal component combined the relative abundances of 231 OTUs (PC1 = 200 and PC2 = 31) and allowed a clear discrimination between samples classified as 3.3 | Identification of ruminal OTUs linked R2 and R3 (Figure S2 and Table S4). The area under the ROC to CH4y emission and associated with the curve corresponding to both R2 (0.97) and R3 (0.92) showed structure of the rumen microbiota high values, suggesting a good ability of the model to cor- In this study, we used an integrative approach combining rectly classify these samples. In addition, the sPLS regression multivariate and clustering analyses to identify microbial model allowed us to identify features that maximize the covar- biomarkers linked to CH4y emission and the structure of the iance between OTU relative abundance and CH4y phenotype rumen microbiota (Figure 2). We are aware that in compari- variation. After tuning parameters (fivefold cross‐validation son with a predictive model (which would take into account repeated 500 times), a single component was identified as the total microbial variation), ruminotype cluster approaches optimal (PC1 Q2h = 0.12), yielding a final selection of 200 may have some limitations for biomarker identification OTUs. We then applied a conservative approach, retaining (Costea et al., 2018; Knights et al., 2014). However, as pre- the 86 OTUs that were common to both sPLS analyses (Table viously reported (Danielsson et al., 2017; Kittelmann et al., S4). We note that 75.5% (65/86) of them were also identified 2014) and confirmed in this work, the link between the struc- as DA among ruminotype clusters, where selected OTUs dis- ture of the bacterial ruminal ecosystem and CH4y emission played consistent differential abundance patterns (Table S4). RAMAYO‐CALDAS et al.    | 55 The taxonomic classification of these OTUs included fami- ruminotype, a finding that agrees with the results of the 16S lies reportedly linked to CH4 such as Succinivibrionaceae, rRNA gene data. Members of the Succinivibrionaceae have Ruminococcaceae, Christensenellaceae, Lachnospiraceae, been previously proposed to be responsible for the lower Gastranaerophilales, Rikenellaceae, Bacteroidales_BS11 and CH4 emissions in wallaby microbiota (Pope et al., 2011). Prevotellaceae (Danielsson et al., 2017; Difford et al., 2018; Moreover Succinivibrionaceae was recently identify among Kamke et al., 2016; Kittelmann et al., 2014). the core of heritable rumen bacteria (Wallace et al., 2019) as well as associated with feed efficiency and levels of pro- pionate in beef cattle (Li et al., 2019). Succinivibrionaceae 3.4 | Rumen microbiome linked to CH4y produce succinate as a main fermentation product (O'Herrin identified through metagenomics & Kenealy, 1993; Santos & Thompson, 2014), which is To gain insight into a more precise taxonomic composition converted to propionate by other members of the micro- and functional processes that differed between ruminotypes, biota and thus less hydrogen might be available for meth- we generated whole‐metagenome data from 30 samples dis- anogens. Co‐exclusion between Succinivibrionaceae and tributed across the two periods and the three ruminotype Methanobrevibacter has been reported in cattle (McCabe et clusters. After removing low‐quality and host contaminant al., 2015). The use of metagenome data allows identifying reads (Supporting Information), an average of 21.3 million other features not revealed by 16S rRNA gene sequences reads per sample was produced. Reads were mapped against including unclassified Bacteroidales and Fibrobacteres the rumen gene reference catalogue (Junhua et al., 2019). A (Table S6) as well as to confirm the patterns of methano- total of 608 MGS were identified by MSPminer (Cervino et genic archaea members of the family Methanobacteriaceae al., 2019) and 264 KEGG functional modules were detected observed using the metataxomic approach. This is in agree- by consolidating the gene abundance matrix into KEGG an- ment with the previously reported negative correlation be- notations (Supporting Information). In a first step, the overall tween Methanosphaera spp. and CH4 production in sheep structure of the ruminal microbial community was explored (Kittelmann et al., 2014). The only MGS taxonomy clas- using non‐metric multidimensional scaling (nMDS) with sified as Methanosphaera (MGS512) showed lower abun- Bray–Curtis and binary Jaccard dissimilarities inferred from dance in samples classified as R2 compared with R1 and the metagenomic gene abundance matrix (Figure S3). nMDS R3 (fold changeR1,R2 = 5.92, corrected p‐value = .03; fold plots showed a clear relationship between ruminotypes and changeR3,R2 = 8.56, corrected p‐value = .008). nMDS patterns for both dissimilarities. After correction for A functional approach using KEGG modules shows a group and period effects, PERMANOVA tests of signifi- classification of samples according to ruminotype cluster cance also showed a significant association between rumi- assignation. The sPLS‐DA analysis revealed 55 modules; notype classification and Bray–Curtis (p < .05) and binary however, the model showed a slightly lower accuracy than Jaccard (p < .01) dissimilarity matrices. The same approach for OTUs and MGS analyses (Table S7 and Figure S5). was followed to explore the relationship between rumen Among the discriminant KEGG modules, there were func- ecosystem and CH4y, with nMDS plots showing a tendency tional modules related to CH4y emission, such as methano- of samples to group according to CH4y values. In this case, genesis (M00567), C4‐dicarboxylic acid cycle, NAD‐malic the PERMANOVA test corrected by group and period was enzyme type (M00171, M00172) and acetyl‐CoA pathway not significant for Bray–Curtis or binary Jaccard (p = .06) (M00422), but also with general functions such as glycol- dissimilarities. ysis (M00001), gluconeogenesis (M00003), formaldehyde The MGS and functional profiles of each sample were assimilation (M00345) and crassulacean acid metabolism explored through sPLS‐DA. Based on MGS, the relative (M00169; Table S7). A detailed examination of KOs in- abundance of 150 MGS enabled the discrimination between volved in the methanogenesis pathway (M00567) revealed samples (Table S5 and Figure S4). In agreement with the tax- that six of them were significantly more abundant in sam- onomic classification of the 86 selected OTUs, members of ples classified as R2 (Figure S6). Interesting, five of them the Succinivibrionaceae, Lachnospiraceae, Prevotellaceae, (K00205, K00672, K00399, K00580 and K03389) includ- Rikenellaceae and Ruminococcaceae families were among ing the methyl‐coenzyme M reductase (which catalyses the the top discriminant MGS. Moreover, 59.3% (89/150) of the rate‐limiting CH4 synthesis (Scheller, Goenrich, Thauer, & MGS were also identified as DA among ruminotype (Table Jaun, 2013; Wongnate et al., 2016) have been suggested as S6). Similarly, to the aforementioned results at the OTU biomarkers for CH4 production across diverse cattle breeds level, a larger number of DA MGS were observed between (Auffret et al., 2018). Overall, our results show that shotgun samples classified as R2 and R3 (Table S6). Interestingly, metagenomics is able to provide additional insights into the shotgun metagenomics shows a lower abundance of a differences observed with 16S rRNA gene data sets, even MGS affiliated to Succinivibrionaceae (MGS116, fold with smaller sample sizes. Interestingly, specific KEGG changeR3,R2 = 1.87, corrected p‐value = .043) in the R2 pathways and enzymes appear to be associated with CH4 and 56 |    RAMAYO‐CALDAS et al. TABLE 3 Estimated heritability (h2) and microbiability (m2) of methane production (CH4) and methane yield (CH4y) of lactating Holstein cows 86 selected OTUs (Bray–Curtis 86 selected OTUs (log‐transformed Whole OTU table (1,198) distance) and standardized OTU table)a Trait h2 (SD) m2 (SD) h2 (SD) m2 (SD) h2 (SD) m2 (SD) CH4 0.144 (0.09) 0.164 (0.10) 0.141 (0.09) 0.192 (0.11) 0.157 (0.09) 0.130 (0.06) CH4y 0.148 (0.10) 0.181 (0.11) 0.143 (0.09) 0.242 (0.14) 0.163 (0.09) 0.174 (0.08) a Method proposed by Ross et al., 2013. Estimated a microbial relationship matrix based on the variance–covariance matrix from the log‐transformed and standardized OTU table. had a moderate discriminant ability of whole‐metagenome ruminal bacterial community structure. OTUs associated data using the sPLS‐DA approach. with CH4y emission were predominantly hydrogen‐produc- ing bacteria and explained up to 24% of the CH4y pheno- typic variance, whereas the host genome contribution was 3.5 | CH4 emission heritability and around 14%. Some discriminant bacterial OTUs identified microbiability by metataxonomic were confirmed by whole metagenome. A mixed model was implemented to estimate the proportion In particular, samples clustered in R2 (high CH4y emis- of CH4 and CH4y phenotypic variance explained by the host sion) showed a lower abundance of Succinivibrionaceae cow genome, the whole OTU matrix (n = 1,198) and the 86 and Methanosphaera spp. as well as a higher abundance of selected OTUs (Figure 2). According to the model based on genes coding for functional modules and enzymes involved the whole OTU matrix, the host genome explained ~14% in methanogenesis. Overall, we report a set of microbial bio- of the CH4 and CH4y, whereas the bacterial community ex- markers that have the potential to be employed for charac- plained around ~16% of CH4 production and 18% of CH4y terizing high‐emitting cattle for targeted management in the (Table 3). The estimated CH4 heritability and microbiability dairy cattle industry. on the whole OTU model were consistent with recent results in dairy cattle (Difford et al., 2018). Meanwhile, in agree- ACKNOWLEDGEMENTS ment with previous studies (Saborío‐Montero, 2018), differ- ences between the estimated microbiability were observed Y. Ramayo‐Caldas was funded by the European Union depending on the method employed to build the microbial H2020 Research and Innovation programme under Marie relationship matrix (Table 3). We acknowledge that the Skłodowska‐Curie grant (P‐Sphere) agreement no 6655919. number of samples in our study may not be large enough to L. Zingareti is recipient of a FPI grant to achieve the PhD accurately estimate parameters, which limit the robustness research from Ministry of Economy and Science (MINECO, of the model. In spite of this limitation, consistent values of Spain), associated with the “Centro de Excelencia Severo CH4 heritability between models were observed (Table 3). Ochoa 2016–2019” award SEV‐2015‐0533 to CRAG. These values are also comparable to those of Difford et al. The authors warmly thank all technical staff at the INRA (2018). It is noted that heritability is similar for the two CH4 Méjusseaume farm for providing high‐quality measure- emission traits but microbiability was always higher when ments. This work was funded by APIS‐GENE and used using CH4y. Our results also suggest a slight improvement production data obtained in the Deffilait project cofounded of CH4y variance component estimation (up to ~24%) after by APIS‐GENE and the Agence Nationale de la Recherche OTU preselection, but no for CH4 alone (Table 3). Further (ANR‐15‐CE20‐0014). studies with larger sample sizes and standardized analyti- cal pipelines could provide more reliable estimates of the CONFLICT OF INTEREST microbiota contribution to this complex CH4 emission trait. The authors declare that there is no conflict of interest. 4 | CO NC LU S ION DATA AVAILABILIT Y STATEMENT Our results confirm the link between the structure of the ru- The data that support the findings of this study are available minal bacterial community and CH4 emission. We identified from APIS‐GENE. Restrictions apply to the availability of 86 OTUs simultaneously linked to CH4y emission and the these data, which were used under license for this study. Data RAMAYO‐CALDAS et al.    | 57 are available from the corresponding author with the permis- Flay, H. E., Kuhn‐Sherlock, B., Macdonald, K. A., Camara, M., sion of APIS‐GENE. Lopez‐Villalobos, N., Donaghy, D. J., & Roche, J. R. (2019). Hot topic: Selecting cattle for low residual feed intake did not affect daily methane production but increased methane yield. Journal ORCID of Dairy Science, 102(3), 2708–2713. https​ ://doi.org/10.3168/ jds.2018-15234​ Yuliaxis Ramayo‐Caldas https://orcid. Gerber, P. J., Hristov, A. N., Henderson, B., Makkar, H., Oh, J., Lee, org/0000-0002-8142-0159 C., … Oosting, S. (2013). Technical options for the mitigation of direct methane and nitrous oxide emissions from livestock: A re- Núria Mach https://orcid.org/0000-0002-8001-6314 view. Animal, 7(s2), 220–234. https​://doi.org/10.1017/S1751​73111​ Miguel Perez‐Enciso https://orcid. 3000876 org/0000-0003-3524-995X Hennig, C. (2007). Cluster‐wise assessment of cluster stability. Computational Statistics & Data Analysis, 52(1), 258–271. https​:// doi.org/10.1016/j.csda.2006.11.025 R E F E R E NC E S Huws, S. A., Creevey, C. J., Oyama, L. B., Mizrahi, I., Denman, S. E., Popova, M., … Morgavi, D. P. (2018). Addressing global ruminant Arumugam, M., Raes, J., Pelletier, E., Le Paslier, D., Yamada, T., agricultural challenges through understanding the rumen microbi- Mende, D. R., … Bork, P. (2011). Enterotypes of the human gut ome: Past, present, and future. Frontiers in Microbiology, 9, 2161. microbiome. Nature, 473, 174. https​://doi.org/10.1038/natur​e09944 https​://doi.org/10.3389/fmicb.2018.02161​ Auffret, M. D., Stewart, R., Dewhurst, R. J., Duthie, C.‐A., Rooke, J. Jami, E., White, B. A., & Mizrahi, I. (2014). Potential role of the bo- A., Wallace, R. J., … Roehe, R. (2018). Identification, comparison, vine rumen microbiome in modulating milk composition and feed and validation of robust rumen microbial biomarkers for methane efficiency. PLoS ONE, 9(1), e85423. https​://doi.org/10.1371/journ​ emissions using diverse bos taurus breeds and basal diets. Frontiers al.pone.0085423 in Microbiology, 8, 2642. Jari Oksanen, F., Blanchet, G., Friendly, M., Kindt, R., Legendre, P., Barchia, I. M., Arthur, P. F., Weber, C., Herd, R. M., Donoghue, K. A., McGlinn, D., … Wagner, H. (2018). vegan: Community ecology Bird‐Gardiner, T., & Hegarty, R. S. (2017). Optimizing test proce- package. R package version 2.5‐3. Retrieved from https​://CRANR- dures for estimating daily methane and carbon dioxide emissions in proje​ctorg/​packa​ge=vegan​ cattle using short‐term breath measures. Journal of Animal Science, Jeraldo, P., Kalari, K., Chen, X., Bhavsar, J., Mangalam, A., White, B., 95(2), 645–656. … Chia, N. (2014). IM‐TORNADO: A tool for comparison of 16S Caliński, T., & Harabasz, J. (1974). A dendrite method for cluster anal- reads from paired‐end libraries. PLoS ONE, 9(12), e114804. https​:// ysis. Communications in Statistics, 3(1), 1–27. doi.org/10.1371/journ​al.pone.0114804 Cervino, A. C. L., Pichaud, M., Plaza Oñate, F., Le Chatelier, E., Jewell, K. A., McCormick, C. A., Odt, C. L., Weimer, P. J., & Suen, G. Gauthier, F., Almeida, M., … Magoulès, F. (2019). MSPminer: (2015). Ruminal Bacterial community composition in dairy cows is abundance‐based reconstitution of microbial pan‐genomes from dynamic over the course of two lactations and correlates with feed shotgun metagenomic data. Bioinformatics, 35(9), 1544–1552. https​ efficiency. Applied and Environmental Microbiology, 81(14), 4697. ://doi.org/10.1093/bioin​forma​tics/bty830 https​://doi.org/10.1128/AEM.00720-15 Cole, N. A. (1991). Effects of animal‐to‐animal exchange of ruminal Junhua, L., Zhong, H., Ramayo‐Caldas, Y., Terrapon, N., Lombard, V., contents on the feed intake and ruminal characteristics of fed and Potocki‐Veronese, G., … Morgavi, D. P. (2019). A catalog of micro- fasted lambs. Journal of Animal Science, 69(4), 1795–1803. https​:// bial genes from the bovine rumen unveils a specialized and diverse doi.org/10.2527/1991.6941795x biomass‐degrading environment. bioRxiv:272690. Costea, P. I., Hildebrand, F., Arumugam, M., Bäckhed, F., Blaser, M. J., Kamke, J., Kittelmann, S., Soni, P., Li, Y., Tavendale, M., Ganesh, S., … Bushman, F. D., … Bork, P. (2018). Enterotypes in the landscape of Attwood, G. T. (2016). Rumen metagenome and metatranscriptome gut microbial community composition. Nature Microbiology, 3(1), analyses of low methane yield sheep reveals a Sharpea‐enriched 8–16. https​://doi.org/10.1038/s41564-017-0072-8 microbiome characterised by lactic acid formation and utilisation. Danielsson, R., Dicksved, J., Sun, L., Gonda, H., Müller, B., Schnürer, Microbiome, 4(1), 56. https​://doi.org/10.1186/s40168-016-0201-2 A., & Bertilsson, J. (2017). Methane production in dairy cows cor- Kittelmann, S., Pinares‐Patiño, C. S., Seedorf, H., Kirk, M. R., Ganesh, relates with rumen methanogenic and bacterial community struc- S., McEwan, J. C., & Janssen, P. H. (2014). Two different bacterial ture. Frontiers in Microbiology, 8, 226. https​ ://doi.org/10.3389/ community types are linked with the low‐methane emission trait fmicb.2017.00226​ in sheep. PLoS ONE, 9(7), e103171. https​://doi.org/10.1371/journ​ Delgado, B., Bach, A., Guasch, I., González, C., Elcoso, G., Pryce, J. E., al.pone.0103171 & Gonzalez‐Recio, O. (2019). Whole rumen metagenome sequenc- Knights, D., Ward, T. L., McKinlay, C. E., Miller, H., Gonzalez, A., ing allows classifying and predicting feed efficiency and intake McDonald, D., & Knight, R. (2014). Rethinking "enterotypes". levels in cattle. Scientific Reports, 9(1), 11. https​://doi.org/10.1038/ Cell Host & Microbe, 16(4), 433–437. https​://doi.org/10.1016/j. s41598-018-36673-w chom.2014.09.013 Difford, G. F., Plichta, D. R., Løvendahl, P., Lassen, J., Noel, S. J., Langmead, B., & Salzberg, S. L. (2012). Fast gapped‐read alignment Højberg, O., … Sahana, G. (2018). Host genetics and the rumen with Bowtie 2. Nature Methods, 9, 357. https​://doi.org/10.1038/ microbiome jointly associate with methane emissions in dairy cows. nmeth.1923 PLOS Genetics, 14(10), e1007580. https​ ://doi.org/10.1371/journ​ Lê Cao, K.‐A., González, I., & Déjean, S. (2009). integrOmics: An al.pgen.1007580 R package to unravel relationships between two omics datasets. 58 |    RAMAYO‐CALDAS et al. Bioinformatics, 25(21), 2855–2856. https​ ://doi.org/10.1093/bioin​ Ross, E. M., Moate, P. J., Marett, L. C., Cocks, B. G., & Hayes, forma​tics/btp515 B. J. (2013). Metagenomic predictions: From microbiome to Leahy, S. C., Kelly, W. J., Ronimus, R. S., Wedlock, N., Altermann, E., complex health and environmental phenotypes in humans and & Attwood, G. T. (2013). Genome sequencing of rumen bacteria and cattle. PLoS ONE, 8(9), e73056. https​://doi.org/10.1371/journ​ archaea and its application to methane mitigation strategies. Animal, al.pone.0073056 7(s2), 235–243. https​://doi.org/10.1017/S1751​73111​3000700 Rousseeuw, P. J. (1987). Silhouettes: A graphical aid to the in- Li, F., Li, C., Chen, Y., Liu, J., Zhang, C., Irving, B., … Guan, L. L. terpretation and validation of cluster analysis. Journal of (2019). Host genetics influence the rumen microbiota and heritable Computational and Applied Mathematics, 20, 53–65. https​://doi. rumen microbial features associate with feed efficiency in cattle. org/10.1016/0377-0427(87)90125-7 Microbiome, 7(1), 92. https​://doi.org/10.1186/s40168-019-0699-1 Saborío‐Montero, A. (2018). Variance components estimation of com- Manafiazar, G., Zimmerman, S., & Basarab, J. A. (2016). Repeatability plex traits including microbiota information. Master of Science. and variability of short‐term spot measurement of methane and car- Universidad Politecnica de Valencia. bon dioxide emissions from beef cattle using GreenFeed emissions Santos, E. dO., & Thompson, F. (2014). Other major lineages of bac- monitoring system. Canadian Journal of Animal Science, 97(1), teria and the archaea. In E. F. D. Eugene Rosenberg, S. Lory, E. 118–126. https​://doi.org/10.1139/CJAS-2015-0190 Stackebrandt & F. Thompson (Eds.), The prokaryotes (pp. 363– Martin, C., Morgavi, D. P., & Doreau, M. (2009). Methane mitigation in 383). Berlin, Heidelberg, Germany: Springer, Berlin Heidelberg. ruminants: From microbe to the farm scale. Animal, 4(3), 351–365. Scheller, S., Goenrich, M., Thauer, R. K., & Jaun, B. (2013). Methyl‐ https​://doi.org/10.1017/S1751​73110​9990620 coenzyme M reductase from methanogenic archaea: Isotope effects McAllister, T. A., Cheng, K. J., Okine, E. K., & Mathison, G. W. (1996). on the formation and anaerobic oxidation of methane. Journal of Dietary, environmental and microbiological aspects of methane pro- the American Chemical Society, 135(40), 14975–14984. https​://doi. duction in ruminants. Canadian Journal of Animal Science, 76(2), org/10.1021/ja406​485z 231–243. https​://doi.org/10.4141/cjas96-035 Schloss, P. D., Westcott, S. L., Ryabin, T., Hall, J. R., Hartmann, M., McCabe, M. S., Cormican, P., Keogh, K., O’Connor, A., O’Hara, E., Hollister, E. B., … Weber, C. F. (2009). Introducing mothur: Open‐ Palladino, R. A., … Waters, S. M. (2015). Illumina MiSeq phylo- Source, platform‐independent, community‐supported software for genetic amplicon sequencing shows a large reduction of an unchar- describing and comparing microbial communities. Applied and acterised succinivibrionaceae and an increase of the methanobrevi- Environmental Microbiology, 75(23), 7537. https​://doi.org/10.1128/ bacter gottschalkii clade in feed restricted cattle. PLoS ONE, 10(7), AEM.01541-09 e0133234. https​://doi.org/10.1371/journ​al.pone.0133234 Schubert, M., Lindgreen, S., & Orlando, L. (2016). AdapterRemoval Myer, P. R., Smith, T. P. L., Wells, J. E., Kuehn, L. A., & Freetly, H. v2: Rapid adapter trimming, identification, and read merg- C. (2015). Rumen microbiome from steers differing in feed effi- ing. BMC Research Notes, 9, 88–88. https​ ://doi.org/10.1186/ ciency. PLoS ONE, 10(6), e0129174. https​://doi.org/10.1371/journ​ s13104-016-1900-2 al.pone.0129174 Seedorf, H., Kittelmann, S., Henderson, G., & Janssen, P. H. (2014). O'Herrin, S. M., & Kenealy, W. R. (1993). Glucose and carbon diox- RIM‐DB: A taxonomic framework for community structure anal- ide metabolism by Succinivibrio dextrinosolvens. Applied and ysis of methanogenic archaea from the rumen and other intesti- Environmental Microbiology, 59(3), 748–755. nal environments. PeerJ, 2, e494–e494. https​://doi.org/10.7717/ Paulson, J. N., Stine, O. C., Bravo, H. C., & Pop, M. (2013). Differential peerj.494 abundance analysis for microbial marker‐gene surveys. Nature Shannon, C. (1984). A mathematical theory of communication. Methods, 10, 1200. https​://doi.org/10.1038/nmeth.2658 Tapio, I., Snelling, T. J., Strozzi, F., & Wallace, R. J. (2017). The rumi- Pérez, P., & de los Campos, G. (2014). Genome‐wide regression and nal microbiome associated with methane emissions from ruminant prediction with the BGLR statistical package. Genetics, 198(2), livestock. Journal of Animal Science and Biotechnology, 8, 7–7. 483–495. https​://doi.org/10.1534/genet​ics.114.164442 https​://doi.org/10.1186/s40104-017-0141-0 Pope, P. B., Smith, W., Denman, S. E., Tringe, S. G., Barry, K., Hugenholtz, Wallace, R. J., Rooke, J. A., McKain, N., Duthie, C.‐A., Hyslop, J. J., P., … Morrison, M. (2011). Isolation of Succinivibrionaceae impli- Ross, D. W., … Roehe, R. (2015). The rumen microbial metagenome cated in low methane emissions from tammar wallabies. Science, associated with high methane production in cattle. BMC Genomics, 333(6042), 646. https​://doi.org/10.1126/scien​ce.1205760 16(1), 839. https​://doi.org/10.1186/s12864-015-2032-0 Quast, C., Pruesse, E., Yilmaz, P., Gerken, J., Schweer, T., Yarza, P., … Wallace, R. J., Sasson, G., Garnsworthy, P. C., Tapio, I., Gregson, E., Glöckner, F. O. (2012). The SILVA ribosomal RNA gene database Bani, P., … Mizrahi, I. (2019). A heritable subset of the core rumen project: Improved data processing and web‐based tools. Nucleic microbiome dictates dairy cow productivity and emissions. Science Acids Research, 41(D1), D590–D596. https​://doi.org/10.1093/nar/ Advances, 5(7), eaav8391. https​://doi.org/10.1126/sciadv.aav8391 gks1219 Weimer, P. J. (2015). Redundancy, resilience, and host specificity of the Ricci, P., Rooke, J. A., Nevison, I., & Waterhouse, A. (2013). Methane ruminal microbiota: Implications for engineering improved ruminal emissions from beef and dairy cattle: Quantifying the effect of phys- fermentations. Frontiers in Microbiology, 6, 296. iological stage and diet characteristics. Journal of Animal Science, Whittaker, R. H. (1972). Evolution and measurement of species diver- 91(11), 5379–5389. https​://doi.org/10.2527/jas.2013-6544 sity. Taxon, 21(2/3), 213–251. https​://doi.org/10.2307/1218190 Rohart, F., Gautier, B., Singh, A., & Lê Cao, K.‐A. (2017). mixOmics: Wongnate, T., Sliwa, D., Ginovska, B., Smith, D., Wolf, M. W., Lehnert, An R package for ‘omics feature selection and multiple data integra- N., … Ragsdale, S. W. (2016). The radical mechanism of biologi- tion. PLOS Computational Biology, 13(11), e1005752. https​://doi. cal methane synthesis by methyl‐coenzyme M reductase. Science, org/10.1371/journ​al.pcbi.1005752 352(6288), 953. https​://doi.org/10.1126/scien​ce.aaf0616 RAMAYO‐CALDAS et al.    | 59 Yu, Z., & Morrison, M. (2004). Improved extraction of PCR‐quality community DNA from digesta and fecal samples. BioTechniques, How to cite this article: Ramayo‐Caldas Y, 36(5), 808–812. https​://doi.org/10.2144/04365​ST04 Zingaretti L, Popova M, et al. Identification of rumen microbial biomarkers linked to methane emission in Holstein dairy cows. J Anim Breed Genet. SUPPORTING INFORMATION 2020;137:49–59. https​://doi.org/10.1111/jbg.12427​ Additional supporting information may be found online in the Supporting Information section at the end of the article.