Associations among Milk Microbiota, Milk Fatty Acids, Milk Glycans, and Inflammation from Lactating Holstein Cows

ABSTRACT Milk oligosaccharides (MOs) can be prebiotic and antiadhesive, while fatty acids (MFAs) can be antimicrobial. Both have been associated with milk microbes or mammary gland inflammation in humans. Relationships between these milk components and milk microbes or inflammation have not been determined for cows and could help elucidate a novel approach for the dairy industry to promote desired milk microbial composition for improvement of milk quality and reduction of milk waste. We aimed to determine relationships among milk microbiota, MFAs, MOs, lactose, and somatic cell counts (SCC) from Holstein cows, using our previously published data. Raw milk samples were collected at three time points, ranging from early to late lactation. Data were analyzed using linear mixed-effects modeling and repeated-measures correlation. Unsaturated MFA and short-chain MFA had mostly negative relationships with potentially pathogenic genera, including Corynebacterium, Pseudomonas, and an unknown Enterobacteriaceae genus but numerous positive relationships with symbionts Bifidobacterium and Bacteroides. Conversely, many MOs were positively correlated with potentially pathogenic genera (e.g., Corynebacterium, Enterococcus, and Pseudomonas), and numerous MOs were negatively correlated with the symbiont Bifidobacterium. The neutral, nonfucosylated MO composed of eight hexoses had a positive relationship with SCC, while lactose had a negative relationship with SCC. One interpretation of these trends might be that in milk, MFAs disrupt primarily pathogenic bacterial cells, causing a relative increase in abundance of beneficial microbial taxa, while MOs respond to and act on pathogenic taxa primarily through antiadhesive methods. Further research is needed to confirm the potential mechanisms driving these correlations. IMPORTANCE Bovine milk can harbor microbes that cause mastitis, milk spoilage, and foodborne illness. Fatty acids found in milk can be antimicrobial and milk oligosaccharides can have antiadhesive, prebiotic, and immune-modulatory effects. Relationships among milk microbes, fatty acids, oligosaccharides, and inflammation have been reported for humans. To our knowledge, associations among the milk microbial composition, fatty acids, oligosaccharides, and lactose have not been reported for healthy lactating cows. Identifying these potential relationships in bovine milk will inform future efforts to characterize direct and indirect interactions of the milk components with the milk microbiota. Since many milk components are associated with herd management practices, determining if these milk components impact milk microbes may provide valuable information for dairy cow management and breeding practices aimed at minimizing harmful and spoilage-causing microbes in raw milk.


RESULTS
Milk fatty acids and microbiota. No relationships between MFA b-diversity and microbial b-diversity, or between MFA a-diversity and microbial a-diversity were identified. Yet, there were a total of 925 relationships between MFAs and microbial genera ( Fig. 1; see also Supplemental File 1). Many of the genera had relationships with 15 or more MFAs, and most had associations with a variety of MFAs: short-chain fatty acids (SCFAs), odd-chain saturated fatty acids (SFAs), even-chain fatty acids, branched-chain saturated fatty acids, monounsaturated fatty acids (MUFAs), and polyunsaturated fatty acids (PUFAs). The relationships tended to be mostly positive or mostly negative for a given genus. Corynebacterium, Prevotella, Chryseobacterium, Pseudomonas, an unknown Enterobacteriaceae genus, and Acinetobacter had mostly negative associations with MFAs. Meanwhile, Bifidobacterium, Staphylococcus, an unknown Muribaculaceae genus, Bacteroides, an unknown Oscillospiraceae genus, Ruminococcus, Mogibacterium, and anunknown Lachnospiraceae genus had mostly positive associations with MFAs. The trends in associations between genera and individual MFAs were mirrored in the associations between genera and total MFA concentration (Fig. 2). Total MFA concentration had the strongest negative associations with an unknown Enterobacteriaceae genus and Acinetobacter and the strongest positive associations with an unknown Muribaculaceae genus, an unknown Oscillospiraceae genus ("uncultured 111"), and Lactobacillus.
Milk oligosaccharides and microbiota. No relationships between MO b-diversity and microbial b-diversity or between MO a-diversity and microbial a-diversity were identified. Still, there were a total of 224 relationships between distinct MOs and microbial genera ( Fig. 3; Supplemental File 2). Most of the genera had associations with a variety of MOs: neutral and fucosylated, neutral and nonfucosylated, and sialylated and nonfucosylated. Corynebacterium, Enterococcus, Aerococcus, Stenotrophomonas, Pseudomonas, and an unknown Enterobacteriaceae genus had many and mostly positive relationships with MOs. Meanwhile, Bifidobacterium, Sphingomonas, a Rhizobiaceae genus, Chryseobacterium, a FIG 1 Associations between distinct milk fatty acids (MFAs) and distinct microbial genera identified with generalized linear mixed-effects modeling. The bar plot depicts the total number of negative or positive relationships between microbial genus counts and square root-transformed MFA concentration (grams per 100 mL of milk) and number by MFA type: monounsaturated fatty acid (MUFA), polyunsaturated fatty acid (PUFA), saturated branched-chain fatty acid (SFA: BCFA), even-chain saturated fatty acid above six carbons (SFA: even chain), odd-chain saturated fatty acid above six carbons (SFA: odd chain), and saturated short-chain fatty acid (SFA: SCFA). Genera are grouped by class-level taxonomy. Genus RF39 belongs to the class Bacilli, genera uncultured.111 and UCG.005 belong to the family Oscillospiraceae, genus UCG.010 belongs to the order Oscillospirales, and genus Family_XIII_AD3011_group belongs to the family Anaerovoracaceae.
Lachnospiraceae genus, and Acinetobacter had many and mostly negative associations with MOs.
Lactose and microbiota. Lactose abundance was not associated with predicted relative abundance of microbial b-galactosidase in the milk microbial community or FIG 2 Associations between total MFAs and distinct microbial genera. Generalized linear mixed-effects modeling was used to identify relationships between square root-transformed MFA (grams per 100 mL of milk) and microbial genus counts. The estimated beta coefficient for total MFA abundance as a fixed effect for microbial genus abundance is depicted. Genera are grouped by class-level taxonomy. Genus RF39 belongs to the class Bacilli, genera uncultured.111 and UCG.005 belong to the family Oscillospiraceae, genus UCG.010 belongs to the order Oscillospirales, and genus Family_XIII_AD3011_group belongs to the family Anaerovoracaceae.

FIG 3
Associations between distinct MOs and distinct microbial genera identified with generalized linear mixed-effects modeling. The bar plot depicts the total number of negative or positive relationships between microbial genus counts and MO relative abundance and numbers by MO type: neutral and fucosylated, neutral and nonfucosylated, and sialylated and nonfucosylated. Genera are grouped by class-level taxonomy. Genus RF39 belongs to the class Bacilli, genus uncultured.96 belongs to the family Peptococcaceae, and genus UCG.005 belongs to the family Oscillospiraceae.
combined predicted relative abundances of microbial b-galactosidase and 6-phospho-b-galactosidase. However, lactose did have relationships with 14 microbial genera (Fig. 4), half of which were negative and half positive. Lactose had pronounced negative estimated beta coefficients as a fixed effect for Corynebacterium, Stenotrophomonas, and Pseudomonas abundances. Lactose also had a large positive beta coefficient as a fixed effect for Romboutsia and Acinetobacter and especially for Paracoccus abundance.
Milk microbiota and SCC. The median SCC in milk samples was 55,000/mL, with 67.2% of milk samples below 100,000/mL and 21.3% of milk samples above 200,000/ mL (Fig. S1 in Supplemental File 3). Microbial a-diversity was not a predictor of SCC. Also, none of the microbial genera were predictors of SCC, after multiple-test correction. Before multiple-test correction, Staphylococcus (Fig. S2) and Aerococcus (Fig. S3) abundances had significant positive relationships with log SCC (P value , 0.05).
SCC and milk oligosaccharides and fatty acids. Log SCC had a significant negative relationship with log lactose as determined by linear mixed-effects modeling (LMM) (beta coefficient = 20.096; P value , 0.001) (Fig. 5) and repeated-measures correlation (RMCORR) (r = 20.228; P value = 0.010). Log SCC did not have a significant relationship with MO a-diversity. However, a positive relationship between log SCC and the neutral, nonfucosylated MO comprised of eight hexoses (8_0_0_0_0) was identified with LMM (beta coefficient = 0.137; false-discovery rate [FDR]-adjusted P value = 0.011) (Fig. 6). No relationship was identified between log SCC and MFA a-diversity or between log SCC and individual MFA abundance.

DISCUSSION
Associations among bovine milk microbiota, milk components, and SCC were investigated to inform our hypotheses that milk components and mammary inflammation work both independently and in concert to limit pathogenic bacteria and promote beneficial microbes in milk. While relationships between diversity measures of MFAs or MOs were not found to correlate with microbial aor b-diversity measures, lactose and one MO were found to correlate with SCC. Additionally, this study identified many associations among individual MFAs, MOs, and microbes that might be indicative of anti-infective functions for MFAs and MOs in milk. FIG 4 Associations between lactose and distinct microbial genera. Generalized linear mixed-effects modeling was used to identify the relationship between lactose (grams per 100 mL of milk) and microbial genus counts. The estimated beta coefficient for lactose abundance as a fixed effect for microbial genus abundance is depicted. Genera are grouped by class-level taxonomy.

Relating Milk Microbes, Components, and Inflammation
Microbiology Spectrum Similar to our study, Moossavi et al. (33) did not find an association of MFA profile, MO profile, or MO a-diversity with milk microbial composition or a-diversity in human milk, yet in the same cohort, they found associations between distinct MFAs or MOs and milk microbial diversity and between distinct MFAs or MOs and individual microbial taxa (14). These findings suggest that certain MFAs or MOs have unique relationships with certain milk microbes such that a more diverse MFA or MO composition does not necessarily support a more diverse milk microbial composition.
As reviewed by Desbois and Smith, antimicrobial activity of fatty acids has been observed for a variety of microbes, including methanogens, viruses, and fungi, with fatty acid structure and shape being important to efficacy (34). The mechanisms underlying the antibacterial activity of fatty acids are varied and require further elucidation, but many appear to relate to destabilization of the cell membrane and may explain why polyunsaturated fatty acids (PUFAs) and monounsaturated fatty acids (MUFAs) tend to be more toxic to microbes than saturated fatty acids (SFAs). Fatty acids tend to have more microbial toxicity with greater numbers of carbon double bonds and cis configuration. However, establishing an acidic environment may be another mode of antibacterial activity as well. The human skin exudes free fatty acids, which are critical in limiting Staphylococcus aureus colonization (34). FAs have been shown to have antibacterial activity against potential mastitis-causing species such as S. aureus and Corynebacterium bovis (34)(35)(36). We found that several bacterial genera to which mastitis-associated and psychrotrophic species belong, including Corynebacterium and Pseudomonas (37,38), and an unknown genus of Enterobacteriaceae (39) (a family that contains several mastitis-associated bacteria [40]) had mostly negative associations with numerous different MFAs, including most MUFAs, PUFAs, and SCFAs. The same trends were observed for Acinetobacter, which is frequently associated with healthy udders (41) but is psychrotrophic (37) and clinically relevant for its antibiotic-resistant infections in humans (42). The majority of MFA 23:0, .99%, is incorporated into FIG 5 Association between somatic cell count per mL (SCC) and lactose (grams per 100 mL of milk). Using linear mixed-effects modeling, an estimated beta coefficient of 20.096 was determined for log SCC as a fixed effect for lactose abundance. All three samples for each cow are plotted, and lines depicting the estimated intercept for each cow and beta coefficient are shown. sphingomyelin (43) and was negatively associated with potentially pathogenic genera: Pseudomonas, Enterococcus, Stenotrophomonas, Corynebacterium, and an unknown Enterobacteriaceae genus (Supplemental File 1). This is particularly interesting because sphingomyelins can be binding sites for pathogenic bacteria and their toxins (44).
Surprisingly, though, Staphylococcus had almost exclusively positive relationships with MFAs, even with MUFAs and PUFAs. This is similar to the case with Staphylococcus in human milk, which was reported to have positive associations with several human MFAs, including MFAs 22:0 and 20:2 n-6 (12), like in this study. Meanwhile, we found that the potentially beneficial or symbiotic bacteria, Bacteroides, Bifidobacterium, Ruminococcus, and an unknown Lachnospiraceae genus, had mostly positive associations with individual SFAs, MUFAs, and PUFAs (45,46). Prevotella, however, had mostly negative associations. SCFAs (represented by butyric acid and caproic acid in this data set) especially showed a trend of positive relationships with potentially beneficial or symbiotic microbes (including Lactobacillus [46]) and negative relationships with potentially harmful bacteria (including Enterococcus [47] and Stenotrophomonas [48]) and with Aerococcus, which, in turn, had a trending positive relationship with log SCC. However, none of the SCFAs or other MFAs had relationships with log SCC, so we do not have evidence to suggest that there may have been a link between MFAs and inflammation. Interestingly, many of the potentially symbiotic microbial genera (e.g., Bifidobacterium, Bacteroides, Lactobacillus, and certain Clostridia genera) that were positively associated with MFAs also tended to correlate positively with each other and negatively with potentially harmful bacteria, Enterococcus, unknown Enterobacteriaceae, Pseudomonas, and Acinetobacter (Fig. S4). These trends may suggest that MFAs are important in promoting a potentially less harmful milk microbiota, but how they may do so in vivo has yet to be determined. According to previous research, one could speculate that MFAs may affect milk microbes in vivo directly through antimicrobial activity and indirectly through modulation of immune responses via binding to free fatty acid receptors which are expressed in a large range of tissues and immune cells (49).
Interestingly, some of the genera that had mostly positive or mostly negative associations with MFAs had mostly opposite associations with MOs. These included Corynebacterium, Pseudomonas, and an unknown Enterobacteriaceae genus which had mostly positive FIG 6 Association between SCC and a neutral, nonfucosylated eight-hexose milk oligosaccharide (MO) 8_0_0_0_0. Using linear mixed-effects modeling, an estimated beta coefficient of 0.137 was determined for log SCC as a fixed effect for log-transformed relative abundance of MO 8_0_0_0_0. All three samples for each cow are plotted, and lines depicting the estimated intercept for each cow and the beta coefficient are shown. associations with MOs, while Bifidobacterium and an unknown Lachnospiraceae genus had mostly negative associations with MOs. These differences may relate to the differing modes of action of MOs versus MFAs on potentially pathogenic milk microbes. For example, MOs might act through antiadhesive effects, which would then retain bacterial cell viability and detection using culture-independent methods. This in contrast to MFAs, which might act primarily through antimicrobial effects such as cell membrane disruption, which would destroy the bacterial cell and expose the DNA to environmental degradation, preventing detection via sequencing.
We might further speculate that within the milk matrix, and perhaps even within the mammary tissue, the beneficial effects of MOs are exhibited through antiadhesive properties against potential pathogens, rather than prebiotic activities to promote the growth of probiotics, as is so often observed in the gut. Notably, neutral, fucosylated MOs and sialylated, nonfucosylated MOs usually had positive associations with potentially harmful or psychrotrophic bacteria, including Staphylococcus, Corynebacterium, Enterococcus, and the unknown Enterobacteriaceae genus. Relatedly, Williams et al. (12) found a positive association between 29-fucosyllactose and Staphylococcus. The neutral, nonfucosylated MO with eight hexoses (8_0_0_0_0) was positively associated with log SCC and negatively associated with Staphylococcus (Supplemental File 2), which trended toward positive association with log SCC. It is interesting to consider that this MO may be connected to inflammation and involved in limiting Staphylococcus to improve mammary health.
Similar to the case with MFAs and MOs, we found several relationships between lactose abundance and milk microbial genera. The most pronounced positive relationship was between lactose abundance and Paracoccus abundance. Although Paracoccus is not frequently detected in bovine milk, this genus can contain species with b-galactosidase that degrade lactose in milk (50). On the other hand, the most pronounced negative relationship of lactose was with Corynebacterium. Interestingly, lactose had negative relationships with a few other potentially harmful bacteria-Enterococcus, Stenotrophomonas, and Pseudomonas-and a negative relationship with Aerococcus, which showed a trend toward positive association with log SCC. Lactose also had a negative association with log SCC, which has previously been observed in human milk (12). Follow-up experiments are needed to determine whether lactose abundance impacts milk microbes that play a role in mammary inflammation or whether lactose is simply produced in higher quantities when mammary inflammation is lower. In support of a negative correlation with inflammation not mediated by prebiotic effects, lactose abundance was not associated with predicted relative abundances of microbial b-galactosidases and 6-phospho-b-galactosidase. While this may indicate that lactose, at least in the concentration range analyzed in this study, did not have a prebiotic effect on milk microbiota, we cannot exclude the possibility that microbes altered expression of lactose-degrading enzymes in response to the lactose.
In this study, we found many associations between milk microbes and MFAs, MOs, and lactose, but we found little evidence of associations between the milk microbes, or milk components, and mammary inflammation measured by SCC. However, if there were relationships between milk components and SCC or milk microbes and SCC in the milk, they may have been difficult to detect in our samples since we did not collect quarter-level milk samples and there were up to 4 days' difference in collection days for SCC measurements versus for measurements of milk components and microbiota. Additional analysis of quarter-level relationships between SCC and milk components is needed.
The findings reported here, along with previous related findings in humans (12-14, 16, 33), suggest that future research should determine if MFAs, MOs, or lactose directly impacts microbes in raw milk or within the udder. This work identifies candidate individual milk components and milk microbes to interrogate specifically in relation to each other. It is unknown if there are viable, metabolically active microbes within healthy udders. However, there does appear to be potential for microbial entry through the teat canal during and following milking and for microbial viability in milk stored in the teat cistern between milkings (51,52). Identifying antimicrobial effects of milk components upon viable milk microbes in vivo could lead to feeding and breeding objectives designed to modify milk composition in order to select for beneficial bacteria over potentially pathogenic or spoilage-causing bacteria.
Limitations. There were several aspects of this study that may have limited the ability to detect relationships among abundances of milk components, milk microbiota, and SCC. First, the milk samples from which lactose, MOs, MFAs, and microbiota were measured were collected on different days than the milk samples from which SCC were measured. Similarly, MOs were measured in milk samples resulting from combining milk from the same cow on two consecutive days within the same period, while milk samples from consecutive days were kept separate for lactose, MFA, and microbiota measurements. A potential implication of these differing collection dates is a loss of sensitivity and accuracy in detecting relationships between milk components/microbiota and SCC and in detecting relationships between MOs and microbiota. Additionally, correlations between absolute abundances of milk microbiota and milk components, or SCC, are not necessarily detectable when using compositional measurement of the milk microbiota. Additional research measuring the absolute abundance of the microbes could improve our understanding of the correlations between milk microbiota and milk components. While numerous relationships were identified in this study between milk microbes and milk components, these relationships may not necessarily apply to viable milk microbes, as a culture-independent approach was taken to detect milk microbes. Accurate detection of viable bacteria at the time of collection was not possible since milk samples were frozen during transport and storage at 280°C, which leads to significant reductions in viable cells and alterations in viable microbial community structure (53). However, there is potential insight gained from measuring microbes that were present in the teat regardless of viability. For example, potentially harmful microbes could have lost viability within the teat prior to sample collection due to recent interactions with the host immune response, milk components, or other bacteria. On the other hand, microbial viability might not be critical to inducing changes in production of some milk components that are involved in pathogen defense. Research on intramammary inoculation of heat-killed mastitis pathogens (54,55) suggests that dead bacteria can elicit immune responses in the mammary tissue that are effective in defending against viable pathogens. Therefore, measuring the abundance of bacteria, regardless of viability, may be critical in understanding potential relationships between potential pathogens and milk components.

MATERIALS AND METHODS
Bovine housing and feed. As described previously (56) and under protocols approved by the University of Wisconsin-Madison Institutional Animal Care and Use Committee (protocol number A005945), 76 lactating Holstein cows were housed in tie stalls at the U.S. Department of Agriculture-Agricultural Research Service U.S. Dairy Forage Research Center Dairy Farm in Prairie du Sac, WI, as part of a 20-week crossover feed intervention trial to determine differences in feed efficiency of cows on a low-starch-high-fiber diet versus on a high-starch-low-fiber diet. Cows were evenly assigned to diets to achieve groups with similar parity, dry matter intake, milk production, and body weight. As a result, there was not a difference in days in milk (DIM) or in parity between cows that began the low-starchhigh-fiber diet in period 1 (median DIM at start of period 1 = 126; median parity = 1) versus the cows that began the high-starch-low-fiber diet in period 1 (median DIM at start of period 1 = 124; median parity = 2). For the entire cohort, parities ranged from 1 to 6, with a median of 2 (see Fig. S1A in the supplemental material). The trial started in October 2017 and ended in April 2018 and was comprised of three diet periods. During period 0, all cows were fed a diet for 31 days composed of 50% low-starchhigh-fiber diet and 50% high-starch-low-fiber diet, and milk samples were collected starting at 3 weeks into the diet. During period 1, cows were fed for 70 days with either the low-starch-high-fiber diet or the high-starch-low-fiber diet, and milk samples were collected starting at 5 weeks into the diet. During period 2, cows were fed the opposite diet for 70 days and milk samples were collected starting at 5 weeks into the diet.
Milk collections. Preceding milk collection, teats were stripped with 3 streams of milk, disinfected with chlorine dioxide-containing Gladiator predip (BouMatic, WI), and towel dried, and milking equipment was disinfected with iodine in water. After milk collection, teats were disinfected with iodine-containing Udderdine postdip (BouMatic). Milk was collected from cows three times daily: at 4 a.m., 10:30 a.m., and 6 p.m. Milk samples for lactose, MO, MFA, microbiota, and SCC measurements were collected at 4 a.m. Milk samples were collected from each of the three diet periods. For lactose, MO, MFA, and microbiota measurements, milk samples were collected on two consecutive days in each period (Fig. 7). For MO measurements, milk samples from consecutive collection dates were combined before processing. For SCC measurements, milk samples were collected before and after the dates of milk collection for lactose, MO, MFA, and microbiota measurements and differed in collection date by 1 to 6 days. On the first day of milk collection, 21 November 2017, cows ranged in DIM from 82 to 193, with a median of 115 (Fig. S1B), and on 19 March 2018 (the last day of milk collection analyzed in this study), cows ranged in DIM from 200 to 311. Milk samples were immediately stored at 210°C, shipped on dry ice, and stored at 270°C while awaiting processing. SCC was measured from fresh milk samples.
Milk analyses. For lactose measurements, milk samples of around 50 mL were warmed to 38°C in a water bath, inverted five times, and measured on a Foss Milkoscan Minor type 78110 to obtain lactose concentration per 100 g of milk. The MOs are from the milk samples that Durham et al. reported previously (57). In short, MOs were measured from 200 mL of milk extracted through skimming, ethanol precipitation of proteins, and C 18 and graphitic carbon solid-phase extraction, prior to isobaric labeling and analysis by nano-liquid chromatography quadrupole time-of-flight tandem mass spectrometry. MOs were identified based on their monosaccharide composition as the number of hexose_N-acetylhexosamine_fucose_N-acetylneuraminic acid_N-glycolylneuraminic acid, followed by an isomer designation, as needed. MOs were reported as abundance relative to the matching oligosaccharide internal standard. The MFAs are from the milk samples that Picklo et al. reported earlier (58). Briefly, MFAs were measured from 0.1 mL of milk in duplicate for each collection day and reported as the mean value for each treatment in grams per 100 mL of milk. MFAs were converted to fatty acid methyl esters using acetyl chloride and were separated and quantified using gas chromatography with flame ionization detection (31). Microbial community analysis was performed on the milk samples that Coates et al. published previously (6). Succinctly, 10 mL of milk were centrifuged and then processed with the ZymoBIOMICS DNA miniprep kit to extract DNA. Particularly low DNA concentrations (,1.9 ng/mL) were obtained from several milk samples, so another round of DNA extraction ("duplicate") was performed on another 10-mL aliquot from the sample milk sample. DNA extracts were sent to the Integrated Microbiome Resource at Dalhousie University (Halifax, NS, Canada) for PCR amplification using primers 515FB and 926R, library preparation, and sequencing of the V4-V5 region of 16S rRNA genes. Amplicons were then sequenced on the Illumina MiSeq platform with 300 paired ends (PE). Demultiplexed reads were processed using QIIME2 (59). Likely contaminant amplicon sequence variants (ASVs) were identified by comparing the prevalence in DNA extracts from blank (negative control) samples versus milk DNA samples using decontam (60) and were removed from analysis. ASVs were assigned to taxonomic groups using the SILVA 138 SSURef NR99 database (61). ASVs classified as mitochondria, chloroplast, and Eukarya were also removed from analysis. Prior to normalization, cows given antibiotics (to treat clinical mastitis and other infections) were removed (see "Sample selection and data filtering" below for additional details). Also, samples with fewer reads than blanks (i.e., 49 or fewer reads) were removed and for the samples that were extracted in duplicate, the replicate with the higher number of reads was retained. Counts were normalized across the data set by rarefying without replacement at the 15th percentile (507 sequences per sample). Fresh milk was preserved with bronopol (2-bromo-2-nitropropane-1,3-diol) and sent to Valley Agricultural Software and AgSource laboratories for SCC measurements from around 30 mL of milk using a FOSSOMATIC 7 flow cytometer, and SCC was reported as counts per milliliter of milk.
Statistical analyses. Relationships between milk microbiota and lactose, milk microbiota and MOs, milk microbiota and MFAs, milk microbiota and SCC, lactose and SCC, MOs and SCC, and MFAs and SCC were explored (Table 1) in R version 4.0.2.
Correlations between microbial (ASV) b-diversity and MO b-diversity or between microbial (ASV) b-diversity and MFA b-diversity were determined within each period by relating the Bray-Curtis dissimilarity using the Mantel test and Procrustes analysis. These analyses were performed with the vegan package v 2.5-7 (62). P values were adjusted for multiple tests using false-discovery rate (FDR), and a cutoff of ,0.05 was used to report significance (Table 1).  Relationships between microbial (ASVs) a-diversity and a-diversity of MOs or MFAs was determined with linear mixed-effects modeling (LMM) and repeated-measures correlation (RMCORR) (63). LMM was performed with lmerTest package v 3.1-3 (64), and RMCORR was performed with the rmcorr package v 0.4.4 (65). LMM and RMCORR were also utilized to identify relationships between microbial a-diversity and single MOs, MFAs, lactose, or SCC. For LMM, the microbial a-diversity was analyzed as the response variable, while the MOs, MFAs, or lactose was analyzed as the predictor variables (fixed effects), and cow was included as a random effect (i.e., random intercept). For RMCORR, cow was set as a participant. Relationships between microbial a-diversity and SCC, and between microbial genera and SCC, were also investigated with LMM, but the microbiota were analyzed as the predictor variable and SCC as the response variable. Shannon and Simpson metrics were used to calculate a-diversity of MOs or MFAs. Shannon, Simpson, and Faith's phylogenetic diversity metrics were used to calculate a-diversity of microbial ASVs.
Relationships between rarefied abundances of microbial genera and abundances of MOs, MFAs, or lactose were identified with generalized linear mixed-effects modeling (GLMM) with Poisson distribution in which the microbial abundance was the response variable, the milk component (fatty acids, oligosaccharides, or lactose) abundance was the predictor variable, and cow was the random effect (i.e., random intercept). GLMM was performed with lme4 package v 1.1-27.1 (66). Relationships between the abundances of milk components and the abundances of milk microbes were investigated for genus-level microbial taxa in order to identify relationships with microbes at the lowest accurate taxonomic level that could be determined from the 16S rRNA V4-V5 gene region (67) since microbial ability to metabolize MOs can vary within taxa even down to species and strain levels (68), to enable interpretation in relation to correlations identified between SCC and milk microbes which have been implicated in mammary health and inflammation at the genus level and lower taxonomic levels (2), and to be comparable to similar studies in humans (12,14,16).
In our previous publication (6), we reported differential abundances of 11 bovine milk microbial genera with diet from this cohort of cows. While many of these diet-associated genera were also associated with abundance of certain MOs or MFAs in these analyses, the reverse was not true. In other words, most of the genera associated with MFA, MO, or lactose abundances were not differentially abundant with diet (see Results). We did not include diet or parity in our models because diet and parity have been found to associate with MFAs (58, 69), MOs (57,70), and lactose (71) in our cohort of cows and others; therefore, if diet and parity impact the milk microbes, then it may be through effects on the milk composition. With this reasoning, it would not be accurate to include diet and parity in our models as fixed effects completely independent of MFAs, MOs, or lactose. Since we have exactly three data points (one from each of the three periods and three diets) for each cow in our models, the effect of diet/period is similar for each cow. We have also included cow as a random effect in our model, which not only accounts for data points not being independent of each other but also allows us to relate the abundance of a milk component to the abundance of a milk microbe within each cow and averaged across cows. Ultimately, our objective was to determine how milk components and milk microbes correlate without necessarily identifying the isolated effect of milk components on milk microbes. Accordingly, we used parsimonious models to identify correlations between milk components and milk microbes.
When identifying relationships between SCC and MOs, MFAs, or lactose, RMCORR and LMM were utilized with SCC as the predictor, milk component as the response, and cow as random effect. The SCC measurements were right skewed and were consequently log 10 transformed before applying LMM or RMCORR (Table 1). Likewise, MO measurements were log 10 transformed before investigating relationships with log SCC but were not transformed for diversity calculations or relationships with microbial genus abundance. MFA concentrations were also transformed prior to LMM, RMCORR, or GLMM to achieve more normal distribution by applying the square root. MFA measurements also were not transformed prior to diversity measures.
Picrust2 (72), as a QIIME2 plugin, was employed to predict relative abundances of milk microbial b-galactosidase and 6-phospho-b-galactosidase from the ASVs. Correlations between lactose and predicted b-galactosidase relative abundance, or lactose and predicted b-galactosidase and 6-phosphob-galactosidase relative abundances combined, were investigated with LMM and RMCORR.
Sample selection and data filtering. Cows 5651 and 6229 were removed from MO and lactose analyses because they had lactose concentrations that were outliers, as determined based on the standard error with a cutoff of 3 (73). Ten cows (5020, 5212, 5297, 5405, 5697, 6076, 6215, 6232, 6233, and 6241) were removed from analyses because they were treated with antibiotics during the study. Eight cows (5053, 5464, 5677, 6206, 6210, 6211, 6214, and 6254) stole food and were retained in all data sets, but milk samples from these cows were not analyzed for MO abundance. Ultimately, 19 MOs, 77 MFAs, and 1,327 genera were measured and analyzed (Supplemental File 4). The parity and DIM of cows before and after these filtering criteria remained similar (Fig. S1). After matching data sets by sample collection date, cows with less than one sample per period were retained for b-diversity analyses since the analyses were performed separately for each period, but these cows were omitted from analyses using LMM, GLMM, and RMCORR.
For analyses of relationships between milk components and SCC, we hypothesized that SCC predicts the milk component abundance; therefore, the SCC samples from 21  Microbial genera, MOs, or MFAs present in less than 50% of samples and less than 90% of cows were not included in analyses, except for b-diversity and a-diversity calculations. Genera that met the prevalence cutoffs differed slightly for each pair of data analyzed (Table S1).
The Spearman rank correlations among all 44 milk microbial genera listed in Table S1 were calculated in R with the function cor.test, and the Spearman rho for each correlation was visualized in a heat map (Fig. S5).
Data availability. The data sets and code used in our analysis are publicly available at the GitHub repository: https://github.com/L-Coates/DGC_study1_MilkMicrobes_MOs_MFAs_lactose_SCC_Associations .git. The milk microbiota data are available at Qiita (study identifier [ID]