Milk microbiome: evaluation study on the differences among cows with a different health status classified by leukocyte pattern

Milk is considered not only a source of nutrient for the offspring but also a font of immunoregulatory compounds capable to predispose the naïve intestinal immune system of the new-born to react to the external environment. In the present study we evaluated the composition of milk microbiome from cows classified according to milk total and differential somatic cell counts. A total of 34, 13 and 13 milk samples of healthy, at risk and subclinical or chronic cows, respectively, were collected during the same milking from a local dairy herd. Through Next Generation Sequencing (NGS) of bacterial 16S rRNA gene, the differences of taxa in terms of relative abundances (RA) and alpha and beta biodiversity were analysed. The RA of several genera were statistically significant in the three groups, such as Arcanobacterium (p=0.001), Rhodococccus (p<0.05) and Rubrobacter (p<0.05), while at species level the presence of Propionibacterium granulosum, Pseudomonas alcaligenes and Prosthecobacter debontii were found. Shannon and Evenness indices computed at the genus level were not significant, while beta biodiversity showed a clear clusterization between groups. The results highlighted that milk microbiome is associated to a different cellular response at udder level, although more specific studies are needed to assess the source of bacteria species identified in milk microbial population of healthy animals.


Introduction
Mastitis is a multifactorial disease with a complex etiology still difficult to understand nowadays. In the past, milk secreted by acino alveolar glands from healthy quarter was thought to be sterile (Tolle, 1980), and the presence of certain bacteria in the udder, and consequently in its secretions, was attributed to an entrance from the teat canal, provoking a colonization of the mammary gland (Zecconi at al., 2000). Recently, the advent of high-throughput sequencing technique has revealed that also wholesome milk contains lower abundances of some bacterial taxa that generally are undetected with culture-based methods (Oikonomou et al., 2014;Lima et al., 2018). In addition, different studies report that milk contains not only mastitis pathogens but also opportunistic and/or commensal strains (Young et al., 2015;Taponen et al., 2019). From this evidence, it is justified to speculate that a gut-mammary axis exists, as the debated origin of these bacteria is probably not only environment enrichment but also internal. In humans and in mice, this endogenous pathway was reported (Rodriguez, 2014) and it was achieved through the migration of viable bacteria from gut to mammary gland via dendritic cells (DCs) and macrophages. These cells, found in gut-associated lymphoid tissue (GALT), are occasionally able to internalize and translocate resident intestinal microbiota to the mammary gland (Roux et al.,1977;Macpherson and Uhr, 2004). Moreover, the alteration of the mammary gland barrier integrity, that occurs in case of infections (e.g. mastitis), might promote translocation of bacteria through the mammary epithelium, even if the migration of gut bacteria through an intact mammary gland epithelium is still debated. Studies on cattle have demonstrated that the transfer of lymphoid cells from intestines to mammary gland is minor, because the majority of lymphocytes supplying local immunity originate from peripheral lymph nodes rather than from the mucosal lumen of the gut (Kehrli and Harp, 2001), advising that this axis is less functional in ruminants. Moreover, the finding of common pathogens causing environmental mastitis (i.e. S. uberis, E. coli) in healthy milk, both from in vivo and in vitro experimental studies (Ganda et al., 2017;Young et al., 2015), could be imputed to an altering of the internal microbial ecology rather than a primary infection (Oikonomou et al., 2014). Considering this, it is possible to assume that mammary gland infections are not simply related to a bacterial pathogen, but also to the consequence of a dysbiosis (Derakhshani et al., 2018). This dysbiotic status decreases biodiversity and alters the composition of the microbial community of the udder, predisposing the bovine to intramammary infections (Taponen et al., 2019). Besides, the progressive minor accuracy of the somatic cell count (SCC) as a mastitis marker due to the decrease of the mean SCC in dairy herd, leads to the need to use another more precise method, as the differential somatic cell count (DSCC), alone or in combination with the previous. DSCC analyses is performed with Fossomatic 7 DC (Foss, Hilleroed, DK), a machinery that allows to identify in milk samples the percentage of polymorphonuclear neutrophils and lymphocytes on the overall count of milk cells . Moreover, DSCC showed to be associated to changes in milk composition; therefore, it could supply information on mammary gland barrier integrity (Stocco et al., 2020;. The purpose of this study was to characterize and compare milk microbiome in terms of relative abundances (RAs) and biodiversity in cows classified according to the DSCC+SCC assessment as healthy, at risk and subclinical or chronic.

Animals and housing
Sixty Holstein Friesian cows from a local farm (46.1134059, 13.2817455;N 46°6'48.261'', E 13°16'54.283'';Italy) were recruited for this study. The cows were preselected based on SCC and DSCC values of the previous monthly official record of the Breeder Association (Associazione Allevatori del Friuli Venezia Giulia, Codroipo, Italy; www.aafvg.it). The cows of herd were divided in 4 groups according to the already existing classification Zecconi et al., 2020a) in healthy (G; <200,000 SCC/mL; DSCC ≤69.3%), at risk (Y; <200,000 SCC/mL; DSCC >69.3%), chronic (O; >200,000 SCC/mL; DSCC ≤69.3%) and with subclinical mastitis (R; >200,000 SCC/mL; DSCC >69.3%) subjects. According to the data of the official record of the day before sampling, the number of cows was 34 in G (11 primiparous and 23 pluriparous; 10 DIM < 70 and 24 DIM > 70), 13 in Y (6 primiparous and 7 pluriparous; 2 DIM < 70 and 11 DIM > 70), 4 in O (4 pluriparous; 1 DIM < 70 and 3 DIM > 70) and 9 in R (2 primiparous and 7 pluriparous; 2 DIM < 70 and 7 DIM > 70). The O group was merged with the R group since they were considered similar at milk microbiome profile having more than 200,000 SCC/mL. The animals were housed in free stalls with cubicles and fed with the same total mixed ration and without having received previous antibiotic treatments. The milking parlour (parallel 12+12) was adjacent to the bedding. All protocols, procedures and the care of the animals complied to the Italian legislation on animal care (DL n.116, 27/1/1992) and the study was approved by the ethical committee of the University of Udine.

Collection of samples
Milk samples were collected during the evening milking. Briefly, the first flow of milk from each quarter was discarded and the teats were dipped in iodine tincture. Thus, each quarter was cleaned with disposable wipes and disinfected (ethanol 70%) for the milking. A sample of 50 mL of milk (a mixture including all the four quarters) was withdrawn from each cow. The samples were kept on ice during the collecting period and once back to the laboratory they were stored at -20°C until analysis.

DNA Extraction, Library Preparation, Sequencing and Taxonomic Annotation
Microbial DNA from milk samples was extracted from 250 µl of starting material using a Faecal DNA Miniprep kit (Zymo Research; Irvine, CA, US), following the manufacturer's instructions, with the addition of a preliminary warming step of 10' at 70°C (Lima et al., 2018) followed by a bead lysis process. DNA concentration was measured with a QubitTM 3 Fluorometer (Thermo Scientific; Waltham, MA, USA) and the 16S rRNA of V3 and V4 regions amplified for library preparation, adding also the Indexes for sequencing, using a Nextera DNA Library Prep kit (Illumina; San Diego, CA, USA), following manufacturer's instruction and primers (Klindworth et al., 2013). The resulting amplicons were sequenced with a NovaSeq6000 (Illumina; San Diego, CA, USA) in 2x300 paired-end mode, following the standard procedures. The Quantitative Insights Into Microbial Ecology (QIIME 2) (Bolyen et al., 2019) was used to process the raw sequences. After demultiplexing, sequenced reads that passed the quality check (Phred score ≥30) were annotated for 16S rRNA against the most recent Greengenes database (version gg.13_8.otus.tar.gz), with 99% identify with reference sequences.
Chimeras were also detected and then filtered from the reads and the remaining sequences were clustered into Operational Taxonomic Units (OTUs) by using an open reference approach in QIIME 2. This procedure is the preferred strategy for OTU picking in QIIME2, which includes taxonomy assignment, sequence alignment, and tree-building steps.

Statistical analysis
Data were imputed on a spreadsheet for analysis, together with the healthy status of the subjects previously explained with the classification based on the SCC and DSCC to allow and facilitate further statistical analysis. The sequences annotated from each sample were normalized to ‰ abundance profiles, known as Relative Abundance (RA), at each taxonomic level. Taxa with RAs lower than 1‰ in more than 50 samples were excluded from the statistical analysis. Alpha biodiversity (the individual richness in terms of microbial diversity observed) was evaluated with Shannon and Evenness indices and the observed RA through a Kruskal-Wallis non-parametric test. Beta biodiversity (difference in species composition between different groups, sites or communities) was evaluated with the phylogeny based on unweighted UniFrac (Lozupone and Kight, 2005) distance metric and visualized using Principal Coordinate Analysis (PCoA) plot. Association between microbiome composition and covariate (SCC and DSCC classification;  was tested using PERMANOVA, a nonparametric test similar to ANOVA but without the requirement of a normally distributed data. Significance of PERMANOVA test was determined using 999 permutations with adjustment for multiple testing. Statistical analyses were performed with XLSTAT (Addinsoft, 2020).
Rhodococcus and Williamsia genera, belonging to the phylum of Actinobacteria, are considered opportunistic bacteria (Keikha, 2017) and interestingly, we found a higher abundance in G group than the other groups. Another contradictory evidence is that we did not find significant differences among taxa commonly known to be responsible for mastitis (Staphylococcus spp., Streptococcus spp., Escherichia spp.).
Some species as Pseudomonas alcaligenes and Prosthecobacter debontii were higher in G group than in Y, and interestingly were absent in R. This finding led to hypothesize that these species, and in a broad sense the belonging genera, might represent the local physiological microbiome of the mammary gland in lactating cows. This concept was supported by the study of Oikonomou et al. (2014) in which they found an iteration of Propionibacterium spp. in every healthy quarter samples and particularly, at species level they found a higher RA of P. acnes. Another explanation for this could be attributed to contamination issues during the sampling, as sustained by Metzger et al. (2018), who have noticed that Pseudomonas spp., known to be an environmental microorganism, was one of the most abundant genera in their negative control samples. It is already known that the sampling techniques is a critical point for obtaining uncontaminated samples from milk, even if sampling in dairy farms is always challenging manly owing to the elevate sources of contamination (Vangroenweghe et al. 2001). The relative poor abundance or absence of Staphylococcus spp., Streptococcus spp. and Escherichia spp. suggest that the level of contamination in our samples should be very low, since these are the most frequently isolated bacteria in milk, and their main reservoirs are environment or teat-skin (Zecconi and Piccinini, 2002). It should be considered that a dilution effect of the mixture of four quarters could have reduced the possibility to detect these pathogens. However, yet another hypothesis was formulated, that is when a dysbiosis occurs, the enrichment of opportunistic and/or environmental bacterial strains in milk could be detected, probably caused by a translocation to the mammary site through blood or lymphatic streams, arousing inflammation process and ultimately, mastitis (Young et al., 2015). Notably, it has to be taken into consideration that here we analysed milk from a mixture of all four quarters, representing the four different ecological niches, and thus the comparison with studies made on quarters could be a little forced. In our case microbial profiles that distinguish healthy from cows with an intramammary infection (IMI) are not univocally consistent with the overall existing studies and a common point could be searched in the association of IMI with dysbiotic microbiota, basing the analyses on the estimation of differences in microbial composition and/or richness (alpha and beta biodiversity). The study of Kuehn et al. (2013), in which microbiota of milk from clinically healthy quarters and those from culture-negative quarters with clinical mastitis were compared, revealed that the latter group had higher proportions of Burkholderia, Sphingomonas and Stenotrophomonas genera, while Pseudomonas, Psychrobacter and Ralstonia genera increased in the microbiota of healthy quarters. In a discriminating study on milk microbiome, Oikonomou et al. (2012) registered the lowest Shannon and Chao 1 index in milk from clinical mastitis cases in comparison to healthy milk. Conversely here, alpha biodiversity indices (Shannon and Evenness) calculated at genus level did not show significant differences between the groups (data not shown), similarly to results of Lima et al. (2018) for Shannon Index in milk fractions from healthy and mastitic cows. Beta biodiversity index calculated on the unweighted UniFrac distance matrices showed a significant differentiation between the three experimental groups (G, Y and R), as shown in PCoA graph (Figure 1). PERMANOVA analysis of the latter data revealed significant differences between G and R (p=0.003) and between Y and R (p=0.002), detecting that R group presented lowest richness in comparison to the other experimental group.

Conclusions
The preliminary results presented in this paper suggest that milk microbiome is associated to a different cellular response at udder level. Indeed, apparently healthy cows (SCC <200,000 cells/mL), but with different proportions of polymorphonuclear neutrophils and lymphocytes (G and Y groups) have a significant different microbial population. Moreover, these latter group of cows have also a significant different microbial population than cows with subclinical or chronic mastitis (O+R group). To confirm the potential immunomodulating activity of milk microbiota or its capability to modulate pathogens growth, more targeted studies are deserved, as well as to define the origin of bacteria species identified in milk microbiota of healthy animals. Considering the potential of the relationship between SCC and microbial population in the milk, further investigation could be of high importance for the characterization of technological quality of milk. If this will be validated, SCC and DSCC can be considered further classification criteria for the destination of milk to dairy industry.

Figure 1
Beta diversity analysis represented by PCoA graph of unweighted UniFrac distance between healthy (G; green dots), subclinical (Y; yellow dots) and mastitic (R; red dots) bovines. PERMANOVA confirmed the differences between the G, Y and R groups of cows at P < 0.05.