Characterization of phenotypes
In this study, previously reported milking traits were obtained from 374 dairy cows [8], and 10 cows with the highest MPY (cows with high milk yield and milk protein content; HH) and 10 cows with the lowest MPY (cows with low milk yield and milk protein content; LL) were selected for metagenome, rumen metabolome, and serum metabolome analyses. Among the phenotypes, milk yield (P < 0.01), milk protein content (P < 0.01), and MPY (P < 0.01) were significantly different between the HH and LL groups (Table S1).
Profiling of the rumen metagenome
Metagenome sequencing generated a total of 1,069,431,480 reads, with 66,839,468 ± 1,168,990 reads (mean ± standard error of the mean [SEM]) per sample (Table S2). After quality control and removing host genes, a total of 1,033,603,420 reads were retained, with 64,600,214 ± 1,165,364 per sample. After de novo assembly, a total of 12,097,293 contigs were generated (the N50 length of 795 ± 28 bp), with 756,081 ± 27,721 per sample. The rumen metagenome consisted of 94.43% bacteria (355,456,488 sequences), 3.80% eukaryotes (14,312,486 sequences), 1.41% archaea (5,292,432 sequences), and 0.16% viruses (601,612 sequences; Figure S1).
The microbial domains were compared between the rumen microbiomes of the two MPY groups, and archaea were significantly different between the two groups (adjusted P < 0.01, Fig. 1a). The permutational multivariate analysis of variance (PERMANOVA) showed that both bacteria and archaea were significantly different (adjusted P < 0.01), while eukaryota and viruses were not different (adjusted P > 0.05) between the two groups (Table S3). The principal coordinate analysis (PCoA) showed separations between the two MPY groups based on bacterial (Fig. 1b) and archaeal species (Fig. 1c), while no separation was found based on eukaryotic or viral species (Figure S2). Thus, the downstream comparison of rumen microbial taxa between the two groups of animals was focused only on bacteria and archaea.
Compositional profiles of the rumen microbiome and taxonomic differences between the HH and LL cows
The dominant bacterial phyla included Bacteroidetes (55.98 ± 1.02%), Firmicutes (27.32 ± 1.14%), and Proteobacteria (7.32 ± 1.57%); the dominant bacterial genus was Prevotella (41.95 ± 0.85%), followed by Bacteroides (7.29 ± 0.31%), unclassified Lachnospiraceae (3.29 ± 0.18%), and Clostridium (2.99 ± 0.19%); and the dominant bacterial species included Prevotella sp. FD3004 (7.01 ± 0.37%), Prevotella ruminicola (4.64 ± 0.21%), Prevotella brevis (3.83 ± 0.21%), Prevotella sp. MA2016 (2.77 ± 0.15%), and Prevotella bryantii (2.57 ± 0.44%). For differential abundance comparison analysis at the phylum level, the abundance of Bacteroidetes was significantly higher in the rumen of LL cows, while that of Proteobacteria was significantly higher in the rumen of HH cows (adjusted P < 0.05, Figure S3). At the species level, 15 species, including 11 Prevotella sp., one Succinimonas sp., one Selenomonas sp. and one unclassified Bacteroidales exhibited significantly higher abundances in the rumen of HH animals (linear discriminant analysis [LDA] > 2, P < 0.05), while 23 species showed significant enrichment in the rumen of LL animals (LDA > 2, P < 0.05; Fig. 2a).
For the differential abundance comparison analysis of archaea, the abundance of the most abundant archaeal phylum, Euryarchaeota (99.01 ± 0.23%), was significantly higher in the rumen of LL cows (adjusted P < 0.01, Figure S4). At the genus level, the abundance of Methanobrevibacter, the most abundant archaeal genus (85.44 ± 2.41%), was significantly higher in the rumen of LL cows, while the abundances of other differential genera were all significantly higher in the rumen of HH cows (adjusted P < 0.05, Figure S4). At the species level, only the abundance of Methanobrevibacter millerae (22.10 ± 2.31%), the most abundant archaeal species, was significantly higher in the rumen of LL cows, while the abundances of the other differential species were all higher in the rumen of HH cows (LDA > 2, P < 0.05, Fig. 2b).
Functional profiles of the rumen microbiome and differential functions between the HH and LL cows
The functions of the rumen microbiome were determined by the Kyoto Encyclopedia of Genes and Genomes (KEGG) profiles and genes encoding CAZymes. For KEGG profiles, 158 endogenous third-level pathways were considered as rumen microbial metabolic pathways (Table S4). These pathways belonged to four first-level categories, including “Metabolism” (72.26 ± 0.46%), “Genetic information processing” (19.08 ± 0.12%), “Environment information processing” (4.42 ± 0.03%), and “Cellular processes” (4.24 ± 0.04%). At the second level, 20 categories were observed, with “Carbohydrate metabolism” (17.33 ± 0.10%), “Amino acid metabolism” (15.96 ± 0.11%), “Nucleotide metabolism” (9.82 ± 0.06%), “Replication and repair” (8.71 ± 0.06%), and “Energy metabolism” (8.07 ± 0.05%) being the most abundant. When the identified KEGG pathways were compared, a total of 13 third-level pathways, including two “Cellular processes” pathways, two “Genetic information processing” pathways, two “Environmental information processing” pathways, and seven “Metabolism” pathways, were significantly enriched in the rumen microbiomes of HH cows, while 18 pathways, including one “Genetic information processing pathway”, two “Cellular processes pathways” and 15 “Metabolism” pathways, were significantly enriched in the rumen of LL animals (LDA > 2 and P < 0.05; Fig. 3a). When the KEGG modules involved in the above differential third-level pathways were compared, 24 HH-enriched and 19 LL-enriched modules were identified (Fig. 3b). Regarding carbohydrate metabolism and energy metabolism, only two downstream functions (ko00290 and M00019, converting pyruvate to valine and isoleucine) were enriched in the rumen of HH cows (Fig. 4a). Four pathways and two modules were significantly enriched in the rumen of LL animals (LDA > 2, P < 0.05). The four pathways included “Glycolysis” (ko00010), “Starch and sucrose metabolism” (ko00500), “Galactose metabolism” (ko00052), and “Methane metabolism” (ko00680). The two modules were “Glycolysis” (M00001) and “Galactose degradation” (M00632). The downstream function of “Valine, leucine and isoleucine degradation” (ko00280) was also enriched in the rumen of the LL cows (Fig. 4b).
For CAZyme profiles, a total of 313 genes encoding CAZymes were identified (Table S5), including 8 auxiliary activities (AAs), 79 carbohydrate-binding modules (CBMs), 16 carbohydrate esterases (CEs), 115 glycoside hydrolases (GHs), 74 glycosyltransferase (GTs), and 21 polysaccharide lyases (PLs). Among them, genes encoding GT2 (8.64 ± 0.04%) were the most dominant, followed by those encoding CE1 (4.66 ± 0.02%), GT4 (4.34 ± 0.02%), GH2 (4.30 ± 0.02%), and GH3 (4.16 ± 0.02%). Among the genes encoding CAZymes involved in deconstructing carbohydrates (including cellulose, hemicellulose, starch, protein, and lignin), 18 were enriched in the rumen of HH cows (15 GH, 1 CE, 1 PL, and 1 AA), while 34 were enriched in the rumen of LL cows (27 GH, 4 CE, 2 PL, and 1 AA; Figure S5). Among the GTs (involved in carbohydrate synthesis), 11 were enriched in the rumen of HH cows, while two were enriched in the rumen of LL cows. Regarding the CBMs, the noncatalytic CAZymes that are involved in the degradation of complex carbohydrates, three were enriched in the rumen of HH cows, while 19 were enriched in the rumen of LL cows.
Associations between microbial species and microbial functions
As protein content is one of the determining measurements of MPY, we further focused on the functions of amino acid metabolism in the rumen microbiome. We found two important pathways involved in branched-chain amino acid (BCAA) metabolism (Fig. 5a), which were “valine, leucine and isoleucine biosynthesis” (ko00290, enriched in the rumen of HH cows) and “valine, leucine and isoleucine degradation” (ko00280, enriched in the rumen of LL cows), and these pathways showed a converse enrichment between the HH and LL groups (LDA > 2, P < 0.05; Figure S6). The abundances of genes encoding enzymes involved in these two pathways were also compared, showing that the abundances of genes encoding enzymes involved in BCAA biosynthesis were all significantly enriched in the rumen of HH cows, while the abundances of genes encoding enzymes involved in BCAA degradation were all significantly higher in the rumen of LL cows (adjusted P < 0.05; Fig. 5a and Figure S7). A Spearman’s rank correlation network between bacterial species and those two BCAA pathways was then created to explore how rumen bacterial species could affect the microbial BCAA functions. A total of 24 species showed significant relationships with two BCAA pathways (R > 0.50 and P < 0.05), 13 showing positive relationships with a BCAA biosynthesis pathway (ko00290). Among those 13 positive relationships between bacterial species and ko00290, the strongest (R > 0.65 and P < 0.05) were detected for five Prevotella species, including P. multisaccharivorax, P. histicola, P. maculosa, P. buccae, and P. albensis (Fig. 5b).
Rumen metabolome and serum metabolome
A total of 263 compounds were identified in the rumen metabolome. After t test and variable importance in projection (VIP) filtering for the relative concentrations of rumen metabolites, 25 metabolites were significantly different between the two MPY groups, all of which were significantly higher in the rumen of HH cows (P < 0.05, VIP > 1; Fig. 6a). Metabolic pathway analysis (MetPA) based on these 25 significantly different rumen metabolites revealed the enrichment of 10 pathways (Fig. 6b), with “vitamin B6 metabolism”, “glycerolipid metabolism”, and “beta-alanine metabolism” being the significantly different pathways (Benjamini-Hochberg false discovery rate [FDR] < 0.01, pathway impact > 0.1). The rumen metabolome was also used for phenotype (MPY) association analysis, and 126 MPY-associated metabotypes (metabolites that were significantly associated with MPY) were detected (see details in Methods, Table S6). The 126 MPY-metabotypes were used for PERMANOVA analysis; 106 of the MPY-metabotypes (all were MPY-positive metabotypes) were correlated with alterations in the rumen microbiome (adjusted P < 0.05; Table S6). These 106 MPY-positive metabotypes were considered as rumen microbiome-responsive metabotypes, which were then found to be significantly associated with 43 microbial modules (P < 0.05; Figure S8). In addition to the relative concentrations of ruminal small molecules that were identified by metabolomics, the absolute concentrations of the total volatile fatty acids (VFAs), propionate, valerate, and isovalerate (Fig. 6c, d) were quantified and were significantly higher in the HH cows (P < 0.05).
For the serum metabolome, we analyzed the 176 compounds identified in our previous study [6]. The comparison analysis revealed that the relative concentrations of 19 metabolites were significantly higher in the serum of HH cows, and the relative concentrations of 12 metabolites were significantly higher in the serum of LL cows (P < 0.05, VIP > 1; Fig. 7a). These 31 significantly different concentrations of metabolites were then used for MetPA analysis, revealing the enrichment of 12 pathways (Fig. 7b), with “glycine, serine, and threonine metabolism”, “nicotinate and nicotinamide metabolism”, and “sphingolipid metabolism”, “aminoacyl-tRNA biosynthesis” and “valine, leucine and isoleucine degradation” being the significantly different pathways (FDR < 0.01, pathway impact > 0.1). The serum metabolome was then identified as MPY-positive metabotypes (21 metabolites) or MPY-negative metabotypes (14 metabolites) using the phenotype (MPY) association analysis as stated above (see details in Methods, Table S7).
To identify whether the MPY-associated metabolites in rumen could be related to those in the serum, we compared the rumen and serum metabolites, including the significantly different metabolites between two MPY groups, MPY-positive metabotypes and MPY-negative metabotypes (Figure S9). A Venn diagram of differential metabolites revealed that a fatty acid, named lauric acid, was shared by the rumen and serum. For the differential metabolite-enriched pathways, three pathways were common in both the rumen and serum of HH cows, including “pyrimidine metabolism”, “glycerolipid metabolism”, and “starch and sucrose metabolism”. The Venn diagram of MPY-associated metabotypes showed that “Arginine and proline metabolism”, “Aminoacyl tRNA biosynthesis”, and “Purine metabolism” were shared by both rumen and serum MPY-positive metabotypes.
Relationships between the rumen microbiome, rumen metabolome and serum metabolome, and their explainabilities for MPY
Spearman’s rank correlations between the rumen microbiota and rumen metabolites were assessed, with the results revealing 65 significant correlations (R > 0.50, P < 0.05; Fig. 8a). Among the 65 correlations, positive correlations existed between mainly 11 Prevotella species (P. albensis, P. maculosa, P. timonensis, P. histicola, P. denticola, P. buccae, P. paludivivens, P. multisaccharivorax, P. corporis, P. bryantii, and P. oralis) and amino acids, peptides, proteins and organic chemicals (0.50 < R < 0.82, P < 0.05). Spearman’s rank correlation network showed 22 relationships between the rumen microbiota and MPY-associated metabotypes (Fig. 8c). Among the 22 correlations, nine Prevotella species (P. maculosa, P. histicola, P. denticola, P. buccae, P. paludivivens, P. multisaccharivorax, P. corporis, P. bryantii, P. oralis) also exihibited correlations with most of the MPY-associated metabotypes, including metabolites involved in glutathione, phenylalanine, starch, sucrose, and galactose metabolism.
To identify the potential rumen microbiome-host metabolic interactions, Spearman’s rank correlations between the rumen microbiota and serum metabolites were performed (Fig. 8b). Fewer relationships existed compared to the relationships identified between the rumen microbiota and rumen metabolites. The relationships between the rumen microbiota and serum MPY-associated metabotypes showed that seven Prevotella species were positively correlated with metabotypes involved in the metabolism of several amino acids, including glycine, serine, threonine, alanine, aspartate, glutamate, cysteine, and methionine (Fig. 8d).
The proportions of variation in MPY due to rumen microbial composition, microbial functions, rumen metabolites, and serum metabolites were estimated using linear mixed effect model (see Methods). The MPY variation explained by the rumen microbial composition, microbial functions, rumen metabolome, and host serum metabolome were 17.81, 21.56, 29.76, and 26.78%, respectively (Fig. 9).