Lactation and rumen fermentation
As the degree of mastitis increased, milk yield, milk fat and lactose contents were significantly reduced (P < 0.01), while SCC and milk protein were significantly increased (P < 0.01). There were no significant differences in pH-values (P = 0.214), the concentrations of RUN (P = 0.851), isobutyrate (P = 0.325), isovalerate (P = 0.213) and acetate to propionate ratio (A/P) (P = 0.111) in rumen fluid among the 3 groups. However, the LA, TVFAs (P < 0.01), acetate, propionate, butyrate and valerate (P < 0.001) contents in the CM group were significantly lower than those of the H and SM groups. In addition, the NH3-N concentration in rumen fluid tended to decline as the degree of mastitis increased (P = 0.061) (Table 1).
Rumen microbial diversity
A total of 3,072,603 high-quality sequences were generated from the 60 rumen fluid samples. We obtained 2,973 OTUs ≥ 97% identity, with the Good’s coverage of H, SM and CM groups all above 98% (98.9% ± 0.0013, 98.7% ± 0.0020 and 98.6% ± 0.0011; mean ± standard deviation). Pan OTU was used to observe the increase of the total number of OTU with the increasing number of samples. The pan species curve exhibited that the curve gradually flattened as the number of samples increased, indicating that the sample size of this sequencing was sufficient (Additional file 1: Fig. S1). Similarly, the rarefaction curve showed that the rate of increase in OTU number slowed down with the increasing reads per sample, illustrating that the sequencing depth was adequate (Additional file 1: Fig. S2). Alpha diversity indices analysis suggested that there were significant differences in all the indexes at the OTU level, except for the Simpson index (P = 0.10). The Shannon (P = 0.02), Ace (P = 0.01) and Chao (P = 0.01) indexes in the SM and CM groups were significantly lower than those in the H group (Table 2). The PCoA based on unweighted Unifrac distance (P = 0.003) (Fig. 1a), weighted normalized Unifrac distance (P = 0.022) (Fig. 1b) and OTU classification level showed that ruminal microbiota in the CM group was distinguishable from the H and SM groups. The ruminal microbiota could also be slightly distinguished between H and SM groups but there were still some intersections.
Rumen microorganism composition and structure
In total, 21 phyla and 356 genera of bacteria were obtained (Additional file 1: Table S4 and S5). Bacteria with relative abundance < 1% of all samples were classified as others. The dominating bacteria at phylum level in the 3 groups were Firmicutes (49.69 ± 0.57%), Proteobacteria (39.14 ± 0.42%), Bacteroidetes (4.65 ± 0.28%), Actinobacteria (3.22 ± 0.06%) and Cyanobacteria (1.29 ± 0.04%) (Additional file 1: Fig. S3a). At the genus level, Prevotella_1 (23.89 ± 0.43%), Succiniclasticum (10.74 ± 0.11%), Mollicutes_RF39 (8.77 ± 0.10%), Pseudobutyrivibrio (4.67 ± 0.27%), Succinivibrionaceae_UCG-001 (3.82 ± 0.16%) and Ruminococcaceae_NK4A214_group (3.08 ± 0.25%) were the top 6 in all 3 groups (Additional file 1: Fig. S3b).
Hierarchical cluster analysis conducted for taxa at genus level (top 55 species in total relative abundance) further demonstrated the distribution of ruminal bacteria identified in different udder health status. Through average-linkage clustering techniques, the rumen microorganisms in the 60 dairy cows were aggregated into 3 clusters. Consistent with the PCoA results, ruminal bacteria in the CM group gathered together and separated from the H and SM groups. Although the H and the SM groups were generally differentiated, they still had some intersections, such as SM13, SM4, SM20, SM18 and SM9. There were 25 major bacteria in H group, including Ruminococcus_2, Pseudoscardovia, Acetitomaculum, Bacteroidales, Bifidobacterium, etc. Twelve bacteria, including Succiniclasticum, Ruminiclostridium_9, vadinBE97, Prevotellaceae_UCG-001 etc., were detected in SM group. The CM group contained 17 genera of bacteria, such as Herbaspirillum, Moraxella, Neisseriaceae, etc. (Additional file 1: Fig. S4).
Differentially abundant ruminal microbiota
Kruskal-Wallis H test was used to perform the analysis of significant differences of ruminal bacteria among the 3 groups. The significantly diverse microbes at genus level were assembled in 6 phylum, incorporating in Firmicutes, Bacteroidetes, Proteobacteria, and Cyanobacteria, etc. (Table 3 and Additional file 1: Fig. S5). The ruminal bacteria with high abundance in the CM group were mainly Pseudobutyrivibrio, Gastranaerophilales, Moraxella, Neisseriaceae, etc. In SM group, Ruminiclostridium_9, Bacteroidales_UCG-001 and Enterorhabdus, etc., were most abundant. Comparing with the other two groups, Lachnospiraceae, Prevotella_1, Mollicutes_RF39, Bifidobacterium, etc., differed significantly in H group.
Prediction of rumen microbial function
To further speculate the contribution of these key flora, we attempted to predict the function of rumen microbiota in H (average NSTI = 0.33 ± 0.005), SM (average NSTI = 0.33 ± 0.007) and CM (average NSTI = 0.35 ± 0.005) groups by PICRUSt. At superclass level, 42.97% of the genes participated in metabolism, 21.79% belonged to genetic information processing, and 11.55% were involved with environment information processing. No significant difference was detected in the functions of rumen microorganisms among the 3 groups at level 1 (Fig. 2a). While at class level, a total of 39 KEGG pathways were noted, and 26 of them were different among H, SM and CM samples (FDR-adjusted P < 0.05). The abundance of flora related to amino acids coenzymes (FDR-adjusted P = 0.008) and vitamins (FDR-adjusted P = 0.004), nucleotides (FDR-adjusted P = 0.001) and lipid metabolism (FDR-adjusted P = 0.014), etc., were prominent in H group. Conversely, flora correlated with energy metabolism (FDR-adjusted P = 0.001), infection diseases (FDR-adjusted P = 0.001), immunity system diseases (FDR-adjusted P = 0.001), etc., accounted for a greater proportion in the CM group. In the SM group, the significantly different functional flora was less than H and CM groups, mainly focusing on translation (FDR-adjusted P = 0.003) and genetic information processing (FDR-adjusted P = 0.002) (Fig. 2b and Additional file 1: Table S6). At subclass level, among the 268 KEGG characteristics, 33 were different in abundance among the 3 groups (FDR-adjusted P < 0.05). Flora related to methane metabolism (FDR-adjusted P = 0.003), bacterial invasion and infection (FDR-adjusted P = 0.04), primary immunodeficiency (FDR-adjusted P = 0.037) and lipopolysaccharide biosynthesis (FDR-adjusted P = 0.043), etc., were significantly enriched in the CM group (Fig. 2c and Additional file 1: Table S7).
Correlation between ruminal flora, and rumen fermentation and lactation performance
According to Spearman correlation coefficients, the rumen microbiota was associated with several rumen fermentation parameters, and milk yield and compositions (Fig. 3). Lachnospiraceae_UCG-006 (r = 0.409, FDR-adjusted P = 0.039) was positively correlated with milk protein. Neisseriaceae (r = 0.434, FDR-adjusted P = 0.008), Moraxella (r = 0.402, FDR-adjusted P = 0.005) and Gastranaerophilales (r = 0.402, FDR-adjusted P = 0.004) were positively correlated with milk SCC. Bifidobacterium (r = 0.349, FDR-adjusted P = 0.036) was positively associated with lactose, while Pasteurellaceae (r = -0.354, FDR-adjusted P = 0.039), Moraxella (r = -0.360, FDR-adjusted P = 0.045) and Gastranaerophilales (r = -0.360, FDR-adjusted P = 0.043) were negatively associated with lactose. Prevotella_1 (r = 0.384, FDR-adjusted P = 0.047) was positively correlated with acetate, but Moraxella (r = -0.396, FDR-adjusted P = 0.045) and Gastranaerophilales (r = -0.396, FDR-adjusted P = 0.044) were negatively correlated to acetate. Herbaspirillum (r = -0.370, FDR-adjusted P = 0.043) and Gastranaerophilales (r = -0.346, FDR-adjusted P = 0.043) were negatively associated with butyrate. In addition, Pseudobutyrivibrio (r = -0.400, FDR-adjusted P = 0.02), Herbaspirillum (r = -0.347, FDR-adjusted P = 0.04) and Lachnobacterium (r = -0.333, FDR-adjusted P = 0.042) were negatively associated with LA (Additional file 1: Table S8).
Comparative analysis of ruminal metabolites
The total ion chromatogram (TIC) of QC samples in positive and negative ion mode showed an overlap in the response intensity and retention time of each chromatographic peak (Additional file 1: Fig. S6), which indicated the high reliability of data quality. To exhibit the overall difference and the degree of the variability of rumen samples among and within the groups, unsupervised PCA was performed under positive and negative ionization modes, respectively (Fig. 4). Symbols representing ruminal metabolites in H, SM and CM groups could be basically separated. Further supervised OPLS-DA analysis of the data was performed in positive (Additional file 1: Fig. S7a, c and e) and negative ion modes (Additional file 1: Fig. S8a, c and e). Significant differences were detected between CM and H groups, SM and H groups, and CM and SM groups. Model parameter R2Y (cum) were all above 0.7 and Q2 (cum) were all above 0.5, indicating a good ability of model prediction. Model evaluation parameters of response permutation testing, Q2, were all below 0, suggesting that there was no overfitting of the model. The results illustrated that there were significant differences in the rumen metabolism under different udder health conditions of dairy cows (Additional file 1: Fig. S7b, d and f and Fig. S8b, d and f).
Differential metabolites identification and evaluation
A total of 642 differential metabolites were selected including 384 and 258 in positive and negative ion mode, respectively, from the 60 rumen fluid samples. Under positive and negative ionization modes, the significantly differential metabolites (VIP > 1 and P-value < 0.05) among the 3 groups were shown in Additional file 1: Table S9-S14. Significantly differential metabolites (top 15) in the rumen between CM and H groups, SM and H groups, and CM and SM groups are shown in Fig. 5 (a), (b) and (c), respectively. Comparing to H group, 10beta-Hydroxy-6beta-isobutyrylfuranoeremophilane (FDR-adjusted P < 0.001) and 12-oxo-20-dihydroxy-leukotriene B4 (FDR-adjusted P < 0.001) in CM group increased 12.23 and 16.36-fold, respectively. Meanwhile, these two in CM group increased 2.88 and 3.64-fold, respectively, compared with SM group (Table 4). Methenamine (FDR-adjusted P < 0.001), 5-Hydroxymethyl-2-furancarboxaldehyde (5-HMF) (FDR-adjusted P < 0.001) and 6-Methoxymellein (FDR-adjusted P < 0.001) in CM and SM groups were both higher than those in H group. Besides, xestoaminol C (FDR-adjusted P = 0.016), lentialexin (FDR-adjusted P = 0.031) and cinnamic acid (FDR-adjusted P = 0.042) in SM group were increased compared to H group. Less abundance of 2-Phenylbutyric acid (2-PBA) (FDR-adjusted P < 0.001) was detected in CM group compared with H group. In CM group, N-Acetylcadaverine (FDR-adjusted P = 0.020) and (3R, 5Z)-5-Octene-1, 3-diol (FDR-adjusted P = 0.011) were reduced compared with SM group (Table 4 and Fig. 5).
To evaluate the diagnostic value for BM of these key differential metabolites in the rumen, receiver operator characteristic (ROC) analysis was performed, which reflected the relationship between false positive rate (x-axis) and true positive rate (y-axis) (Fig. 6). In the current study, the AUC values of the 5 key differential metabolites, including 10beta-Hydroxy-6beta-isobutyrylfuranoeremophilane, 12-oxo-20-dihydroxy-leukotriene B4, methenamine, 5-HMF and 6-Methoxymellein, were all above 0.9, and the ROC curves were all close to the upper left corner. This indicated that these 5 differential metabolites in the rumen screened in the present study were of high value to be served as the biomarkers for bovine mastitis.
Differential metabolites classification and metabolic pathways enrichment analysis.
Compound classification of differential metabolites in the rumen fluid samples were processed based on comparisons with the HMDB. Differentially clustered metabolites belonging to lipid-like molecules, organ heterocyclic compounds, and organic acids and derivatives compounds accounted for the top 3 among the 3 groups (Table 4 and Additional file 1: Fig. S9).
Furthermore, the enrichment of metabolic pathway analysis showed that, 10beta-Hydroxy-6beta-isobutyrylfuranoeremophilane was enriched in limonene and pinene degradation pathway (FDR-adjusted P = 0.008). 12-oxo-20-dihydroxy-leukotriene B4 was enriched in arachidonic acid metabolism pathway (FDR-adjusted P = 0.026). Besides, methenamine was the intermediate product of formaldehyde biosynthesis (FDR-adjusted P = 0.036). 5-Hydroxymethyl-2-furancarboxaldehyde participated in the furfural degradation (FDR-adjusted P = 0.032). 6-Methoxymellein was enriched in the biosynthesis of plant secondary metabolites pathway (FDR-adjusted P = 0.045). And 2-Phenylbutyric acid was enriched in the butanoate metabolism pathway (FDR-adjusted P = 0.042) (Table 5).
Correlations between ruminal metabolites, and lactation and rumen fermentation parameters.
Based on Spearman correlation coefficients and Euclidean distance matrix, we observed significant correlations between significantly differential metabolites, and lactation as well as rumen fermentation parameters (Fig. 7). 10beta-Hydroxy-6beta-isobutyrylfuranoeremophilane (r = -0.584, FDR-adjusted P < 0.001) and 12-oxo-20-dihydroxy-leukotriene B4 (r = 0.601, FDR-adjusted P < 0.001) were negatively correlated to acetate, whereas the 2-PBA (r = 0.520, FDR-adjusted P < 0.001) was positively correlated to acetate. Meanwhile, 10beta-Hydroxy-6beta-isobutyrylfuranoeremophilane (r = -0.472, FDR-adjusted P = 0.046) and 12-oxo-20-dihydroxy-leukotriene B4 (r = -0.504, FDR-adjusted P = 0.002) were also negatively associated with butyrate. Besides, (3R, 5Z)-5-Octene-1, 3-diol (r = 0.612, FDR-adjusted P < 0.001) and N-Acetylcadaverine (r = 0.670, FDR-adjusted P < 0.001) were positively associated with LA. 10beta-Hydroxy-6beta-isobutyrylfuranoeremophilane (r = -0.700, FDR-adjusted P < 0.001), 12-oxo-20-dihydroxy-leukotriene B4 (r = -0.688, FDR-adjusted P < 0.001), 5-HMF (r = -0.537, FDR-adjusted P < 0.001) and 6-Methoxymellein (r = -0.542, FDR-adjusted P < 0.001) were negatively associated with lactose, while 2-PBA was positively associated with lactose (r = 0.547, FDR-adjusted P < 0.001). Furthermore, 5-HMF was negatively associated with milk fat (r = -0.604, FDR-adjusted P < 0.001). (3R, 5Z)-5-Octene-1, 3-diol (r = -0.606, FDR adjusted P < 0.001) and N-Acetylcadaverine (r = -0.582, FDR-adjusted P < 0.001) were negatively correlated to milk protein. 10beta-Hydroxy-6beta-isobutyrylfuranoeremophilane (r = 0.716, FDR-adjusted P < 0.001), 12-oxo-20-dihydroxy-leukotriene B4 (r = 0.712, FDR-adjusted P < 0.001), methenamine (r = 0.603, FDR-adjusted P < 0.001), 5-HMF (r = 0.550, FDR-adjusted P < 0.001) and 6-Methoxymellein (r = 0.600, FDR adjusted P < 0.001) were positively associated with milk SCC, but 2-PBA was negatively associated with milk SCC (r = -0.573, FDR-adjusted P < 0.001) (Additional file 1: Table S15).
Associations between BM-linked rumen microbes and metabolites
The Spearman correlation coefficients were calculated to reflect the correlations between significantly differential bacteria and metabolites via heat map (Fig. 8). In CM group, Lachnospiraceae_UCG-006 (r = 0.493, FDR-adjusted P = 0.038), Pasteurellaceae (r = 0.565, FDR-adjusted P = 0.043), Moraxella (r = 0.599, FDR-adjusted P = 0.033), and Gastranaerophilales (r = 0.589, FDR-adjusted P = 0.031) were positively related to 10beta-Hydroxy-6beta-isobutyrylfuranoeremophilane. Lachnospiraceae_UCG-006 (r =0.542, FDR-adjusted P = 0.042), Moraxella (r = 0.574, FDR-adjusted P = 0.038) and Gastranaerophilales (r = 0.561, FDR-adjusted P = 0.037) were also positively related to 12-oxo-20-dihydroxy-leukotriene B4. Pasteurellaceae (r = -0.556, FDR-adjusted P = 0.036), Pseudobutyrivibrio (r = -0.588, FDR-adjusted P = 0.01), Moraxella (r = -0.595, FDR-adjusted P = 0.025) and Gastranaerophilales (r = -0.585, FDR-adjusted P = 0.025) were negatively correlated with N-Acetylcadaverine. In SM group, Enterorhabdus (r = 0.527, FDR-adjusted P = 0.047) were positively associated with 5-HMF. The preponderant strains in H group, including Prevotella_1 (r = 0.572, FDR-adjusted P = 0.040), Lachnospiraceae (r = 0.578, FDR-adjusted P = 0.042), and Izimaplasmatales (r = 0.545, FDR-adjusted P = 0.042) were positively correlated with 2-PBA, while Izimaplasmatales was negatively associated with methenamine (r = -0.577, FDR-adjusted P = 0.041) and 6-Methoxymellein (r = -0.499, FDR-adjusted P = 0.002) (Additional file 1:Table S16).