Nisin A-producing Lactococcus cremoris formulations for pre- and post-milking teat disinfection modulate the bovine milk microbiota

Abstract

Background

Bovine mastitis is a major challenge in dairy farms. Since the agents commonly used for pre- and post-dipping can affect the udder health by modifying milk microbiota, alternative products are needed. This study aimed to evaluate the effect of the use of pre- and post-dipping formulations containing the fermented broth of Nisin A-producing Lactococcus cremoris FT27 strain (treated group, TR) on the abundance and biodiversity of milk microbiota as compared to iodine-based commercial disinfectants (control group, CTR) during a three-month trial. The experiment was conducted on 20 dairy cows, divided into two groups (CTR and TR) of 10 lactating cows each. Milk samples were collected from two selected healthy quarters of each cow at 3 time-points. Microbial communities were investigated by cultural and sequence-based methods, and analyzed through bioinformatic and statistical approaches.

Results

Clear differences in bacterial community composition were observed among groups, with higher species richness in TR, especially of Staphylococcus, Enterococcus, Lactococcus, and Streptococcus genera. The microbiota was dominated by Firmicutes, followed by Actinobacteriota, Proteobacteria, and Bacteroidota. Staphylococcaceae family was significantly higher in TR (p < 0.009), whereas Carnobacteriaceae, Mycobacteriaceae, and Pseudomonadaceae were significantly lower (p = 0.005, p = 0.001, and p = 0.040, respectively). CTR had considerably higher abundances of the genera Alkalibacterium (p = 0.011), Pseudomonas_E (p = 0.045), Corynebacterium (p = 0.004), and Alloiococcus (p = 0.004), and lower abundances of Staphylococcus (p < 0.009). Milk microbiota changed noticeably during the experimental period, regardless of treatment. A significant decrease was observed in both groups for Firmicutes_A phylum, with an increment in Actinobacteriota phylum, Propionibacteriaceae family, and Cutibacterium genus. Streptococcaceae significantly decreased in CTR (p = 0.013) and rose in TR (p = 0.001). Several differences were observed between the two groups during the experimental period. Streptococcus genus almost disappeared in CTR (p = 0.013), whereas it significantly increased in TR (p = 0.001). Three and twelve enriched groups were significantly identified respectively in CTR and TR using LEfSe.

Conclusions

The use of Nisin A-based teat dip formulations could be linked to greater microbial diversity compared to commercial products. Despite the influence of seasonality, the experimental formulations maintained higher milk biodiversity, suggesting that lactic acid bacteria metabolites prevent alterations in the milk microbiota.

Background

Bovine mastitis is one of the most important challenges in dairy farms, as it can result in significant economic losses due to reduced milk yield and quality and to the high costs of the consequent antimicrobial treatments [1].

The overuse of antibiotics in past years has promoted the emergence of multidrug-resistant pathogens, which pose a threat to public health [2]. Furthermore, antibiotics can have a negative impact on the intra-mammary ecosystem, frequently leading to reduced abundance and diversity of commensal microbial species and alterations in their metabolic activities. This condition known as dysbiosis can lead to the development of intra-mammary infections due to the loss of both the resistance to exogenous pathogens and the competition between endogenous microbial species [3, 4]. Thus, the maintenance of milk microbial diversity is important for the cow’s health. Pre- and post-milking teat disinfection are important mastitis management tools to reduce the incidence of new intra-mammary infections in dairy cows. However, these practices, which are routinely performed using commercial solutions containing iodine, hydrogen peroxide, or chlorine, may modify the composition of the milk microbial population [5]. Indeed, earlier research demonstrated that the microbial population and the balance between spoilage organisms and beneficial cheese-making microbes can be affected by various milking practices, including equipment use and udder preparation both before and after milking [6].

Since lactic acid bacteria (LAB) are key players in the production of raw milk cheeses such as Grana Padano and Parmigiano Reggiano, an imbalance in the milk microbiota or the presence of disinfectant or antimicrobial residues in milk can significantly affect the dairy industry [7]. In this context, alternative teat dip products are needed to prevent mastitis, reduce teat skin irritation, and maintain a balanced microbial population. Previous studies evaluated the effects of various milking products, hydrocolloidal water-in-oil emulsion [8] or a novel probiotic lactobacilli-based [9, 10] teat disinfectant, on the microbial communities on the teat skin of cows and in their milk.

Antimicrobial peptides, such as bacteriocins, are small molecules produced by bacteria for defense or niche competition since they are biologically active against other microorganisms [11]. These compounds are commonly produced by LAB, a group of Gram-positive bacteria which comprises several genera including Lactococcus [12]. Interestingly, during the milking process, these peptides have been proposed as disinfectants or teat dip products, either in addition to or instead of the current commercial teat dip solutions [9, 13].

The present study aimed at evaluating the effect of the use of pre- and post-dipping products containing the fermented broth of Nisin A-producing Lc. cremoris FT27 strain on the abundance and biodiversity of the milk microbiota as compared to commercial disinfectants.

Results

Strain typing

The Total Mesophilic aerobic Bacteria (TMB) values of the CTR and TR samples did not differ significantly, as shown in Table 1; the TMB load of all single quarter milk samples was around 2 log at all time-points examined. Additionally, there were no discernible variations in the LAB content between the experimental groups (Table 1).

Table 1 Microbial counts of raw milk samples collected from cows treated with a commercial pre-milking and a post-milking disinfectants (CTR) and natural formulations containing the fermented broth of Nisin A-producing Lc. cremoris FT27 strain (TR)

A total of 282 strains (135 from CTR samples and 147 from TR samples) were isolated from the milk samples at different time-points. Forty isolates (19 from CTR and 21 from TR samples) did not grow after plate isolation. Clear differences in bacterial community composition were observed among the samples (Fig. 1). The CTR samples were dominated by Staphylococcus epidermidis (56.9%) and S. haemolyticus (34.5%), followed by low percentages of Enterococcus hirae (3.4%), E. villorum (1.7%), Lc. cremoris (1.7%), E. pseudoavium (0.9%) and Lacticaseibacillus paracasei (0.9%). On the other hand, milk samples of the TR group showed a higher level of biodiversity (12 bacterial species). These samples mainly consisted of S. epidermidis (55.6%), S. haemolyticus (18.3%), and S. warneri (9.5%), followed by E. faecium (4.8%), E. hirae, L. paracasei, and Streptococcus uberis (2.4%), Lc. garviae (1.6%), and E. casseliflavus, E. pseudoavium, Lc. lactis, and St. mitis (0.8%). Interestingly, E. villorum was detected only in the CTR samples, while S. warneri, E. faecium, E. casseliflavus, E. hirae, Lc. garviae, St. uberis, and St. mitis were peculiar to the TR samples. Lc. cremoris used to produce the alternative disinfectants was not detected in the experimental samples. In fact, the Lc. cremoris strains isolated from CTR and TR samples had a different RAPD pattern (similarity level < 85%) (data not shown). These findings confirm that the FT27 strain was not alive in the pre- and post-dipping products.

Fig. 1

Distribution of the bacterial species isolates in CTR and TR raw milk samples

Sequencing results

Out of the 120 quarter milk samples collected during the trial (40 samples for each time-point), 3 samples were discarded following experimental and sequencing issues (i.e. <100 reads), resulting in a final database of 117 samples that produced a total of 65,862 sequence variants (average non-chimeric ASV per sample, after filtering and denoising: 10916.0 ± 5725.7).

Microbiota diversity analysis

Multiple biodiversity indices were used to assess the biodiversity of the milk samples for treatment and at various time-points. While the PD whole tree metric showed significantly different biodiversity among groups (p < 0.001), alpha-diversity analysis did not reveal any significant differences between the CTR and TR groups for any of the metrics used (Observed species, p = 0.159; Chao1, p = 0.132; Shannon, p = 0.172), as shown in Fig. 2A.

Fig. 2

Microbiota diversity. Alpha-diversity analysis for treatment (A) and time-points (B) with Chao1, Observed species, Shannon index, and Faith’s phylogenetic tree (“PD whole tree”) metrics. (C) Principal coordinates analysis showing the unweighted Unifrac beta-diversity over time; the first and third principal coordinates (PCoA1 and PCoA3) are reported

No significant differences were observed between CTR and TR when paired by sampling time, whereas bacterial biodiversity was observed to be significantly changed within the groups over the experimental period (Fig. 2B). The CTR group only reported a change in the biodiversity between T1 and T2 (PD whole tree, p = 0.017); the TR group, on the other hand, significantly differed from the baseline both at T1 (Chao1, p = 0.037; Observed species, p = 0.010; Shannon, p = 0.006) and at T2 (p = 0.026, Chao1; p = 0.003, Observed species; p = 0.002, Shannon) with an increasing trend.

The beta-diversity analysis reported that the structure of the CTR milk microbiota differed significantly from that of the TR group according to both unweighted and weighted (p = 0.031 for both comparisons) Unifrac distances (data not shown).

In particular, regardless of the treatment, the principal coordinates analysis (PCoA) revealed that the milk microbiota changed noticeably during the experimental period, revealing sample groups that were more comparable to one another (higher similarity between T0 and T1 for both CTR and TR, while CTR T2 and TR T2 are similar and separated from the others), as shown in Fig. 2C. Within the CTR group, significant differences were observed between T0-T2 (adjusted p-values = 0.006 and 0.011, respectively unweighted and weighted Unifrac) and T1-T2 (adj. p-values < 0.02 for both Unifrac distances). As for the TR group, fermented broth containing nisin-A changed the milk microbiota significantly compared to the baseline composition (T0-T1, adj. p-value < 0.02 for both distances; T0-T2, adj. p-values = 0.002 and 0.007, respectively unweighted and weighted Unifrac).

Milk microbiota taxonomy profiles

For both CTR and TR groups, most of the reads belonged to Firmicutes, typically dominants in cow milk core microbiota, accounting for almost half of the average microbial composition (relative abundance: 44.6% CTR vs. 41.8% TR), followed by Actinobacteriota (21.9% CTR vs. 17.4% TR), Proteobacteria (18.8% CTR vs. 19.1% TR), and Bacteroidota (5.8% CTR vs. 7.6% TR); among these, the only significant differences were observed between the Actinobacteriota (p = 0.017) and the Bacteroidota phyla (p = 0.027).

Out of the most relatively abundant families, the Staphylococcaceae family was significantly higher in the nisin-A treated group (1.4% CTR vs. 4.7% TR; p < 0.009); in contrast, we found that Carnobacteriaceae, Mycobacteriaceae and Pseudomonadaceae were significantly lower in the TR group (p = 0.005, p = 0.001, and p = 0.040, respectively). Although not significant, the same trend was observed for the Propionibacteriaceae and the Oscillospiraceae families; other noteworthy (p < 0.050) families with an overall relative abundance less than 1.5% were Bifidobacteriaceae, significantly lower in the TR group, while Xanthomonadaceae and Desulfovibrionaceae were significantly higher in the TR group. The microbiota composition at the phylum and family levels is reported in Table S1.

At the genus level (Fig. 3A), the CTR group had considerably higher abundances of Alkalibacterium (p = 0.011), Pseudomonas_E (p = 0.045), Corynebacterium (p = 0.004), and Alloiococcus (p = 0.004). In contrast, the CTR group had lower abundances of Staphylococcus (p < 0.009), with Jeotgalibaca showing a similar but non-significant trend.

Fig. 3

Genus-level analysis. (A) Taxonomy classification for CTR and TR groups. (B) LAB relative abundance, reported as mean and SD. (C) Mirror plot of genera relative abundances grouped according to the Gram stain. (D) Relative abundance differences of T1 and T2 vs. T0 within the CTR and (E) the TR groups; significant differences refer to the comparisons T0-T1 and T0-T2. For all plots in the panel, stars indicate a statistically significant (i.e. adjusted p < 0.05) diversity between groups or times

While the previously detailed benchmark investigation did not highlight differences in the lactic acid bacteria, the 16S sequencing analysis identified 5 bacterial genera accounting for the LAB group (Fig. 3B): Alkalibacterium, Alloiococcus, members of the Carnobacteriaceae family, Lactobacillus, and Streptococcus; the first 3 were statistically different between the experimental groups (respectively: p = 0.011, p = 0.004, p = 0.048).

Milk microbiota taxonomy profiles (time-points)

The relative abundance of the bacterial distribution was further investigated.

At the trial baseline (T0), the two groups were highly comparable, net of inter-sample variations, as we observed p-values = 1 for 45 out of the first 50 relatively abundant genera.

Throughout the time-points, considering the average microbial composition at the different phylogenetic levels, an overall tendency towards an increment of the Gram-negative bacteria was observed within the TR group, while the CTR group reported higher Gram-positive relative abundances at each experimental sampling point (Fig. 3C).

Detailing the taxonomy phylogenetic levels over time, a decrease was observed in both CTR and TR groups for the Firmicutes_A (CTR-T0 vs. -T2, from 27.5 to 13.3%, p = 0.006; CTR-T1 vs. -T2, from 27.0 to 13.3%, p = 0.005) phylum, while an increment in the Actinobacteriota phylum was noticed in both TR and CTR groups between baseline and the last sampling (TR-T0 vs. -T2, from 13.8 to 21.0%; p = 0.034; CTR-T0 vs. -T2, from 17.2 to 27.0%; p = 0.034). The Oscillospiraceae family relative abundance decreased for both TR and CTR groups over time (from 14.5% CTR-T0 and 13.8% TR-T0 to 6.3% CTR-T2 and 7.7% TR-T2), with a significant difference within the CTR (T0-T2 p = 0.001; T1-T2 p = 0.032). On the other hand, the Propionibacteriaceae family significantly rose from baseline to the T2 sampling for both groups, while bacteria belonging to the Streptococcaceae family significantly decreased their relative abundance in the CTR samples (from T0 2.4% to T2 0.2%, p = 0.013) and rose in the TR group (T0 0.7% vs. T2 1.8%, p = 0.001). Detailed phylum and family abundances are reported in Table S2.

As for the genus level, several differences have been observed between the two groups during the experimental period (Figure S1). In particular, when comparing the experimental samplings to the trial baseline (Fig. 3D and E), the Clostridiales bacterium group “UBA1774” significantly decreased its relative abundance from baseline within the CTR group (T0 12.5% vs. 5.3% T2, p = 0.002), while the TR group did not report such variations. Notably, the Cutibacterium genus significantly rose from T0 to T2 among both the CTR and the TR groups, respectively from baseline relative abundances of 10.6% (CTR-T0) and 8.6% (TR-T0) to 19.9% (CTR-T2, p = 0.007) and 16.1% (TR-T2, p = 0.006). Contrariwise, Phocaeicola showed a significant difference from T0 to T2 within the TR with a decreasing trend (4.4% vs. 1.8%, p = 0.043). As observed at the corresponding family level, over time Streptococcus almost disappeared among the CTR samples (from T0 2.4% to T2 0.2%, p = 0.013) and significantly increased within the TR group (T0 0.7% vs. T2 1.8%, p = 0.001).

Co-abundant bacterial groups

When analyzing the entire microbiota for each experimental group, a few joint bacterial variations were highlighted. Among the CTR samples (Fig. 4A), two main and highly correlated Co-Abundant Groups (CAG) were observed as opposed: one comprising the Beijerinckiaceae bacterium RH-AL1, Cutibacterium, bacterial members of the Burkholderiaceae family, Moraxella_A, and Pseudomonas_E, all with strong positive interactions (“Cutibacterium CAG”); and one made up of the Clostridium, Romboutsia, Phocaeicola, and UBA-1774 genera, all with negative correlations among one another (“UBA-1774 CAG”). The same interactions were observed in terms of positive/negative direction in the TR group (Fig. 4B), but not with the same intensity and significance; the two CAGs comprised other genera that in the CTR group did not report correlations, forming a wider bacterial relationship variety: in the TR group, the “UBA-1774 CAG’’ also included the LAB bacteria Alloiococcus, Alkalibacterium, Corynebacterium, and members of the Carnobacteriaceae family; the “Cutibacterium CAG” incorporated Staphylococcus and Akkermansia.

Fig. 4

Co-abundance analysis. Matrices and statistical significance (i.e. p < 0.05) for both (A) CTR and (B) TR experimental groups. Spearman correlation with Ward’s clustering has been employed for the identification of bacterial co-present (positive correlation, red tiles) or co-absent (negative correlation, blue tiles) genera

The LEfSe analysis

The search for presumed distinctive microbial biomarkers was performed using the LEfSe method. Linear discriminant analysis scores computed for differentially abundant marker genera between CTR and TR groups revealed that, at all phylogenetic levels, respectively 3 and 12 enriched groups were significantly identified. As reported in Fig. 5, the Paludibacteraceae family and the Xanthomonadaceae bacterium RF16, at both genus and species levels, have been identified as enriched within the CTR group, while the TR group was mainly enriched in Epilithonimonas (both the genus and the E. hispanica species), the Flavobacteriales order and its Weeksellaceae family, and in members of the Carnobacteriaceae family.

Fig. 5

LEfSe analysis. (A) Cladogram and (B) bacterial markers identified as belonging to the CTR (red) and TR (green) groups

Discussion

Nisin is the most studied bacteriocin since its discovery in 1928, and it is largely used for biomedical applications thanks to its antimicrobial and immunomodulatory properties [14]. In the veterinary field, the use of nisin-containing products for the prevention and treatment of bovine mastitis, as well as its potential as an alternative to commercial teat dip products, has been well explored and studied [9, 10, 13, 15]. A dip product for the disinfection of the cows’ teats before and after milking must not unbalance the milk microbiota, since it might impact the animal health and the subsequent raw milk processing in the dairy industry [6, 10]. We evaluated this crucial aspect by comparing the changes in cows’ milk microbiota composition when using a Nisin A-based teat dip product or a commonly used iodine-based disinfectant.

Especially over time and at later time-points, we found that the use of the Nisin A-based teat dip products was linked to a higher microbial diversity than using the commercial one. The microbiota of the TR samples displayed a greater species richness, according to the relative abundance exploration, the correlation profiles, and the linear discriminant analysis. Several distinct biomarkers characterized both groups, as it appeared that the microbial populations in the TR samples varied more throughout the experimental period, showing a higher species richness, while the overall composition of the CTR samples remained quite constant during the study. Similar trends were observed among both groups, but the interactions between bacterial elements differed: the TR group exhibited a wider variety of bacterial relationships, while the CTR was characterized by stronger but narrower correlations.

The greater biodiversity in the TR milk samples was highlighted also by the cultivable microbiota analysis: notably, S. warneri, E. faecium, E. casseliflavus, E. hirae, Lc. garviae, St. uberis, and St. mitis were peculiar to the TR samples, while E. villorum was detected only in the CTR samples.

Detailing the bacterial composition, we observed that the relative abundance of the most prevalent phyla was similar in the CTR and the TR groups: they were dominated by Firmicutes, Actinobacteriota, and Bacteroidota, typically reported phyla in bovine milk [3]. Firmicutes alone accounted for almost half of the average microbial relative abundance, while the microbial composition of the other two groups changed similarly throughout the experimental period. The increase of Actinobacteriota, in particular, was mostly due to the increase of the corresponding Propionibacteriaceae family and the Cutibacterium genus, which we observed to be significantly increased for both CTR and TR groups at T2. The seasonality may have influenced these changes in the taxonomic abundances, as there is a known effect of temperature and humidity on the raw milk microbiota composition [16,17,18]. Of particular interest is the increase of Actinobacteriota, and its Propionibacteriaceae family, previously reported from winter to spring [19, 20]; our findings place themselves in line with these studies, as our T0 was set in January and T2 in April. Intriguingly, we observed for the experimental group a consistent Cutibacterium increment along the study period, with the TR abundances being higher at T1 and T2 than in the CTR group.

Previous studies also reported an association between high temperature and higher percentages of both Streptococcus and Staphylococcus genera [19, 20]. Throughout our experimental period, we observed a net increase of the Streptococcus genus in the TR group in April (T2), while its abundance decreased in the CTR group. The average relative abundances of the Staphylococcaceae family and the Staphylococcus genus were as well significantly higher in the TR group, but with no clear trend over time.

Accordingly, the analyses of the cultivable microbiota regarding the Staphylococcaceae family showed that both CTR and TR groups were dominated by S. epidermidis (56.9% and 55.6%, respectively), followed by S. haemolyticus, especially abundant in CTR (34.5% vs. 18.3%). The third most abundant species in the TR group was S. warneri (9.5%), though not detected in the CTR group. Apart from the seasonality association, the role of non-aureus Staphylococci (NAS) in udder health is still controversial [21, 22], although it is known that their ability to produce bacteriocins may influence the composition of the microbial community of raw milk [23].

The use of Nisin A, which is a bacteriocin naturally produced by Lc. cremoris, as the active principle of both the experimental formulations could explain the observed lower alteration of the udder microbial composition and, as a consequence, higher respect of the resident microbiome as compared to the commercial iodine-based formulations. It is well known that the use of iodine-based commercial teat dip products affects the milk microbial composition, with consequences on the cheesemaking processes [24]. The results of both the bacteriological analysis, which is addressed only to the cultivable bacterial species, and the NGS analysis showed a greater diversity of the microbial composition in the TR group compared to the CTR group. The co-abundance analysis consolidated this observation as the TR group reported wider positive and negative associations between the bacterial genera, compared to the same but narrower CAGs observed for the CTR group.

In addition, we observed an overall tendency towards an increment of the Gram-negative bacteria within the TR group, while the CTR group reported higher Gram-positive relative abundances at each experimental sampling point. This was expected considering that Nisin-A had previously shown in vitro antimicrobial activity against Gram-positive pathogens and, contrarily, low activity against Gram-negative bacteria. This tendency confirms the efficacy of Nisin, which is reported as characterized by a wide spectrum of antibacterial activity against Gram-positive bacteria [25].

Lastly, the bacteriocins produced by LAB seem to protect the microbial species associated with a healthy udder. In fact, in our study, at the beginning of the trial, the CTR group showed higher relative abundances of LAB; throughout the time-points we recorded an increase of these species in the TR group, reaching similar abundances in both the experimental groups. Lactic acid bacteria are natural inhabitants of udder and teats and exogenous bacteria in milk. A recent study highlighted their protective activity against major mastitis pathogens such as Escherichia spp. and Serratia spp. (Addis et al. 2016)

Conclusion

The use of pre- and post-milking teat formulations containing the metabolites including Nisin A produced by Lactococcus cremoris FT27 appears to modulate the bovine milk microbiota by maintaining a greater microbial diversity as compared to the iodine-based commercial disinfectants. Further studies are required to better explain the role of these compounds on the udder microbial resident species and the role of these bacteria in maintaining the udder health.

Methods

Experimental design and sampling

The experimental design of the study has been previously detailed [26] and reported in Fig. 6. Briefly, the 210 lactating cows of the herd had been divided into two groups for a 3-month trial: (1) control group (CTR) – subjected to a pre-milking and a post-milking using commercial disinfectants containing lactic acid or iodine (3,000 ppm), respectively; (2) treated group (TR) – with pre- and post-milking disinfection performed using natural formulations including the fermented broth of Lc. cremoris FT27 strain, containing Nisin A. The fresh culture used in the preparation of pre‑ and post‑dipping products contained 9.3 log₁₀ CFU/mL of the FT27 strain and about 950 IU/mL of Nisin A.

Fig. 6

Experimental design reporting the sampling process, the formulations, and the time-points involved

The basic components of the two experimental formulations were: lactic acid for the pre-dipping, and emollients, skin protective substances, and film-forming elements for the post-dipping, as those used in commercial disinfectants. Based on the previous in vitro results of the antimicrobial activity of the Nisin A produced by the Lc. cremoris FT27 strain against the main mastitis pathogens [26], the dilution of the fermented broth was 1:4 in both the experimental formulations. The addition of the basic components caused the complete loss of vitality by Lc. cremoris FT27 [26].

The udder health was monitored by bacteriological analysis and somatic cell counts (SCC) every two weeks, sampling all the cows at quarter level. At the laboratory, 10 µL of milk was plated on blood agar plates containing 5% defibrinated bovine blood and incubated aerobically at 30 °C with evaluation after 24 and 48 h according to the guidelines of the National Mastitis Council [27] NMC, 2017). The colonies were identified using the matrix-assisted laser desorption/ionization time-of-flight mass spectrometry (MALDI-TOF MS). Moreover, SCC was determined by an automated fluorescent microscopic somatic cell counter (Bentley Somacount 150, Bentley Instrument, Chaska, MN, USA).

For the microbiota analysis, 10 cows from each experimental group were chosen, based on negative bacteriological results and low SCC (< 100,000 SCC/ml), and milk samples were collected at T0 (January), after one (T1, February) and three months (T2, April) of treatment. In order to consider a minimum number of animals, and to limit the costs for genomic analyses, two quarters were selected from for each cow, with the attention of avoiding those quarters with the homolateral one affected by high SCC.

Bacterial enumeration and isolation

One milliliter of each raw milk sample was serially diluted in quarter-strength Ringer’s solution (Scharlau Microbiology, Barcelona, Spain) for bacterial enumeration. Total mesophilic aerobic bacteria (TMB) were counted on Petrifilm Aerobic Count Plate (3 M, Minneapolis, MN, USA) after incubation at 30 °C for 72 h (ISO 4833–1:2022). Enterobacteriaceae were enumerated on Petrifilm Enterobacteriaceae Count Plate (3 M) at 37 °C for 24 h (ISO 21528–1:2017) and lactic acid bacteria (LAB) were determined on Petrifilm Lactic Acid Bacteria Count Plate (3 M) at 30 °C for 48 h.

As bacterial cells in quarter milk samples were very low, milk samples were also incubated at 37 °C overnight and then striked on Homofermentative Heterofermentative Differential (HHD) agar (Biolife Italiana), to increase the possibility of bacterial isolation.

All the colonies with different morphology were picked up from countable plates and incubated milk and were sub-cultured in de Man Rogosa and Sharpe (MRS) broth (Biolife Italiana, Milan, Italy) and in Brain Heart Infusion (BHI) broth (Biolife Italiana) at 30 °C. The purity of the isolates was checked by streaking repeatedly on HHD agar. After purification, the isolates were examined for cell morphology and catalase activity and then stored at − 20 °C in Litmus milk (Biolife Italiana). The statistical analysis of microbiological data was presented as means ± standard deviation (SD). Significant differences (P < 0.05) among the data were calculated by one-way ANOVA using Minitab ver. 14.13 (Minitab Inc.).

Strain typing and identification

Randomly Amplified Polymorphic DNA-PCR (RAPD-PCR) method was applied to explore the biodiversity and genetic relatedness within the isolated strains. The genomic DNA of the selected isolates was extracted from overnight cultures using the Microlysis kit (Aurogene, Rome, Italy) following the manufacturer’s instructions. Yield and purity of DNA were evaluated using the Infinite F200 PRO microplate reader (Tecan, Mannedorf, Switzerland). RAPD-PCR analysis was performed with 3 primers (M13, D11344, and D8635) as described by Morandi et al. [28]. The resulting fingerprints were compared with the BioNumeric 5.0 software package (Applied Maths, Sint-Martens-Latem, Belgium), using the UPGMA (unweighted pair group method with arithmetic averages) cluster analysis. Strains with a similarity coefficient ≥ 90% were considered to belong to the same biotype [28]. One representative strain of each cluster obtained by RAPD analysis was afterwards identified by matrix-assisted laser desorption ionization-time of flight mass spectrometry (MALDI-TOF MS); in detail, overnight cultures were plated on MRS and BHI agar plates and cell material from an isolated colony was deposited on the target plate using a toothpick. Samples were overlaid with 1 µL of -cyano-4-hydroxycinnamic acid in 50% acetonitrile with 2.5% trifluoroacetic acid (Bruker Daltonik GmbH, Bremen, Germany). The spectra were acquired with a microFlex™ mass spectrometer (Bruker Daltonik GmbH) in the positive mode. Bacterial Test Standard (Bruker Daltonik GmbH) was used for the Instrument Calibration. Spectra were automatically interpreted by the database MBT Compass^® 4.1. A log (score) of 1.7 was the threshold for the genus-level identification and a log (score) of 2.0 was the threshold for the species-level identification.

DNA extraction and purification

For each quarter, 5 mL of milk was analyzed by using a DNA extraction method based on the combination of a chaotropic agent, guanidium thiocyanate, with silica particles, to obtain bacterial cell lysis and nuclease inactivation as previously described [29, 30]. All the DNA samples were assessed using a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) to check the quality and quantity. The isolated DNA was stored at − 20 °C until use.

16S rRNA gene library preparation and sequencing

Bacterial DNA was amplified using the primers described in the literature [31] which target the V3–V4 hypervariable regions of the 16S rRNA gene. All PCR amplifications were performed in 25 µL volume per sample.

A total of 12.5 µL of Phusion High-Fidelity Master Mix 2x (ThermoFisher Scientific, Waltham, MA, USA) and 0.2 µL of each primer (100 µM) were added to 2 µL of genomic DNA (5 ng/µL). A first amplification step was performed in an Applied Biosystem 2,700 thermal cycler (ThermoFisher Scientific). Samples were denatured at 98 °C for 3 min, followed by 25 cycles with a denaturing step at 98 °C for 30 s, annealing at 56 °C for 1 min, and extension at 72 °C for 1 min, with a final extension at 72 °C for 7 min. Amplicons were cleaned with Agencourt AMPure XP (Beckman, Coulter Brea, CA, USA) and libraries were prepared following the 16S Metagenomic Sequencing Library Preparation Protocol (Illumina, San Diego, CA, USA). The libraries obtained were quantified by RealTime PCR with KAPA Library Quantification Kits (Kapa Biosystems, Inc., MA, USA), pooled in equimolar proportion, and sequenced in a MiSeq (Illumina) run with 2 × 250-base paired-end reads.

Bioinformatic pipeline and statistical analysis

Demultiplexed paired-end reads from 16S rRNA-gene sequencing were first checked for quality using FastQC [32] for an initial assessment. Amplicon sequence variants (ASVs) were identified from 16S paired-end sequencing using the Divisive Amplicon Denoising Algorithm (DADA2) pipeline, including filtering and trimming of the reads (version 1.18.0) [33]. Reads per sample were trimmed to 4195 to compensate for the sequencing unevenness of the samples and to provide a consistent minimum amount for the downstream analysis, carried out through the “phyloseq” package (version 1.34.0) [34].

Alpha-diversity evaluation was performed according to several microbial diversity metrics (i.e., chao1, Shannon Index, Observed species, and Faith’s phylogenetic distance/PD whole tree). Beta‐diversity analysis was conducted using both weighted and unweighted Unifrac metrics [35], and through the principal coordinates analysis (PCoA). For both alpha- and beta-diversity, statistical significance was corrected using the Benjamini-Hochberg procedure. Taxonomy was assigned to the ASVs using the 8‐mer‐based classifier from the 11.5 release of the RDP database [36] and using the GTDB 16S rRNA database (release r207) [37]. Statistical comparisons between groups were assessed through the Mann-Whitney U-test for unpaired data (i.e. the experimental groups) and the Dunn’s test for multiple comparisons (i.e. time-points). The LEfSe analysis was conducted through the MicrobiomeMarker R package (version 1.8.0) [38], identifying the group markers via the Limma-Voom differential analysis and the Benjamini-Hochberg correction. Throughout the paper, all p-values reported are adjusted following the methods listed above, and a significance threshold of p = 0.05 was applied.