The experimental procedures were approved by the Animal Care Committee, Beijing University of Agriculture (Beijing, China). Raw milk samples (n = 266) were obtained directly from commercial dairy farms in Beijing in the summer, when S. agalactiae numbers are elevated. The selected animals with similar parity were between 115 and 230 days in milk (DIM) and had an average SCC of 250,000 cells/ml, ranging from 110,000 to 440,000 cells/ml, measured 3 days prior to transportation to the barn. No samples were collected from cows with GBS or high SCCs, in accordance with previously published methods (Jørgensen et al., 2016). Bacteriological analyses of the 266 milk samples were performed, which revealed those cows with GBS. The experimental animals were distributed to the control [CON; healthy animals, SCC (85.13 ± 39.83 × 10³ cells/ml)] and GBS [SCC (329.03 ± 84.98 × 10³ cells/ml)] groups. There were no initial differences among the six Chinese Holstein dairy cows in each group in terms of milk yield (24.42 ± 1.68 kg/d; p = 0.16), DIM (154.60 ± 7.58 days; p = 0.71), parity (2.63 ± 0.38; p = 0.36), or body weight (BW; 570.25 ± 24.87 kg; p = 0.90). The animals were housed individually in stalls bedded with sawdust and were fed the same basal diet. The ingredients and nutrient composition of the total mixed ration (TMR) are presented in Table 1. A total of 15 ml of composite milk from each animal, with approximately equal volumes from each lactating udder quarter, was transferred to a sterile plastic bottle (Corning, Inc., Corning, NY, United States) and kept on ice until transport to the laboratory and storage at −80°C for further analysis.

For the isolation and detection of S. agalactiae, all raw milk samples were first processed as described by Jørgensen et al. (2016) and Gao et al. (2017). Briefly, the milk samples were shaken, and 100 μl aliquots were removed under sterile conditions, streaked onto 5% rabbit blood agar and incubated at 37°C for 24 h. Colony morphology was observed, and several colonies were picked for Gram staining and microscopic examination. Biochemical profiling was conducted by testing for catalase activity and lactose, esculin, mannitol, salicin, raffinose, sorbitol, monosaccharide, xylose, and cellobiose fermentation. In addition, isolates were confirmed as S. agalactiae using the slide agglutination test for Lancefield GBS (Oxoid). The milk bacteriological analysis results are shown in Supplementary Table 1.

DNA was extracted from the milk bacteria using an InstaGene Matrix DNA Extraction Kit (Bio-Rad Laboratories, Hercules, CA, United States) according to the manufacturer’s instructions. S. agalactiae were then identified by PCR (Bio-Rad Laboratories) amplification of the 16S rRNA locus of the DNA. The primers were synthesized by Sangon Biotech, Co., Ltd. (Shanghai, China). A negative control (a sample without genomic DNA) and a positive control (a sample containing DNA of S. agalactiae, CVCC 3940, China Institute of Veterinary Drug Control, Beijing, China) were included in all PCR assays. The following primers were used: 27F: 5′-AGAGTTTGATCCTGGCTCAG-3′ and 1492R: 5′-GGTTACCTTGTTACGACTT-3′ (Saffary et al., 2002). The amplification programs were performed as described by Rubens et al. (1991). The sequence data of the isolated S. agalactiae strains were analyzed using the Basic Local Alignment Search Tool (BLAST) program available from the National Center for Biotechnology Information (NCBI). During universal 16S rRNA fragment sequence analysis, the cut-off level for pairwise comparisons was 99%.

Microbial DNA was extracted from the milk samples using a Power Soil DNA Isolation Kit (Qiagen, Crawley, United Kingdom) according to the manufacturer’s protocols. The yield and purity of the extracted DNA were assessed with a NanoDrop 1000 spectrophotometer (Nanodrop Technologies, Wilmington, DE, United States), and DNA was stored at −80°C until use. The V1–V2 region of the bacterial 16S rRNA gene was amplified with a common primer pair (forward, 5′-CGTATCGCCT-CCCTCGCGCCATCAG-3′; reverse, 5′-CTA-TGCGCCTTGCCAGCCCGCTCAG-3′) combined with adapter sequences and barcode sequences, as previously described (Oikonomou et al., 2014). PCR amplification was performed in a total volume of 50 μl containing 10 μl of buffer, 0.2 μl of Q5 high-fidelity DNA polymerase, 10 μl of high-GC enhancer, 1 μl of dNTPs, 10 μM of each primer and 60 ng of genomic DNA. The thermocycling conditions were as follows: initial denaturation at 95°C for 5 min, followed by 15 cycles of 95°C for 1 min, 50°C for 1 min and 72°C for 1 min and a final extension at 72°C for 7 min. The products from the first round of PCR were purified through VAHTS^TM DNA Clean Beads. A second round of PCR was then performed in a 40 μl reaction containing 20 μl of 2 × Phusion HF MM, 8 μl of ddH₂O, 10 μM of each primer and 10 μl of PCR products from the first step. The thermocycling conditions were as follows: initial denaturation at 98°C for 30 s, followed by 10 cycles of 98°C for 10 s, 65°C for 30 s, and 72°C for 30 s and a final extension at 72°C for 5 min. Samples with a bright main band of approximately 450 bp were chosen and mixed in equal concentration ratios as previously described (Li et al., 2018). The mixture of PCR products was purified using a GeneJET Gel Extraction Kit (Thermo Fisher Scientific, Waltham, MA, United States). The sequencing libraries were validated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States) and quantified with a Qubit 2.0 fluorometer (Thermo Fisher). Finally, high-throughput sequencing analysis of bacterial rRNA genes was performed on the purified pooled samples using the Illumina HiSeq 2500 platform (2 × 250 paired ends) at Biomarker Technologies Corporation, Beijing, China.

Bioinformatics analysis was performed using both FLASH version 1.2.11 and Quantitative Insights Into Microbial Ecology (QIIME) version 1.9.1. These programs were considered to give similar results, as previously reported (Li et al., 2018). Briefly, the sequences were classified into taxa with BLAST searches in the Ribosomal Database Project (RDP) database at a 97% similarity threshold. Operational taxonomic units (OTUs) with counts greater than three in at least one of the samples were retained for further analysis. The selected OTUs were normalized to the relative abundance for each sample. The diversity of the microbial communities was estimated by using the R package phyloseq on the online BioCloud Platform¹. Within-sample diversity (alpha diversity) was assessed through the Shannon diversity index calculated in a randomly selected subset of the OTU database. Between-sample microbial diversity (beta diversity) was assessed through phylogenetic-based weighted UniFrac distances (Catherine and Rob, 2005). For a deeper analysis of the diversity of the major evolutionary clades in the ruminal microbiota, the OTUs were filtered to acquire those with a relative abundance of at least 1% in at least one sample.

Milk samples were analyzed for specific components using an automated ultrasonic milk composition analyzer (Bentley Instruments, Chaska, MN, United States) and a method similar to our previously published procedure (Wang et al., 2017). The untargeted metabolic profiling analysis was performed by means of GC-TOFMS, and sample preparation was performed as previously described (Yan et al., 2012). The data were analyzed using an XploreMET v3.0 system (Metabo-Profile, Shanghai, China) (Chen W. L. et al., 2016). MetaboAnalyst was used to analyze range/scale data and to perform principal component analysis (PCA; De Barros et al., 2010) and KEGG pathway analysis of the altered metabolites².

The data on milk yield, milk composition and alpha diversity indices were analyzed using the PROC MIXED program in SAS 9.4 (SAS Institute, Inc., Cary, NC, United States), as shown in the following model: Yij = μ + Ti + ej, where Yij is the dependent variable, μ is the overall mean, Ti is the effect of treatment (no disease or GBS; considered fixed), and ej is the residual. p < 0.05 was considered statistically significant, and a trend was indicated by p < 0.10. The alpha diversity indices are presented as the mean ± SD. PCA and orthogonal partial least-squares-discriminant analysis (OPLS-DA) were performed to visualize the metabolic alterations between the experimental groups after mean centering and unit variance scaling. Variables with variable importance in the projection (VIP) values exceeding 1.0 were considered relevant for group discrimination. In this study, the OPLS-DA model was validated with sevenfold permutation tests. Significant differences in metabolites between groups were assessed using Wilcoxon rank-sum tests. A correlation matrix between metabolites and ruminal bacterial species was generated using Spearman correlation coefficients and visualized by using the R language.

All raw sequences were submitted to the NCBI Sequence Read Archive (SRA³) under accession number SRP192494.

The milk yield, percentage of milk fat and protein, and energy-corrected milk (ECM) content showed no difference between the two groups (Table 2). However, the yield of lactose was lower (p = 0.02) in the cows with GBS than in the CON cows. The percentages of lactose (p = 0.05) and fat (p = 0.08) tended to be lower in the GBS group than in CON cows. The milk from the GBS group also showed a significantly higher SCC than milk from the CON group (p = 0.01). Although subclinical mastitis did not significantly affect the milk yield, the percentages of lactose and fat were lower in the GBS group, indicating that the quality of the milk may have been reduced. Whether the metabolite alterations in milk from GBS cows are associated with diminished milk quality needs to be investigated further.

A total of 910,984 merged sequences passed quality control, and 802,166 high-quality sequences for 12 samples were available for downstream analysis. The samples had an average length of 422 bp, and there was > 99% depth coverage, which indicated that the data volume was reasonable and could reflect changes in most bacterial flora. The alpha diversity indices of the bacterial communities in the two groups are shown in Table 3. The ACE and Chao 1 estimates of richness were not different between the two groups (p > 0.05). The Shannon index revealed that the diversity of microbial communities in the GBS group was significantly lower than that in the CON group (p = 0.05). The Shannon index and the Simpson index both show bacterial diversity. However, these values can vary depending on whether the method of calculation is based on abundance or biomass. Therefore, no significant difference was found between groups for the other alpha diversity indices. Principal coordinate analysis (PCoA) was performed to reveal distinct microbiota between the two groups, and the weighted UniFrac metrics were also determined (Figure 1). PCoA showed that principal coordinates 1 and 2 accounted for 73.33 and 12.35% of the total variance, respectively. This indicated that the microbiota in milk from cows with GBS and CON cows were markedly distinct.

Fifty-four bacterial phyla were identified among all the samples. Firmicutes, Proteobacteria, and Actinobacteria were the three predominant phyla, representing 51.45, 32.48, and 5.34% of the total sequences, respectively (Figure 2). The remaining phyla represented less than 1% of all the sequences. We found that the GBS group had a significantly lower relative abundance of Proteobacteria, Actinobacteria, and Acidobacteria and higher relative abundance of Firmicutes at the phylum level than did the CON group (p < 0.01; Figure 3). The 10 predominant genera were identified: Streptococcus (24.22%), Acinetobacter (21.37%), Romboutsia (4.99%), Turicibacter (2.64%), Stenotrophomonas (2.33%), Enterococcus (1.86%), Microbacterium (1.66%), Aerococcus (1.60%), Corynebacterium-1 (1.49%), and Bacteroides (1.44%). At the genus level, taxa with a relative abundance of 1% in at least one sample were further analyzed, and the relevant genera are presented in Figure 4. Compared to the CON group, the GBS group had a significantly increased relative abundance of Streptococcus and tended to show increased abundances of Turicibacter (p = 0.07) and Enterococcus sp. (p = 0.07). In contrast, the relative abundances of Acinetobacter, Stenotrophomonas, Microbacterium, and Corynebacterium-1 were significantly decreased in the GBS group (p < 0.05; Figure 4).

In total, 254 measurable and reproducible metabolite signals were obtained from the milk samples. After strict quality control, 155 metabolites were identified using the JiaLib^TM metabolomics library. They consisted mainly of glycerophosphoethanolamines, fatty acids, amino acids, peptides, carbohydrates, and pyrimidines. To characterize the differences in metabolic profiles between the CON and GBS groups, PCA and OPLS-DA were conducted. The PCA score plots (Figure 5A) revealed that the milk samples from the CON group were readily separable from those of the GBS group. Principal components 1 and 2 accounted for 23.41 and 18.69% of the variation, respectively. The parameters for assessment of OPLS-DA model quality could be represented by validation plots (Figure 5B). The corresponding Q2Y value of the OPLS-DA model in milk was 0.715, and the R2Y value was 0.995. The intercept of Q2Y with a threshold < zero indicated a valid model (Q2Y = −0.012). The OPLS-DA results showed that the two groups had distinct, significantly different metabolite compositions (Figure 5C). Thus, PCA and OPLS-DA were both good indicators of the differences in metabolites in milk from control cows compared to cows with GBS.

Through a combination of the statistical analysis results and the VIP values obtained from the OPLS-DA, 22 different metabolites (p < 0.05 and VIP > 1) were found to be significantly different between the control group and the mastitis group (Table 4). These different metabolites were diverse, including organic acids, carbohydrates, nucleotides, amino acids, and alcohols. Milk from cows with GBS contained higher levels of pyrrole metabolites, phenylpyruvic acid, uridine, and glycerol and had a higher ratio of homogentisic acid to 4-hydroxyphenylpyruvic acid and a higher ratio of xanthine to guanine (14.29-fold higher) than did CON cows.

To better evaluate how multiple pathways differed between the CON group and the GBS group, a Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis of the pathways related to the different metabolites was conducted (Table 5). The altered metabolites included organic acids, carbohydrates, nucleotides, amino acids and alcohols. Pathway topology analysis (Figure 6) with a criterion of p < 0.05 showed that eight main metabolic pathways were enriched between the groups, including galactose metabolism; pentose and glucuronate interconversion; starch and sucrose metabolism; alanine, aspartate, and glutamate metabolism; arginine biosynthesis; the citrate cycle (TCA, tricarboxylic acid cycle); D-glutamine and D-glutamate metabolism; and neomycin, kanamycin, and gentamicin biosynthesis pathways. The top four altered metabolic pathways between the two groups of cows were galactose metabolism, pentose and glucuronate interconversion, starch and sucrose metabolism, and the alanine, aspartate and glutamate metabolism pathways.

Spearman correlation coefficients were calculated, and clear correlations were observed between the perturbed milk microbiome and the altered metabolite profiles (r > 0.5 or < −0.5, p < 0.05). Figure 7 shows that the significantly different types of metabolites between the groups were highly correlated with specific bacteria, demonstrating functional correlations between the milk microbiome and associated metabolites. We also found that Streptococcaceae, Lachnospiraceae, Lactobacillaceae, and Corynebacteriaceae were highly correlated with most metabolites. Only Streptococcaceae were significantly correlated with D-mannose levels. The most significantly changed metabolite ratio, the xanthine: guanine ratio, was positively correlated with Streptococcaceae and negatively correlated with Lachnospiraceae, Lactobacillaceae, Xanthomonadaceae, Microbacteriaceae, and Brevibacteriaceae.

The microbial diversity in milk was dramatically decreased in the GBS group. Next-generation DNA sequencing has been widely used to evaluate microbial communities for in-depth characterization of the microbial profiles of specific ecosystems (Baltasar et al., 2014). Previous studies have indicated that there is a close relationship between reduced microbial diversity and mastitis in both human milk (Fernández et al., 2013) and bovine milk (Bhatt et al., 2012; Falentin et al., 2016). We found that the GBS group had a significantly lower Shannon index than the CON group (Table 3). These findings are consistent with those of previous studies indicating that there were more bacterial genera in non-infected quarters than in infected quarters (showing that diversity was related to health status; Braem et al., 2012) and that there was lower diversity in the uterine microbiota of cows suffering from metritis than in that of CON cows (Santos and Bicalho, 2012). However, several factors should be considered in accounting for the discrepancies between the groups, including the bacterial lysis and genomic DNA extraction methods, the different quarters sampled, etc. In a previous study, Henderson et al. (2013) clearly demonstrated that sampling techniques and DNA extraction methods can significantly affect the microbial compositions of cow and sheep rumens.

It has been shown that spoilage by Acinetobacter (Gurung et al., 2013), a serious problem in the dairy industry, as well as spoilage by Streptococcus (Zanardi et al., 2014), is strongly correlated with high SCCs (greatest size effect) as well as risk of GBS. Similar to previous studies that have described milk microbial communities, our study revealed the presence of Acinetobacter and Streptococcus in the milk samples (Figure 4). Streptococcus thermophilus generally produces bacteriocins that protect the final product from microbial spoilage (Kabuki et al., 2009). In our study, the bacterial population was extremely high in milk from cows with GBS, indicating that the metabolic end-products and enzymes of the bacteria had a considerable influence on milk quality.

One interesting observation from our study was the presence of Firmicutes and Proteobacteria in our milk samples (Figure 3). These species are common in the rumen microbiome (Tong et al., 2018; Wang et al., 2018), which may indicate a transfer between the rumen and the mammary gland, similar to the gut-breast axis observed in humans that allows transfer from mothers to neonates (Ted et al., 2014). Anne et al. (2010) reported that immune cells might be involved in the transport of microorganisms to different sites in the body, which may explain the high SCCs in milk from cows with GBS. Alternatively, the alterations in bacterial diversity may have resulted from fecal contamination through ascending colonization of the teat canal. In our study, we also found that Lactobacillaceae were highly negatively correlated with most metabolites associated with the GBS group (Figure 7). It has been reported that Lactobacillus lactic acid bacteria can degrade the disaccharide lactose, the major carbohydrate in milk (Huys et al., 2006). The proteolytic and lipolytic capacities of Lactobacillus in carbohydrate metabolism have considerable effects on the texture, taste, and aroma of fermented milk (Kim, 2014). These findings suggest that changes in the milk microbiota may be related to changes in milk yield or metabolic processes involved in milk synthesis.

Metabolomics, which integrates both technological and analytical approaches, can enable a comprehensive understanding of an organism’s physiological and biochemical status (Zhao et al., 2014). In this study, GC-TOFMS was used to compare metabolite levels and to determine differences in milk metabolites between healthy cows and those with GBS to discover potential biomarkers for milk quality responses to GBS. We identified 155 compounds in milk, more than previously identified in studies using nuclear magnetic resonance (NMR) (Hailemariam et al., 2014), liquid chromatography-mass spectrometry (LC-MS/MS) (Boudonck et al., 2009) and gas chromatography-mass spectrometry (GC–MS) approaches. The higher number of metabolites identified in this study shows that more sensitive detection and identification methods enable better qualitative identification of metabolites. Among the 22 significantly different metabolite parameters, five were more prevalent in the GBS group than in the control group, including phenylpyruvic acid, uridine, glycerol, the homogentisic acid: 4-hydroxyphenylpyruvic acid ratio, and the xanthine: guanine ratio (Table 4). It is known that aromatic compounds are derived from aromatic amino acids (Sébastien et al., 2012) such as tyrosine, which could account for the increases in phenylpyruvic acid (by 2.63-fold) and the homogentisic acid: 4-hydroxyphenylpyruvic acid ratio (by 2.92-fold) in the GBS group. Furthermore, phenylalanine is an essential amino acid and the precursor of catecholamines, which are neurotransmitters and adrenaline-like substances (Waisbren et al., 2007). Phenylalanine is transformed into tyrosine through the action of phenylalanine hydroxylase and a biopterin cofactor (Gilbert, 1992). Tyrosine is an important amino acid in many proteins, peptides, and even enkephalins and is also the precursor for hormones such as thyroxin and catecholoestrogens (Lemmon and Schlessinger, 2000). These characteristics may partially explain the variations in milk yield and quality between the groups.

It has previously been demonstrated that glycerol oxidation produces glyceric acid, which is mediated by the transfer of glucose carbon to serine in cows (Black et al., 1955). When glyceric acid is secreted excessively in milk, the lactating animal may suffer from D-glyceric aciduria and D-glycerate anemia, resulting in metabolic acidosis, progressive neurological impairment, seizures, hypotonia, and other adverse effects (JöRn Oliver et al., 2010). The current study revealed that glycerol was dramatically increased in the GBS group compared with the CON group, which could explain why GBS not only reduces milk yield and quality but also may lead to poor health. Thus, our data indicate that glycerol could be a potential biomarker for the diagnosis of GBS.

The current study not only identified differences in the metabolite profiles between the groups but also pinpointed the likely pathways that generated these metabolite differences. Based on an integrated analysis of the key metabolic pathways producing the 22 identified milk metabolites, carbohydrate metabolism appears to be the most important pathway associated with mastitis-induced metabolic changes in cows (Figure 6 and Table 5). It is well-known that carbohydrate metabolism performed by microorganisms in ruminant animals normally involves fatty acyl-CoA synthetases, which produce fatty acids, lactate and other products that can be absorbed in the intestine and transported to the mammary gland (Friggens et al., 1998). It has also been reported that alterations in the starch and sucrose pathway cause disorders of carbohydrate metabolism accompanied by significant increases in lactate and propionate and decreases in acetate (Ametaj et al., 2010; Pan et al., 2016). Moreover, since propionate is the predominant glucose precursor in ruminants (Ungerfeld, 2015), elevations in glucose stimulate milk protein synthesis (Rulquin et al., 2004). Interestingly, our correlation results revealed that the presence of Streptococcaceae was negatively correlated with the levels of D-mannose and glucose 1-phosphate (Figure 7). Thus, it is reasonable to conclude that S. agalactiae causes mammary gland metabolic dysfunction and milk composition changes. These data suggest a potential mechanism by which GBS affects milk composition.

In our pathway analysis, the TCA cycle was significantly stimulated in the GBS group compared to the CON group. The TCA cycle is the primary metabolic pathway in mitochondria and is critically important for the health of dairy cows. ATP is produced from a surplus of acetyl-CoA in the mitochondria (Wellen and Thompson, 2012) that originates from glucose, fatty acid, and lactate metabolism (Hardie and Carling, 2010). Previous studies have reported that microorganisms convert carbohydrates to pyruvate and acetyl-CoA through the glycolytic pathway and the pentose phosphate pathway (Chen L. et al., 2016). Due to the different needs of cells, ATP allocation is flexible. Pyrophosphate is a high-energy phosphate precursor that acts as a primary energy source for mitochondria in cells. Following creatinine formation, the pyrophosphate pathway plays an important role in monitoring energy status in tissues (Wyss and Kaddurahdaouk, 2000; Antunes-Fernandes et al., 2016). In our study, pyrophosphate levels were significantly lower in the milk of cows with GBS than in that of CON cows, which could have resulted in lower energy allocation to the mitochondria for milk synthesis in the infected cows (Table 4). Pyrophosphate and phosphoenolpyruvic acid were positively correlated with Corynebacteriaceae and with Bacteroidaceae, Clostridiaceae, Ruminococcaceae, and Rikenellaceae, respectively (Figure 7). Taken together, our results support the hypothesis that GBS alters the metabolite balance in cow milk by disrupting the TCA cycle in the mammary gland. This effect appears to be closely associated with milk yield and composition via changes in mitochondrial production of ATP as a fuel for metabolic activity.

In the present study, alanine, aspartate, and glutamate metabolism and arginine biosynthesis were significantly decreased in GBS cows compared to CON cows (Table 5). Previous reports have indicated that arginine, glycine, and methionine are the predominant amino acids for creatine synthesis in the liver (Poortmans et al., 2005). Creatine, an intermediate metabolite in energy reactions, plays an important role in regulating negative energy balance in dairy cows (Wang et al., 2016) and could also be a potential diagnostic biomarker for heat stress (Wang et al., 2016). Aspartate is a precursor for the biosynthesis of many essential compounds and participates in some cellular signaling pathways, such as gluconeogenesis, in dairy cows (Piccioli-Cappelli et al., 2014). Therefore, our results provide the potential molecular mechanisms underlying the observed differences in the cow milk metabolomes between the GBS and CON groups.