Microbiome mapping in dairy industry reveals new species and genes for probiotic and bioprotective activities

Abstract

The resident microbiome in food industries may impact on food quality and safety. In particular, microbes residing on surfaces in dairy industries may actively participate in cheese fermentation and ripening and contribute to the typical flavor and texture. In this work, we carried out an extensive microbiome mapping in 73 cheese-making industries producing different types of cheeses (fresh, medium and long ripened) and located in 4 European countries. We sequenced and analyzed metagenomes from cheese samples, raw materials and environmental swabs collected from both food contact and non-food contact surfaces, as well as operators' hands and aprons. Dairy plants were shown to harbor a very complex microbiome, characterized by high prevalence of genes potentially involved in flavor development, probiotic activities, and resistance to gastro-intestinal transit, suggesting that these microbes may potentially be transferred to the human gut microbiome. More than 6100 high-quality Metagenome Assembled Genomes (MAGs) were reconstructed, including MAGs from several Lactic Acid Bacteria species and putative new species. Although microbial pathogens were not prevalent, we found several MAGs harboring genes related to antibiotic resistance, highlighting that dairy industry surfaces represent a potential hotspot for antimicrobial resistance (AR) spreading along the food chain. Finally, we identified facility-specific strains that can represent clear microbial signatures of different cheesemaking facilities, suggesting an interesting potential of microbiome tracking for the traceability of cheese origin.

Fig. 1. Different taxonomic profiles are present in cheese and environmental swabs.

A. Taxonomic composition of the 30 most abundant species in all samples, with the average relative abundances of each species in each category. B. Alpha diversity, in terms of observed species, Shannon and Simpson analysis of samples by category, and (C) beta diversity analysis represented by a Principal-Coordinate Analysis (PCoA) plot of Bray-Curtis distance, with ellipses representing clustering by category. Boxes represent the interquartile range (IQR) between the first and third quartiles, and the line inside represents the median (2nd quartile). Whiskers denote the lowest and the highest values within 1.5 x IQR from the first and third quartiles, respectively. The significance was tested by applying pairwise Wilcoxon test. Average values are obtained from n = 308, 199, 37, 524, and 182 biologically independent samples from food contact, non food contact, operators’ swabs, final products and cheese-related materials, respectively. The category “cheese-related materials” groups together milk, brine and whey culture.

Fig. 2. Different functional profiles are present in cheese and environmental swabs.

A Non-metric Multidimensional Scaling (NMDS) plot based on Jaccard’s distance of functional profiles obtained by HUMAnN. Samples are colored according to the sample type. B Boxplots showing the abundance (log values) of UniRef50 genes detected in the different sample groups. Boxes represent the interquartile range (IQR) between the first and third quartiles, and the line inside represents the median (2nd quartile). Whiskers denote the lowest and the highest values within 1.5 × IQR from the first and third quartiles, respectively. The significance was tested by applying pairwise Wilcoxon test. The category “cheese-related materials” groups together milk, brine and whey culture. Average values are obtained from n = 308, 199, 37, 524, and 182 biologically independent samples from food contact, non food contact, operators’ swabs, final products and cheese-related materials, respectively.

Fig. 3. Genes related to gastrointestinal tract transit stress resistance and engraftment encoded in the metagenomes.

A Heatplot showing the abundance (Reads per Kilobase per Million, RPKM) of predicted genes related to gastrointestinal tract transit stress resistance and engraftment. B–F Boxplots showing the abundance (log RPKM) in the different sample groups of specific genes. Boxes represent the interquartile range (IQR) between the first and third quartiles, and the line inside represents the median (2nd quartile). Whiskers denote the lowest and the highest values within 1.5 × IQR from the first and third quartiles, respectively. The significance was tested by applying pairwise Wilcoxon test. The category “cheese-related materials” groups together milk, brine and whey culture. Average values are obtained from n = 308, 199, 37, 524, and 182 biologically independent samples from food contact, non food contact, operators’ swabs, final products, and cheese-related materials, respectively.

Fig. 4. Analysis of the bacteriocins encoded in the metagenomes.

The category “cheese-related materials” groups together milk, brine and whey culture. A Boxplot reporting the number of observed bacteriocin genes for each sample. B Genes observed with a significatively higher frequency in one of the sample groups. Percentages report the proportion of samples from each category testing positive for the selected genes. Boxes represent the interquartile range (IQR) between the first and third quartiles, and the line inside represents the median (2nd quartile). Whiskers denote the lowest and the highest values within 1.5 × IQR from the first and third quartiles, respectively. The significance was tested by applying pairwise Wilcoxon test. Average values are obtained from n = 308, 199, 37, 524, and 182 biologically independent samples from food contact, non food contact, operators’ swabs, final products and cheese-related materials, respectively.

Fig. 5. Analysis of antimicrobial-resistance and virulence-associated genes.

The category “cheese-related materials” groups together milk, brine and whey culture. A Connected barplot showing the variation in the average Copies Per Million (CPM) abundance of each Antimicrobial Resistance Gene (AMRG) family between the classes of samples. B Taxonomic assignment of contigs encoding for AMRG. For each AMRG - taxon pair, each slice is proportional to the percentage of contigs reconstructed from a category of samples. The size of each pie is proportional to the number of contigs linked to an AMRG - taxon pair. C Sum of the CPM abundance of genes linked with Adherence, Biofilm formation and Motility virulence traits. Statistical differences between the groups were calculated through the Wilcoxon’s rank-sum test. Boxes represent the interquartile range (IQR) between the first and third quartiles, and the line inside represents the median (2nd quartile). Whiskers denote the lowest and the highest values within 1.5 x IQR from the first and third quartiles, respectively. The category “cheese-related materials” groups together milk, brine and whey culture. Average values are obtained from n = 308, 199, 37, 524, and 182 biologically independent samples from food contact, non food contact, operators’ swabs, final products and cheese-related materials, respectively.

Fig. 6. Phylogenetic tree of Metagenome-Assembled Genomes (MAGs) reconstructed from cheeses and dairy environment.

Phylogenetic tree of all the MAGs reconstructed in this study, spanning 2111 Species-level Genome Bins (SGBs). From outer to inner, rings are colored according to phylum-level taxonomic assignment, sample type and identification of SGBs at species level (as reported in the Methods). The category “cheese-related materials” groups together milk, brine and whey culture. MAGs have been reconstructed from n = 308, 199, 37, 524, and 182 biologically independent samples from food contact, non food contact, operators’ swabs, final products, and cheese-related materials, respectively.

Fig. 7. MAGs distribution in cheeses and dairy environment.

A Bar chart showing the number of Metagenome-Assembled Genomes (MAGs) for the top 30 Species-level Genome Bins (SGBs) reconstructed for each sample group. B Heatplot reporting, for each of the top 30 SGBs, the proportion (%) of MAGs reconstructed from each sample group. FC, food-contact surfaces; NFC, non food-contact surfaces; Ripened, ripened cheeses (>30 days); Unripened, unripened cheeses (<30 days); CRM (cheese-related materials), that groups together milk, brine, and whey culture.

Fig. 8. MAGs harbor Antimicrobial Resistance and bacteriocin coding genes.

For each group of samples, the heatmaps show the proportion (%) of Metagenome-Assembled Genomes (MAGs) within a Species-level Genome Bin (SGB) containing at least 1 Antimicrobial Resistance (AMR; A) or bacteriocin (B) coding gene. The side color bar is colored according to order-level taxonomic assignment of the SGB. The bar charts represent the proportion (%) of each AMR or bacteriocin class within the SGB. The category “cheese-related materials” groups together milk, brine, and whey culture.

Fig. 9. MAGs of Lb. delbrueckii can be clustered in different putative strains.

Non-metric Multidimensional Scaling (NMDS) based on ANI distance matrix of MAGs belonging to SGB_0 (Lb. delbrueckii). A All MAGs in SGB_0 are shown. Points are colored according to the cheese type, while different shapes indicate the sample types. B Only MAGs reconstructed from Caciocavallo cheese facilities are included and points are colored according to the facility code. C Phylogenetic tree of MAGs belonging to SGB_0 reconstructed from Caciocavallo cheese facilities. The category “cheese-related materials” groups together milk, brine and whey culture.