- Open Access
Association of bovine major histocompatibility complex (BoLA) gene polymorphism with colostrum and milk microbiota of dairy cows during the first week of lactation
Microbiomevolume 6, Article number: 203 (2018)
The interplay between host genotype and commensal microbiota at different body sites can have important implications for health and disease. In dairy cows, polymorphism of bovine major histocompatibility complex (BoLA) gene has been associated with susceptibility to several infectious diseases, most importantly mastitis. However, mechanisms underlying this association are yet poorly understood. In the present study, we sought to explore the association of BoLA gene polymorphism with the dynamics of mammary microbiota during the first week of lactation.
Colostrum and milk samples were collected from multiparous Holstein dairy cows at the day of calving and days 1 and 6 after calving. Microbiota profiling was performed using high-throughput sequencing of the V1-V2 regions of the bacterial 16S rRNA genes and ITS2 region of the fungal ribosomal DNA. Polymorphism of BoLA genes was determined using PCR-RFLP of exon 2 of the BoLA-DRB3. In general, transition from colostrum to milk resulted in increased species richness and diversity of both bacterial and fungal communities. The most dominant members of intramammary microbiota included Staphylococcus, Ruminococcaceae, and Clostridiales within the bacterial community and Alternaria, Aspergillus, Candida, and Cryptococcus within the fungal community. Comparing the composition of intramammary microbiota between identified BoLA-DRB3.2 variants (n = 2) revealed distinct clustering pattern on day 0, whereas this effect was not significant on the microbiota of milk samples collected on subsequent days. On day 0, proportions of several non-aureus Staphylococcus (NAS) OTUs, including those aligned to Staphylococcus equorum, Staphylococcus gallinarum, Staphylococcus sciuri, and Staphylococcus haemolyticus, were enriched within the microbiota of one of the BoLA-DRB3.2 variants, whereas lactic acid bacteria (LAB) including Lactobacillus and Enterococcus were enriched within the colostrum microbiota of the other variant.
Our results suggest a potential role for BoLA-gene polymorphism in modulating the composition of colostrum microbiota in dairy cows. Determining whether BoLA-mediated shifts in the composition of colostrum microbiota are regulated directly by immune system or indirectly by microbiota-derived colonization resistant can have important implications for future development of preventive/therapeutic strategies for controlling mastitis.
Different body sites of vertebrates provide stable and nutrient-rich ecosystems for a diverse range of microbial symbionts to thrive [1, 2]. Microbial symbionts, comprised of different species of bacteria, fungi, and archaea, in turn perform a number of functions that are essential for host physiology, including digestion of complex carbohydrates, biosynthesis of vitamins, and modulation of immune homeostasis [3, 4]. The relationship between host animal and microbiota, however, is not always beneficial. Indeed, host-associated microbiota usually contain a number of opportunistic members (i.e., pathobionts), which under certain circumstances, such as immunosuppression of the host, are capable of dominating the microbiota and disrupting the immune homeostasis at different body sites . Moreover, at the interface with host immune system, commensal microbes can stimulate the activation of a number of pro-inflammatory responses . In order to cope with this constant immunological challenge while maintaining the equilibrium between pro- and anti-inflammatory responses, the immune system of mammals has adapted specific mechanisms for discriminating between commensals, pathobionts, and strict pathogens .
An example of this evolutionary adaptation is the polymorphism of the class II region of the major histocompatibility complex (MHC) gene . This gene encodes highly polymorphic glycoproteins on the surface of the professional antigen-presenting cells (APCs) of the immune system, including macrophages, dendritic cells, and B lymphocytes . The main function of the MHC class II glycoproteins is to facilitate the watchdog role of the immune system in recognition of various antigens and triggering pro- and anti-inflammatory responses accordingly . Following internalization by APCs, exogenous antigens are degraded into smaller peptides that can bind to the antigen-binding groove of the MHC class II glycoproteins. The resulting peptide-MHC complex is then transported to the plasma membrane of the APCs and presented to the receptors of the CD4+ helper T lymphocytes . Depending on the nature of the antigens presented, activated T lymphocytes will be differentiated into effector TH1, TH2, or TH17 cells, which can in turn trigger a number of pro- and anti-inflammatory responses . MHC heterozygosity and/or certain allelic variations within MHC gene can influence the diversity of the T lymphocyte receptors and therefore confer resistance to microbial colonization by virtue of their ability to trigger immune responses against a wider range of antigens . Indeed, several studies in humans , mice [12, 13], and fish  have demonstrated that MHC polymorphism can influence the composition of the microbiota that inhabit different niches of the gastrointestinal tract.
In cattle, similar to other species, the physiologic relevance of MHC gene polymorphism has been traditionally viewed from the perspective of host-pathogen interaction and infectious tolerance. In particular, several studies have focused on exploring the associations of allelic variants of MHC gene family, also known as bovine leukocyte antigens (BoLA) genes, with susceptibility to infectious diseases caused by pathogenic bacteria, viruses, and parasites, among others [15,16,17]. Particularly, in the context of bovine mastitis, which is the most prevalent and costly disease of dairy cows [18, 19], several investigations have detected associations between certain allelic variants of exon 2 of the BoLA-DRB3.2 gene and resistance to intramammary infection (IMI) by major mastitis pathogens and/or severity of the inflammatory responses of the mammary gland (MG) [20,21,22,23,24]. Nonetheless, direct and indirect mechanisms by which BoLA gene polymorphism contributes to the immune homeostasis of the MG have not yet been described in detail. Recent metagenomics studies have suggested that in addition to major mastitis pathogens, which usually possess a variety of virulence factors evolved to resist the defense mechanisms of the MG [25, 26], intramammary secretions of dairy cows can harbor a broad range of non-pathogenic bacterial groups [27,28,29,30,31,32]. However, whether BoLA gene polymorphism can modulate the composition of the intramammary microbiota remains unknown. Indeed, the paucity of evidence for crosstalk between MG microbiota and cows’ genotype has raised concerns regarding the validity of the concept of intramammary commensal microbiota . Hence, the main objective of the present study was to evaluate the potential association of BoLA-DRB3 polymorphism with the composition and dynamics of bacterial and fungal components of the intramammary microbiota of dairy cows during the first week of lactation, coinciding with the transition from colostrum to milk.
Identification of BoLA-DRB3.2 variants and sequencing results
Based on the digestion activity of the BstYI enzyme, three distinct BoLA-DRB3.2 variants were identified within the studied population of cows (Additional file 1: Figure S1), named herein BstYa (n = 24), BstYb (n = 25), and BstYc (n = 5). Due to the limited number of cows that were categorized under the BstYc variant, and in order to maintain statistical power of between-variants comparisons, all the colostrum and milk samples belonging to this variant were excluded from downstream microbial analyses. For BstYa and BstYb variants, after quality filtering of the sequencing reads and removing samples with low sequencing depth (< 5000 sequence per sample), a total of 132 and 82 samples were kept for downstream analyses of bacterial and fungal communities. On average, 28,727 (6889–100,740; SD = 13,921) and 35,328 (8657–119,756; SD = 17,175) high-quality 16S rRNA and ITS2 sequences were obtained per sample, resulting in identification of 590 (83–1155; SD = 290) and 24 (6–80; SD = 13) representative bacterial and fungal OTUs, respectively.
Biodiversity dynamics of intramammary secretions during the first week of lactation: associations with BoLA-DRB3.2 polymorphism
Overall, diversity of intramammary microbiota (Shannon’s index) underwent significant changes during the first week of lactation (Fig. 1a, b). Diversity of bacterial communities increased over time (p < 0.001), with day 6 samples harboring the most diverse communities across all three sampling time points. On the other hand, diversity of fungal communities was higher on day 1 compared to day 0 (p = 0.002) and day 6 (p < 0.001). The association between BoLA-DRB3.2 polymorphism and the diversity of intramammary microbiota was only significant on day 0, when colostrum samples of BstYb variant harbored bacterial communities that were more diverse than those of the samples belonging to BstYa variant (p = 0.002). There was no association between BoLA-DRB3.2 polymorphism and Shannon’s index of diversity of fungal communities (p = 0.082; summary statistics for alpha-diversity comparisons available in Additional file 2: Table S1).
PCoA and permutational multivariate analysis of variance (PERMANOVA) of Bray-Curtis dissimilarities revealed that the composition of both bacterial (ADONIS R2 = 0.19, p(PERMANOVA) < 0.001; Fig. 2a) and fungal (ADONIS R2 = 0.11, p(PERMANOVA) < 0.001; Fig. 3a) communities of intramammary secretions differed among all sampling time points. The association between BoLA-DRB3.2 polymorphism and the composition of bacterial communities was only significant on day 0, when colostrum samples belonging to BstYa and BstYb variants showed distinct clustering pattern (ADONIS R2 = 0.22, p(PERMANOVA) < 0.001; Fig. 2b–d). BoLA-DRB3.2 polymorphism was also associated with the composition of fungal communities of intramammary secretions on day 0 (ADONIS R2 = 0.06, p(PERMANOVA) = 0.070) and day 1 (ADONIS R2 = 0.06, p(PERMANOVA) = 0.032), but not on day 6 (ADONIS R2 = 0.04, p(PERMANOVA) = 0.226; Fig. 3b–d). PCoA and PERMANOVA based on Jaccard binary distances of bacterial and fungal communities also revealed similar clustering patterns based on sampling time points and BoLA-DRB3.2 polymorphism (Additional file 1: Figures S2 and S3). Summary statistics for PERMANOVA analyses are presented in Additional file 2: Table S2.
Dynamics of intramammary microbiota in relation to days in milk and BoLA-DRB3.2 polymorphism
The bacterial communities of colostrum and milk samples were mainly composed of the members of phyla Firmicutes, Proteobacteria, Bacteroidetes, Actinobacteria, and Fusobacteria. Dominant bacterial genera of mammary microbiota had different trajectory patterns over the first week of lactation and in relation to BoLA variants (Fig. 4a, b). The proportion of Staphylococcus, the predominant genus within the colostrum microbiota of both BstY variants, decreased during transition from colostrum to milk. Similarly, proportions of Streptococcus, Fusobacterium, Enterococcus, and Bacteroides, which were all overrepresented within the microbiota of BstYa variant on day 0, decreased significantly during the first week of lactation and became proportionally similar between the microbiota of the two BoLA variants on day 6. On the other hand, proportions of Sphingobacterium and Acinetobacter increased during transition from colostrum to milk within the microbiota of the two BoLA variants. Summary statistics for all pairwise contrasts are presented in Additional file 2: Table S3.
The fungal communities of colostrum and milk samples were mainly composed of phylum Ascomycota, followed by considerably lower proportions of genera belonging to phyla Basidiomycota and Zygomycota. Similar to bacterial community, dominant fungal genera of mammary secretions showed different trajectory patterns over the first week of lactation and in relation to BoLA variants (Fig. 5a, b). For example, Alternaria, the predominant genus within day 0 colostrum microbiota of both BstY variants, showed a decreasing pattern during the first week of lactation in the microbiota of BstYa variant, whereas the proportion of this genus remained unchanged in the microbiota of BstYb variant during the same time period. On the other hand, proportions of genera Aspergillus and Cryptococcus showed similar decreasing pattern following transition from colostrum to day 6 milk in BstYb variant. The proportion of Candida increased from day 1 to day 6 in the microbiota of both BoLA variants. Summary statistics for all pairwise contrasts are presented in Additional file 2: Table S4.
Association of BoLA-DRB3.2 polymorphism with the composition of day 0 colostrum microbiota
Unsupervised hierarchical cluster analysis (UPGMA) based on the proportions of most abundant bacterial OTUs (> 0.1% of the community) was performed to test whether samples belonging to the two BoLA variants tend to cluster distinctly based only on the proportion of abundant bacterial OTUs. This resulted in generation of two main clusters of samples corresponding to the two BoLA variants. In order to test the significance of clustering pattern of the two BoLA variants based on the proportion of abundant OTUs, the Bray_Curtis dissimilarity matrix of samples was calculated and subjected to PERMANOVA, resulting in identification of distinct clustering pattern (p(PERMANOVA) = 0.002; Fig. 6). Statistical comparison (LEfSe) at the OTU level revealed that the proportions of several non-aureus staphylococci (NAS) OTUs, including OTU109 (S. equorum), OTU3 (S. gallinarum), OTU9 (S. sciuri), and OTU21 (S. haemolyticus), Streptococcus OTU124 (S. equinus), actinobacterial OTUs belonging to Corynebacterium and Mycobacterium, and proteobacterial OTUs belonging to udder-associated opportunistic/pathogenic bacteria such as Stenotrophomonas, Pseudomonas, Acinetobacter, and Alcaligenaceae were all enriched within the colostrum microbiota of day 0 colostrum samples of BstYb variant (Fig. 6; NCBI BLASTN bit-scores for species level classification are presented in Additional file 2: Table S5). On the other hand, day 0 colostrum samples of BstYa variant were enriched with Streptococcus OTU1 (S. dysgalactiae), Staphylococcus OTU2 (S. chromogenes), several Fusobacterium OTUs, and OTUs belonging to lactic acid bacteria (LAB) including Lactobacillus and Enterococcus. Other than abovementioned bacterial OTUs, the proportion of several gut-associated bacterial OTUs were also enriched within the colostrum microbiota of each BoLA-DRB3.2 variant; Firmicutes OTUs belonging to Butyrivibrio and Clostridiales were enriched within the colostrum microbiota of BstYb variant, whereas Bacteroidetes OTUs belonging to Bacteroides and Prevotella were overrepresented within the microbiota of BstYa variant (Fig. 5).
UPGMA based on the proportion of abundant fungal OTUs (> 0.1% of community) revealed no distinct clustering pattern between the mycobiota of different BoLA-DRB3.2 variants (Additional file 1: Figures S6 and S7). This was further confirmed by complementary PERMANOVA test on the Bray-Curtis dissimilarity matrix of samples based on the proportion of abundant OTUs (day 0, p(PERMANOVA) = 0.098; day 1, p(PERMANOVA) = 0.25). Nonetheless, statistical comparison (LEfSe) revealed enrichment of several fungal OTUs within the mycobiota of each variant. On day 0, the proportions of several Penicillium OTUs (including OTU32: Penicillium ilerdanum, OTU22: Penicillium psychrosexualis, and OTU556: Penicillium chrysogenum) were enriched within the colostrum microbiota of BstYa variant (Additional file 1: Figure S6). On the other hand, on day 1, proportions of Alternaria OTUs (including OTU6: Alternaria infectoria and OTU177: Unclassified Alternaria) were higher within the colostrum microbiota of BstYb variant (Additional file 1: Figure S7).
Association of BoLA-DRB3.2 polymorphism with co-occurrence patterns of intramammary bacterial communities
Correlation network analysis was further performed to evaluate potential association of BoLA-DRB3.2 polymorphism with co-occurrence patterns of intramammary bacterial communities (Fig. 7a–c, and Additional file 1: Figures S4 and S5a–f). By using an ensemble method (a combination of diverse measures of correlation, i.e., Pearson and Spearman, and dissimilarity, i.e., Bray-Curtis and Kullback-Leibler), we were able to track changes in the co-occurrence patterns of bacterial OTUs within the microbiota of each BoLA variant at different time points. Similar to the trends observed for alpha- and beta-diversity analyses, the association of BoLA-DRB3.2 with the structure of microbe-microbe relationships was more pronounced within the microbiota of day 0 colostrum samples than day 1 or 6. On day 0, Firmicutes OTUs played central role in the overall structure of the bacterial community of BstYb variant, while the contributions of OTUs from other bacterial phyla were negligible. Indeed, with the exception of two OTUs belonging to Fusobacterium and Bacteroides, the entire hub OTUs (i.e., OTUs showing the highest number of positive/negative connections with other members of the community) within day 0 microbiota of BstYb variant belonged to Firmicutes, more specifically, genera Streptococcus, Butyrivibrio, Enterococcus, Helcococcus, and unclassified Clostridiaceae. On the other hand, within day 0 microbiota of BstYa variant, Firmicutes and Bacteroidetes OTUs showed equal number of connections with other members of the community. Hub OTUs within this group belonged to Firmicutes (including genera Oscillospira, unclassified Clostridiaceae, and unclassified Ruminococcaceae), Bacteroidetes (including genera Bacteroides, Prevotella, S24–7, and unclassified Bacteroidales), and Fusobacterium.
Following the transition from colostrum to milk, the proportions of main bacterial phyla within the microbiota of the two BoLA-DRB3.2 variants became similar to each other, and so did the overall co-occurrence pattern of different bacterial OTUs belonging to these phyla (Fig. 7a–c). Although Firmicutes continued to remain the most dominant bacterial phylum within intramammary microbiota of both BoLA-DRB3.2 variants, overall connectedness of Bacteroidetes and Proteobacteria OTUs took over those of the Firmicutes OTUs by progression of DIM. In addition, during this transitional period, the total number of negative connections (mutual exclusion relationships) of intramammary bacterial communities was generally reduced. On day 6, the hub OTUs within the bacterial communities of both variants became more or less similar to each other, with several Sphingobacterium OTUs, including OTU75 (S. daejeonense) and OTU45 (S. alimentarium) showing the highest number of connections within both communities. Nonetheless, day 6 microbiota of each variant also contained few exclusive hub OTUs, including Staphylococcus (OTU9) within the microbiota of BstYa variant and Microbacterium (OTU103) and Acinetobacter (OTU25) within the microbiota of BstYb variant (Additional file 1: Figure S4).
Although the impacts of environmental drivers, such as housing conditions and antimicrobial use, on the composition of bovine milk microbiota have been described [29, 34,35,36], it remains unknown whether cows’ genotype can also play a role in modulating the diversity of microbial communities of intramammary secretions. In the present study, we demonstrated that BoLA-DRB3.2 polymorphism is strongly associated with the composition of both bacterial and fungal communities of intramammary secretions during the first week of lactation. Particularly, we observed that the distribution of dominant bacterial OTUs were widely different between the colostrum microbiota of the two identified BoLA-DRB3.2 variants, suggesting that variations in BoLA class II molecules may result in colonization resistance against certain microbial lineages within the MG ecosystem. In addition, by characterizing the microbiota profile of intramammary secretions collected at multiple sampling time points during the first week of lactation, our study provides novel insights into the dynamics of bacterial and fungal communities of intramammary secretions during transition from colostrum to milk.
Association between BoLA polymorphism and intramammary microbiota
The contribution of local immunity of MG in modifying intramammary microbiota is unknown. MHC polymorphism can influence the composition of mucosa-associated microbiota via inducing biases in the selective activity of secretory IgAs against bacterial surface epitopes [13, 37]. Secretory IgA is the predominant immune component of the mucus layer responsible for exclusion of exogenous antigens . However, due to the very low concentrations of IgA in bovine colostrum and milk (approximately 3.5 and 0.13 mg/mL, respectively ), it is less likely that this antibody can serve as the main driver of BoLA-mediated shifts in the composition of mammary microbiota. Conversely, bovine colostrum contains high concentrations of IgG (32–212 mg/mL) . IgG can participate in MHC-mediated selective immune responses against exogenous antigens via inactivation of bacterial antigens by binding to specific surface sites, bacterial opsonization, and complement activation [40, 41]. IgG1, the predominant subclass of IgG in bovine colostrum (~ 80% of total Igs), is mainly derived from blood, whereas IgG2, as well as IgA and IgM, are mainly synthesized by plasma cells or epithelial cells of the MG . In the present work, the influence of BoLA-DRB3.2 polymorphism on the composition of intramammary microbiota was more evident on day 0 colostrum samples and faded abruptly thereafter. This trend is very similar to peripartal fluctuations in the concentration of intramammary IgG1. During the prepartum period, IgG1 is gradually accumulated in colostrum in order to prepare for passive transfer of maternal immunity to newborn calves ; however, transition from colostrum to milk is accompanied by a rapid drop in intramammary concentrations of IgG1 [39, 44]. An important limitation of the present study was the lack of measurement of immunomodulatory components of intramammary secretions, most importantly IgGs concentrations. Nonetheless, we speculate that during the prepartum period, prolonged isolation of mammary gland from environmental sources of microbes (i.e., due to teat canal closure) provides an ideal ecosystem for IgG1 to exert BoLA-mediated shifts in the composition of persistent intramammary microbiota. Nonetheless, BoLA class II expression by epithelial cells of the MG and/or intramammary lymphocytes also has the potential to influence local immune responses and therefore modify the composition of microbiota. To what extent each of the abovementioned systemic and local immune responses contribute to regulation of intramammary microbiota remains largely unknown.
The main biologic relevance of BoLA polymorphism for udder health is modulation of mastitis susceptibility. In dairy cows, allelic variations of BoLA-DRB3.2 have been associated with susceptibility to IMI by mastitis pathogens and/or severity of mastitis [20,21,22,23,24]. Yet, underlying mechanisms remain unknown. Recently, Kubinak et al.  demonstrated that MHC-mediated shifts in the composition of mucosa-associated microbiota can result in colonization resistant against Salmonella enterica serovar Typhimurium infection in mice. The authors suggested an indirect mechanism by which MHC gene polymorphism can modulate infection susceptibility by altering the robustness of the intestinal microbiota rather than direct involvement of the immune system. In the present study, increased proportions of LAB, especially Lactobacillus and Enterococcus, within the colostrum microbiota of BstYa variant was concomitant with decreased proportions of several important bacterial groups that are commonly regarded as mastitis pathogens or opportunists, including Staphylococcus, Streptococcus, Mycobacterium, and Corynebacterium. With the exception of few pathogenic strains, Lactobacillus and Enterococcus strains isolated from bovine milk are commonly regarded as beneficial bacteria capable to inhibit the growth of mastitis pathogens via the production of inhibitory substances such as bacteriocins, organic acids, and hydrogen peroxide [45,46,47]. On the other hand, corynebacterial and staphylococcal strains isolated from bovine milk, particularly NAS, are also capable to produce a wide range of antimicrobial molecules [48, 49]. Our results suggest that each of the BoLA-DRB3.2 variants may promote the growth of different NAS species; S. equorum, S. gallinarum, S. sciuri, and S. haemolyticus were all enriched within the colostrum microbiota of BstYb variant, whereas S. chromogenes was enriched within the colostrum microbiota of a sub-group of cows within the BstYa variant. Altogether, these observations may be indicative of development of specific patterns of microbiota-derived colonization resistant that are mediated by BoLA-DRB3.2 polymorphism. Nonetheless, it is also possible that variations in MHC class II expression directly influence the composition of microbiota by activating different profiles of pro- and anti-inflammatory responses and therefore exerting positive (tolerance) or negative (clearance) effect on different bacterial groups . Given the important implications that microbiota-derived resistance against colonization by mastitis pathogens can have for development of alternative preventive therapeutics, this hypothesis warrants further mechanistic investigations of potential links among BoLA gene polymorphisms, immunomodulatory components of the immune system, and intramammary microbiota.
An important caveat of the present study was the inability to assess potential impacts of BoLA-DRB3.2 polymorphism on the stability and resilience of the intramammary microbiota beyond the first week of lactation. Indeed, the main objective of the present study was to provide evidence for potential links between cows’ genotype and intramammary microbiota in the first place. Otherwise, assessment of the impact of BoLA-DRB3.2 polymorphism on the resilience of intramammary microbiota would have required quarter-based longitudinal milk sampling over a relatively long period of time during early and mid stage of lactation. Notwithstanding this limitation, comparison of the interrelationship networks of the microbial communities during the first week of lactation suggests that BoLA-mediated effects on the structure of intramammary microbiota are not long-lasting. It appears that during the colostrogenesis period, when accumulation of immune components results in a high selective pressure on the microbial ecosystem of the MG, community assembly process follows a “habitat filtering” pattern in which phylogenetically related species that possess specific traits required for survival under certain environmental conditions (i.e., the selective pressure of the immune system) tends to co-occur . This might be the reason why day 0 microbiota of BoLA-DRB3.2 variants showed distinct clustering patterns, each orchestrated by a completely different set of hub OTUs. It is also noteworthy to mention that other than potential contribution of selective pressure (i.e., immunomodulatory components of the milk), other stochastic ecological processes, such as drift , might have also contributed to shaping the composition of the isolated microbial ecosystem of the mammary gland during the dry period.
Dynamics of intramammary microbiota during transition from colostrum to milk
During the early stages of lactation, the changeover from colostrum to transitional and then mature milk is accompanied by considerable shifts not only in the concentrations of immune components [39, 44] but also other bioactive compounds such as milk oligosaccharides . Many of these factors have the ability to impact the composition of mammary and intestinal microbiota [53, 54] via several mechanisms such as the ability of milk oligosaccharides to compete with microbes over attachment to epithelial cell surface receptors and therefore confer colonization resistance . In addition, lactation increases the exposure of intramammary ecosystem to environmental sources of microbes including, among others, milking equipment and bedding materials. These environmental factors have also proven to be influential in modifying the composition of intramammary microbiota [34, 56]. Prior to entering the sampling procedure of the present study, selected cows were subjected to blanket dry cow therapy using a combination of internal teat sealant and long-lasting antimicrobial. Not surprisingly, we observed that transition from colostrum to milk, which coincided with the exposure of isolated mammary gland to environmental sources of microbes (i.e., milking equipment, bedding material, etc.), resulted in significant increases in the diversity of both bacterial and fungal components of the intramammary microbiota. Regarding the bacterial communities of colostrum microbiota, our results are in partial agreement with the observations of Lima et al.  who also reported the predominance of Staphylococcus and gut-associated bacteria including Ruminococcaceae and Clostridiales in the colostrum microbiota of both primiparous and multiparous cows. Staphylococcus spp., in particular NAS, are predominant colonizer of bovine MG commonly considered as the main causes of subclinical mastitis . In the present study, we observed that Staphylococcus was the most abundant bacterial genera throughout the first week of lactation within the microbiota of both BoLA variants. Fusobacterium was another dominant member of the colostrum microbiota the proportion of which remained unchanged during the first week of lactation. Lima et al.  reported that the colostrum microbiota of MG quarters that had mastitis during the early stages of lactation contained a higher proportion of Fusobacterium and were significantly less diverse than the microbiota of quarters that did not experience clinical mastitis. In our study, the colostrum microbiota of the BstYa variant, which contained a higher proportion of phylum Fusobacteria, was significantly less diverse than BstYb variant. However, the composite milk sampling procedure of our study precluded us form relating this observation with future incidences of clinical mastitis. Another notable dynamic in the bacterial communities during transition from colostrum to milk was enrichment of the genus Sphingobacterium. Indeed, our correlation network analysis revealed that on day 6, OTUs belonging to this genus were among the most influential (hub) OTUs in shaping the overall structure of the milk microbiota. This is in agreement with the result of our parallel study , where we observed positive associations between the proportions of hub Sphingobacterial OTUs and teat-end hyperkeratosis or future incidences of clinical mastitis. Oikonomou et al.  also reported a strong correlation between the proportion of Sphingobacterium and increased milk somatic cell count. Others [59, 60] have also detected members of this bacterial genus in the milk samples obtained from clinical and subclinical cases of mastitis. All together, these observations suggest that Sphingobacterium is a potent colonizer of the MG and perhaps an emerging mastitis pathogen in different countries.
Our study provides the first report on the composition of fungal microbiota of intramammary secretions. Fungi are ubiquitously present in dairy environment and therefore cross-contamination of samples with environmental fungi through air-born spores and/or aerosols during sampling might be inevitable under field conditions. Nonetheless, our results indicated significant impact of BoLA-DRB3.2 polymorphism on the composition of fungal component of colostrum microbiota. This observation suggests the existence of a resident, perhaps commensal, community of fungi within the MG ecosystem. The prevalence of mycotic mastitis is relatively low (~ 7–12% [61, 62]); however, strains belonging to genera Candida, Cryptococcus, Alternaria, and Aspergillus have been commonly isolated from bovine milk [62,63,64,65]. Here, we also observed that Alternaria, Aspergillus, Candida, and Cryptococcus were the most abundant fungal components of intramammary microbiota during the first week of lactation. Compared to the bacterial community, the composition of fungal community was overall more stable during transition from colostrum to milk. Scedosporium was the only abundant fungal genera the proportion of which increased during this transitional phase. Scedosporium spp. have been characterized as opportunistic pathogens with the ability to cause infection in immunocompromised animals and humans [66, 67]. Our observation suggests that this genus might be an environmental opportunistic that is capable to colonize the MG after the dry period and following exposure of the intramammary ecosystem to the environmental sources.
The present study provides novel insights into the influence of BoLA-DRB3.2 polymorphism on the composition of both bacterial and fungal components of intramammary secretions. Overall, our results indicated that BoLA-mediated shifts were more pronounced on the microbiota composition of colostrum samples rather than milk. Particularly, we observed that within colostrum microbiota of each BoLA-DRB3.2 variant, specific groups of bacteria that are known to exert colonization resistance via production of antimicrobial peptides were enriched. Determining whether BoLA-mediated shifts in the composition of intramammary microbiota are regulated directly by immune system or indirectly by microbiota-derived colonization resistance can have important implications for future development of preventive/therapeutic strategies for controlling mastitis. Even if this modulatory impact is restricted to the colostrogenesis period, it can still have great implications for udder health and mastitis susceptibility as a high incidence rate of peripartal mastitis remains a major problem for dairy industry worldwide.
Selection of animals and sample collection
This study was conducted in a commercial 500-head milking dairy farm located in Manitoba, Canada, from December 2014 to February 2015. Late pregnant multiparous cows were identified during their last 3 weeks of pregnancy and deemed eligible for gradual enrollment in the sampling procedure of the present study based on the following criteria: presence of four quarters with healthy appearance (no visible sign of clinical mastitis such as swelling or redness, and devoid of anatomically damaged teat ends), and no incidence of clinical mastitis during the last 90 days of the previous lactation. At the end of the previous lactation cycle, all cows were subjected to blanket dry cow therapy using internal teat sealant (Orbeseal, Zoetis, CA, USA) in combination with long-acting antimicrobial product containing 200,000 units penicillin G and 400 mg novobiocin per 10 mL tube (NovaDry; Pfizer Animal Health, QC, Canada). During the late-pregnancy period till 3 days post-calving, cows were housed in designated transition pens, bedded with fresh and dry straw. Afterwards, all cows were transferred to free-stall pens bedded every other day with recycled bedding material. Following parturition, cows were milked three times a day at 4 am, noon, and 8 pm. Cows selected for the purpose of the present study (n = 60) were subjected to three sampling time points during the first week of lactation:( a) fresh colostrum collected within 6 h of calving (day 0), (b) colostrum collected at noon of the first day post-calving (day 1), and (c) fresh milk collected at noon of the sixth day post-calving (day 6). On the day of sampling, all four quarters of the MG of each cow were checked for clinical signs of mastitis (swelling or redness). Prior to sampling, four streams of foremilk from all quarters were discarded to (a) minimize chances of sample contamination from bacteria colonizing the teat canal, and (b) to check for abnormal appearance of colostrum/milk (i.e., watery secretions, presence of blood, flake, or abnormal color). Cows diagnosed with clinical mastitis at any sampling time point were excluded from the study (n = 6). For the remaining cows (n = 54), pre-milking teat disinfection was performed using 0.5% iodine pre-dip solution, teats were thoroughly dried using individual paper towels and then scrubbed for 15 s using cotton pads moistened in 70% alcohol. Composite colostrum and milk samples (~ 10 mL from each quarter combined into a 50 mL sterile containers) were then collected and placed on ice in the milking parlor, temporary stored at − 20 °C on the farm, and transferred over dry-ice to the laboratory for further processing.
DNA extraction from colostrum and milk samples
Colostrum and milk samples collected from 54 multiparous cows were subjected to genomic DNA extraction using the ZR-96 Fungal/Bacterial DNA Kit (Zymo Research, Irvine, CA, USA) following modified protocols of the manufacturer as follows: 1.5 mL of each sample was centrifuged at 12,000×g for 20 min at 4 °C. Supernatants were carefully removed and 200 μL of TE buffer and 300 μL of 0.5 M EDTA (pH = 8) were added to the pellet. The mixture was incubated for 20 min at room temperature and again centrifuged at 12,000×g for 10 min. Supernatants were removed and pellets were resuspended by adding 200 μL of PBS buffer and vortexing for 30 s. Next, 1 g of autoclaved 0.5 mm silica beads, 400 μL of Lysis Solution (Zymo Research, CA, USA), and 18 μL of 20 mg/mL Proteinase K (Zymo Research, CA, USA) were added to each tube and incubated in a heated shaker at 45 °C for 45 min. Tubes were then vortexed for 2 min using a 2010 GenoGrinder (SPEX SamplePrep, NJ, USA) at 1700 stroke per minutes. Then, 400 μL of the resulting mixture was then transferred to the deep-well plate of the Fungal/Bacterial DNA Kit and extraction process continued following manufacturer protocol. Negative controls (n = 3) were included in the extraction process by replacing 1 mL of sterile water instead of colostrum/milk samples.
PCR-RFLP identification of the BoLA-DRB3.2 variants
Genomic variations within the exon 2 of the BoLA-DRB3 were determined according to the hemi-nested PCR-RFLP method described by van Eijk et al.  with minor modifications. For the first round of amplification, PCR reactions contained 3.0 μL of extracted DNA from day 0 colostrum samples, 1.0 μL (5 μM) of each of the primers HLO30 (5′-ATC CTC TCT CTG CAG CAC ATT TCC-3′) and HLO31 (5′-TTT AAA TTC GCG CTC ACC TCG CCG CT-3′), 0.4 μL (20 mg/mL) BSA (ThermoFisher Scientific, Mississauga, ON, Canada), 9.6 μL nuclease-free water (ThermoFisher Scientific), and 10 μL of 5 Prime Hot MasterMix (5 Prime Inc., Gaithersburg, MD, USA). PCR conditions were set up as follows: an initial denaturing step at 94 °C for 3 min followed by ten amplification cycles at 94 °C for 25 s, 60 °C for 30 s, and 72 °C for 30 s with a final extension step at 72 °C for 5 min in an Eppendorf Mastercycler pro (Eppendorf, Hamburg, Germany). A second hemi-nested amplification round was performed using 3.0 μL of the PCR product of the first round as input genomic material. Reaction ingredients were similar to the first round except that primer HLO32 (5′-TCG CCG CTG CAC AGT GAA ACT CTC-3′) was substituted for HLO31. Reaction conditions consisted of 25 amplification cycles at 94 °C for 40 s and 65 °C for 30 s with a single final extension of 72 °C for 5 min. Further, 10 μL of the PCR products of the second amplification round were incubated with 5 units of BstYI restriction enzyme (New England Biolabs Inc., MA, USA) at 60 °C for 70 min. After incubation, reactions were incubated at 85 °C for 5 min to denature the enzyme. DNA fragments were visually detected by conducting vertical electrophoresis on Mini-Protein III unit (Bio-Rad Laboratories, Mississauga, ON, Canada) using 10% Mini-PROTEAN TBE gels (Bio-Rad Laboratories) at 300 mA for 45 min.
PCR amplification and construction of sequencing libraries
Samples were subjected to two separate PCR rounds targeted to amplify the V1-V2 regions of the bacterial 16S rRNA genes and the internal transcribed spacer 2 (ITS2) of the eukaryotic ribosomal genes (the sequence of primers used for PCR amplification and sequencing reactions can be found in Additional file 2: Tables S6 and S7). The forward PCR primers of the V1-V2 primer set and the reverse PCR primer of the ITS2 primer set were indexed with 12-base Golay barcodes, allowing for multiplexing of samples. For each sample, the PCR reaction targeting each region was performed in duplicate and contained 3.0 μL of extracted metagenomic DNA, 1.0 μL of each forward and reverse primer (5 μM), 0.4 μL of 20 mg/mL BSA (ThermoFisher Scientific), 11.6 μL nuclease-free water (ThermoFisher Scientific), and 10 μL of 5 Prime Hot MasterMix (5 Prime Inc.). A total of three negative controls were also included in PCR reactions by replacing metagenomic DNA with 3 μL of nuclease-free water. The PCR conditions for the V1-V2 primer set consisted of an initial denaturing step at 94 °C for 3 min followed by 35 amplification cycles at 94 °C for 30 s, 55 °C for 20 s, and 72 °C for 20 s, with a final extension step at 72 °C for 5 min in an Eppendorf Mastercycler pro (Eppendorf, Hamburg, Germany). The PCR conditions for the ITS2 primer set were identical to those used for V1-V2 primer set with the exception of annealing temperature being set at 52 °C for 30 s. The sequencing libraries were then generated as explained by Derakhshani et al.  and sequenced using MiSeq Reagent Kits V3 (600-cycle; Illumina, San Diego, CA, USA) at the Gut Microbiome Laboratory at the University of Manitoba.
Bioinformatics and statistical analyses
The default settings of FLASH assembler ver. 1.2.11  were used to merge the overlapping paired-end Illumina fastq files. UPARSE algorithm (version 8.0.1623)  was used for (a) quality filtering of the reads based on maximum expected error value = 1.0, (b) identification of unique sequences, (c) abundance sorting and removal of singletons, (d) clustering the reads into operational taxonomic units (OTUs) based on 97% identity threshold, (e) de novo and reference-based chimera checking (against GOLD database  for V1-V2 sequences and UNITE database for ITS2 sequences ), and (f) construction of OTU table. For V1-V2 sequences, taxonomies were assigned using the QIIME  implementation of the UCLUST consensus taxonomy assigner (version = 1.2.22)  against the Greengenes database (release May 2013). Due to the importance of staphylococci and streptococci species to udder health and pathogenesis of mastitis, sequences of representative OTUs within these two genera were further assigned taxonomy at the species level based on top BLASTN bit-scores with minimum sequence identities of 98% to representative sequences deposited in the NCBI-16S rRNA gene database. For ITS2 sequences, taxonomies were assigned using the SortMeRNA classifier (version 2.0, 29/11/2014)  against the UNITE database (dynamic version 7.0) .
Prior to performing downstream analyses, the resulting bacterial and fungal OTU tables were filtered to remove all the samples with low sequencing depths (< 5000 sequences per sample). Community richness (Chao1) and diversity (Shannon) indices were then calculated using QIIME default scripts at even depths of 6000 and 5000 sequences per sample for bacterial and fungal OTU tables, respectively (see Additional file 1: Figure S8: rarefaction curves of richness and diversity indices of bacterial and fungal communities). Bray-Curtis dissimilarities of microbial communities were calculated following normalization of the final OTU tables using the cumulative sum scaling (CSS) transformation . Jaccard binary distances of microbial communities were calculated using rarefied OTU tables at even depths of 6000 and 5000 sequences per sample for bacterial and fungal OTU tables, respectively. Normalization of the final OTU tables using principal coordinate analysis (PCoA) was applied on the resulting distance matrices to generate two-dimensional plots using default settings of the PRIMER-E software ver. 6.1.18 .
Unsupervised clustering analysis was performed to relate clustering patterns of samples to the proportion of the most abundant OTUs. The relative abundances of selected OTUs (> 0.1% of the community) were normalized across samples (values divided by the Euclidean length of the row vector). Bray-Curtis dissimilarities were calculated using R “vegan” package  and the resulting matrix was subjected to unsupervised agglomerative hierarchical clustering method (unweighted pair group method with arithmetic mean (UPGMA)) using R “dendextend” package  and visualized over a heatmap of the abundance matrix using R “complexheatmap” package .
The UNIVARIATE procedure of SAS (SAS 9.3, 2012) was used for testing the normality of residuals for α-diversity measurements. Non-normally distributed data were Box-Cox power transformed and then fitted in a repeated measurement model using MIXED procedure of SAS with DIM, BstYI variants, and the interaction between DIM and BstYI variants being treated as fixed factors and the effect of individual cows treated as the random effect nested within the BstYI variants. Permutational multivariate analysis of variance (PERMANOVA ; implemented in Primer6 software) was used to test for significant differences in the β-diversity of microbial communities belonging to different DIM (n = 3) and BstYI variants (n = 2). Label permutations (n = 9999) were used in PERMANOVA to estimate the distribution of test statistics under the null hypothesis that within-group Bray-Curtis and binary Jaccard distances were not significantly different from between-group ones. A repeated measurement PERMANOVA model was defined in which DIM, BstYI variants, and the interaction between DIM and BstYI variants were treated as fixed factors while the effect of individual cows was treated as a random factor nested within the BstYI variants. Associations of bacterial genera with BoLA variants and DIM were analyzed with a generalized linear mixed effect model (repeated measurement) using package glmmTMB . The total count of OTUs assigned to each genus were offset to the library size (total count of OTUs detected in each sample) and then used as the response variable in a negative binomial model where BoLA variants, DIM, and their interaction were included as fixed effects, whereas the effect of individual cows were included as the random effect. Estimated group means, confidence intervals, and pairwise comparisons for effects of BstY variants and DIM were derived using the package emmeans  and visualized using the package ggplot2 . Multiple hypotheses were adjusted by Benjamini and Hochberg false discovery rate (FDR ). Statistical comparisons between the relative abundances of bacterial/fungal OTUs and BstYI variants at selected time points were performed using linear discriminant analysis effective size (LEfSe; Segata et al. ). In this approach, the nonparametric factorial Kruskal-Wallis sum rank test was used to test whether the relative abundances of OTUs were differentially distributed between different classes (variants) of BstYI. This was followed by linear discriminant analysis (LDA) to estimate the effect size of each differentially abundant OTU. The threshold on the logarithmic LDA score was set at 2.0, so that only OTUs with at least 2.0 log fold change were considered to be differentially distributed between classes.
Finally, correlation network analysis (CoNet, ) was used to explore microbial co-occurrence/mutual exclusion relationships and identify hub OTUs that show the highest number of positive/negative correlations with other OTUs. In this ensemble method, combinations of different correlation (Pearson’s and Spearman’s rank correlation coefficient) and dissimilarity (Bray-Curtis and Kullback-Leibler) measures were used to infer co-occurrence networks. In brief, for each measure, distributions of all pairwise scores between the nodes (a node representing the relative abundance of a non-singleton OTU that was found in at least 20% of the samples) were computed. For each measure and edge (an edge representing a positive or negative correlation between two nodes), 1000 permutation was conducted (this included a renormalization step for Pearson and Spearman measures in order to address the issue of compositionality introduced by different sequencing depths for each sample). For all correlation and dissimilarity measures, p values were computed as the probability of the null value (represented by the mean of the null distribution) under a Gauss curve generated from the mean and standard deviation of the bootstrap distribution. Measure-specific p values were then merged using Brown’s method  and subjected to Benjamini-Hochberg’s FDR correction. An edge was considered significant and kept in the final network if (a) it was supported by at least three measures, (b) it had a merged p value below 0.05, and (c) it was within the 95% confidence interval defined by the bootstrap distribution.
Cho I, Blaser MJ. The human microbiome: at the interface of health and disease. Nat Rev Genet. 2012;13:260–70.
Ley RE, Lozupone CA, Hamady M, Knight R, Gordon JI. Worlds within worlds: evolution of the vertebrate gut microbiota. Nature Rev Microbiol. 2008;6:776–88.
Nicholson JK, Holmes E, Kinross J, Burcelin R, Gibson G, Jia W, Pettersson S. Host-gut microbiota metabolic interactions. Science. 2012;336:1262–7.
Hooper LV, Littman DR, Macpherson AJ. Interactions between the microbiota and the immune system. Science. 2012;336:1268–73.
Round JL, Mazmanian SK. The gut microbiota shapes intestinal immune responses during health and disease. Nat Rev Immunol. 2009;9:313–23.
Lee YK, Mazmanian SK. Has the microbiota played a critical role in the evolution of the adaptive immune system? Science. 2010;330:1768–73.
Sommer S. The importance of immune gene variability (MHC) in evolutionary ecology and conservation. Front Zool. 2005;2:16.
Roche PA, Furuta K. The ins and outs of MHC class II-mediated antigen processing and presentation. Nat Rev Immunol. 2015;15:203.
Klein J. Natural history of the major histocompatibility complex. New York: Wiley; 1986.
Teh HS, Kisielow P, Scott B, Kishi H, Uematsu Y, Blüthmann H, von Boehmer H. Thymic major histocompatibility complex antigens and the αβ T-cell receptor determine the CD4/CD8 phenotype of T cells. Nature. 1988;335:229–33.
De Palma G, Capilla A, Nadal I, Nova E, Pozo T, Varea V, Polanco I, Castillejo G, López A, Garrote JA. Interplay between human leukocyte antigen genes and the microbial colonization process of the newborn intestine. Curr Issues Mol Biol. 2010;12:1.
Toivanen P, Vaahtovuo J, Eerola E. Influence of major histocompatibility complex on bacterial composition of fecal flora. Infect Immun. 2001;69:2372–7.
Kubinak JL, Stephens WZ, Soto R, Petersen C, Chiaro T, Gogokhia L, Bell R, Ajami NJ, Petrosino JF, Morrison L. MHC variation sculpts individualized microbial communities that control susceptibility to enteric infection. Nat Commun. 2015;6:8642.
Bolnick DI, Snowberg LK, Caporaso JG, Lauber C, Knight R, Stutz WE. Major histocompatibility complex class IIb polymorphism influences gut microbiota composition and diversity. Mol Ecol. 2014;23:4831–45.
Behl JD, Verma N, Tyagi N, Mishra P, Behl R, Joshi B. The major histocompatibility complex in bovines: a review. ISRN Vet Sci. 2012;872710:1–12.
Rupp R, Boichard D. Genetics of resistance to mastitis in dairy cattle. Vet Res. 2003;34:671–88.
Takeshima SN, Aida Y. Structure, function and disease susceptibility of the bovine major histocompatibility complex. Animal Sci J. 2006;77:138–50.
Bradley AJ. Bovine mastitis: an evolving disease. Vet J. 2002;164:116–28.
USDA. Dairy 2014, milk quality, milking procedures, and mastitis in the United States. Fort Collins: USDA–APHIS–VS–CEAH–NAHMS; 2016.
Sharif S, Mallard B, Wilkie B, Sargeant J, Scott H, Dekkers J, Leslie K. Associations of the bovine major histocompatibility complex DRB3 (BoLA-DRB3) alleles with occurrence of disease and milk somatic cell score in Canadian dairy cattle. Anim Genet. 1998;29:185–93.
Rupp R, Hernandez A, Mallard B. Association of bovine leukocyte antigen (BoLA) DRB3. 2 with immune response, mastitis, and production and type traits in Canadian Holsteins. J Dairy Sci. 2007;90:1029–38.
Dietz AB, Cohen ND, Timms L, Kehrli ME. Bovine lymphocyte antigen class II alleles as risk factors for high somatic cell counts in milk of lactating dairy cows. J Dairy Sci. 1997;80:406–12.
Kelm SC, Dettilleux JC, Freeman AE, Kehrli ME, Dietz AB, Fox LK, Butler JE, Kasckovics I, Kelley DH. Genetic association between parameters of innate immunity and measures of mastitis in periparturient Holstein cattle. J Dairy Sci. 1997;80:1767–75.
Yoshida T, Mukoyama H, Furuta H, Kondo Y, Takeshima SN, Aida Y, Kosugiyama M, Tomogane H. Association of BoLA-DRB3 alleles identified by a sequence-based typing method with mastitis pathogens in Japanese Holstein cows. Animal Sci J. 2009;80:498–509.
Melchior MB, Vaarkamp H, Fink-Gremmels J. Biofilms: a role in recurrent mastitis infections? Vet J. 2006;171:398–407.
Barkema HW, Schukken YH, Zadoks RN. Invited review: the role of cow, pathogen, and treatment regimen in the therapeutic success of bovine Staphylococcus aureus mastitis. J Dairy Sci. 2006;89:1877–95.
Kuehn JS, Gorden PJ, Munro D, Rong R, Dong Q, Plummer PJ, Wang C, Phillips GJ. Bacterial community profiling of milk samples as a means to understand culture-negative bovine clinical mastitis. PLoS One. 2013;8:e61959.
Oikonomou G, Bicalho ML, Meira E, Rossi RE, Foditsch C, Machado VS, Teixeira AGV, Santisteban C, Schukken YH, Bicalho RC. Microbiota of cow’s milk; distinguishing healthy, sub-clinically and clinically diseased quarters. PLoS One. 2014;9:e85904.
Ganda EK, Gaeta N, Sipka A, Pomeroy B, Oikonomou G, Schukken YH, Bicalho RC. Normal milk microbiome is reestablished following experimental infection with Escherichia coli independent of intramammary antibiotic treatment with a third-generation cephalosporin in bovines. Microbiome. 2017;5:74.
Lima SF, Teixeira AG, Lima FS, Ganda EK, Higgins CH, Oikonomou G, Bicalho RC. The bovine colostrum microbiome and its association with clinical mastitis. J Dairy Sci. 2017;100:3031–42.
Derakhshani H, Plaizier JC, De Buck J, Barkema HW, Khafipour E. Composition of teat canal and intramammary microbiota of dairy cows subjected to antimicrobial dry cow therapy and internal teat sealant. J Dairy Sci. 2018;101:10191–10205. In Press
Derakhshani H, Fehr KB, Sepehri S, Francoz D, De Buck J, Barkema HW, Plaizier JC, Khafipour E. Invited review: microbiota of the bovine udder: contributing factors and potential implications for udder health and mastitis susceptibility. J Dairy Sci. 2018;101:1–21. https://doi.org/10.3168/jds.2018-14860.
Rainard P. Mammary microbiota of dairy ruminants: fact or fiction? Vet Res. 2017;48:25.
Doyle CJ, Gleeson D, O'Toole PW, Cotter PD. Impacts of seasonal housing and teat preparation on raw milk microbiota: a high-throughput sequencing study. Appl Environ Microbiol. 2017;83:e02694–16.
Bonsaglia EC, Gomes MS, Canisso IF, Zhou Z, Lima SF, Rall VL, Oikonomou G, Bicalho RC, Lima FS. Milk microbiome and bacterial load following dry cow therapy without antibiotics in dairy cows with healthy mammary gland. Sci Rep. 2017;7:8067.
Ganda EK, Bisinotto RS, Lima SF, Kronauer K, Decter DH, Oikonomou G, Schukken YH, Bicalho RC. Longitudinal metagenomic profiling of bovine milk to assess the impact of intramammary treatment using a third-generation cephalosporin. Sci Rep. 2016;6:37565.
Peterson DA, McNulty NP, Guruge JL, Gordon JI. IgA response to symbiotic bacteria as a mediator of gut homeostasis. Cell Host Microbe. 2007;2:328–39.
Brandtzaeg P. The mucosal immune system and its integration with the mammary glands. J Pediatr. 2010;156:S8–S15.
Gapper LW, Copestake DE, Otter DE, Indyk HE. Analysis of bovine immunoglobulin G in milk, colostrum and dietary supplements: a review. Anal Bioanal Chem. 2007;389:93–109.
Schwab I, Nimmerjahn F. Intravenous immunoglobulin therapy: how does IgG modulate the immune system? Nat Rev Immunol. 2013;13:176.
Amigorena S, Bonnerot C. Fc receptors for IgG and antigen presentation on MHC class I and class II molecules. Semin Immunol. 1999;11:385–90.
Fox P. Milk proteins: general and historical aspects. In: Advanced dairy chemistry—proteins: Springer; 2003. p. 1–48.
Barrington GM, McFadden TB, Huyler MT, Besser TE. Regulation of colostrogenesis in cattle. Livest Prod Sci. 2001;70:95–104.
Levieux D, Ollier A. Bovine immunoglobulin G, β-lactoglobulin, α-lactalbumin and serum albumin in colostrum and milk during the early post partum period. J Dairy Res. 1999;66:421–30.
Espeche MC, Otero MC, Sesma F, Nader-Macias MEF. Screening of surface properties and antagonistic substances production by lactic acid bacteria isolated from the mammary gland of healthy and mastitic cows. Vet Microbiol. 2009;135:346–57.
Espeche MC, Pellegrino M, Frola I, Larriestra A, Bogni C, Nader-Macías MF. Lactic acid bacteria from raw milk as potentially beneficial strains to prevent bovine mastitis. Anaerobe. 2012;18:103–9.
Bouchard DS, Seridan B, Saraoui T, Rault L, Germon P, Gonzalez-Moreno C, Nader-Macias FM, Baud D, François P, Chuat V. Lactic acid bacteria isolated from bovine mammary microbiota: potential allies against bovine mastitis. PLoS One. 2015;10:e0144831.
Hogan J, Smith K, Todhunter D, Schoenberger P. Rate of environmental mastitis in quarters infected with Corynebacterium bovis and Staphylococcus species. J Dairy Sci. 1988;71:2520–5.
Carson DA, Barkema HW, Naushad S, De Buck J. Bacteriocins of non-aureus staphylococci isolated from bovine milk. Appl Environ Microbiol. 2017;83:e01015–7.
Gerhold P, Cahill JF, Winter M, Bartish IV, Prinzing A. Phylogenetic patterns are not proxies of community assembly mechanisms (they are far better). Funct Ecol. 2015;29:600–14.
Costello EK, Stagaman K, Dethlefsen L, Bohannan BJ, Relman DA. The application of ecological theory toward an understanding of the human microbiome. Science. 2012;336:1255–62.
Gopal PK, Gill H. Oligosaccharides and glycoconjugates in bovine milk and colostrum. Br J Nutr. 2000;84:69–74.
Aakko J, Kumar H, Rautava S, Wise A, Autran C, Bode L, Isolauri E, Salminen S. Human milk oligosaccharide categories define the microbiota composition in human colostrum. Benef Microbes. 2017;8:563–7.
Addis MF, Tanca A, Uzzau S, Oikonomou G, Bicalho RC, Moroni P. The bovine milk microbiota: insights and perspectives from-omics studies. Mol BioSyst. 2016;12:2359–72.
Ruiz-Palacios GM, Cervantes LE, Ramos P, Chavez-Munguia B, Newburg DS. Campylobacter jejuni binds intestinal H (O) antigen (Fuca1, 2Galb1, 4GlcNAc), and fucosyloligosaccharides of human milk ihibit its binding and infection. J Biol Chem. 2003;278:14112–20.
Fehr KB, Derakhshani H, Sepehri S, Plaizier JC, Khafipour E. Effects of dairy environment on milk microbiota and mammary inflammation. J Dairy Sci. 2017;100(Suppl. 2):142.
Condas LAZ, De Buck J, Nobrega DB, Carson DA, Roy J-P, Keefe GP, DeVries TJ, Middleton JR, Dufour S, Barkema HW. Distribution of non-aureus staphylococci species in udder quarters with low and high somatic cell count, and clinical mastitis. J Dairy Sci. 2017;100:5613–27.
Derakhshani H, Plaizier JC, Khafipour E. Bacterial ecosystem of the bovine mammary gland: potential role of foundation taxa in shaping mammary gland microbiota and modulating udder homeostasis. J Dairy Sci. 2017;100(Suppl. 2):232.
Jiang Q, Yang Y, Xin K, Bian X, Li M, Zhang B, Li C, Yao Z, Hu J, Sun D. Microbial diversity analysis of subclinical mastitis in dairy cattle in Northeast China. Afr J Microbiol Res. 2015;9:687–94.
Petrovski KR, Williamson NB, Lopez-Villalobos N, Parkinson TJ, Tucker IG. Culture results from milk samples submitted to veterinary diagnostic laboratories from august 2003 to December 2006 in New Zealand. N Z Vet J. 2011;59:317–22.
Wawron W, Bochniarz M, Piech T. Yeast mastitis in dairy cows in the middle-eastern part of Poland. Bull Vet Inst Pulawy. 2010;54:201–4.
Costa E, Gandra C, Pires M, Coutinho S, Castilho W, Teixeira C. Survey of bovine mycotic mastitis in dairy herds in the state of São Paulo, Brazil. Mycopathologia. 1993;124:13–7.
AalbæK B, Stenderup J, Jensen HE, Valbak J, Nylin B, Huda A. Mycotic and algal bovine mastitis in Denmark. APMIS. 1994;102:451–6.
Pachauri S, Varshney P, Dash SK, Gupta MK. Involvement of fungal species in bovine mastitis in and around Mathura, India. Vet World. 2013;6:393–5.
Farnsworth R, Sorensen D. Prevalence and species distribution of yeast in mammary glands of dairy cows in Minnesota. Can J Comp Med. 1972;36:329.
Bouza E, Muñoz P. Invasive infections caused by Blastoschizomyces capitatus and Scedosporium spp. Clin Microbiol Infect. 2004;10:76–85.
Maslen M, Peel M. Human and animal isolates of Pseudallescheria boydii and Scedosporium species, from Melbourne, Australia, 1977–1995. Mycoses. 2011;54:442–9.
Eijk M, Stewart-Haynes J, Lewin H. Extensive polymorphism of the BOLA-DRB3 gene distinguished by PCR-RFLP. Anim Genet. 1992;23:483–96.
Derakhshani H, Tun HM, Khafipour E. An extended single-index multiplexed 16S rRNA sequencing for microbial community analysis on MiSeq illumina platforms. J Basic Microbiol. 2016;56:1–6.
Magoc T, Salzberg SL. FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics. 2011;27:2957–63.
Edgar RC. UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat Methods. 2013;10:996–8.
Nilsson RH, Tedersoo L, Ryberg M, Kristiansson E, Hartmann M, Unterseher M, Porter TM, Bengtsson-Palme J, Walker DM, De Sousa F. A comprehensive, automatically updated fungal ITS sequence dataset for reference-based chimera control in environmental sequencing efforts. Microbes Environ. 2015;30:145–50.
Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, Fierer N, Pena AG, Goodrich JK, Gordon JI. QIIME allows analysis of high-throughput community sequencing data. Nat Methods. 2010;7:335–6.
Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010;26:2460–1.
Kopylova E, Noé L, Touzet H. SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data. Bioinformatics. 2012;28:3211–7.
Paulson JN, Stine OC, Bravo HC, Pop M. Differential abundance analysis for microbial marker-gene surveys. Nat Methods. 2013;10:1200–2.
Warwick R, Clarke K. PRIMER 6. Plymouth: PRIMER-E Ltd; 2006.
Oksanen J, Kindt R, Legendre P, O’Hara B, Stevens MHH, Oksanen MJ, Suggests M. The vegan package. Community ecology package; 2007. p. 10.
Galili T. Dendextend: an R package for visualizing, adjusting and comparing trees of hierarchical clustering. Bioinformatics. 2015;31:3718–20.
Zuguang G. ComplexHeatmap: making complex Heatmaps. 2015.
Anderson MJ. A new method for nonparametric multivariate analysis of variance. Austral Ecol. 2001;26:32–46.
Magnusson A, Skaug H, Nielsen A, Berg C, Kristensen K, Maechler M, van Bentham K, Bolker B, Brooks M, Brooks MM: Package ‘glmmTMB’. 2017.
Lenth R, Love J, Hervé M. Package ‘emmeans’. Underst Stat. 2017;34:216–21.
Wickham H. ggplot2: elegant graphics for data analysis. J Stat Softw. 2010;35:65–88.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B Stat Methodol. 1995;1:289–300.
Segata N, Izard J, Waldron L, Gevers D, Miropolsky L, Garrett WS, Huttenhower C. Metagenomic biomarker discovery and explanation. Genome Biol. 2011;12:R60.
Faust K, Raes J. CoNet app: inference of biological association networks using Cytoscape [version 2; referees: 2 approved]. F1000Res. 2016;5:1519.
Brown MB. 400: a method for combining non-independent, one-sided tests of significance. Biometrics. 1975:987–92.
The authors thank the staff of the Rosser Holsteins Ltd. for their assistance with conduct of study and sample collection. Special thanks to Jake Cowper Szamosi (McMaster University) for her assistance with bioinformatics and statistical modelling.
The present research was supported by grants from Growing Forward 2—Agricultural Rural Development Initiative Program of the Province of Manitoba, Canada; and Dairy Farmers of Manitoba, Canada.
Availability of data and materials
Paired-end Illumina fastq files were merged, demultiplexed, and deposited into the Sequence Read Archive (SRA) of NCBI (http://www.ncbi.nlm.nih.gov/sra) and can be accessed via accession numbers SRR7262561-SRR7262777. Metadata used to perform bioinformatics and statistical analyses for bacterial and fungal communities are presented in Additional file 2: Tables S8 and S9, respectively.
The sampling protocol of the present study was reviewed and approved by the Animal Care Committee of the University of Manitoba, Winnipeg, Canada (protocol number F14-027).
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. PCR-RFLP Identification of BoLA-DRB3.2 variants. Figure S2. Beta-diversity of the bacterial communities of intramammary secretions during the first week of lactation. Figure S3. Beta-diversity of the fungal communities of intramammary secretions during the first week of lactation. Figure S4. Association of lactation stage and BoLA-DRB3.2 polymorphism with interrelationship patterns and hub species of intramammary bacterial communities. Figure S5. Bacterial co-occurrence and co-exclusion networks. Figure S6. Unsupervised cluster analysis of day 0 colostrum samples based on the distribution of fungal OTUs. Figure S7. Unsupervised cluster analysis of day 1 colostrum samples based on the distribution of fungal OTUs. Figure S8. Rarefaction curves of richness (Observed-OTUs and Chao1) and diversity (Shannon) indices for bacterial (a-c) and fungal (d-f) communities. (PDF 6512 kb)
Table S1. Summary statistics for comparison of alpha-diversity of microbial communities. Table S2. Summary statistics for comparison of beta-diversity of microbial communities (PERMANOVA test statistics). Table S3. Summary statistics for repeated measurement analysis of bacterial communities. Table S4. Summary statistics for repeated measurement analysis of fungal communities. Table S5. NCBI blastn search result for selected bacterial OTUs. Table S6. Primer design and equences for barcoded PCR amplification of the V1-V2 regions of bacterial 16S rRNA genes. Table S7. Primer design and sequences for barcoded PCR amplification of the ITS2 of the eukaryotic ribosomal genes. Table S8. Metadata used for processing of 16S rRNA gene raw sequences, bioinformatics and statistical analyses of bacterial communities. Table S9. Metadata used for processing of ITS2 raw sequences, bioinformatics and statistical analyses of fungal communities. (XLSX 48 kb)