Variety of rumen microbial populations involved in biohydrogenation related to individual milk fat percentage of dairy cows

Our objective was to investigate the contribution of the rumen microbiome on the individual milk fat percentage (MFP) of Holstein dairy cows under the same nutritional and management conditions. From 92 early lactation dairy cows, the top 10 with the highest MFP (HF; n = 10) and the last 10 with the lowest MFP (LF; n = 10) were selected for the study. As a result, the milk trans-10, cis-12 C18:2 content was significant lower in the HF group than that in the LF group (P < 0.001). The rumen acetate to propionate ratio was significant higher in the HF group than that in the LF group (P = 0.035). According to the results of 16S rRNA gene sequencing, a minor but significant difference existed between the groups (P = 0.040). Three genera of the family Lachnospiraceae and four genera of the order Bacteroidales were identified to be the biomarkers for the LF group and HF group in the LEfSe analysis, respectively. Three microbial modules enriched by the family Lachnospiraceae were positively related to the milk trans-10, cis-12 C18:2 content (r_s > 0.60, P < 0.05). According to the results of shotgun metagenome sequencing, three kinds of linoleic acid (LA) isomerase genes were present in the gene pools of the rumen microbiome. Among them, the relative abundance of Bifidobacterium LA isomerase (BBI) was higher in the HF group than that in the LF group (P = 0.007). Three metagenome-assembled genomes (MAGs) with LA isomerase genes were positively correlated to the milk trans-10, cis-12 C18:2 content (r_s> 0.40, P < 0.05). Furthermore, all of these three MAGs were found to be able to produce lactate. Taken together, these results indicate that the increased relative abundance of microbial population with the trans-10 biohydrogenation pathway within the rumen microbiome contributes to the decrease of MFP via the increase of rumen trans-10, cis-12 C18:2 production. This study provides a new perspective for the development of measures for improving the milking performance of dairy cows.

1. Introduction

Milk fat percentage (MFP) is an important indicator of the milk quality and performance of dairy cows. Compared with other nutritional components of the milk, MFP vary to a greater extent (1). Many factors have been found to affect the MFP, such as the breed, diet formula, lactation period, and feeding system (2–4). Even under the same diet and feeding system, the MFP of Holstein dairy cows can vary between 3.0 and 5.0% (5). Recent studies have shown that the individual differences in the composition and metabolic gene profile of the rumen microbiome contribute to variances in MFP (5, 6). However, the underlying roles of rumen microbiome in the modulation of MFP are not as yet understood.

Milk fatty acids (MFA) are the basic units of milk fat. MFAs are derived in two ways: (1) short and medium chain FA (SMCFA, C4-C14) and 50% of C16 FA are de novo synthesized from mainly acetate by mammary epithelial cells, and (2) long chain fatty acids (LCFA, >C16) and another 50% of C16 FA are absorbed from blood (7). In ruminants, the rumen microbiome is responsible for the fermentation of dietary carbohydrate into acetate and also for the lipolysis of dietary fat into free FA. Accordingly, the supply amounts of the above MFA precursors partly depend on the fermentation and lipolysis efficiency of the rumen microbiome. In addition, dietary fats for ruminants mostly consist of unsaturated FA (UFA), especially linoleic acid (LA; cis-9,12 C18:2) and a-linolenic acid (ALA; cis-9,12,15 C18:3). As a detoxifying adaptation, the isomerization from a cis-to a trans-geometric configuration and the reduction of UFA to saturated FA (SFA), i.e., biohydrogenation (BH), is extensively implemented by ruminal microbes (8). Studies have shown that the composition of BH intermediates, such as the contents of trans-10, cis-12 C18:2 and trans-10 C18:1, and the ratio of trans-10 to trans-11 C18:2, is associated with the variances of MFP (7, 9, 10). Accordingly, rumen BH pathways and levels depend on the compositions and metabolic properties of the BH microbial population and are also considered to be the important factors affecting MFP.

In this study, we investigated the differences in the fermentation performance, species composition, and lipids metabolism process of the rumen microbiome, with special attention being paid to the composition and metabolism of the microbial population involved in BH, in dairy cows with extremely high and low MFP. We used 16S rRNA gene sequencing and shotgun metagenome sequencing methods in order to expand our understanding of the contribution of the rumen microbiome to the individual performance of MFP.

2. Materials and methods

2.1. Animals and diets

All experimental procedures involving animals were approved by the Animal Care and Use Committee of Nanjing Agricultural University, in compliance with the Regulations for the Care and Use of Animals (Nanjing Agricultural University, 2018). In total, 92 early lactation Holstein cows of the same parity reared at a commercial dairy farm (Shanghai, China) were used in this study, where they were housed in tie stalls and had free access to water. The total mixed ration (TMR) was provided three times per day at 0,300, 1,200, and 1,700 (forge/concentrate, 40:60, Supplementary Table 1). The chemical analysis of feed samples was based on dry matter (DM) content after drying at 105°C for 2 h (AOAC) (11).

The period of sampling was 3 days. The milk-yield-record and sample-collections of the 92 cows were conducted on d 1–3 at 0,200, 1,100, and 1,600 during experiment. The daily milk samples were mixed in 4:3:3 corresponding to 0,200, 1,100, and 1,600 milking time for the measurement of compositions by a spectrophotometer in Bright Dairy & Food Co., Ltd (Shanghai). The power analysis conducted by G^*Power v.3.1.9.2 revealed that the sample size of five of each group could receive more than 80% power-value. Subsequently, the top 10 cows with the highest MFP (HF group: MFP > 4.11%, Parity = 2, days in milking, DIM = 63 ± 4), and the last 10 cows with the lowest MFP (LF group: MFP <3.38%, Parity = 2, DIM = 59 ± 2) were selected for the following study (Supplementary Figure 1A). The milk samples of HF and LF cows on d 3 were mixed as before and used for the measurement of LCFA profile. The lactating performance of both groups was recorded during the following 5 weeks.

2.2. Rumen contents collection and fermentation index measurement

Rumen contents from 20 cows of two groups were collected at 4 h after a morning feed via oral stomach tubes on d 3 and were then divided into two portions (12). The first portion used to extract microbial DNA for 16S rRNA gene and metagenomic sequencing were transferred into sterile tubes, which were immediately placed into liquid nitrogen. The second portion was first squeezed through 4-layers-gauze for pH measurement and then stored at −20°C for the measurement of the fermentation index and LCFA profile. The methods of determination of NH₃-N and volatile fatty acids (VFA) concentrations were described by Weatherburn et al. (13) and Qin et al. (14), respectively.

2.3. Milk and rumen FA profile analysis

The extractions of LCFA from milk and rumen fluid were based on fatty acid (FA) methylation (15). Non-adecanoic acid methyl ester (M102326, Aladdin, Shanghai, China) was used as an internal standard. The GC system comprised an Agilent 8,890 instrument fitted with an autosampler (Agilent 7693A) and equipped with a CP-Sil 88 capillary column. The standards for quantification were FA methyl ester (FAME) mixtures (18919-1AMP, Supleco, Germany), Methyl trans-11 C18:1 (CDAA-253185M, ANPEL, China) and Methyl trans-10, c12 C18:2 (CDAA-258061M, ANPEL, China). The initial column oven temperature of 150°C for 5 minutes was increased at 2°C/min to 175°C, which was held for 15 min, followed by an increase at a rate of 7°C/min to 200°C, which was held for 20 min, and a final increase at 5°C/min to 220°C, which was held for 25 min. The carrier gas and fuel gas were nitrogen and hydrogen at a flow rate of 1.1 and 40 mL/min, respectively. Meanwhile, the temperature of the injector and detector were 260 and 280°C, and the split ratio was set at 20:1.

2.4. Total DNA extraction

The rumen contents from 20 cows of LF and HF groups were thawed at 4°C, and then, their microbial DNA were extracted using hexadecyltrimethylammonium bromide (16). The concentration and quality of extracted DNA were evaluated by a spectrophotometer and gel electrophoresis, respectively. Subsequently, the extracted DNAs were used in 16S rRNA gene and shotgun metagenomic sequencing.

2.5. 16S rRNA gene sequencing and analysis

The 341F/806R primers were used to amplify the V3-V4 region of the bacterial 16S rRNA genes, and the amplicons were sequenced on an Illumina MiSeq PE250 platform (17). The QIIME2 software suite (version: 2021.2) was applied for the later analysis (18). Briefly, raw reads were qualified, denoised, classified, and counted to generate an amplicon sequence variants (ASVs) table via the DADA2 pipeline (19). The SILVA database (SSU Ref NR 99, release 138) was used to annotated ASVs (20). After the random rarefication of the ASV counts of all samples to 18,984 (the lowest number), the α-diversity and β-diversity of bacterial communities were calculated based on the plugin diversity. Non-parametric permutational multivariate ANOVA (PERMANOVA) was conducted by the scikit-bio package (http://scikit-bio.org/) in Python. The differentially abundant species with the biological relevance to the groups, referred to as the biomarker, were assigned to those ASVs whose Linear discriminant analysis (LDA) scores were more than 2 or <-2 in the LDA Effect Size (LEfSe) test (21).

The co-occurrence network of ASVs was constructed according to Spearman's correlations analysis, where only significantly positive correlations (Spearman's R > 0.6) were visualized via the Fruchterman-Reingold layout with 9,999 permutations by igraph in R (22). The significantly different ASVs within the network were calculated by Deseq2 (23). We used fast-greedy modularity optimization algorithm (24) to identify modules. The correlations of the detected modules and C18 FAs were identified by Spearman's rank correlation analysis.

2.6. Shotgun metagenomic sequencing and analysis

The metagenome libraries were sequenced on an Illumina HiSeq X Ten platform (PE150 mode). The quality of raw reads was assessed by Trimmomatic v.0.33 (25). The qualified reads were compared with the NCBI reference genomes of both host (Bos Tarus, GCA_002263795.2) and plant-originated components in the diets, which were maize (Zea mays, GCA_003185045.1), soybean (Glycine max, GCA_000004515.4), medicago (Medicago truncatula, GCA_000219495.2), oat (Avena sativa, GCA_022788535.1) and cotton (Gossypium arboretum, GCA_025698485.2), to eliminate the data contamination by BWA-MEM v.0.07 (26). The contigs were assembled by MEGAHIT v.1.1.1 (27) from the remaining clean reads, and then used to predicted the open reading frames (ORFs) by Prodigal v.2.6.3 (28). Based on these ORFs, a non-redundant gene set was generated by CD-HIT (29) with the criteria of identity more than 95% and overlap more than 90%. The relative abundance was calculated in transcripts per million (TPM) (30).

The available gene sequences of lipase were downloaded from LED databases (31). So far, LA/ALA isomerase has been divided into Bifidobacterium isomerase (BBI, produce trans-11 isomers), Lactobacillus isomerase (LAI, produce trans-10 or trans-11 isomers) and Propionibacterium isomerase (PAI, produce trans-10 isomers), based on the species and isomeric production (32). In this study, the available gene sequences of PAI (ADE00997), BBI (PCT/CN2019/1218), and LAI (ADD22720.1; CBY89653.1; QTP12276.1) were downloaded from GenBank database (33), respectively. BLASTp (34) was used to detected the above genes from the non-redundant gene set, with the parameters of E-value <1e-5 and minimum sequence identity of 40%. The sequence of isomerase was first aligned with mafft (35) and then used to construct a phylogeny tree with fasttree2 (36). The differences in the relative abundances of the genes were analyzed by the Wilcoxon rank-sum test. Differences were considered as significance at a false discovery rate (FDR) <0.05.

MetaBAT2 v.2.11.1 was used to assemble metagenomic bins (MAGs) (37). CheckM v.1.0.7 (38) was used to evaluated the completeness and contamination of aseembled MAGs. dRep v.3.0.1 (39) was used to remove the replicated MAGs. The taxonomy of dereplicated MAGs whose completeness more than 50% and contamination <10%, referred to as high-quality MAGs, were assigned by GTDB-Tk v.1.5.0 (40). We used metaWRAP v.1.1.0 (41) to calculate the relative abundance of each MAG, and used PhyloPhlAn v.3.0 (42) to analyze the phylogenetic relationships of the high-quality MAGs. Finally, the genes were predicted from the sequences of MAGs. The phylogeny tree of MAGs which was visualized by ggtree package (43) in R.

2.7. Statistical analysis

To make sure the statistical efficiency of this study, the differences in the rumen fermentation index, milk FAs, and rumen FAs concentrations between the groups were analyzed by a two-sided Student t-test, where P < 0.05 was defined as significantly different.

3. Results

3.1. Lactation performance

The lactation performances of 92 cows are listed in Supplementary Table 2. With regard to the lactation performance of the 20 selected dairy cows, no significant differences (P > 0.05), except for those in MFP and milk fat yield (P < 0.001), were found in body weight, DIM, milk yield, milk protein percentage, or lactose percentage between the groups (Table 1). Moreover, the differences in MFP remained significant (P < 0.001) during the following 5 weeks after sampling (Supplementary Figure 1B).

TABLE 1

Table 1. The lactation performances of Holstein cows with highest milk fat percentage (HF) or lowest milk fat percentage (LF).

3.2. Rumen fermentation index and LCFA compositions

Among the rumen fermentation index, the propionate proportion was significant lower (P = 0.032, Student t-test), whereas the ratio of acetate to propionate (Acetate/Propionate) was significant higher (P = 0.035, Student t-test) in the HF cows when compared with those in the LF cows (Table 2). However, no significance (P > 0.05, Student t-test) was found for other aspects of the fermentation index or the FA (C > 8) composition of the rumen between the groups (Table 2).

TABLE 2

Table 2. The rumen fermentation index of Holstein cows with highest milk fat percentage (HF) or lowest milk fat percentage (LF).

3.3. Milk fatty acid compositions

The compositions of milk FA (C > 8) in the two groups were shown in Table 3. Compared with the HF cows, the levels of C8:0, C10:0, C14:0, C15:0, trans-10, cis-12 C18:2, and cis-8,11,14 C20:3 were higher (P < 0.05, Student t-test), and the concentration of cis-9 C16:1 was lower (P = 0.003, Student t-test) in LF cows.

TABLE 3

Table 3. Fatty acid (longer than 8C) compositions in milk and rumen fluid of Holstein cows with highest milk fat percentage (HF) or lowest milk fat percentage (LF) (g/100g of TFA).

3.4. Microbial diversity and taxonomic distribution of two groups

The rarefaction curves of amplicon sequencing were shown in Supplementary Figure 2. The α-diversity of bacterial community showed no significant difference (P > 0.05) between the groups (Supplementary Table 3). However, β-diversity exhibited a significant difference (P = 0.040, PERMANOVA) between two groups (Figure 1A). There were no significant differences (P > 0.05, Wilcoxon rank-sum test) in taxa between the two groups in the phylum or genus levels (Supplementary Figure 3). Nevertheless, LEfSe analysis revealed that Prevotellaceae NK3B31 group, Bacteroidales p-251-o5, Prevotellaceae UCG-003, and Prevotellaceae UCG-001 contributed to the greatest extent (biomarkers in the LFfSe analysis) to the construction of the rumen bacterial community in the HF groups (Figure 1B). On the other hand, Butyrivibrio, Lachnospiraceae AC2044 group, and Lachnoclostridium, all of which belong to the family Lachnospiraceae, contributed the greatest extent (biomarkers in the LFfSe analysis) to the construction of the rumen bacterial community in the LF group.

FIGURE 1

Figure 1. (A) Comparisons of bacterial community composition. (B) LEfSe identification of biomarker genera of the groups. HF, the top 10 cows with the highest milk fat percentage; LF, the last 10 cows with the lowest milk fat percentage; LDA, linear discriminant analysis.

3.5. Co-occurrence patterns of rumen microbes

To detect the ecological interactions of the ruminal bacteria, the co-occurrence network was constructed, with the ASVs having significantly different abundances falling into four modules (Figure 2A). Calculated as the total relative abundance of ASVs in modules, the relative abundances of modules 1, 2, and 4 were significant lower (P < 0.05), and the relative abundance of module 3 was significant higher (P < 0.05) in the HF cows when compared with those of the LF cows (Figure 2B). Furthermore, in the analysis of the compositions of these modules, the relative abundance of Lachnospiraceae was within the top 3 in modules 1, 2, and 4 at the family level (Figure 2C). In the analysis of correlations between the concentration of C18 FA and the relative abundances of the modules, a positive correlation of trans-10, cis-12 C18:2 with module 1 (Spearman's correlation coefficient r_s = 0.63, P = 0.004) and module 4 (r_s = 0.41, P = 0.073), and a negative correlation between trans-10, cis-12 C18:2 and module 3 (r_s = −0.45, P = 0.047) were detected (Figure 2D).

FIGURE 2

Figure 2. (A) The cooccurrence network of the significantly different amplicon sequence variants (ASVs). (B) The relative abundance of the predominant modules who contained the largest number of significantly different ASVs. (C) The compositions of bacterial family of the predominant modules. The ranking of family was based on the relative abundance. The ranking of Lachnospiraceae was marked with the box. (D) The spearman correlations between the relative abundance of predominant modules and the relative content of milk C18 fatty acids (FA). HF, the top 10 cows with the highest milk fat percentage; LF, the last 10 cows with the lowest milk fat percentage. *indicated P < 0.05 in the Wilcoxon rank-sum test.

3.6. Variance in biohydrogenation

Lipases and three types of LA/ALA isomerases have been detected in our metagenome sequencing data. The roles of these enzymes in the rumen BH process, based on previous reports (8), are shown in Figure 3A. Herein, the relative abundance of lipase showed no significant difference (P > 0.05) between the two groups (Figure 3B). However, the relative abundance of BBI was higher (P = 0.007, Wilcoxon rank-sum test), whereas the relative abundance of LAI tended to be lower (P = 0.070, Wilcoxon rank-sum test) in the HF cows when compared with the LF cows (Figure 3C). By analyzing the species with lipase genes, there was no significant difference (P > 0.05, Wilcoxon rank-sum test) of the predominant genera between the two groups (Supplementary Figure 4). In the analysis of the species with these isomerase genes, the relative abundance of Butyrivibrio with LAI tended to be lower (P = 0.052, Wilcoxon rank-sum test), and the relative abundance of Ruminococcus with LAI was lower (P = 0.037, Wilcoxon rank-sum test) in the HF cows, when compared with those in the LF cows (Figure 3D).

FIGURE 3

Figure 3. (A) Major biohydrogenation pathways in the production of C18 fatty acids in the rumen. (B) The relative abundance of lipase and isomerase genes identified in the metagenome data. (C) The classification and relative abundance of three types of isomerase genes identified in the metagenome data. (D) The taxonomic annotation and relative abundance of the isomerase genes. HF, the top 10 cows with the highest milk fat percentage; LF, the last 10 cows with the lowest milk fat percentage; BBI, Bifidobacterium linoleic acid (LA) isomerase; LAI, Lactobacillus LA isomerase; PAI Propionibacterium LA isomerase. *indicated that P < 0.05 in Wilcoxon rank-sum test.

3.7. Correlations between MAGs and trans-10, cis-12 C18:2

As shown in the phylogeny tree, a total of 228 MAGs was obtained in this study (Figure 4A). Among them, a total of 23 MAGs carrying at least one kind of LA/ALA isomerase genes were identified. According to the taxonomy annotation results, 7 MAGs with PAI were assigned to Lachnospiraceae (two), Oscillospiraceae (one), and unclassified bacteria (four), 6 MAGs with LAI were Pseudoscardovia radai (three), Oscillospiraceae (one), Lachnospiraceae (one), and unclassified bacterium (one), and 10 MAGs with BBI were Paludibacteraceae (five), Lachnospiraceae (two), and unclassified bacteria (three).

FIGURE 4

Figure 4. (A) The phylogenetic tree of 228 metagenome-assembled genomes (MAGs). The color orange and blue on the heatmap next to the tree represented significantly positive correlations and negative correlations between the relative abundance of MAGs and the relative content of milk trans-10, cis-12 C18:2, respectively. The MAGs who contained the linoleic acid (LA) isomerase were marked by the black arrow. The MAGs in black frame were selected to do the genomic analysis. (B) The glucose metabolism pathways of the selected MAGs. The detailed information of the enzyme genes was listed in Supplementary Table 4. (C) The hypothesis concerning the relationship of the major suppliers of reducing equivalents and the production of C18:2 isomers for the biohydrogenation bacteria in rumen. HF, the top 10 cows with the highest milk fat percentage; LF, the last 10 cows with the lowest milk fat percentage; BBI, Bifidobacterium linoleic acid (LA) isomerase; LAI, Lactobacillus LA isomerase; PAI Propionibacterium LA isomerase; NAD, nicotinamide adenine dinucleotide; FAD, flavin adenine dinucleotide.

Spearman's correlation analysis was performed between the relative abundance of MAGs with LA/ALA isomerase genes and the relative content of trans-10, cis-12 C18:2 in milk fat. The results showed that MAG with PAI (MAG198; Lachnospiraceae bacterium) and two MAGs with LAI (MAG17 & MAG18, Pseudoscardovia radai) were positively correlated with milk trans-10, cis-12 C18:2 (r_s > 0.4, P < 0.05, Figure 4A). On the other hand, one MAG with BBI (MAG130; Paludibacteraceae bacterium) was negatively correlated with milk trans-10, cis-12 C18:2 (r_s <−0.4, P < 0.05, Figure 4A). A reconstruction of the metabolic pathways of the above four MAGs (Supplementary Table 4) revealed that the three MAGs (MAG17, 18, and 198) positively correlated with milk trans-10, cis-12 C18:2 were able to produce lactate, whereas the MAG130 negatively correlated with milk trans-10, cis-12 C18:2 was not able to produced lactate (Figure 4B).

4. Discussion

The present study aimed to investigate the relationships between rumen microbiome and MFP, where two groups with no significant differences in their body weight, DIM, parity, milk yield, milk protein, and lactose percentage were selected for experiments. Results showed a significant difference in MFP between two groups for a total of 5 weeks, indicating that the disparity was not attributed to a random error of sampling. According to our results, the rumen pH, rumen VFA concentrations, and rumen FAs concentrations exhibited no significant differences between the groups. However, the rumen acetate/propionate, being an important index reflecting microbial community structure, and the β-diversity of ASVs showed the minor but significant differences between the groups (44). Accordingly, we inferred that a minor difference involving a small number of microbes existed between the composition of the rumen microbiome in the two groups. Notably, the concentration of milk trans-10, cis-12 C18:2, being a potent inhibitor of milk fat synthesis (6), significantly increased in the LF cows compared with the HF cows. Because of the absence of Δ12 desaturase in animal tissue, milk trans-10, cis-12 C18:2 should be derived from the BH process of ruminal microbes (45), which suggested that the rumen BH pathways or levels would be different between the groups. In the present study, the inconsistency of rumen trans-10, cis-12 C18:2 might have been induced by the difference in its ruminal escape or in the metabolism efficiency between the groups. Overall, these results implied that the variance in the microbial population undergoing BH within the rumen microbiome contribute to the variance of MFP.

As shown in Figure 3A, two key enzymes, the lipases for lipolysis and the isomerases for LA/ALA isomerization, were involved in BH process (8). There was no difference in the relative abundance of lipases and their taxonomy, which indicated the variance of rumen BH pathways might be in LA/ALA isomerization between two groups. According to the available studies (6, 8), the trans-10 BH pathway that produces trans-10, cis-12 C18:2 intermediates and the tran-11 BH pathway that produces cis-9, trans-11 C18:2 intermediates are two major pathways for LA/ALA isomerization in the rumen (Figure 3A). However, knowledge concerning the species involved in these pathways is limited. A recent study has found that the enrichment of the trans-10 BH pathway is related to an increase in the relative abundance of Lachnospiraceae (6). In our study, Butyrivibrio, Lachnospiraceae AC2044 group and Lachnoclostridium, all of which belong to the family Lachnospiraceae, have been found to be the biomarker genera for the LF group, suggesting a dominant role of the trans-10 BH pathway in the LF group. Our co-occurrence network analysis has also revealed that Lachnospiraceae is the predominant family in the modules 1, 2 and 4 who were significantly enriched in the LF group. Furthermore, modules 1 and 4 are positively correlated with trans-10, cis-12 C18:2. Bauman and Griinari (7) have proposed that the production of trans-10, cis-12 C18:2, rather than cis-9, trans-11 C18:2, from LA isomerization, is one reason for the occurrence of low milk fat syndrome in dairy cows. Taking these data together, we inferred that the dominance of bacteria with trans-10 BH pathway in rumen microbiome, especially Lachnospiraceae, leads to an increased production efficiency of trans-10, cis-12 C18:2, contributing to the decrease of MFP in the LF group.

So far, three types of LA/ALA isomerase with different isomerization production and derived species have been reported: (1) BBI, which was identified from Bifidobacterium and was found to produce trans-11 isomers (32), (2) PAI that was identified from Propionibacterium and was found to produce trans-10 isomers (46), and (3) LAI that was identified from Lactobacillus and was found to produce trans-11 or trans-10 isomers (47). Here, the above three isomerase genes have been detected in our results, where we have found a significantly increased relative abundance of BBI and a decreased trend regarding the relative abundance of LAI in the HF cows compared with the LF cows (Figure 3C). Our results further support our previous hypothesis at the gene level that changes in the BH microbial populations of the rumen leads to changes in the production efficiency of trans-10, cis-12 C18:2 in the LF group.

For a greater understanding of the possible reasons for the shaping of the BH microbial populations, we reconstructed the glucose metabolism pathways of the MAGs that were highly correlated to the concentration of milk trans-10, cis-12 C18:2 (Figure 4B). As a result, in the metabolism pathways of three MAGs (MAG17, 18 and 198) that were positively correlated to trans-10, cis-12 C18:2, the reduction of pyruvate was found to be associated with the oxidation of reduced nicotinamide adenine dinucleotide (NADH). On the other hand, in the metabolism pathways of MAG130 that was negatively correlated to trans-10, cis-12 C18:2, the reduction of pyruvate was associated with the oxidation of reduced flavin adenine dinucleotide (FADH₂). Previous studies showed that reduction of cis-9, trans-11 C18:2 required the oxidation of NADH, and the production of trans-10, cis-12 C18:2 required the oxidation of FAD H₂ (48, 49). Hence, we infer that the usage of NADH or FADH₂ as the major supplier of reducing equivalents affected the choice of trans-11 or trans-10 BH pathway, and thereafter, the production of cis-9, trans-11 C18:2 or trans-10, cis-12 C18:2 for the BH bacteria in the rumen (Figure 4C). However, the hypothesis needs further investigation.

5. Conclusion

The increased relative abundance of microbial population with the trans-10 BH pathway within the rumen microbiome, especially Lachnospiraceae, contributes to the decrease of MFP via the increase of rumen trans-10, cis-12 C18:2 production. Our study provides a new perspective for the development of measures for improving the milking performance of dairy cows.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/sra/PRJNA883555, https://www.ncbi.nlm.nih.gov/sra/PRJNA883576, https://doi.org/10.6084/m9.figshare.21211136.

Ethics statement

The animal study was reviewed and approved by Regulations for the Care and Use of Animals (Nanjing Agricultural University, 2018).

Author contributions

LZ and JZ conducted animal experiments. The statistical analysis and manuscript were performed by LZ and HS. The experiment design and funding were from SM. All authors have read the manuscript.

Funding

This research was funded by the National Key R&D Program of China (2021YFF1000703-01).

Acknowledgments

Our project was supported by the High-Performance Computing Platform of Bioinformatics Center, Nanjing Agricultural University.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fvets.2023.1106834/full#supplementary-material

References

1. Jenkins TC, McGuire MA. Major advances in nutrition: impact on milk composition. J Dairy Sci. (2006) 89:1302–10. doi: 10.3168/jds.S0022-0302(06)72198-1

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Soyeurt H, Dardenne P, Gillon A, Croquet C, Vanderick S, Mayeres P, et al. Variation in fatty acid contents of milk and milk fat within and across breeds. J Dairy Sci. (2006) 89:4858–65. doi: 10.3168/jds.S0022-0302(06)72534-6

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Schingoethe DJ. A 100-year review: total mixed ration feeding of dairy cows. J Dairy Sci. (2017) 100:10143–50. doi: 10.3168/jds.2017-12967

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Tse C, Barkema HW, DeVries TJ, Rushen J, Pajor EA. Impact of automatic milking systems on dairy cattle producers' reports of milking labour management, milk production and milk quality. Animal. (2018) 12:2649–56. doi: 10.1017/S1751731118000654

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Xue M, Sun H, Wu X, Guan LL, Liu J. Assessment of rumen microbiota from a large dairy cattle cohort reveals the pan and core bacteriomes contributing to varied phenotypes. Appl Environ Microbiol. (2018) 84:e00970–18. doi: 10.1128/AEM.00970-18

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Dewanckele L, Toral PG, Vlaeminck B, Fievez V. Invited review: role of rumen biohydrogenation intermediates and rumen microbes in diet-induced milk fat depression: an update. J Dairy Sci. (2020) 103:7655–81. doi: 10.3168/jds.2019-17662

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Bauman DE, Griinari JM. Nutritional regulation of milk fat synthesis. Annu Rev Nutr. (2003) 23:203–27. doi: 10.1146/annurev.nutr.23.011702.073408

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Buccioni A, Decandia M, Minieri S, Molle G, Cabiddu A. Lipid metabolism in the rumen: new insights on lipolysis and biohydrogenation with an emphasis on the role of endogenous plant factors. Anim Feed Sci Technol. (2012) 174:1–25. doi: 10.1016/j.anifeedsci.2012.02.009

CrossRef Full Text | Google Scholar

9. Ventto L, Leskinen H, Kairenius P, Stefański T, Bayat AR, Vilkki J, et al. Diet-Induced milk fat depression is associated with alterations in ruminal biohydrogenation pathways and formation of novel fatty acid intermediates in lactating cows. Br J Nutr. (2017) 117:364–76. doi: 10.1017/S0007114517000010

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Perfield Ii JW, Lock AL, Griinari JM, Sæbø A, Delmonte P, Dwyer DA, et al. Trans-9, Cis-11 conjugated linoleic acid reduces milk fat synthesis in lactating dairy cows. J Dairy Sci. (2007) 90:2211–8. doi: 10.3168/jds.2006-745

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Chemists AA, Horwitz W. Official Methods of Analysis. Arlington, VA: AOAC (1990), p. 489.

Google Scholar

12. Shen JS, Chai Z, Song LJ, Liu JX, Wu YM. Insertion depth of oral stomach tubes may affect the fermentation parameters of ruminal fluid collected in dairy cows. J Dairy Sci. (2012) 95:5978–84. doi: 10.3168/jds.2012-5499

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Weatherburn MW. Phenol-hypochlorite reaction for determination of ammonia. Anal Chem. (1967) 39:971–4. doi: 10.1021/ac60252a045

CrossRef Full Text | Google Scholar

14. Qin WL. Determination of rumen volatile fatty acids by means of gas chromatography. J Nanjing Agric Coll. (1982) 4:110–6.

Google Scholar

15. Shen X, Dannenberger D, Nuernberg K, Nuernberg G, Zhao R. Trans-18:1 and Cla isomers in rumen and duodenal digesta of bulls fed N-3 and N-6 pufa-based diets. Lipids. (2011) 46:831–41. doi: 10.1007/s11745-011-3586-5

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Khafipour E, Li S, Plaizier Jan C, Krause Denis O. Rumen Microbiome composition determined using two nutritional models of subacute ruminal acidosis. Appl Environ Microbiol. (2009) 75:7115–24. doi: 10.1128/AEM.00739-09

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Mu YY, Qi WP, Zhang T, Zhang JY, Mao SY. Gene function adjustment for carbohydrate metabolism and enrichment of rumen microbiota with antibiotic resistance genes during subacute rumen acidosis induced by a high-grain diet in lactating dairy cows. J Dairy Sci. (2021) 104:2087–105. doi: 10.3168/jds.2020-19118

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, et al. Reproducible, interactive, scalable and extensible microbiome data science using qiime 2. Nat Biotechnol. (2019) 37:852–7. doi: 10.1038/s41587-019-0209-9

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. Dada2: high-resolution sample inference from illumina amplicon data. Nat Methods. (2016) 13:581–3. doi: 10.1038/nmeth.3869

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The silva ribosomal rna gene database project: improved data processing and web-based tools. Nucleic Acids Res. (2013) 41:D590–D6. doi: 10.1093/nar/gks1219

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Segata N, Izard J, Waldron L, Gevers D, Miropolsky L, Garrett WS, et al. Metagenomic biomarker discovery and explanation. Genome Biol. (2011) 12:R60. doi: 10.1186/gb-2011-12-6-r60

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Csardi G, Nepusz T. The igraph software package for complex network research. Int J Complex Syst. (2006) 1695:1–9. Available online at: https://igraph.org

Google Scholar

23. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for Rna-Seq data with Deseq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Clauset A, Newman ME, Moore C. Finding community structure in very large networks. Phys Rev E Stat Nonlin Soft Matter Phys. (2004) 70:066111. doi: 10.1103/PhysRevE.70.066111

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for illumina sequence data. Bioinformatics. (2014) 30:2114–20. doi: 10.1093/bioinformatics/btu170

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Li H, Durbin R. Fast and accurate long-read alignment with burrows–wheeler transform. Bioinformatics. (2010) 26:589–95. doi: 10.1093/bioinformatics/btp698

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Li D, Liu C-M, Luo R, Sadakane K, Lam T-W. Megahit: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de bruijn graph. Bioinformatics. (2015) 31:1674–6. doi: 10.1093/bioinformatics/btv033

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Hyatt D, Chen G-L, LoCascio PF, Land ML, Larimer FW, Hauser LJ. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC bioinformatics. (2010) 11:119. doi: 10.1186/1471-2105-11-119

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Fu L, Niu B, Zhu Z, Wu S, Li W. Cd-Hit: accelerated for clustering the next-generation sequencing data. Bioinformatics. (2012) 28:3150–2. doi: 10.1093/bioinformatics/bts565

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Wagner GP, Kin K, Lynch VJ. Measurement of Mrna abundance using rna-seq data: rpkm measure is inconsistent among samples. Theory Biosci. (2012) 131:281–5. doi: 10.1007/s12064-012-0162-3

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Fischer M, Pleiss J. The lipase engineering database: a navigation and analysis tool for protein families. Nucleic Acids Res. (2003) 31:319–21. doi: 10.1093/nar/gkg015

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Mei Y, Chen H, Yang B, Zhao J, Zhang H, Chen W. Research Progress on Conjugated Linoleic Acid Bio-Conversion in Bifidobacterium. Int J Food Microbiol. (2022) 369:109593. doi: 10.1016/j.ijfoodmicro.2022.109593

PubMed Abstract | CrossRef Full Text | Google Scholar

33. O'Leary NA, Wright MW, Brister JR, Ciufo S, Haddad D, McVeigh R, et al. Reference sequence (refseq) database at ncbi: current status, taxonomic expansion, and functional annotation. Nucleic Acids Res. (2016) 44:D733–D45. doi: 10.1093/nar/gkv1189

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. Blast+: architecture and applications. BMC bioinformatics. (2009) 10:421. doi: 10.1186/1471-2105-10-421

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Katoh K, Misawa K, Kuma K, Miyata T. Mafft: a novel method for rapid multiple sequence alignment based on fast fourier transform. Nucleic Acids Res. (2002) 30:3059–66. doi: 10.1093/nar/gkf436

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Price MN, Dehal PS, Arkin AP. Fasttree 2–Approximately Maximum-Likelihood Trees for Large Alignments. PloS One. (2010) 5:e9490. doi: 10.1371/journal.pone.0009490

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Kang DD, Li F, Kirton E, Thomas A, Egan R, An H, et al. Metabat 2: an adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. PeerJ. (2019) 7:e7359. doi: 10.7717/peerj.7359

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW. Checkm: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. (2015) 25:1043–55. doi: 10.1101/gr.186072.114

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Olm MR, Brown CT, Brooks B, Banfield JF. Drep: a tool for fast and accurate genomic comparisons that enables improved genome recovery from metagenomes through de-replication. ISME J. (2017) 11:2864–8. doi: 10.1038/ismej.2017.126

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Chaumeil P-A, Mussig AJ, Hugenholtz P, Parks DH. Gtdb-Tk: a toolkit to classify genomes with the genome taxonomy database. Bioinformatics. (2020) 36:1925–7. doi: 10.1093/bioinformatics/btz848

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Uritskiy GV, DiRuggiero J, Taylor J. Metawrap—a flexible pipeline for genome-resolved metagenomic data analysis. Microbiome. (2018) 6:158. doi: 10.1186/s40168-018-0541-1

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Asnicar F, Thomas AM, Beghini F, Mengoni C, Manara S, Manghi P, et al. Precise phylogenetic analysis of microbial isolates and genomes from metagenomes using phylophlan 3.0. Nat Commun. (2020) 11:2500. doi: 10.1038/s41467-020-16366-7

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Yu G, Smith DK, Zhu H, Guan Y, Lam TTY. Ggtree: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods Ecol Evol. (2017) 8:28–36. doi: 10.1111/2041-210X.12628

CrossRef Full Text | Google Scholar

44. Baker L, Kraft J, Karnezos T, Greenwood S. The effects of dietary yeast and yeast-derived extracts on rumen microbiota and their function. Anim Feed Sci Technol. (2022) 22:115476. doi: 10.1016/j.anifeedsci.2022.115476

CrossRef Full Text | Google Scholar

45. Pereira SL, Leonard AE, Mukerji P. Recent advances in the study of fatty acid desaturases from animals and lower eukaryotes. Prostaglandins Leukot Essent Fatty Acids. (2003) 68:97–106. doi: 10.1016/S0952-3278(02)00259-4

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Liavonchanka A, Hornung E, Feussner I, Rudolph MG. Structure and mechanism of the propionibacterium acnes polyunsaturated fatty acid isomerase. Proc Natl Acad Sci U.S.A. (2006) 103:2576–81. doi: 10.1073/pnas.0510144103

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Kishino S, Takeuchi M, Park S-B, Hirata A, Kitamura N, Kunisawa J, et al. Polyunsaturated fatty acid saturation by gut lactic acid bacteria affecting host lipid composition. Proc Natl Acad Sci U.S.A. (2013) 110:17808–13. doi: 10.1073/pnas.1312937110

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Enjalbert F, Combes S, Zened A, Meynadier A. Rumen microbiota and dietary fat: a mutual shaping. J Appl Microbiol. (2017) 123:782–97. doi: 10.1111/jam.13501

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Hackmann TJ, Firkins JL. Electron transport phosphorylation in rumen butyrivibrios: unprecedented atp yield for glucose fermentation to butyrate. Front Microbiol. (2015) 6:22. doi: 10.3389/fmicb.2015.00622

PubMed Abstract | CrossRef Full Text | Google Scholar