Potential Role of the Rumen Microbiome in Modulating Milk Protein and Fat in Dairy Cow Using Microgenomic Sequencing


 Background

The rumen is the main digestive and absorption organ of dairy cows that contains abundant microorganisms and effectively utilizes human-indigestible plant mass. Investigation on microbiome in the rumen from lactating dairy cows using metagenomic sequencing is reasonable for identifying ruminal microorganisms that contribute to milk composition traits.
Results

We used the Illumina HiSeq platform to generate the rumen microbiome of the six lactating Holstein cows with extremely high and low milk protein and fat percentages (high and low groups of PP and FP). In total, 6977 microorganism species were detected in which Bacteroidetes (51.4%) and Prevotella (38.48%) was the most predominant phylum and genus, respectively. Between high and low groups, we observed significantly differential microorganism abundances in genus and species levels. By performing LEfSe and Metastats analyses, we identified 38 top abundant species displaying differential richness between two groups in common (LDA > 3, p < 0.05, q = 0.037 ~ 0.048), in which Prevotella accounted for 68.8% of the species with higher abundance in high group. Function annotation with KEGG, eggnog and CAZy databases showed the species with significantly higher abundance in high group were enriched in carbohydrate, amino acid, pyruvate, insulin and lipid metabolism and transportation, indicating their higher capability of digesting feed and subsequently providing substrate for milk composition synthesis in mammary gland. In addition, a kind of anaerobic fungi, Neocallimastix californiae, was identified in high group that could coexist with rumen microbes and promote cellulose digestion.
Conclusion

This study investigated the rumen microbiome in lactating Holstein cows using metagenomic sequencing. Significant differential bacterial richness were observed between the cows with extremely high and low PP and FP. Function annotation showed the abundant species in high group were involved in carbohydrate, amino acid, pyruvate, insulin and lipid metabolism and transportation, indicating the significant correlation between rumen microbiota and milk compositions formation in dairy cattle.


Introduction
Cow's milk acts as an abundant source of high-quality proteins, minerals and trace elements and is generally considered balanced and nutritive foods, being frequently included as important components of a healthy diet [1]. The rumen is the main digestive and absorption organ of dairy cows and it can effectively use grass, food and non-protein nitrogen to produce milk [2][3][4]. The rumen contains abundant microorganisms, comprising bacteria, protozoa and fungi, and is like a large anaerobic fermenter [2,5,6], which digested complex brous substrates into fermentable sugars, ultimately, were fermented by rumen bacteria and converted primarily into volatile fatty acids (VFA) [7]. Approximately 90% of the VFA produced in rumen are acetate, butyrate, propionate that are absorbed into the blood through the rumen wall and transported to the liver. Subsequently, the liver transports cholesterol to the mammary gland for lipid synthesis in the form of lipoproteins [8][9][10][11][12][13][14]. As for lactose synthesis, glucose from carbohydrate digestion in rumen and gluconeogenesis in liver was transported into the mammary gland through blood circulation and converted into lactose [18]. Amino acids produced by the breakdown of rumen microbiota are absorbed into the blood in small intestine and converted into free amino acids (FAAs) that are subsequently owed into mammary gland where FAAs were synthesized into milk protein [19].
Metagenomics is an analysis of microbial communities at particular habitat by using high-throughput sequencing without necessary requirement of laboratory isolation and culture of individual strains [15,16]. It has been widely used to study microbial diversity and metabolic capabilities of microbiome in different ecological niches, fermented food, waste-water treatment facilities, and gastrointestinal tracts in human and animals [17][18][19][20][21][22][23]. In dairy cattle, rumen microbes played a vital role in the decomposition of plant lignocellulosic matter [24]. Jewell et al. found the ruminal bacterial community is dynamic in terms of membership and diversity and speci c members are associated with milk production e ciency over two lactation cycles [25]. Jami et al. reported that rumen microbial diversity were closely related to milk yield and composition in which the ratio of the Firmicutes to the Bacteroides was correlated with milk fat yield [26]. So far, most researches on rumen microbiota in dairy cows mainly focused on the phylum and genus of rumen bacteria and their effects on milk production, feed conversion ratio and methane emissions mainly through 16S rRNA sequencing. From species level, regulation roles of rumen microbiota on milk composition traits has been still under investigated in dairy cattle.
In the present study, the purpose was to investigate the rumen microorganism constitutions of the lactating Holstein cows who exhibited extremely high and low milk PP and FP with high-throughput metagenomics technology, and identify the microbiota with signi cantly differential abundance between the high and low groups of PP and FP. Further, functional enrichment analysis were performed to elucidate the potential correlation of rumen microbiome with milk protein and fat biosynthesis.

Animals and sample collection
We selected six healthy lactating Chinese Holstein cows from the Beijing Sanyuanlvhe Dairy Farming Center, which were in the almost same lactation period (230, 272, 230, 281, 233, and 196 days in milk of the secondlactation, respectively) and were fed the same total mixed rations. The cows were chosen among more than 4000 lactating Holstein cows based on their monthly test-day milk PP and FP records for the rst and second lactations, which were provided by the Dairy Data Center of China (https://www.holstein.org.cn/). The six cows were divided into two groups with extremes of the average phenotypic values for PP and FP: three cows (high group) had high PP (3.7%, 3.6 and 3.8%) and FP (3.9%, 4.3% and 4.5%); the other three cows (low group) showed low PP (3.0%, 3.1% and 3.2%) and FP (3.1%, 3.2% and 2.9%).
Rumen uid sample was collected from the rumen of each cow within 15 minutes after slaughter in the same slaughterhouse. All samples were collected in the cryogenic vials and stored in liquid nitrogen for further analysis. All animal works were carried out in accordance with the approved guidelines for the care and use of experimental animals from the Ministry of Agriculture of China. The protocol was specially approved by the Animal Welfare Committee of China Agricultural University (permit number: DK996).

Metagenomic data processing
Quality control and metagenome assembly The raw data from Illumina HiSeq sequencing platform were preprocessed using Readfq (https://github.com/cj elds/readfq) to obtain clean data for subsequent analysis. The reads which contain low quality bases less than length of 40 bp, the N base that reached 10 bp and overlap above length of 15 bp were removed. In order to avoid the possibility of host pollution, the clean data was compared against the bovine reference genome UMD3.1.69 to lter out the reads that may be from the host with the parameters [27,28] as follows: --end-to-end, --sensitive, -I 200 and -X 400. The ltered reads were assembled using SOAPdenovo 2.04 software (http://soap.genomics.org.cn/soapdenovo). The scaffolds longer than 500 bp were extracted for subsequently genetic prediction and annotation [29][30][31].
Comparative analysis on microorganism abundance between high and low groups Krona analysis was performed to display relative abundance and abundance cluster heat map. Principal component analysis (PCA) [38] (R ade4 package, Version 2.15.3) and NMDS [39] (R vegan package, Version 2.15.3) decrease-dimension analysis were performed based on the abundance of each taxonomic hierarchy. Differences of microorganism abundance within groups and between groups were tested by Anosim analysis (R vegan package, Version 2.15.3). LEfSe and Metastats analysis were used to look for the species with different abundance between groups. LEfSe software was conducted with the default threshold LDA score of 3 [40]. For Metastats analysis, permutation test between groups was used for each taxonomy and get the raw P value, then Benjamini and Hochberg False Discovery Rate were used to correct P value and acquire q value [41].
Functional annotation analysis DIAMOND software (V0.9.9) was adopted to blast the unigenes to functional databases to analyze the functions of rumen microorganisms, including Kyoto Encyclopedia of Gene and Genome (KEGG, http://www.kegg.jp/kegg/), Evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG, http://eggnogdb.embl.de/#/app/home) and Carbohydrate-Active enzymes Database (CAZy, http://www.cazy.org/). The best Blast Hit was used for subsequent analysis. LEfSe and Metastats analysis were used to look for the different functions between the high and low groups.

Sequencing of the rumen microbiota
With Illumina HiSeq PE-150 platform, we sequenced the microbiota DNA of each rumen uid sample from the six lactating Holstein cows with extremely high and low milk PP and FP with three individuals in each group (high and low groups), and 76 Gbp clean data in total were generated after quality control (Supplementary Table 1). Through prediction and de-redundancy with MetaGeneMark and CD-HIT softwares, we identi ed a total of 1,078,009 ORFs with a total length of 667.24 Mbp in the scaffolds longer than 500 bp, from which 278,274 complete genes were annotated, accounting for 25.81% (Supplementary Table 2).
Differential abundance of the rumen microbiome between high and low groups with PP and FP With PCA analysis, we observed that the six cows with extremely high and low milk PP and FP were obviously separated into two clusters based on the rumen microbiome abundance at genus and species levels( Fig. 1A-B), in whichprincipal coordinate 1 accounted for 55.31%, 60.68% and principal coordinate 2 accounted for 14.21%, 12.81% of the total variation in two groups, respectively. Differential abundance comparison analysis showed there were signi cant differences on the abundances of rumen microorganisms between the cows in high and low groups (R = 0.889, P = 0.1; R = 1, P = 0.1; Fig. 2). Among the top ten abundant genera, Prevotella was the most genus with proportion of 42.8% and 34.2%% in high and low groups, respectively; Fibrobacter accounted for 3.13% and 2.69%, and Bacteroides accounted for 3.02% and 2.25% in two groups, respectively; the abundances of another seven genera were 0.44% ~ 0.8% (Supplementary Fig. 2E).To identify differentially abundant species, we performed LEfSe analysis and detected 40 kinds of microbiota that had signi cantly differential abundances in the rumen of cows between high and low groups (LDA > 3, p < 0.05, Fig. 3 and Fig. 4), including 8Prevotella sp., 2 unclassi ed_Bacteroidales, 2 unclassi ed_Bacteria and 1 eukaryota of Neocallimastix californiae exhibited signi cantly higher abundances in high group, while 27 species showed signi cant enrichment in low group. Simultaneously, through Metastats analysis, 2797 microorganisms showed signi cantly differential abundances between two groups (p = 1.48E-05 ~ 0.049; q = 0.037 ~ 0.048), in which 38 species were commonly identi ed with both methods (Supplementary Table 4). Out of these, the top 35 abundant microorganisms(abundance >0.02%; Fig. 5) included 16 species from phylum Bacteroidetes in high group and 19 bacteria from phylumBacteroidetes, Firmicutes, Fibrobacteres and Proteobacteriain low group; of note, 11 of 16 species in high group werePrevotella which functions as carbohydrate degradation.
Functional enrichment of the rumen microbiome with differential abundances between high and low groups with PP and FP To understand the potential correlation of rumen microbiome with milk protein and fat traits, the functions of the differentially abundant rumen microbiome were determined by the Kyoto Encyclopedia of Gene and Genome (KEGG), Evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG) and Carbohydrate-Active enzymes Database (CAZy) functional annotation.
Regarding KEGG, 108 pathways were signi cantly enriched in third-level pathways (LDA > 2.0 , P < 0.05, Fig. 6 and Supplementary Table 5), including 26 involved for the species with signi cantly high abundance in high group such as ve "Metabolism", six "Cellular Processes", six "Environmental Information Processing", seven "Organismal Systems", and two "Human Diseases" pathways, while 42 pathways were enriched in low group. At the second level, 28 categories were observed, in which the most abundant pathways included "Glycan biosynthesis and metabolism", "Carbohydrate metabolism", "Amino acid metabolism", "Metabolism of cofactors and vitamins" and "Signal transduction". For KEGG orthologies, 42 signi cant KOs were determined (LDA >2,P < 0.05, Fig. 7 and Supplementary Table 6) including 6 KOs (K01190, K05349, K01811, K03737, K01006 and K00688) related to carbohydrate, energy and pyruvate metabolisms were enriched in the rumen of cows in high group, and another 17 KOs with higher abundances in low group such as DNA-damage-inducible protein J (K07473), ATP-dependent DNA helicase RecG (K03655), and agellin (K02406).
Based on the eggNOG database, seven signi cantly enriched pathways for the rumen microorganisms in high group were involved in carbohydrate transport and metabolism, signal transduction metabolism, energy production and conversion, translation ribosomal structure and biogenesis, lipid transport and metabolism, secondary metabolites biosynthesis and transport and catabolism, and posttranslational modi cation; four signi cantly enriched pathways in the low group included replication recombination and repair, animo acid transport and metabolism, coenzyme transport and metabolism, and cell motility (LDA >3,P < 0.05, Fig. 8 and Supplementary Table 7).
Based on the CAZyme database, a total of 215 genes encoding CAZymes were identi ed (Supplementary Table 8), including 96 glycoside hydrolases (GHs), 46 carbohydrate-binding modules (CBMs), 45 glycosyltransferase (GTs), 12 carbohydrate esterases (CEs), 11 polysaccharide lyases (PLs), and four auxiliary activities (AAs). Among these, GH51, GH97, GH31, GH2, GH3 and GH43 family glycoside hydrolases, and CE1 were signi cantly enriched in high group that were involved in deconstructing carbohydrates such as cellulose, hemicellulose and starch (LDA > 3, P < 0.05, Fig. 9 and Supplementary  Table 9) indicating higher digestive ability for feedings of the rumen of cows in high group. In summary, Prevotella bacteria was found signi cantly higher enriched in high group, which was capable of digesting brous hemicellulose starch, produced lactose and galactose, and then converted them to pyruvate.
Pyruvate converted into VFA and amino acid, then transported through the blood to the mammary epithelium to make milk fat and protein (Fig. 10).

Discussion
In this study, we investigated the rumen microbiome from the milking Holstein cows in mid-lactation with extremely high and low milk protein percentage and fat percentages using metagenomics technology, and identi ed differential abundances of the microorganisms at phylum, genus and species levels. It was found that signi cantly differential species were signi cantly associated with milk protein and fat.
Here, we detected 6977 microorganisms in the rumen of cows, the top three most abundant phylum out of them were Bacteroidetes, Firmicutes and Proteobacteria, which was consistent with the discovers in previous studies [42][43][44]. In addition, we found the predominant anaerobic fungi was Neocallimastigomycota that was consistent with the result reported by Zhang et al. [43]. However, Huang et al. found Actinobacteria was one of the most abundant predominant phylum besides Bacteroidetes and Firmicutes in the rumen of lactating Chinese Holstein cows with 16S rRNA sequencing [45]. This was likely due to the different feeding management and environment of the farms where cows were collected from that in this study. At the genera level, we identi ed that Prevotella was the most abundant genus in the rumen of lactating cows with high protein trait. This is coincident with the previous report by Xue et al. [44].
We observed 2798 microbial species that showed signi cantly differential abundances in the rumen uid between the lactating cows with extremely high and low milk protein percentage and fat percentages. Based on KEGG, eggNOG and CAZy databases, it was found that such species were mainly involved in amino acid, carbohydrate, pyruvate and lipid metabolism and transportation. Especially, majority of the most abundant species in high group belonged to the Prevotella genus. It was well known that Prevotella genus couldmetabolize cellulose and starch, its fermentation products include acetate, butyrate and propionate [46]. The previous studies indicated that higher Prevotella abundance in the rumen accelerated the decomposition of fructose and starch [47] and could produce propionate or succinate and acetate [48].
Besides, some studies found that Prevotella was also polysaccharide-degrading genus, had the ability to degrade mucin and plant carbohydrates [49,50]. Previous research had shown that the P. bryantii 25A in Prevotella could increase ruminal fermentation products and milk fat concentration after supplied into early lactating cows [51]. Our results were consistent with these reports. Thus, since rumen uid contained more Prevotella species in the high group, we speculated it could digest the cellulose, hemicellulose and starch well in the feed and degrade them into galactosidase and glucose. Under the action of microbiotas in the rumen, galactosidase and glucose converted to pyruvate. Then, part of pyruvate were broken down into VFA, which provided a su cient material basis for lactation in cows. The rumen VFA in dairy cows were basically acetate, propionate and butyrate. After passing through the rumen wall, acetate was transported to the liver via the blood [52] . Butyrate was conversion to β-hydroxybutyric acid (BHBA) and absorbed by the liver. In the liver, these substances were converted into cholesterol. Finally, cholesterol was transported to the mammary gland in the form of lipoproteins, involved in the synthesis of milk fat [53][54][55] . Propionate was the main precursor of glucose synthesis in ruminants, which was bene cial to the supply and transformation of energy in the body and provided energy for protein synthesis. The increase of propionate could also stimulate the secretion of insulin, the blood ow of the breast, to promote the synthesis of milk protein [56][57][58] . Amino acids produced in the rumen were used by the small intestine and converted to FAA. FAA was absorbed as it owed through the the mammary gland to synthesize milk protein. Additionally, Neocallimastix californiae, a kind of anaerobic fungi, was highly enriched in our high group which had been reported that it could reduce the inner tension of plant ber and make it easier to be degraded by rumen microbiomes [59]. In subsequent experiments, we made bacteria with signi cant differences into enzyme preparations and fed them to cows to test whether the bacteria could affect milking traits. This might warrant further studies in the future.
Interestingly, in top 35 abundant microorganisms, 16 species with high abundance in the high group belong to Bacteroidetes genera, and 11 of 16 species in high group werePrevotella which functions as carbohydrate degradation. While 19 species with high abundance in the low group were distributed in four genera. The high group, with high milk fat and protein, had the characteristics of lower diversity and higher advantage in genus content. As we all known, microbial diversity had an important effect on cow metabolism. And Shabat et al.' study found that lower richness of microbiome gene content and taxa was tightly linked to higher feed e ciency [60].The diversity and abundance of microbiomes, and mutualism between microbiomes and anaerobic bacteria were signi cantly associated with milking traits in dairy cow.

Conclusion
The present study obtained the rumen microbiome pro le of lactating Holstein cows using microgenomic sequencing, and identi ed 38 top abundant species that showed signi cantly differential richness between the cows with extremely high and low milk PP and FP. Function annotation indicated that the differential species with higher abundance in high group were associated with carbohydrate, amino acid, pyruvate, insulin and lipid metabolism and transportation, implying their higher capability of digesting feed and subsequently providing substrate for milk composition synthesis in mammary gland.