Testing Two Somatic Cell Count Cutoff Values for Bovine Subclinical Mastitis Detection Based on Milk Microbiota and Peripheral Blood Leukocyte Transcriptome Profile

Abstract

Somatic cell count (SCC) is an important indicator of the health state of bovine udders. However, the exact cut-off value used for differentiating the cows with healthy quarters from the cows with subclinical mastitis remains controversial. Here, we collected composite milk (milk from four udder quarters) and peripheral blood samples from individual cows in two different dairy farms and used 16S rRNA gene sequencing combined with RNA-seq to explore the differences in the milk microbial composition and transcriptome of cows with three different SCC levels (LSCC: <100,000 cells/mL, MSCC: 100,000−200,000 cells/mL, HSCC: >200,000 cells/mL). Results showed that the milk microbial profiles and gene expression profiles of samples derived from cows in the MSCC group were indeed relatively easily discriminated from those from cows in the LSCC group. Discriminative analysis also uncovered some differentially abundant microbiota at the genus level, such as Bifidobacterium and Lachnospiraceae_AC2044_group, which were more abundant in milk samples from cows with SCC below 100,000 cells/mL. As for the transcriptome profiling, 79 differentially expressed genes (DEGs) were found to have the same direction of regulation in two sites, and functional analyses also showed that biological processes involved in inflammatory responses were more active in MSCC and HSCC cows. Overall, these results showed a similarity between the milk microbiota and gene expression profiles of MSCC and HSCC cows, which presented further evidence that 100,000 cells/ml is a more optimal cut-off value than 200,000 cells/mL for intramammary infection detection at the cow level.

Keywords: dairy cows; milk microbiota; somatic cell count; subclinical mastitis; transcriptome.

Figure 1

Microbiota composition at the phylum and genus level. (a): Top 5 most dominant phyla identified based on the mean relative abundance of all samples. (b): Top 10 most dominant genera identified based on the mean relative abundance of all samples.

Figure 2

Boxplot showing different α diversity indices in each group. Estimated marginal means of each diversity index was assessed using mixed linear model.

Figure 3

Principal coordinate analysis (PCoA) of Bray−Curtis distances among all milk samples.

Figure 4

Differentially abundant genera and metabolic pathways obtained by pairwise comparison (a) Differentially abundant genera identified in each comparison. (b,c) Boxplot showing log10-transformed relative abundance of two most significantly differentially abundant genera between LSCC and MSCC (d) Venn diagram showing the number of the predicted metabolic pathways identified as significantly differentially abundant in each comparison and log10-transformed relative abundance of pathways that were shared in all comparisons. H−M: HSCC vs. MSCC, L−M: LSCC vs. MSCC, L−H: LSCC vs. HSCC.

Figure 4

Differentially abundant genera and metabolic pathways obtained by pairwise comparison (a) Differentially abundant genera identified in each comparison. (b,c) Boxplot showing log10-transformed relative abundance of two most significantly differentially abundant genera between LSCC and MSCC (d) Venn diagram showing the number of the predicted metabolic pathways identified as significantly differentially abundant in each comparison and log10-transformed relative abundance of pathways that were shared in all comparisons. H−M: HSCC vs. MSCC, L−M: LSCC vs. MSCC, L−H: LSCC vs. HSCC.

Figure 5

Co-occurrence networks of the top 100 predominant bacteria in different groups (spearman correlation > 0.8). Each node represents one genus. The size of the node represents mean relative abundance in the corresponding group. Red edges represent positive correlation and blue edges represent negative correlation.

Figure 6

Differentially expressed genes identification. (a,c): Venn diagram showing the number of differentially expressed genes (DEGs) identified in each comparison in each site. (b,d): Heatmap showing the expression profile of top DEGs (FDR < 0.05). Left panel: site1, Right panel: site2. H−M: HSCC vs. MSCC, L−M: LSCC vs. MSCC, L−H: LSCC vs. HSCC.

Figure 7

Gene functional analyses. (a,b) Bubble diagram and treemap illustrating the top 7 KEGG pathways and GO terms enriched by DEGs with the same regulation direction in two sites. (c) Heatmap showing the result of GSVA using biological process−related GO terms.

Figure 8

Correlation analysis between DEGs and differentially abundant taxa. (a) Associations between differentially abundant genera and the top discriminating DEGs identified in both sites using linear mixed model (‘*’stands for p < 0.05, ‘**’stands for p < 0.01). (b,c) Rlog-transformed read counts of S100A9 and GADD45G in different groups.