This study was conducted in strict adherence to the Animal Protection Law as stipulated in the Guide for the Care and Use of Laboratory Animals and approved by the Ethics Committee of Nanjing Agricultural University, Nanjing, China.

A total of 18 lactating Holstein cows on second parity and 56 (±14.81 SD) DIM were selected from 94 cows on a commercial dairy farm (Shanghai Yi Nan Dairy farm, Shanghai, China). All the cows were sired using artificial insemination and the semen used was purchased from the same company. The cows have been on the farm and fed the same diet for more than 60 days. The nutritional composition of the diets was based on the N.R.C. (2001) requirements for lactating Holstein cows (Table 1) and the average dry matter intake (DMI) for all the cows within the same row was recorded. All the animals were healthy and had access to clean drinking water ad libitum.

Milk samples were collected from all the cows three times daily (2 a.m., 10 a.m., and 5 p.m.) into a 50 ml tube containing a preservative (K₂Cr₂O₂) for three consecutive days. The milk sample for each collection day was mixed at a ratio of 4:3:3 for the samples collected at 2 a.m., 10 a.m., and 5 p.m., respectively. The milk samples were then stored at 4°C while waiting to be sent for milk composition analysis. The milk samples were sent to the laboratory (Shanghai DHI company, Shanghai, China) within 2 days of collection for the determination of milk composition using near-infrared spectrometry (MilkoScan Minor FT120, Denmark). The determination of milk composition and somatic cell count was carried out separately for each collection day. Milk yield was recorded with the aid of a milk sampler (Waikato Milking Systems NZ Ltd., New Zealand). Another portion of the milk sample were collected from each cow into a 2 ml tube for the determination of milk metabolites. The milk samples were collected for three consecutive days and mixed to obtain a composite sample (day 1, day 2, and day 3) before analysis of milk metabolites. The milk samples were then stored at −20°C until required for analysis of metabolites. Based on the milk yield of the cows, 18 cows were selected from the 94 cows and divided into 2 groups of 9 cows each. To ensure that cows in the two groups have significantly different milk yields, nine cows with the highest milk yield (44.57 ± 2.11 kg/day) were placed in the high-producing (HP) group while nine cows with the lowest milk yield (26.71 ± 0.70 kg/day) were allotted to the low-producing (LP) group. The selected cows have an average weight of 639.17 ± 17.49 kg.

Rumen content was collected from each of the 18 cows before morning feeding into 50 ml tubes (×2) using an oral rumen tube. To maintain the purity of rumen content and avoid mixing with saliva, the first 150 ml of rumen content collected was discarded. The rumen content (10 ml) to be used for metagenomics and MCP determination was immediately stored in liquid nitrogen. Another portion of the rumen content was filtered through four layers of sterile cheesecloth into 5 ml tubes (three tubes) to be used for the detection of VFAs, ammonia nitrogen (NH₃-N), and rumen fluid metabolites. For the rumen fluid sample meant for VFA analysis, 1 ml of 250 g/L metaphosphoric acid was added to the 5 ml rumen fluid filtrate. For rumen fluid filtrate meant for the determination of NH₃-N, exactly 1 ml of 20 g/L (w/v) H₂SO₄ was added to the 5 ml and NH₃-N was detected using the phenol-hypochlorite reaction (Weatherburn, 1967). The samples were stored at −20° until needed for the determination of rumen fluid metabolites and rumen fermentation parameters. The analysis of VFAs of the ruminal fluid samples was performed with gas chromatography (GC-2014B, Shimadzu, Japan) equipped with a capillary column (film thickness: 30 m × 0.32 mm × 0.25 μm, column temperature: 110°C). The temperatures of both the injector and detector were preset at 180°C.

Microbial protein was determined using the Bradford protein assay. Prior to the protein determination, the rumen fluid bacterial cells were extracted as described by Makkar et al. (1982). Briefly, about 10 ml of rumen fluid was stirred with a magnetic stirrer at 400 rpm for 45 s to remove the microbes attached to the feed particles and then centrifuged for 5 min at 408 × g to remove the protozoa. Five milliliters aliquot was taken from each of the samples and centrifuged again for 20 min at 25,000 × g. The supernatants were discarded, and pellets were washed with distilled water and centrifuged again at 25,000 × g for 20 min. The supernatants were discarded again, and the cells were suspended in 15 ml of 0.25 N NaOH and boiled for 10 min in a water bath before centrifuging at 25,000 × g for 30 min. Finally, the supernatant was collected and MCP was determined from the supernatant using the Bradford protein assay (Bradford, 1976).

Exactly 5 μl of the supernatant for each sample and 5 μl of a standard (Bovine serum albumin) were taken and transferred into microplate wells. Exactly 250 μl of 1× dye reagent (dye solution + methanol and phosphoric acid) was added to each of the wells and mixed thoroughly using a pipette by pressing the plunger repeatedly. A new pipette tip was used for each sample. The mixture was allowed to incubate at room temperature for 1 h and the absorbance was measured with a microplate reader (BioTek, Synergy H1, United States) set at 595 nm. The MCP concentration of the rumen fluid is estimated using a standard curve by plotting the 595 nm values on the y-axis versus their concentration in μg/ml on the x-axis.

Blood samples to be used for detection of serum metabolites were collected once from each cow via the tail vein 1 h before morning feeding into a 5 ml tube containing heparin. Another set of 5 ml tubes without anticoagulant was used to collect blood samples for detection of serum biochemical indices. A warm water bath (37°C) was used to keep the clotted blood for at least 20 min and then all the blood samples collected were centrifuged at 3000 × g for 15 min. The serum was collected with a pipette into 2 ml tubes and stored in a refrigerator at −80° until required for detection of serum biochemical parameters and serum metabolites. Serum biochemical indices were detected by injecting the serum samples into the tubes of an automated clinical chemistry analyzer (model CA-800) at the Jiangsu Hospital of Traditional Chinese Medicine, Nanjing, China.

The analysis of metabolites in the rumen fluid, serum, and milk was carried out at Shanghai cluster Biotechnology Co., Ltd., Shanghai, China. The stored samples (rumen fluid, serum, and milk) were slowly thawed on ice, and 100 μl of each sample was transferred into a 1.5 ml centrifuge tube. Exactly 300 μl of methanol was added and the solution was vortexed for 30 s. The solution was incubated for 1 h at 40°C and vortexed again for 30 s before centrifuging (12,000 rpm) for 15 min at 4°C. Exactly 200 μl of the supernatant was taken and mixed with 5 μl of internal standard (0.14 mg/ml dichlorophenylalanine) before transferring the solution into an injection vial.

The detection of metabolites was performed using an LC–MS (Thermo, ultimate 3000lc). The machine (LC–MS) is equipped with a C18 column [hyper gold C18 (100 × 2.1 mm, 1.9 μm)] preset at a flow rate of 0.3 ml/min and a temperature of 40°C. The mobile phase is made up of A (water, 5% acetonitrile, and 0.1% formic acid) and B (acetonitrile and 0.1% formic acid). The elution process of the mobile phase is shown in Supplementary Table 1. Exactly 4 μl of each sample was injected into the autosampler, which was preset at 4°C. The spray voltage was set at 3.0 kV for the positive and 3.2 kV for the negative modes while other parameters were the same for both the positive and negative modes (Sheath gas flow rate 15 arb heater temperature 300°C, and capillary temperature 350°C). The extraction of the raw peak, filtration and calibration of the data baseline, identification of peak, alignment of the peak, integration of the peak area, and deconvolution analysis were performed on the LECO-Fiehn Rtx5 database (Kind et al., 2009). Metabolite profiling and identification were performed as described by Wang et al. (2018).

Feed samples were randomly collected into a plastic bag at different points during feeding and stored at −80°C until required for analysis of chemical composition. The crude protein, crude fiber, crude fat, and ash were analyzed according to the official analytical methods of the A.O.A.C. (1991) whereas the acid detergent fiber (ADF) and neutral detergent fiber (NDF) were determined following the procedure of Van Soest et al. (1991). The feed ingredients and their nutrient composition are provided in Table 1.

Microbial DNA extraction was carried out using the E.Z.N.A.^® stool DNA Kit while strictly adhering to the procedures stipulated by the company (Omega Bio-Tek, Norcross, GA, United States). In brief, the rumen samples were previously obtained from each cow in duplicates via oral rumen tube and were mixed to get a composite sample before storage in liquid nitrogen. Frozen rumen samples were thawed at 4°C, centrifuged at low speed to obtain the particulate matter and about 200 mg of the composite sample from each cow was collected into a disruptor tube. Exactly 725 μl of a lysis buffer was added and the samples were vortexed on a bead-ruptor for 3--5 min to lyse the samples. DNA was extracted from each sample following the guidelines of the DNA Kit’s manufacturer (Omega Bio-Tek, Norcross, GA, United States) and DNA concentration was measured with a NanoDrop ND-1000 Spectrophotometer (Thermo Scientific, United States). The construction of the metagenomic shotgun sequencing libraries and the DNA sequencing were done at Shanghai Biozeron Biological Technology Co., Ltd. (Shanghai, China). A TruSeq DNA PCR-Free library preparation kit (Illumina, San Diego, CA, United States) was used to construct sequencing libraries (length of about 350 bp) from the extracted DNA. An Illumina NovaSeq instrument was used for sequencing all the samples. Thereafter, raw sequence reads were subjected to quality trimming using Trimmomatic¹ following the procedures described by Bolger et al. (2014). Good quality reads that passed the quality trimming were mapped against Bos taurus ARS-UCD1.2, using the BWA mem algorithm (parameters: -M -k 32 -t 16²). Good quality reads devoid of host genome contaminations (clean reads) were used for phylogenetic annotation.

Clean reads obtained from each sample were further subjected to taxonomic classification with Kraken2 software (Wood and Salzberg, 2014), which uses the National Centre for Biotechnology Institute (NCBI) RefSeq database represented by a series of clade-specific k-mers in which the clean reads are matched against. The classification was done at seven phylogenetic levels which are domain, phylum, class, order, family, genus, species, or unclassified. Bracken³ was used to calculate the abundance of each taxonomy. The functional annotations of the microbiota were obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database⁴ using the KEGG Orthology Based Annotation System (KOBAS version 2.0).⁵

A set of contigs was obtained from the clean reads using MegaHit (parameters: –min-contig-len 500) (Li et al., 2015). The prediction of the open reading frames (ORFs) of assembled contigs was performed with Prodigal software (v2.6.3) (Hyatt et al., 2010), while CD-HIT (parameters: -n 9 -c 0.95 -G 0 -M 0 -d 0 -aS 0.9 -r 1) (Fu et al., 2012) was used for clustering to obtain a unique-gene set. Within a unique gene set, the longest sequence in each of the clusters is the representative sequence of each gene. Salmon software (Patro et al., 2017) was used to obtain the number of reads for each gene, and the gene abundance of all the samples was calculated using the formula below (Li et al., 2014):

where, Ab(S), gene abundance; Ab(U), single-mapping reads abundance; Ab(M), multi-mapping reads abundance; l, length of gene sequence.

To obtain the list and functional annotations of Carbohydrate-Active Enzymes (CAZymes), the gene sequence was converted into an amino acid sequence, then DIAMOND software (Buchfink et al., 2015) was used to align the sequence to the CAZymes database (Lombard et al., 2014). The sequences of rumen microbial metagenome were submitted to the Sequence Read Archive (SRA) of the NCBI with the accession number PRJNA748874.⁶

At the commencement of the study, a power analysis was conducted with SAS (9.4) to determine the number of cows per group, and the results revealed that a minimum of seven cows per group is required to obtain a statistical power of 0.84.

The results obtained from rumen fermentation parameters, DMI, milk composition, and yield were analyzed with the Student’s T-test using the IBM SPSS (v20). The results were considered significant at P < 0.05. The results of the metabolites from milk, rumen fluid, and serum were first subjected to principal component analysis (PCA) and partial least square discriminant analysis (PLS-DA) using SIMCA (13.0). Those metabolites with variable important (VIP) scores greater than 1 were further analyzed using the Student’s T-test (P < 0.05) in IBM SPSS (v20). Differential metabolites between the HP and LP group were matched to the KEGG pathway (B. taurus library) on MetaboAnalyst 5.0.⁷ The rumen microbiota for all the domains (Archaea, Bacteria, Eukaryota, and Virus) at phylum and genus levels, the abundance of genes encoding CAZymes, KEGG genes, and functions were filtered by relative abundance greater than 0.1% in at least 1 sample before subjecting them to the linear discriminant analysis effect size (LEfSe) (Segata et al., 2011) using web-based software.⁸ LEfSe uses the Kruskal–Wallis test to identify significantly different parameters before subjecting them to linear discriminant analysis (LDA) to estimate the effect size of each differential parameter. The taxa, CAZymes, KEGG genes, and functions were declared significant at P < 0.05 and LDA > 2.

The PCA, principal coordinate analysis (PCoA), and the permutational multivariate analysis of variance (PERMANOVA) of the rumen microbiota were conducted using the Vegan package in R software (4.05). Spearman’s correlation analysis was conducted using IBM SPSS (v20) to identify the relationship among the metabolites (milk, rumen fluid, and serum), predominant microbiota (relative abundance >1%), fermentation parameters, and milk yield. The correlation between metabolites (milk, rumen fluid, and serum) and milk yield was conducted and those metabolites with a strong positive association with milk yield (R > 0.6, P < 0.05) were tagged MP-metabolites (Supplementary Table 10). Furthermore, the predominant rumen microbiota were correlated (R > 0.6, P < 0.05) with the MP-metabolites (milk, rumen fluid, serum) (Supplementary Table 12), and the fermentation parameters in the rumen (R > 0.6, P < 0.05) (Supplementary Table 13). All the significant correlation networks were visualized using Cytoscape (3.8.2).

The results obtained from milk composition analysis and somatic cell count showed that the milk yield (kg/day), milk fat yield (kg/day), milk protein yield (kg/day), lactose yield (kg/day), and milk solid yield (kg/day) were significantly higher (P < 0.01) in the HP group as compared to the LP group. The average milk yield of the HP group is 44.57 kg/cow per day whereas that of the LP group is 26.71 kg/cow per day. The milk protein yield of the HP group is 1.33 kg/cow per day while that of the LP group is 0.81 kg/cow per day (Table 2 and Figure 1).

The results of rumen fluid analysis revealed that the percentages of the rumen VFAs (acetate, propionate, butyrate, isobutyrate) and the acetate: propionate ratio did not differ significantly (P > 0.05) between the HP and LP groups. However, the percentages of valerate and isovalerate were found to be significantly higher (P < 0.05) in the HP group whereas, the LP group recorded a significantly higher (P = 0.029) proportion of ruminal NH₃-N (mg/dl). The average concentration of NH₃-N observed in the HP group is 12.53 mg/dl while that of the LP group is 19.60 mg/dl (Table 3).

The metagenomic sequencing of the 18 samples generated a total of 1,299,004,492 raw reads with an average of 72,166,916.22 reads ± 2,167,434.65 per sample. After filtering and removing host genes, the total clean reads recorded were 1,178,340,750 with an average of 65,463,375 ± 1,922,941.62 reads per sample (Supplementary Table 2).

The PCoA was carried out based on Bray Curtis distance metrics and the plot represents 70% of the variation (Axis 1 and Axis 2) yet no clear separation between the HP and the LP group. Similarly, the PCA showed 82% of the variation (PC 1 and PC 2) yet no clear distinction between the HP and the LP group. However, the samples in the LP as observed in both the PCoA (Supplementary Figure 3) and PCA (Supplementary Figure 4) plots tend to cluster closer to each other than those of the HP group. Furthermore, PERMANOVA was performed on the abundance of Archaea, Bacteria, Eukaryota, and Viruses. The results revealed that none of the domains differed significantly (P > 0.05) between the HP and the LP groups (Supplementary Table 3 and Supplementary Figure 1).

After the taxonomic classification, a total of 3982 genera and 212 phyla were obtained in all 4 domains. The Archaea had 22 and 147 phyla and genera, respectively. In the Bacteria, 146 phyla and 3033 genera were found whereas the Eukaryota had 43 and 595 phyla and genera, respectively. However, the Viruses were unclassified at the phylum level but there were 207 genera identified. In the Archaea, the Phylum Euryarchaeota predominated with a relative abundance of 96.74 and 96.68% in the HP and LP groups, respectively (Supplementary Table 5A and Supplementary Figure 5B). At the genus level, Methanobrevibacter was the most abundant genera with a relative abundance of 83.70% in the HP group and 82.02% in the LP group (Supplementary Table 5B and Supplementary Figure 5A). At the genus level, Methanothermobacter was significantly higher (P = 0.012, LDA > 2) in the HP group while Thermoprotei_norank was found to be significantly higher (P = 0.047, LDA > 2) in the LP group (Figure 2D).

Among the Bacteria, Bacteroidetes are the most abundant phylum in both groups with HP and LP occupying 68.34 and 71.82%, respectively Supplementary Figure 5D. Firmicutes also represent a large portion of the rumen Bacterial phyla with the HP group occupying 21.43% while the LP group representing 17.73% (Supplementary Table 6A). For the Bacterial genus, Prevotella topped the list of predominant genera with the HP and LP groups having a relative abundance of 50.63 and 53.67%, respectively Supplementary Figure 5C. None of the Bacterial Phyla was found to differ significantly between the HP and LP groups. However, at the genus level, the relative abundance of the Sharpea was significantly higher (P = 0.047, LDA > 2) in the HP group when compared to the LP group (Figure 2E).

In the Eukaryota, the Ciliophora is the most abundant Phylum which represents about 62–63% in both the HP and the LP group (Supplementary Table 7A and Supplementary Figures 6A,B). At the phylum level, the relative abundance of Euglenozoa was significantly higher (P = 0.015, LDA > 2) in the LP group (Figure 2B and Supplementary Figure 7). As for the genera, eight genera were found to differ significantly between the HP and LP groups (Figure 2A). The HP group recorded a significantly higher abundance of the genera Trichomonas (P = 0.047, LDA > 2), Spironucleus (P = 0.038, LDA > 2), Allomyces (P = 0.019, LDA > 2), Saitoella (0.003, LDA > 2), and Neospora (P = 0.005, LDA > 2). On the other hand, the genera Kluyveromyces (P = 0.047, LDA > 2), Trypanosoma (P = 0.007, LDA > 2), and Bipolaris (0.046, LDA > 2) were found to be significantly higher in the LP group as compared to the HP group.

In the Viruses, the most abundant genus is the Podoviridae_norank which represents 37.77 and 36.46% in the HP and LP group, respectively (Supplementary Table 8 and Supplementary Figure 6C). The next most abundant genus is the Siphoviridae_norank and this genus was found to be significantly higher (P = 0.031, LDA > 2) in the HP group (26.68%) when compared to the LP group (21.49%). However, the LP group recorded a significantly higher abundance (P = 0.024, LDA > 2) of Teseptimavirus (Figures 2C, 3).

The functional annotation of the rumen microbiome was defined by the KEGG profiles. A total of 368 level-3 KEGG pathways were identified. After subjecting the KEGG profiles to LEfSe, ko05416 (Viral myocarditis) was significantly upregulated (P = 0.012, LDA > 2) in the LP group (Figure 4C).

A total of 342 genes associated with CAZymes belonging to the following families were identified. These are; 132 glycoside hydrolases (GHs), 30 polysaccharide lyases (PL), 80 carbohydrate-binding modules (CBMs), 16 carbohydrate esterases (CEs), 74 glycosyltransferases (GTs), and 10 auxiliary activities (AAs). However, only 15 CAZymes were found to differ significantly (P < 0.05, LDA > 2) between the HP and LP groups. The HP group recorded a significantly higher (P < 0.05, LDA > 2) expression of CAZymes belonging to GH149, GT21, PL22, and PL5. However, in the LP group, several CAZymes (CBM56, CBM62, GH22, GH46, GT25, GT32, GT45, GT6, PL11, PL12, and PL27) were significantly enriched (P < 0.05, LDA > 2) (Supplementary Table 4 and Figure 4A).

A total of 2720 unique KEGG Orthology (KO) genes were found and subjected to LEfSe. Two genes (K02316 and K01673) were identified to differ significantly (P < 0.05, LDA > 2) between the HP and LP groups. After matching the identified genes to the KO database, K02316 (dnaG) was identified to belong to the class of transferases while K01673 (cynT) is found to be involved in nitrogen metabolism. Interestingly, both K02316 and K01673 were enriched in the LP group (Figure 4B).

Predominant rumen microbiota were correlated with milk yield and all the 45 MP-metabolites in the rumen, serum, and milk (Supplementary Table 12). However, only 1 MP-metabolite in the milk, 3 MP-metabolites in the rumen, and 17 serum metabolites were found to strongly associate (R > 0.6/R < −0.6, P < 0.05) with the rumen microbiota. The rumen microbes with high and significant correlation (R > 0.6/R < −0.6, P < 0.05) with the MP-metabolites were considered as being associated with milk yield. A total of one Archaeon (Methanomassiliicoccales_norank), five Bacteria genera (Eubacterium, Butyrivibrio, Lachnospiraceae_norank, Ruminococcus, and Ruminococcaceae_norank), five Eukaryota genera (Kluyveromyces, Tritrichomonas, Neocallimastix, Anaeromyces, and Vitrella), and two Viruses (Myoviridae_norank, and Cequinquevirus) were identified to have a strong association (R > 0.6/R < −0.6, P < 0.05) with MP-metabolites (Figure 6). In addition, the correlation between predominant microbiota in the rumen and rumen fermentation parameters revealed a strong and positive association (R > 0.6, P < 0.05) between fermentation parameters and some genera from Bacteria, Archaea, Eukaryota, and viruses (Supplementary Table 13 and Supplementary Figure 2).