The study was conducted from September 2018 to March 2020 on 17 dairy farms in the Netherlands. In this time frame, farms were visited one to three times ( Supplementary Table 1 ), each time for a period of at least four weeks; a minimum of two weeks of adaptation followed by two weeks of measurement period. The farms were selected from the project “Koeien & Kansen”, a multi-annual research and demonstration project in the Netherlands that aims to be representative of the Dutch dairy sector. The study was conducted in accordance with the Dutch Animal Experiments Act in compliance with European Union Directive 2010/63 and approved by the Central Committee for Animal Experiments (The Hague, The Netherlands, 2016.D-0066.001).

In total 212 dairy cows were sampled in this study, of which 15 cows were measured multiple times over the study period with a maximum of three times. Because there was considerable time between each measurement period, including changes in diet, lactation stage and parity, observations of the same cows were considered to be independent. All cows were housed in free-stall barns with deep-litter sand or rubber matrass with sawdust or in a (composted) bedded pack barn. All barns had open sides for natural ventilation and free access to clean drinking water. All cows were fed ad libitum and received a commercial dairy diet that differed per farm and season. The diets were formulated according to commercial practices complied in collaboration with a dairy feed advisor. Practices and diets per farm were maintained unchanged during the adaptation period and measurement period as the objective was to have a representative overview of Dutch farming on CH₄ emission and rumen microbiota. On average (mean ± SD) the diets consisted of 32 ± 13.6% grass silage, 25 ± 11.8% maize silage, 12 ± 13.8% fresh grass, 25 ± 7.9% concentrates, and 6 ± 10.7% by-products, on a DM basis. Supplementary Table 1 presents the diet composition provided on each farm during each measurement period. Fifteen farms applied grazing during part of the year (spring, summer and/or autumn) and thus also during some of the measurement periods. Grazing time was not controlled and differed per farm, ranging from a few hours a day to unrestricted grazing with free choice to stay outside the barn. Farmers were advised to keep the grazing management as homogenous as possible during the measurement period.

Depending on farm, cows were milked either twice daily in a milking parlor or two to three times daily by a milking robot. Milk composition was determined at the individual cow level using milk samples collected on a single day in the second measurement week. Milk samples were analyzed by Fourier transform mid-infrared spectrometry (MIRS) as part of the routine milk recording programs, for 15 farms performed by Qlip B.V. (Zutphen, the Netherlands) and for 2 farms by VVB Veluwe-Ijsselstreek (Nunspeet, the Netherlands). Milk yield and the percentage of protein and fat in the milk of a sample day were used to calculate fat- and protein-corrected milk (FPCM) according to the following equation (CVB, 2016).

Enteric CH₄ production was measured non-invasively using the GreenFeed system (GF, C-lock Inc. Rapid City, SD, USA). The GF is an adapted feeding station that measures individual CH₄ and CO₂ production in grams per day during each visit, as described in detail by van Gastelen et al. (2022). For the present study, the average recovery of CO₂ was 99.1% (for individual GF systems between 97.6 and 102.7%).

Methane production was measured over a 14-day period, with a minimum 7-day acclimation to the GF prior to the measurement period and visit times per cow ranging from 3 to 6 min. Cows received 2 to 6 kg of compound feed per day via the GF, depending on milk yield and lactation stage. Administration of the compound feed was divided into 4 to 6 feeding periods per day, with at least 3 to 4 h between each feeding period. Per feeding period, 0.5 to 1 kg of compound feed was administered in 12–25 drops of approximately 40 g of feed (ranging between 31 and 51 g per drop, depending on the type of compound feed used). There were 10–30 s between each drop (depending on the maximum number of drops) to ensure a minimum visit time of 3 min, but no longer than 6 min.

On each farm, one GF was installed in the barn. In case grazing was applied at farm level, an additional GF on a pasture trailer was used to ensure adequate CH₄ measurements throughout the day. Due to the high diurnal variation of enteric CH₄ production and using short-term breath measurements, multiple records are needed to provide a representative average. Therefore, based on what was found in the study of Manafiazar et al. (2017), CH₄ measurements from the 14-day measurement period were averaged per cow and 32 cows with fewer than 20 valid records were excluded from the analysis. Methane intensity (g CH₄.kg FPCM⁻¹.d⁻¹) was used as a metric for CH₄ emission rather than CH₄ yield (g CH₄.kg DM⁻¹.d⁻¹). Methane production (g.d⁻¹) is strongly related to total feed intake, which makes it necessary to take measures of feed intake into account when comparing CH₄ emission. However, feed intake was not available on individual cow level, thus CH₄ production was corrected for FPCM yield, as it is the closest available measure related to feed intake. Note that CH₄ production was averaged over two weeks while FPCM yield resulted from a single day of milk samples.

Rumen fluid samples were collected in the second week of the measurement period ( Supplementary Table 1 ) using the oral stomach tube (OST) (Muizelaar et al., 2020). Briefly, the OST consisted of a manual pump and a 190 cm long spiral probe with a perforated suction head at the end, inserted through the esophagus into the dorsal cranial part of the rumen. Rumen fluid was collected by the manual pump. The first 500 ml of the rumen fluid was discarded to minimize contamination with saliva. Subsequently, samples of 3 ml were collected and immediately frozen in dry ice. Samples were stored at –80°C until analysis.

Extraction of microbial DNA from rumen samples, library construction of hypervariable region V4 (from 16S rRNA gene) and subsequent sequencing on an Illumina HiSeq platform were performed at Genotypic Technology Pvt. Ltd. In Bangalore, India.

DNA from rumen fluid samples was purified using the Qiagen Dneasy Blood and tissue Kit (Qiagen, Hilden, Germany). Prior to processing the samples using this method, approximately 200μl of ruminal fluid was taken in sterile Tomy tubes containing 3–4 beads and homogenization was carried out at 4,000 rpm for 120 s. A volume of 200 μl (10mg ml⁻¹) lysozyme (MilliporeSigma, St. Louis, USA) was added to the homogenates. The homogenate tube was invert mixed and incubated for 30 min at 37°C. A volume of 200μl AL buffer was added to the samples and vortexed briefly. The samples were subjected to Proteinase K treatment at 56°C for 2 hours followed by Rnase A treatment (MP Biomedicals, Solon, USA) at 65°C for 20 min. The lysate was mixed well with 100% ethanol and loaded onto Qiagen Dneasy blood and tissue column (Qiagen, Hilden, Germany). The samples were further processed following manufacturer’s instructions. Finally, DNA was eluted in 45μl 10mM Tris.Cl, pH 8.0 (Sigma-Aldrich, St. Louis, USA). The concentration and purity of genomic DNA was quantified using the Nanodrop Spectrophotometer (Thermo Scientific; 2000). The integrity of the DNA was assessed by agarose gel electrophoresis.

Sequencing libraries were prepared by a two-step polymerase chain reaction (PCR)-based workflow based on primers specific to V4 region of the 16S rRNA gene (Caporaso et al., 2011). In the first round of PCR, the V4 region of the 16S rRNA gene was amplified using the region-specific primers V4-515F (GTGCCAGCMGCCGCGGTA) and V4-806R (GACTACHVGGGTATCTAATCC) designed by Genotypic Technology Pvt. Ltd. In Bangalore, India. Using the KAPA HiFiHotStart PCR Kit (KAPA Biosystems Inc., Boston, MA USA) and a primer concentration of 5μM, 50ng of genomic DNA was amplified by following first round of PCR condition: 3 min at 95°C, 26 × (30 s at 95°C, 30 s at 64°C, and 30 s at 72°C), and 5 min at 72°C. The amplicons generated were analyzed on a 1.2% agarose gel. The second round of PCR was performed to index the amplicons generated in the first round of PCR. 1μl of the 1:2 diluted PCR amplicons generated in the first round were amplified by following second round of PCR condition: 3 min at 95°C, 10 × (30 s at 95°C, 30 s at 55°C, and 30 s at 72°C), and 5 min at 72°C to add Illumina barcoded adaptors for sequencing (Nextera XT v2 Index Kit, Illumina, USA). Amplicons (sequencing libraries) generated by the second round of PCR were analyzed on a 1.2% agarose gel. The libraries were normalized and pooled for high-throughput multiplex sequencing. Finally, these pools were quantified using the Qubit dsDNA HS assay and fluorometer (Thermo Fisher Scientific, MA, USA). The normalized sample was denatured and fed into the Illumina HiSeqXTen sequencer, where it was sequenced for 150*2 cycles to generate at least 0.7 million paired-end reads. Upon completion of the sequencing run, the data were demultiplexed using bcl2fastq v2.20 software (Illumina, San Diego, USA).

Reads were preprocessed using QIIME2 suite v2020.8 (Bolyen et al., 2019). Reads were first checked for low quality bases (Q30) with FAST QC v 0.11.9, then trimmed of primers with Cutadapt plug-in before being merged using FASTQ-join. The merged reads were denoised, dereplicated, and chimera sequences were removed using the DADA2 plugin (Callahan et al., 2016) with default parameters from the QIIME2 suite. The resulting Amplicon Sequence Variants (ASV) were classified based on SILVA v.138 database (Quast et al., 2013). To achieve the classification a Naive Bayesian classifier, pre-trained, was used with the ‘feature-classifier classify-sklearn’ command implemented in Qiime2 (Bokulich et al., 2018). The classifier was optimized for 515F/806R region at 99% similarity and can be found here: https://resources.qiime2.org/#naive-bayes-classifiers-2. Confidence level was set at 0.7% by default.

A balanced experimental design was established from this large cohort after data generation to balance the factors influencing CH₄ intensity and rumen microbiota ( Figure 1A ). First, two dairy cows with outlier CH₄ intensity (above third quartile added of three times the interquartile range) were removed from the dataset. Second, dairy cows were grouped into three lactation stages, namely early (0–93 days in milk (DIM)), mid (93–183 DIM) and late (>183 DIM), according to Bainbridge et al. (2016). Third, dairy cows within each lactation stage were split in two dietary categories, a group that had access to fresh grass during the measurement period and a group that had not access to fresh grass during the measurement period. This resulted in six groups. Fourth, the top fifth lowest and top fifth highest CH₄ emitters within each of these six groups were selected ( Figure 1B ), after checking for statistical significance of the difference of their mean CH₄ intensity by a Welch test ( Figure 1A ), and placed in two groups; a low and high CH₄ emitting group (N=30 per group) that was balanced for lactation stage and diet type (i.e., with or without access to fresh grass). All the high and low emitters cows within this experimental design were either Holstein Friesian or Holsten Friesian cross breed.

All statistical analysis were performed in R v4.0.3 (R Core Team, 2020) using phyloseq (McMurdie and Holmes, 2013), vegan (Oksanen et al., 2012) and ropls (Thévenot et al., 2015) packages. Methane intensity was compared between low and high CH₄ emitter groups by a Welch test. The low-high dataset previously described was filtered for low abundant and low prevalent taxa (more than five in at least a third of individuals) to limit the zero-inflation issue (from 90.1% 0 in the dataset prior filtering to 34.0% 0 after filtration). Relative abundance of ASVs were then transform by centered log-ratio (CLR) which allows to overcome differences in sequencing depth without wasting data as in rarefaction (McMurdie and Holmes, 2014; Gloor et al., 2017) and at the same time replace data into a Euclidean space (Quinn et al., 2019). Multivariate dimension reduction approaches, namely principal component analysis (PCA) and orthogonal partial least-square discriminant analysis (o-PLS-DA) were performed following the CLR normalization and a z-score scaling to adjust feature-wise homogeneity. Orthogonal signal correction was used to decorrelate variation that would be unrelated to the discriminant low and high emitter groups. Significance of the predictive performance (Q2) was assessed by random permutation of 10% of the samples in the discriminant groups repeated 200 times. Variables important for projection with a value above two were defined as associated with low or high emitters depending on their contribution on the predictive component. The relative abundance of taxa was used to plot the VIP rather than CLR transformed value for comprehension purpose.

Methane intensity from the broad cohort of 212 cows, depleted of two outliers and 32 cows that visited the GF less than 20 times, ranged from 6.55 to 25.02 g.kg FPCM⁻¹.d⁻¹ ( Figure 2 ). The low CH₄ emitters were predominantly cows in early or mid-lactation, whereas the high CH₄ emitters were mostly cows in late lactation. To isolate individual effects on CH₄ intensity from important confounding factors, we created a balanced experimental design that account for lactation stage (i.e., early, mid, and late) and diet type (i.e., with or without access to fresh grass) for two groups of low and high CH₄ emitters. The average CH₄ intensity in the low CH₄ group (9.5 ± 1.57 g.kg FPCM⁻¹.d⁻¹) was significantly different from average CH₄ intensity in the high CH₄ group (16.8 ± 4.11 g.kg FPCM⁻¹.d⁻¹, P < 0.001). The group composition was not biased for specific farms ( Supplementary Figure 1 ) and only one cow in each group was sampled twice.

The dataset from the balanced experimental design contained 23,578 ASV for 60 samples before filtering for low abundant and low prevalent ASV. After the filtering, 2,250 taxa remained. On average, 98.10 ± 2.10% of the ASV were identified as bacteria and 1.90 ± 2.10% were identified as archaea ( Supplementary Figure 2 ). At the Phylum level, Firmicutes (56.6 ± 9.46%), Bacteroidota (31.90 ± 10.17%) and Proteobacteria (3.10 ± 3.55%) composed most of the taxa ( Supplementary Figure 3 ).

The microbial composition at Phylum level across the 17 farms is presented in Supplementary Figure 4 . In farm “I”, the relative abundance of the most abundant phylum, Firmicutes, varied considerably, ranging from a minimum of 47.12 ± 9.18% to a maximum of 71.04 ± 10.55%. A PCA was performed to have an overview of the main factors that influenced the composition of rumen microbiota of cows in the balanced experimental design ( Figure 3 ). The first two principal components (PC) of the PCA explained 29% of the total variance. The most visible separation on the PCA, if any, was related to diet type, namely, the presence of fresh grass in the diet. No discrimination by the lactation stage or by the grouping of cows as low or high CH₄ emitters was observed.

A multivariate discriminant analysis, o-PLS-DA, was performed to reveal ASV associated with low and high CH₄ emitter groups in the balanced experimental design ( Figure 4 ). The predictive accuracy of the o-PLS-DA was 0.16 (Q2 = 0.16) and significant (P = 0.01). A total of 88 VIP (ASVs used in the model) were above the defined threshold of two, i.e., for VIP highly discriminant for the two groups ( Figure 5 ). Collectively, these VIP represented 14.0 ± 4.9% of the total relative abundance across all the 60 samples of the balanced experimental design ( Supplementary Figure 5 ). Thirty-six taxa were more abundant in the high CH₄ emitter group than in the low emitter group, whereas 52 were more abundant in the low CH₄ emitter group than in the high emitter group. The 36 ASV that were relatively more abundant in the high CH₄ emitter group were identified in the family [Eubacterium]_coprostanoligenes, Absconditabacteriales_(SR1), Christensenellaceae, Clostridia_UCG−014, Gastranaerophilales, Hungateiclostridiaceae, Lachnospiraceae, Oscillospiraceae, Prevotellaceae, Ruminococcaceae, Saccharimonadaceae, Spirochaetaceae and one uncultured family from the Armatimonadota Phylum. The 52 ASV more abundant in low CH₄ emitter group were identified in the family [Eubacterium]_coprostanoligenes, Anaerolineaceae, Christensenellaceae, Clostridia_UCG−014, Desulfobulbaceae, Eubacteriaceae, F082 from the Bacteroidales order, Lachnospiraceae, Methanobacteriaceae, unknown family from the Clostridia class, Oscillospiraceae, Prevotellaceae, Rikenellaceae and Ruminococcaceae.