Abstract The present study was conducted to analyze the correlation between the milk fat content of Binglangjiang buffaloes and their microbial and host metabolites. The 10 buffaloes with the highest milk fat content (HF, 5.60 ± 0.61%) and the 10 with the lowest milk fat content (LF, 1.49 ± 0.13%) were selected. Their rumen fluid and plasma were collected for rumen microbiota and metabolome analysis. The results showed that the rumen bacteria abundance of Synergistota, Quinella, Selenomonas, and Fretibacterium was significantly higher in the HF buffaloes. The abundance of 14 rumen fungi, including Candida, Talaromyces, Cyrenella, and Stilbella, was significantly higher in the HF buffaloes. The analysis of the metabolites in the rumen and plasma showed that several metabolites differed between the HF and LF buffaloes. A total of 68 and 42 differential metabolites were identified in the rumen and plasma, respectively. By clustering these differential metabolites, most of those clustered in the HF group were lipid and lipid-like molecules such as secoeremopetasitolide B, lucidenic acid J LysoPE (0:0/18:2 (9Z, 12Z)), and 5-tetradecenoic acid. Spearman’s rank correlations showed that Quinella, Fretibacterium, Selenomonas, Cyrenella, and Stilbella were significantly positively correlated with the metabolites of the lipids and lipid-like molecules in the rumen and plasma. The results suggest that rumen microbiota such as Quinella, Fretibacterium, Selenomonas, and Cyrenella may regulate milk fat synthesis by influencing the lipid metabolites in the rumen and plasma. In addition, the combined analysis of the rumen microbiota and host metabolites may provide a fundamental understanding of the role of the microbiota and host in regulating milk fat synthesis. 1. Introduction Milk is rich in a variety of essential nutrients, including fat, protein, and various trace elements that are important for healthy human development, especially many kinds of unsaturated fatty acids. Milk fat is a milk component with high nutritional value, which is the most important determinant of milk quality. A high fat content is an important characteristic of buffalo milk [[40]1]. In addition, milk fat is an important source of essential fatty acids, phospholipids, and fat-soluble vitamins, providing the body with a large amount of nutrients and energy [[41]2,[42]3]. Buffalo milk, as a best-selling dairy raw material in domestic and international markets, is milky and thick, with a high content of nutrients. Its content of milk fat, milk protein, total solids, and minerals is higher than that of Holstein milk; therefore, the development of buffalo milk is particularly important [[43]4,[44]5]. Ruminants have been domesticated for more than 10,000 years [[45]6] and rely on fermentation by the rumen microbial community to convert complex polysaccharides, which make up a major portion of the plant, into available nutritive substances [[46]7,[47]8]. Furthermore, the rumen microbiota is one of the key factors affecting the formation of milk fat precursors. Changes in milk fat production and milk fatty acid composition are usually caused by alterations in rumen fermentation and biohydrogenation pathways [[48]9]. Previous research has indicated that there are interactions between rumen microbes and the mammary gland. Gastrointestinal microbial metabolites can enter the bloodstream through the intestinal epithelium into the mammary gland, which then influences the composition of milk secreted by the mammary gland [[49]10]. It seems that rumen metabolites, as important precursors for milk fat production, have a direct impact on milk fat percentage in dairy cows [[50]11]. Several studies have shown that the difference in milk fat percentage may be caused by the interactions between gut microorganisms and their metabolites, especially Firmicutes myristic acid and Proteus cholin [[51]12]. However, the levels of metabolites in rumen fluid do not directly reflect the levels of relevant metabolites in the blood [[52]13]. Blood is a commonly used biofluid in metabolomics analyses, providing metabolites from all organs [[53]14]. Blood metabolites can be used as physiological biomarkers to reflect metabolic performance and dysfunction [[54]15]. Meanwhile, some typical metabolites are highly correlated with specific rumen microbiota, suggesting a functional correlation between the microbiome and its associated metabolites [[55]16]. With the development of biotechnology, more researchers are using macro-genomics, macro-transcriptomics, and proteomics to explore the role of the rumen microbiota and their involvement in milk fat synthesis, milk production, and lipid metabolism [[56]12,[57]17]. However, this relationship has not yet been studied in Binglangjiang buffaloes. This buffalo species can be found across all continents and boasts abundant resources. It is classified into two subspecies: river buffaloes (50 chromosomes), primarily inhabiting India, Pakistan, and Mediterranean coastal countries and regions, and swamp buffaloes (48 chromosomes), which are mainly distributed in China and southeast Asian countries [[58]18]. The Binglangjiang buffalo is the only species of river buffalo that has been discovered in China. A previous study found that the fat and protein contents in the milk of Binglangjiang buffaloes was about 2.1 and 1.2 times higher than that of Holstein milk, suggesting that their milk has excellent qualities [[59]19]. However, the individual Binglangjiang buffalo vary greatly in terms of milk composition and reproductive performance, and their daily milk production is still at a low level compared to other countries. In this study, we speculate that the milk fat content of cows fed the same diet may be influenced by rumen microbiota composition and its association with the host metabolism. Therefore, we investigated the correlation between the rumen microbiota and the metabolomics of rumen and plasma in high fat content and low fat content Binglangjiang buffaloes. 2. Materials and Methods 2.1. Animals and Management All the experimental protocols were approved by the Animal Care Committee at the Faculty of Animal Science and Technology, Yunnan Agricultural University (Kunming, P. R. China, approval No.: YNAU202205020). In this experiment, we randomly selected 75 healthy mid-lactation Binglangjiang buffaloes with similar lactation days (120–140 days) and similar body condition (body condition score from 2.5 to 3). Then, they were kept in a commercial buffalo farm in Tengchong City, Yunnan Province. The buffaloes were kept in individually tethered stalls in a barn with good ventilation and were fed 2 times daily at 08:00 h and 16:00 h. The buffaloes were fed the same diet ([60]Tables S1 and S2) with a concentrate-to-forage ratio of 35:65 (DM basis) and had free access to water. Then, we continuously measured the composition of the milk for 6 weeks once a week, and the 10 highest milk fat content buffaloes (HF, 5.60 ± 0.61%) and 10 lowest milk fat content buffaloes (LF, 1.49 ± 0.13%) were selected accordingly ([61]Figure S1). 2.2. Sample Collection and Measurement Milk samples were collected at 06:30 with a bronopol tablet (milk preservative, D & F Control Systems, San Ramon, CA) and then stored at 4 °C before the infrared analysis of milk components and somatic cells using a spectrophotometer (Foss-4000, Foss, Hillerød, Denmark). Milk production was recorded for 7 consecutive days using a mobile milking machine, and feed intake was recorded prior to formal sample collection. Blood samples were collected from the jugular vein, and plasma was separated and then stored at −80 °C for metabolomic analysis. Rumen contents were sampled using oral stomach tubes and were used to measure microbiota and metabolites. The milk, blood and rumen samples were collected on the same day after milk production was recorded. The pretreatment and determination of rumen fluid VFAs were carried out with reference to Hu et al. [[62]20]. 2.3. 16S rDNA and ITS Sequencing and Data Processing Rumen microorganisms DNA were extracted from 2 mL rumen fluid by the kit method (Omega Bio-tek, Norcross, GA, USA). Before PCR amplification, we used 1.0% agarose gel electrophoresis and NanoDrop2000 to detect the purity and concentration of DNA, and we took an appropriate amount of sample from the centrifuge tube for dilution. Diluted DNA was used as a template. Specific primers were selected to amplify the V3 + V4 region of 16SrDNA. The fungus was subjected to ITS sequencing, and the amplified region was ITS1F-ITS2R. The primer sequences and PCR amplification conditions were referred to Yu et al. [[63]21]. We used FASTP (v0.20.0) software to control the quality of the original sequencing sequence [[64]22] and used FLASH (v1.2.7) software to splice [[65]23]. Using UPARSE software (v7.1), OTU clustering was performed on sequences according to 97% similarity, and chimeras were eliminated [[66]24]. RDP classifier (v2.2) was used to classify and annotate each sequence [[67]25], compared with the Silva 16S rDNA database (v138), and we set the comparison threshold to 70%. Then, we used Mothur software (v1.30.2) for alpha diversity analysis. Beta diversity was analyzed by principal component analysis (PCA) based on the Bray–Curtis distance in R language (v 3.3.1), while the 16S and ITS sequence data were compared between the two groups (HF vs. LF) using the Wilcoxon rank-sum test with the p-value < 0.05 considered significantly different. 2.4. Analysis of Rumen and Plasma Metabolome Rumen and plasma metabolites were extracted according to Guo et al. [[68]26]. First, 100 μL liquid sample was added to a 1.5 mL centrifuge tube with 400 μL solution (acetonitrile: methanol = 1:1 (v:v)) containing 0.02 mg/mL internal standard (L-2-chlorophenylalanine) to extract metabolites. The samples were mixed by vortex for 30 s and low-temperature sonicated for 30 min (5 °C, 40 KHz). The samples were placed at −20 °C for 30 min to precipitate the proteins. Then, the samples were centrifuged for 15 min (4 °C, 13,000× g). The supernatant was removed and blown dry under nitrogen. The sample was then re-solubilized with 100 µL solution (acetonitrile: water = 1:1) and extracted by low-temperature ultrasonication for 5 min (5 °C, 40 KHz), which was followed by centrifugation at 13,000× g and 4 °C for 10 min. The supernatant was transferred to sample vials for LC-MS/MS analysis. As a part of the system conditioning and quality control process, a pooled quality control sample (QC) was prepared by mixing equal volumes of all samples. The LC-MS analysis was performed on the Thermo UHPLC-Q system of Majorbio bio-pharma Technology Co. (Shanghai, China). A Thermo UHPLC-Q Exactive HF-X mass spectrometer was used to collect mass spectrum data, which was divided into two working modes: positive mode and negative mode. The data-related acquisition (DDA) mode was used for data acquisition. It was tested in the quality range of 70–1050 m/z. The LC/MS raw data were pre-processed using Progenesis QI (Waters Corporation, Milford, CT, USA) software (version 2.0). At the same time, we used the metabolite identification of a searchable database, which was the main database for the KEGG database ([69]http://www.genome.jp/kegg/, accessed: 20 February 2023). We uploaded a searchable database of the data matrix to the Majorbio cloud platform ([70]https://cloud.majorbio.com, accessed: 15 March 2020) for data analysis. First, the data matrix was preprocessed: at least 80% of the metabolic signatures detected in each set of samples were retained. After filtration, minimum metabolite values were estimated for specific samples whose metabolite levels were below the lower quantization limit, and various metabolic characteristics were normalized to sum. The response intensity of the sample quality spectrum peak was normalized by the sum normalization method. At the same time, the variables of QC samples with a relative standard deviation (RSD) of >30% were removed, and log10 was used for subsequent analysis. Then, the R software package “ropls” (version 1.6.2) was used for orthogonal least partial binary discriminant analysis (OPLS-DA) and interactive verification to evaluate the stability of the model. According to the p-values generated by the OPLS-DA model and Student’s t-test, metabolites with p < 0.05 were significantly different metabolites. Differential metabolites among the two groups were mapped into their biochemical pathways through metabolic enrichment and pathway analysis based on the KEGG database. Python packages “scipy.stats” ([71]https://docs.scipy.org/doc/scipy/, accessed: 27 March 2023) was used to perform enrichment analysis to obtain the most relevant biological pathways for experimental treatments. Correlation analysis between the rumen microbiota, rumen and plasma metabolome were performed using the Spearman ([72]https://cloud.majorbio.com/page/tools.html, accessed: 5 April 2023), any p-value below 0.05 was regarded as indicating a significant correlation. 3. Results 3.1. Characterization of Phenotypes In this study, 10 buffaloes with the highest milk fat content (HF) and 10 buffaloes with the lowest milk fat content (LF) were selected for 16S, ITS, rumen metabolome, and plasma metabolome analyses. Among the phenotypes, the milk fat and lactose content was significantly different between the HF and LF (p < 0.001) ([73]Table 1). Concentrations of rumen volatile fatty acids were quantified and found to be significantly higher (p < 0.05) in HF buffaloes for total volatile fatty acids (TVFAs) and acetate ([74]Table S3). Table 1. Physiological parameters of HF and LF buffaloes. Items Mean ± SEM p HF (n = 10) LF (n = 10) Milk yield (kg/d) 8.05 ± 0.31 7.03 ± 0.37 0.81 Milk fat (%) 5.60 ± 0.61 1.49 ± 0.13 <0.001 Milk protein (%) 4.85 ± 0.21 4.42 ± 0.10 0.076 Lactose (%) 5.24 ± 0.06 5.57 ± 0.08 0.003 Parity 2.10 ± 0.18 2.50 ± 0.17 0.12 Days in milk (DIM, d) 182.30 ± 13.36 170.40 ± 21.97 0.65 Body weight (kg) 505.00 ± 6.56 508.73 ± 7.02 0.37 Age 4.90 ± 0.23 5.10 ± 0.28 0.60 Dry matter intake (DMI, kg) 10.11 ± 0.16 9.76 ± 0.21 0.21 [75]Open in a new tab 3.2. Rumen Bacteria and Taxonomic Differences Between the HF and LF Buffaloes In 16S rDNA sequencing, a total of 960,813 sequences were generated, with 48,040 ±1523 sequences (mean ± standard error of the mean [SEM]) per sample, and the average sequences in group HF and LF were 48,050 and 48,031, respectively ([76]Table S4). The Chao, Shannon and coverage index of HF and LF buffaloes were not significantly different ([77]Figure 1A–C). The partial least squares discriminant analysis (PLS-DA) showed separations between the groups based on bacterial phyla ([78]Figure 1D) and genus ([79]Figure 1E). Figure 1. [80]Figure 1 [81]Open in a new tab Composition and differential bacteria in the rumen. (A) Shannon index. (B) Chao index. (C) Coverage index. (D) The partial least squares discriminant analysis at phylum level. (E) The partial least squares discriminant analysis at genus level. (F) Community structure at the phylum level. (G) Community structure at the genus level. 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 (21.20 ± 0.45%), which was followed by Rikenellaceae_ RC9_ gut_ group (10.05 ± 0.05%), norank_f__UCG-011 (6.34 ±0.49%) and Succiniclasticum (4.20 ± 0.34%) (p < 0.05, [82]Figure 1F,G). For differential abundance comparison analysis at the phylum level, the abundance of Synergistota was significantly higher in the rumen of HF buffaloes. At the genus level, Quinella, Selenomonas and Fretibacterium exhibited significantly higher abundances in the rumen of HF buffaloes (p < 0.05), while 12 genera showed significant enrichment in the rumen of LF buffaloes (p < 0.05; [83]Figure 2). Figure 2. [84]Figure 2 [85]Open in a new tab Differential bacteria. * means p < 0.05, ** means p < 0.01. (A) At phylum level. (B) At genus level. 3.3. Rumen Fungi and Taxonomic Differences Between the HF and LF Buffaloes In ITS sequencing, a total of 1,209,963 sequences were generated, with 60,498 ±2540 sequences (mean ± standard error of the mean [SEM]) per sample, and the average sequences in the HF and LF groups were 62,041 and 58,965, respectively ([86]Table S5). The Chao, Shannon and coverage indexes of HF and LF buffaloes were not significantly different ([87]Figure 3A–C). The partial least squares discriminant analysis (PLS-DA) showed separations between the groups based on fungi phyla ([88]Figure 3D) and genus ([89]Figure 3E). Figure 3. [90]Figure 3 [91]Open in a new tab Composition and differential fungi the rumen. (A) Shannon index. (B) Chao index. (C) Coverage index. (D) The partial least squares discriminant analysis at phylum level. (E) The partial least squares discriminant analysis at genus level. (F) Community structure at the phylum level. (G) Community structure at the genus level. The dominant fungal phyla included Ascomycota (55.98 ± 1.02%), Basidiomycota (27.32 ± 1.14%), and Neocallimastigomycota (7.32 ± 1.57%); the dominant genus was Orpinomyces (14.83 ± 0.11%), which was followed by Cutaneotrichosporon (8.32 ± 0.41%), norank_f_UCG-011 (6.34 ± 0.49%) and Trichosporon (7.49 ± 2.90%) ([92]Figure 3F,G). For differential abundance comparison analysis at the phylum level, there was no significant difference between HF and LF buffaloes (p > 0.05). At the genus level, the abundance of 14 rumen fungi, included Candida, Talaromyces, Cyrenella and Stilbella, was significantly higher in the rumen of HF buffaloes (p < 0.05), while Anaeromyces showed significant enrichment in the rumen of LF buffaloes (p < 0.05; [93]Figure 4). Figure 4. [94]Figure 4 [95]Open in a new tab Differential fungi at genus level. * means p < 0.05. 3.4. Rumen Metabolome To discriminate the metabolic profiles across groups, we performed clustering analyses based on orthogonal partial least square discriminant analysis (OPLS-DA). The rumen samples from distinct groups were largely separated according to the OPLS-DA plots ([96]Figure 5A). A total of 1451 compounds were identified in the rumen metabolome, and nine of these metabolites were identified in the HF group only, e.g., PE (15:0/16:0), LPC (18:1), and LPC (18:2) ([97]Figure S2). After t-test filtering for the relative concentrations of rumen metabolites, 68 metabolites were significantly different between the two groups, 40 of which were significantly higher in the rumen of HF buffaloes ([98]Figure 5B, [99]Table S6). By clustering these differential metabolites, most of those clustered in the HF group were lipid and lipid-like molecules such as secoeremopetasitolide B and lucidenic acid J ([100]Figure 5C,D). KEGG pathway enrichment analysis based on these 68 significantly different rumen metabolites revealed the enrichment of five pathways, which included “Bile secretion”, “Drug metabolism-cytochrome P450”, “Arginine and proline metabolism”, “Pentose and glucuronate interconversions”, and “ Phenylalanine metabolism”; among these, the “Pentose and glucuronate interconversions” pathway was significantly upregulated in the HF group (p < 0.05, DA score > 0.1) ([101]Figure 5E). Figure 5. [102]Figure 5 [103]Open in a new tab Rumen metabolome of HF and LF buffaloes. (A) PLS-DA analysis of rumen metabolome. (B) Volcanogram of rumen differential metabolome. (C) Comparison of rumen metabolome between HF and LF buffaloes visualized. (D) Cluster analysis of differential rumen metabolites, the left is a dendrogram of metabolite clustering, on the right is the name of the metabolite, and the closer two metabolite branches are to each other, the closer they are in terms of expression. (E) Differential abundance score map of KEGG pathway of rumen metabolome. A positive DA score indicates that the expression trend of all annotated differential metabolites in the pathway was upregulated in the HF group, the dots are distributed to the right of the center axis and the longer the line segment, the more the overall expression of the pathway tends to be upregulated. * means p < 0.05, *** means p < 0.001. 3.5. Plasma Metabolome Similar to the rumen, to discriminate the metabolic profiles across groups, we performed clustering analyses based on orthogonal partial least square discriminant analysis (OPLS-DA). The rumen samples from distinct groups were largely separated according to the OPLS-DA plots ([104]Figure 6A). A total of 770 compounds were identified in the plasma metabolome; among these, three of the metabolites were identified in the HF group only, namely alpha-crocetin glucosyl ester, 3,4,5-trihydroxy-6-[(3-oxo-1,7-diphenylheptan-2-yl)oxy]oxane-2-carboxyl ic acid, and isovalerylglutamic acid ([105]Figure S3). After t-test filtering for the relative concentrations of plasma metabolites, 42 metabolites were significantly different between the two groups, 36 of which were significantly higher in the rumen of HF buffaloes ([106]Figure 6B, [107]Table S7). By clustering these differential metabolites, most of those clustered in the HF group were lipid and lipid-like molecules such as LysoPE (0:0/18:2 (9Z, 12Z)), 5-tetradecenoic acid and hexadecanedioic acid ([108]Figure 6C,D). KEGG pathway enrichment analysis based on these 42 significantly different plasma metabolites revealed the enrichment of six pathways, which included “Pyrimidine metabolism”, “Lysine degradation”, “Pantothenate and CoA biosynthesis”, “Choline metabolism in cancer” and “Glycerophospholipid metabolism”; among these, the “Pantothenate and CoA biosynthesis” pathway was significantly upregulated in the HF buffaloes (p < 0.05, DA score > 0.1) ([109]Figure 6E). Figure 6. [110]Figure 6 [111]Open in a new tab Plasma metabolome of HF and LF buffaloes. (A) PLS-DA analysis of plasma metabolome. (B) Volcanogram of plasma differential metabolome. (C) Comparison of plasma metabolome between HF and LF buffaloes visualized. (D) Cluster analysis of differential plasma metabolites. On the left is a dendrogram of metabolite clustering, on the right is the name of the metabolite, and the closer two metabolite branches are to each other, the closer they are in terms of expression. (E) Differential abundance score map of KEGG pathway of plasma metabolome. A positive DA score indicates that the expression trend of all annotated differential metabolites in the pathway was upregulated in the HF group; the dots are distributed to the right of the center axis, and the longer the line segment, the more the overall expression of the pathway tends to be upregulated. * means p < 0.05. 3.6. Relationships Between the Rumen Microbiome, Rumen Metabolome and Plasma Metabolome Spearman’s rank correlations between the rumen microbiota and rumen metabolites were assessed. The results showed that Quinella, Fretibacterium, Selenomonas, Cyrenella, Stilbella and Candida were significantly enriched in HF buffaloes and showed a significant positive correlation with rumen metabolites such as lucidenic acid J, secoeremopetasitolide B and erinacine P (p < 0.05). Meanwhile, they were negatively correlated with organic compounds such as 3-amino-2-naphthoic acid, N-acetyl-b-glucosaminylamine, and flazine (p < 0.05) ([112]Figure 7A,B). Similarly, Spearman’s rank correlations between the rumen microbiota and plasma metabolites showed that Quinella, Fretibacterium and Selenomonas were significantly enriched in HF buffaloes and showed significant positive correlation with plasma metabolites such as corchorifatty acid F, PC (18:3 (6Z, 9Z, 12Z)/18:3(9Z, 12Z, 15Z)), hexadecanedioic acid, LysoPE (0:0/18:2 (9Z, 12Z)) and 5-tetradecenoic acid (p < 0.05). The rumen fungi such as Cyrenella, Stilbella and Beauveria showed significant positive correlations with dehydroadoreone, 3′,4′,5′-trimethoxycinnamyl alcohol acetate and prostaglandin E1. Meanwhile, Candida, Phaeosphaeriopsis and Microsphaeropsis showed significant negative correlations with LysoPC (22:4 (7Z, 10Z, 13Z, 16Z)) ([113]Figure 7C,D). Spearman’s rank correlations between the differential rumen microbiota and rumen-specific metabolites were performed. The results showed that Quinella, Fretibacterium and Selenomonas were significantly enriched in the HF group and showed significant positive correlations with metabolites such as Bisnorbiotin, PE (15:0/16:0), LPC (18:1), and isovalerylglutamic acid (p < 0.05). Then, Paeniclostridium, Lachnoclostridium, and Monogobus were significantly and negatively correlated with isovalerylglutamic acid (p < 0.05) ([114]Figures S4 and S5). Figure 7. [115]Figure 7 [116]Figure 7 [117]Open in a new tab Spearman’s rank correlations between the rumen microbiota and metabolites. (A) Differential bacteria and rumen differential metabolites. (B) Differential fungi and rumen differential metabolites. (C) Differential bacteria and plasma differential metabolites. (D) Differential fungi and plasma differential metabolites. * means p < 0.05, ** means p ≤ 0.01, *** means p ≤ 0.001. To identify the rumen microbiota–metabolic interactions, Spearman’s rank correlations between the rumen microbiota, rumen and plasma metabolites were performed (R > 0.7/R < –0.7, p < 0.05). Among these correlations, Quienlla, Stilbella, Cyrenella, Microsphaeropsis, and Paraphaeosphaeria showed correlations with rumen and plasma metabolites, which included lucidenic acid J, N-acetyl-b-glucosaminylamine, flazine, corchorifatty acid F, LysoPE (16:1 (9Z)/0:0), hexadecanedioic acid, isovalerylglutamic acid, PC (18:3 (6Z, 9Z, 12Z)/18:3 (9Z, 12Z, 15Z)), and other metabolites ([118]Figure 8). Figure 8. [119]Figure 8 [120]Open in a new tab Network analysis of rumen microbiota, rumen metabolites and plasma metabolites (R > 0.7/R < −0.7, p < 0.05). The size of the node in the graph indicates the degree size of the node, different colors indicate different species; the color of the connecting line indicates positive and negative correlation, red indicates positive correlation, blue indicates negative correlation; the thickness of the line indicates the size of the correlation coefficient, the thicker the line, the higher the correlation between the species; the shorter the line, the closer the connection between the node. 4. Discussion The dominant rumen microbiota plays an important role in the lactation performance of the host. Current studies have shown that the dominant phyla of cows’ rumen microorganisms were Bacteroidetes, Firmicutes and Proteobacteria [[121]26], which is consistent with our study using 16S rDNA gene amplicon sequencing. Among these, Bacteroidetes mainly decompose non-fiber carbohydrates, while Firmicutes mainly decompose fibers, and there was a strong positive correlation between the milk fat yield and the proportion of Firmicute and Bacteroidetes [[122]27]. The proportions of the Firmicutes and Bacteroidetes in the rumen of high yielding cows were significantly higher compared with low yielding cows [[123]28]. Meanwhile, changes in diet or feeding preferences can