Shared characteristics and patterns of succession
The experimental setup included cheeses from three different producers who used the same cave for ripening. In addition, for one producer, cheeses derived from a single production batch were ripened in three different caves, in order to identify microbiome signatures associated with the respective cave microenvironments (Fig. S1).
Regardless of producer or cave used for ripening, cheeses maintained certain shared characteristics and microbial succession patterns. Before cheeses entered the ripening caves (Stage1; 30 days post-manufacture) they showed a mild acidic pH (5.41-6.58), and some degree of drying (water activity of 0.828-0.924). Cabrales cheeses are routinely rind-washed throughout ripening in the caves, where the relative humidity is >90%. Consequently, significant (p<0.05) increases in cheese water activity, relative humidity and pH, and decreases in dry matter and salt content occurred thereafter, marked mostly on cheese rinds at intermediate (Stage2) and final (Stage3) stages of ripening (Fig. S2A). Total fat and protein contents remained fairly stable through ripening (Fig. S2A).
Forty-two volatile compounds, belonging to eight chemical classes (Table S1), and seven biogenic amines were detected. While the profile of volatile compounds and biogenic amines was highly variable depending on the producer and ripening cave (as described in detail below), a global view shows that the total content in volatile compounds increased on cheese rinds as ripening proceeded while no significant changes were observed for the cheese cores (Fig. S2B). The biogenic amines content significantly increased with ripening time for both core and rind samples (Fig. S2B).
Regarding the microbiome composition, counts of total aerobic bacteria, lactic acid bacteria, Enterobacteriaceae, and yeasts and molds reached maximum levels at Stage1 and then gradually declined during ripening, both for cheese cores and rinds (Fig S3A). At Stage1, cheeses were dominated by lactic acid bacteria, mainly Lactococcus and members of the former Lactobacillus genus, and by the yeast Debaryomyces and the fungi Penicillium and Geotrichum (Fig. S3B). Whereas fungal populations remained fairly stable along ripening, a significant decrease was evident in the relative abundance of Lactococcus and Lactobacillus, among others (Fig. S3B), while other taxa, such as Tetragenococcus, Brevibacterium and Corynebacterium emerged with high relative abundances at Stage2 and Stage3, especially at rind level (Fig. S3B).
Metabolome and microbiome differences by producer
A principal component analysis, performed to reflect all of the volatile compounds and biogenic amines, identified that the producer variable had a significant influence on the metabolomic profile of cheeses, both for the core and rind (adonis, p=0.001), with the first two components explaining 80.40 % and 74.80 % of the total variation, respectively (Fig. 1A,B). Focusing on the eight most abundant volatiles, some trends were apparent; higher levels of 2-heptanone, 2-nonanone, 2-heptanol and 2 -pentanone were found in cheeses from Producer3; higher concentrations of ethyl hexanoate and ethyl butanoate were found in Producer2 cheeses; and the abundance of hexanoic and butanoinc acids were lower for Producer3 cheeses (Fig. 1C, Table S2). Regarding amines, tryptamine and tyramine were more abundant on Producer2 and Producer3 cores and cadaverine on Producer1 rinds (Fig. 1D, TableS2).
The cheese producer also had a significant influence on the bacterial and fungal taxonomic profile of samples (Fig 2A,B, Fig S4A,B), both for cheese cores (adonis, p=0.001 for both bacteria and fungi) and rinds (adonis, p=0.001 for bacteria and p=0.018 for fungi). At cheese core level, several bacterial genera within the main ones such as Lactobacillus, Corynebacterium, and Enterococcus, among others, were associated with Producer2, while Lactobacillus and Tetragenococcus were predominantly linked to Producer3 (Fig. 2A,C, Table S2). At cheese rind level, several bacterial genera were clearly associated with Producer2 (e.g., Corynebacterium and Tetragenococcus) while Lactococcus was associated with Producer1 and Producer3 (Fig. 2B,C). Regarding fungi, core samples were highly abundant in Penicillium, but, for Producer2, similar relative abundances of Geotrichum and Penicillium were observed. At rind level Penicillium and Debaryomyces showed high abundances in samples from Producer1 and Producer2, and Debaryomyces was the dominant genus for Producer3 (Fig. 2D, Table S2).
Metabolome and microbiome differences by ripening cave
For Producer2, the ripening cave also had a significant influence on the metabolomic profile of cheese cores and rinds (adonis, p=0.001), with the first two components explaining 71.3 % and 73.2 % of the total variation, respectively (Fig. 3A,B). Focusing on the eight most abundant volatiles, Cave3 cheeses differed most due to lower levels of ethyl hexanoate and ethyl octanoate, and higher levels of 2-heptanone on cores and rinds and hexanoic acid, butanoic acid and octanoic acid on rinds (Fig. 3C, Table S3). Regarding amines, tyramine was more abundant in Cave2 cores, tryptamine in Cave2 cores and Cave3 rinds, and cadaverine in Cave1 cores and rinds (Fig. 3D, TableS3).
The cave also had a significant influence, although lower than that found for the variable producer, on the taxonomic profile of cheese cores (adonis, p=0.007 for bacteria and p=0.049 for fungi) and rinds (adonis, p=0.027 for bacteria and 0.653 for fungi) (Fig. 4A,B, Fig. S4C,D). The most relevant differences among caves in the abundance of some particular taxa were those observed for Lactococcus and Lactobacillus, which were more abundant on core samples from Cave1 and Cave3, respectively, at Stage3; Tetragenococcus, Corynebacterium and Staphylococcus, which were more abundant on rind samples from Cave1, Cave2 and Cave3, respectively (Fig. 4A,C, Table S3); and Debaryomyces, which was significantly more abundant on cheese rinds from Cave1 than on those from the other two caves (Fig. 4D, Table S3).
The cave microbiome strongly shapes the rind cheese microbiome
The characterisation of the microbial communities prevailing in some primary sources that could determine the cheese microbiome (i.e., milk, curd, different factory processing environments and cave environments) showed that milk was dominated by Pseudomonas, followed by Acinetobacter, Staphylococcus or Lactococcus, at different levels depending on the producer (Fig. S5). Similar profiles but with a higher abundance of some lactic acid bacteria, such as Lactococcus, were found in curd samples (Fig. S5). Brevibacterium, Psychrobacter, Pseudomonas, Acinetobacter, Penicillium, Debaryomyces and Kluyveromyces were the dominant genera on food processing environments (both food contact and non food contact) of the three cheese producing plants (Fig. S5). Brevibacterium, Corynebacterium and Debaryomyces were the dominant genera in cave food contact surfaces, while Penicillium and a wide range of non-dominant bacterial and fungal genera prevailed on the caves´ non food contact surfaces (Fig. S5).
A source tracking analysis revealed that the bacterial composition of cheese cores, and of cheese rinds at Stage1, was mainly determined, for Producer1 and Producer2, by the curd microbiota, and for Producer3 by the microbiota of food contact environments from the processing plant (Fig. 5). The main bacterial genera traced back to the curd and dairy food contact surfaces were Lactococcus, Lactobacillus and, in the case of Producer2, Hafnia (Fig. 5). Over subsequent stages of ripening, the curd microbiota still represented an important source for the cheese core microbiome, although the cave environment had an important role in shaping the microbiome of cheese cores on Cave2 at Stage2 and Stage3, mainly being a source of Tetragenococcus and Corynebacterium for Producer1 and Producer3, and Corynebacterium for Producer2 (Fig. 5). Other less abundant genera (e.g., Brachybacterium, Alkalibacterium, Staphylococcus) were also traced back to the cave environments. Regarding cheese rinds, their bacterial microbiota at Stage2 and Stage3 were traced back mainly to the cave environments, which were revealed as a likely source of Brevibacterium, Corynebacterium and Tetragenococcus (on Stage3) and other minority genera (on Stage2), for Producer1; Tetragenococcus in Cave1 and Cave3, together with Corynebacterium in Cave2, for Producer2; and Tetragenococcus, Brevibacterium and other minority genera (on Stage2) and Tetragenococcus (on Stage3) for Producer3 (Fig. 5).
A source tracking analysis undertaken on the fungal communities showed that food contact surfaces from Cave2 were a source of Debaryomyces in rind samples from Producer1 and Producer3, and milk was a source of Geotrichum in core samples from Producer2 (Fig. S6).
MAGs analysis confirmed the relevance of cave environments as a source of bacteria
A total of 110 high-quality metagenome assembled genomes (MAGs), 227 medium quality MAGs with completeness > 90%, and 276 medium-quality MAGs with completeness of 50-90% were obtained. They were classified into 70 genera, with the main genera represented being Lactococcus (103 MAGs), Corynebacterium (74 MAGs), Lactobacillus (68 MAGs), Tetragenococcus (54 MAGs), Bifidobacterium (30 MAGs), Staphylococcus (28 MAGs), Brevibacterium (28 MAGs) and Yaniella (22 MAGs). Cheese samples were the ones with the highest number of MAGs (413 out of 613 MAGs), while raw materials, processing environments at factory level and cave environments yielded 43, 74 and 83 MAGs, respectively (Fig. 6A, Fig. S7).
MAGs from several putative new species were obtained on samples from cave surfaces or cheeses. These included 14 Brevibacterium sp. MAGs most closely related to B. sandarakinum, 8 Corynebacterium sp. MAGs related to C. nuruki, 6 Corynebacterium sp. MAGs related to C. glyciniphilum, 2 Psychrobacter sp. MAGs and one Tetragenococcus sp. MAG. related to T. koreensis. Remarkably, despite the low relative abundance of Yaniella found at read level, a total of 22 MAGs from this genus, closely related to the recently discovered Candidatus Yaniella excrementavium, were obtained (Fig. 6A, Fig. S7). Furthermore, 3 MAGs assigned to the former genus Lactobacillus, obtained from cheese core samples from Producer3 at Stage1, were not assigned at species level, with L. selangorensis and L. camelliae being the most closely related species (Fig. 6A).
It is worth noting that MAGs assigned to Lactococcus, a dominant genus at Stage1, presented phylogenetic differences by producer, with Lactococcus lactis MAGs from Producer2 clustering separately from MAGs from other producers, L. piscium/carnosum MAGs being exclusively found for Producer2 and L. raffinolactis and L. garvieae for Producer3 (Fig. 6A). On the other hand, different closely related MAGs from Corynebacterium casei, Staphylococcus equorum, Brevibacterium sp., Tetragenococcus halophilus, T. koreensis and Yaniella sp. were obtained from both cheese and cave environmental samples. Moreover, in the case of cheese MAGs from Brevibacterium, Tetragenococcus and Yaniella, they were only obtained from cheeses at Stage2 and Stage3 (Fig. 6A, Fig. S7). Similar results were observed at read level, with the absence of these genera at Stage1, while they were within the main ones found at Stage2 and Stage3 (Fig. S8). It was also interesting to see that two different clusters of S. equorum, with differentiated functional potential, were detected colonizing food processing environments within factories and cheeses at Stage1, and cave environments and cheeses at Stage2 and Stage3, respectively (Fig. 6A, Fig. S9).
The analysis of the functional potential of MAGs evidenced that Lactococcus, Lactobacillus, Staphylococcus and Tetragenococcus had a relatively similar functional profile, which was characterised by the high abundance of functions related to the metabolism of carbohydrates, while Corynebacterium, Brevibacterium and Yaniella MAGs had a very different functional profile, characterised by the high abundance of different functions mainly related to protein and fatty acids metabolism (Fig. 6B).
The cheese microbiome is associated with cheese quality markers
Correlation analyses undertaken with respect to cheese samples from stages 2 and 3 of ripening revealed that, among the species for which MAGs were generated, at core level, Lactococcus (L. lactis, L. piscium) Lactobacillus (mainly L. paracasei) and Bifidobacterium (B. crudilactis, B. mongoliense) were those that had the strongest positive correlations with the abundance of various volatile compounds and biogenic amines. The most remarkable associations found were for L. lactis with pentanoic and hexanoic acids; L. piscium with 3-methyl-butanal, propanoic acid, ethyl-butanoate, propyl-butanoate, trimethylpyrazine, ethyl-octanoate, putrescine and cadaverine; L. paracasei with 3-methyl-butanal, propyl-butanoate, trimethylpyrazine and isobutyl-decanoate; B. crudilactis with 3-methyl-butanal and trimethylpyrazine; and B. mongoliense with trimethylpyrazine, isoamyll-hexanoate, isopentyl-octanoate, isobutyl-decanoate and tryptamine (Fig. 7A,B). At rind level, these species also showed positive correlations with the abundance of certain volatiles and amines, but other representative taxa, traced back to the cave environment, such as Tetragenococcus (mainly T. halophilus), Staphylococcus equorum, and Corynebacterium (C. flavescens, C. variabile) were also strongly associated with some cheese volatile and amine markers. The most remarkable associations were found for T. halophilus with propyl-butanoate, 2-heptanol, trimethylpyrazine, 2-pentylbutyrate, 2-nonanone, 2-nonanol, ethyl-octanoate, 2-undecanone, ethyl-decanoate, isobutyl-decanoate, ethyl dodecanoate, putrescine and tyramine; S. equorum with 2-heptanone, 2-heptanol, 2-octanone, 2-nonanol, γ-hexalactone, 2-undecanone, isopenthyl-octanoate, isoamyl-hexanoate and putrescine; and C. variabile and C. flavescens with pentanoic and hexanoic acids (Fig.7B). Interestingly, both at cheese core and cheese rind level, with the exception of Brevibacterium and the former Lactobacillus genus, species from the same genera showed highly different microbiome-volatilome correlation patterns while species belonging to the same genera, with the exception of Corynebacterium and Tetragenococcus, had similar correlations with amine concentrations (Fig. 7).
The cheese microbiome is rich in Horizontal Gene Transfer events
In order to identify among the main bacterial taxa cues of microbiome acclimatisation to cave and cheese microenvironments, signs of horizontal gene transfer (HGT) events were seeked in the available MAGs. A total of 23,001 HGT events, containing 67,411 coding regions, were detected between 56 taxa above species level, with Lactococcus, Tetragenococcus and Staphylococcus being the main genera involved (59.3, 52.0 and 23.3 % of total HGT events, respectively), and Lactococcus-Tetragenococcus, Staphylococcus-Tetragenococcus and Lactococcus-Staphylococcus as the main HGT pairs (29.6, 10.9 and 9.4 % of total HGT events, respectively) (Fig. 8A). HGT events were identified from MAGs belonging to all surfaces sampled, but cheese core and rind samples at stages 2 and 3 of ripening contributed the most (56.6 % of HGT events) (Fig. 8B).
Up to 35.3 % of the HGT events were associated with plasmid sequences, according to plasflow analysis, while only 0.5 % carried relaxase encoding genes, associated with plasmid mobilization. Moreover, 11.7, 8.1 and 0.9 % harboured prophages, transposases and integrases, respectively. Although 55.5 % of the coding regions could not be assigned to ko codes of the KEGG Orthology database, 10.6, 7.2 and 5.4 % were assigned at level2 to the groups Protein families: signaling and cellular processes, Carbohydrate metabolism and Protein families: genetic information processing, respectively. At level3 of KO classification, Prokaryotic defense system, Galactose metabolism, Glycolysis/Gluconeogenesis, Purine metabolism and Fructose and mannose metabolism were within the dominant functions (Fig. 8C). Finally, 68 HGT-associated contigs were longer than 20kb, containing from 12 to 45 coding regions (CDS) (Fig. 8D), and clustered into 26 groups. Lactococcus-Tetragenococcus and Alkalibacterium-Tetragenococcus were the dominant pairs within these longest HGT events, with 14 and 4 clusters, respectively (Table S4). Half of the obtained HGT-clusters contained genes related to mobile genetic elements such as plasmids, transposases, phages and integrons; 3 clusters harboured genes related to β-lactams resistance, and up to 9 clusters contained genes related to protein, carbohydrate or lipid metabolism (Table S4).