Shotgun Metagenomics of Traditional Bulgarian Green Cheese Reveals Key Roles of Brevibacterium aurantiacum and Aspergillus puulaauensis in the Ripening Process

1. Introduction

The diverse microbial communities responsible for cheese maturation (ripening) represent the major factors that determine the attributes of the final product, following the specifics of production [1]. These attributes include flavor, odor, texture and color, which can represent the signature characteristics of a particular type of cheese. Therefore, to better understand the relationship between the composition of microbial communities and the corresponding characteristics of the cheese, shotgun metagenomic sequencing could be used as an approach. This method not only identifies the present species but also allows for complex functional analyses based on their genomic material.

Traditional Bulgarian green cheese, produced from raw ewe or goat milk and relying on spontaneous fermentation, is produced only in Cherni Vit village (near the town of Tetevem, Bulgaria). It was recently subject to amplicon sequencing in a prior study [2]. In the current study, we rely on the more detailed approach of short-read shotgun metagenomic sequencing to further elucidate the function and better catalogue the profile of the microbial community within the cheese, as well as to establish the presence of any antibiotic resistance or virulence genes. While amplicon sequencing, which targets a specific marker gene (most commonly the 16S rRNA for bacteria) and uses the abundance of that marker as a proxy for species abundance, is as very fast and detailed approach to microbial profiling, the presence of a PCR amplification step, which is required for isolating the marker gene, introduces a level of uncertainty to the final abundance estimates for a particular clade. An alternative approach is shotgun metagenomic sequencing, in which the entire collection of genomes in the sample (the metagenome) is sequenced after random fragmentation without the need for amplification. Shotgun sequencing also provides additional information, as it encompasses all parts of the genome, not just a single marker gene. This additional information can be used to identify the presence of virulence and antibiotic resistance genes, as well as to extract metagenome-assembled genomes (MAGs) and to map out the metabolic pathways present in the community, which, in the case of cheese, are the pathways that ultimately contribute to the product’s attributes.

The process of production of traditional Bulgarian green cheese was recently described in detail in the previous amplicon-based analysis [2]. In brief, it is produced from unpasturazed sheep milk, and follows a process of brining. The final, ripening step occurs at 16-20 °C for 35-40 days, after which the surface molds are brushed, and the cheese is conserved for up to three years. We believe the unique production steps and environment the cheese is maturated under sets the stage for a specific combination of microorganisms that give the product its sensory properties. This, along with the rarity of the food product, and the interest to assure its consumption is safe, prompted us to carry out the current analysis.

2. Materials and Methods

2.1. Sample Collection, Storage and Preparation

A batch of 12-month-ripened Bulgarian green cheese, prepared from ewe’s milk, was collected directly from the dairy plant on the first days of June 2024. The cheese was transported to the laboratory at 4 °C, where a single 250 mg sample was taken and homogenized using a microcentrifuge tube homoginizer under sterile conditions. The sample included both the crust and the interior, and was cut out using a pyramid technique. Total DNA from the sample was isolated using the Quick-DNA Fecal/Soil Microbe Miniprep Kit, cat. # D6010 (Zymo Research, USA). The DNA concentration was determined by the Quantus fluorimeter (Promega, USA). The obtained DNA was stored at -20 °C. Part of the DNA sample was shipped in dry ice for commercial metagenomic sequencing (BGI, China), and the raw sequencing data were downloaded for downstream processing and analyses.

2.2. Sample Sequencing and Quality Control

Reads were sequenced using Illumina’s NovaSeq 6000 150 bp pair-end sequencing platform. cutadapt (v.5.1) [3] was used to quality filiter the reads using a minimum read length of 150 bp, a maximum error rate of 1, and assuring no left over adapters from sequencing. The reference human genome (GRCh38 - hg38) was used to construct a Kraken2 database, against which the reads were aligned, and human reads were filtered from the data.

2.3. Taxonomic Profiling

For all software, default settings were used unless specified otherwise. The taxonomy was assigned using Kraken2 (v. 2.1.6) (23) against NCBI’s complete nt database (accessed 27.09.2025). The reference genomes of prevailing species and subspecies (>0.1% of classified read abundance) were used to create a KrakenUniq (v.1.0.4) database [3] to assess the presence of false species positives using the unique k-mer counting approach. At 73 million reads, a threshold of 146,000 unique k-mers was chosen to filter out species false positives, which was based on the inflection point of the cummulitive plot of unique k-mers, effectively separating true positives from spurious assignments. False species and subspecies counts were redistributed to the above taxonomic level, and reads were then redistributed to the species level using Bracken (v.3.0.1)[4]. To get the final species abundance, read counts were normalized by the genome size of the species reference genomes.

2.4. Functional Profiling

For protein assignment, forward and reverse reads were independently aligned against NCBI’s nr_clustered database (accessed 27.09.2025), which contains representative sequences of proteins grouped within a 90% similarity threshold. To accelerate alignment, the database and reads were divided into prokaryotic and eukaryotic categories, and only reads with assigned taxonomy were considered. Diamond (v.2.1.13)[5] was used for the alignment, and 61.59% of reads were assigned an NCBI protein ID. Hits in the database were then assigned a KEGG ortholog (KO) using kofamscan (v. 1.3.0) (accessed on 25.03.2025), which was able to assign a KO to 34.2% of proteins [6]. The number of hits per gene was normalized to gene length (kb) to obtain reads per kilobase (RPK). KO counts were pooled at the genus or species level, and module counts were calculated for the genera using the module definitions without the use of a one-off strategy. Module counts were then further summed based on KEGG’s BRITE hierarchy to the pathways to which they belong.

2.5. Metagenome Assembled Genome (MAG) Assembly, Binning and Annotation

Metagenomic scaffolds were assembled from the quality-filtered reads using metaSPAdes (v. 4.0.0)[7]. Scaffolds with a length <1,500 bp were discarded. Assembled scaffolds were assigned a taxonomy using Kraken2 and the same NCBI nt database. Genomes were then binned using TaxVAMB (v. 4.1.4)[9], which incorporates taxonomic data with sequence characteristics for binning. Bin completeness and contamination were assessed using CheckM2 (v.1.0.2).

Bakta (v.1.11.0)[8] was used to annotate bacterial scaffolds, while metaeuk (7.bba0d80)[9] was used for the fungal ones. For annotating resistance and virulence genes, we used NCBI’s tool AMRFinderPlus (v. 4.0.23) (accessed 25.03.2025) [10]. PlasmidFinder (v.2.1.6)[11] was used to annotate plasmid origins of replication within the contigs, along with additional origins assigned by bakta.

3. Results

3.1. Sequencing

A total of 73,074,454 forward and reverse reads of 150 bp in length were sequenced on the Illumina platform from a sample of green cheese, with an average Phred score of 36. Reads were subjected to quality filtering and removal of the human genome before downstream analysis.

3.2. Taxonomic Profiling

78% of reads were classified to at least the domain level. Twenty bacterial species were confidently assigned using a minimal unique k-mer approach, as well as three fungal species. Species abundance is represented in Figure 1 with a threshold of >1% relative abundance. The dominant groups of bacteria are the genera Brevibacterium, followed by Lactobacillus, Streptococcus, Lactococcus and Brachybacterium. These were further stratified into the dominant species Brevibacterium aurantiacum and Brachybacterium sp. P6-10-X1, Lactobacillus delbrueckii subsp. bulgaricus, Lactiplantibacillus plantarum, and Streptococcus thermophilus, among other lactic acid bacteria. For the fungi, Aspergillus puulaauensis accounted for 95% of fungal reads, while the remaining reads were shared between A. versicolor and A. sydowii.

Before Bracken read re-distribution, Brevibacterium aurantiacum reads represented 46% of the genus, followed by B. sandarakinum, B. limosum, B. spongiae, B. BDJS002 and others, but all with <5% of reads. Additionally, a third of reads in this genus could not be identified to the species level. 93% of these genus-level reads successfully mapped back to the B. aurantiacum reference genome, meaning that they represent sequences shared between the member species of Brevibacterium and could not be confidently assigned by Kraken2.

A fourth of Lactobacillus delbrueckii reads were identified as the subspecies bulgaricus, and 10% as the subspecies lactis. Still, the latter did not meet the unique k-mer threshold to be kept as a valid assignment.

The only identified Brachybacterium species that wasn’t a false positive was Brachybacterium sp. P6-10-X1, but represented only 20% of Brachybacterium assignments.

In summary, the dominance of Brevibacterium aurantiacum is consistent with its known role in smear-ripened cheeses [12]. Brachybacterium has a significant role in carbohydrate and amino acid metabolism during cheese ripening, which is common for cheeses using raw milk [13,14]. Lactic acid bacteria are present in significant quantities even during ripening. Aspergillus puulaauensis appears to be a major player in lipid and amino acid metabolism in green cheese.

3.3. Functional Profiling

Considering that a significant number of reads could not be assigned to a species, especially for the dominant genus Brevibacterium, the genus level was chosen to represent the metabolic profile of the community as a whole, with only the dominant genera considered. At the genus level, a total of 165 KEGG modules were reconstructed across the chosen genera using the module definitions without a gap fill strategy. To better summarize the function of the community, the modules were collapsed to their higher brite taxonomy in Figure 2, with a minimal pathway count threshold of 25. To further interrogate amino acid metabolism, the separate modules of the amino acid metabolism pathways are also shown in Figure 3. Brevibacterium, Brachybacterium and Aspergillus are the genera that contributed most heavily to the metabolic potential of the community, and have all previously been implicated as playing a crucial role in cheese ripening [12,13,14,15,16,17].

3.4. Metagenome Assembly and Binning

MetaSPAdes was used for scaffold assembly, and TaxVAMB was used to group the scaffolds into bins by incorporating taxonomic data from Kraken2 and innate sequence features, with a minimum scaffold length of 1.5 kb. The features of the resulting scaffolds are summarized in the blob plot of Figure 4. The scaffolds with the highest coverage belonged to a species of fungal virus from the Genomoviridae family, which was excluded from the figure for clarity.

At a threshold of 50% completeness, as assigned by CheckM2, a total of 17 bins were recovered (Table 1). All of them had contamination <3%, except for the bin of Brevibacterium aurantiacum, for which the contamination was nearly 20%, indicating that the assembly is most likely a hybrid from different species of the genus, which again relates to the high level of reads only assigned to the genus level for that group.

3.5. Bacterial Resistance Annotation

Genes that relate to antimicrobial resistance, the taxonomy of the scaffold they’re located within, and the presence or absence of a point mutation are summarized in Table 2. While several resistance genes were identified, none of the potential pathogen species to which they belong were detected within confidence thresholds.

4. Discussion

4.1. Role of Actinomycetes in Green Cheese Ripening

Brevibacterium species play a crucial role in the ripening of surface-ripened cheeses, owing to their ability to break down lipids and proteins, and produce volatile sulfur compounds, which are essential for the cheese’s aroma profile. Additionally, they provide a typically red-orange pigment to the cheese surface [12]. The genus is also implicated in utilizing histamine as a carbon source [18]. The primary species of Brevibacterium involved in cheese ripening, and most often found on the cheese surface, are B. linens, B. aurantiacum, and B. antiquum [19]. The source of Brevibacterium in cheese is often from the use of raw milk, which serves as the inoculant for the starting culture of microorganisms in the cheese [13,14].

Similarly, our study identifies B. aurantiacum as the dominant species, but the presence of other species from the genus is difficult to discern. While taxonomic annotation reads show a relatively even distribution of a few reads between other species in the genus, only the genomes of B. limosum and B.spongiae could be binned with a completeness above 50%, which makes them the most likely candidates for other species present in the sample.

Since low-contamination and high-completeness genomes could not be acquired for every prevalent species, reads were used for functional annotation, rather than assembled scaffolds. To reduce ambiguity arising from the genus and species-level assignment of the reads by Kraken2, functional profiling was also performed at the genus level.

In our study, Brevibacterium seems to be the only group capable of metabolizing the present galactose. Galactose is produced in cheese as a result of the activity of Streptococcus thermophilus, and can thus remain in the ripening phase of cheese production [20]. Additionally, fatty acid degradation pathways were only detected in this genus, suggesting that remaining fats in the milk are metabolized mainly by these organisms.

Multiple pathways for secondary metabolism of cofactors and vitamins were assigned to Brevibacterium as well. These include lipoic acid biosynthesis, pyridoxalphosphate (vitamin B6) biosynthesis, pantothenate biosynthesis (vitamin B5), the thiamine salvage pathway, porphyrin (heme) biosynthesis, and others. Lipoic acid is crucial for anaerobic respiration in mammals and acts as an antioxidant [21]. It has also been proposed as a treatment for various human diseases due to its ability to quench free radicals [22]. Another beneficial secondary metabolite produced in high abundance by the microbial community is pyridoxal phosphate (vitamin B6), which is indispensable for numerous biological functions in the human body. Its deficiency can lead to diabetes, neuropathy, immune disorders, and other health complications [23]. Pantothenate (Vitamin B5) was also detected and is important for energy generation, hormone synthesis, and other regular bodily functions [24]. Taken together, we believe our study shows the abundance of these secondary metabolites could have a positive effect on digestion and overall health.

Brachybacterium is commonly found in cheese rings, and its origins can be attributed to the use of raw milk, in which different Brachybacterium species are identified [13,14]. In the current study, despite its lower abundance relative to Brevibacterium, we found that this genus exhibits many of the same metabolic functions as Brevibacterium during cheese ripening.

The unique roles Brachybacterium might fill in the community could be related to the pathways for taurine, histidine, and phenylalanine metabolism, as well as its very high activity in the pentose phosphate pathway. PPP is an alternative pathway for glycose that ultimately serves cell anabolism. Upon closer inspection, the pathway for the oxidative phase of the PPP is missing in the Brevibacteria KOs. This observation suggests that glycolysis is a significant source of NADH for the cell, which would explain the disparity in carbohydrate metabolism between the two genera. Our study additionally showed that porphyrin (more specifically, heme) biosynthesis was detected in Brevibacterium at high counts, implying that the genus is capable of efficient aerobic respiration.

We found that the Brevibacterium genus is a major contributor to the biosynthesis of amino acids within green cheese. This observation encompasses the biosynthesis of Pro, Arg, Met, Cys, Thr, Ser, Lys, Ile, and Val, as well as ornithine and ectoine. Interestingly, no pathways for the degradation of amino acids were detected in our study. Catabolism of amino acids during ripening is very significant to the flavor profile of cheese [28]. It has been shown that the potential for amino acid catabolism varies within the genus Brevibacterium [29]. It’s most likely that the pathways used by Brevibacterium aurantiacum for amino acid degradation include the synthesis of aromatic sulfur compounds which may be missing from the KEGG module definitions, and are thus not accounted for in our study.

Despite the lower abundance of Brachybacterium, the genus still appears to play a significant role in amino acid metabolism in our study, primarily through the synthesis of Pro, Arg, Cys, Thr, Ile, Ser, and Gly, as well as ornithine, ectoine, and betaine. The genus also produces Lys, which is lacking from the amino acid repertoire of Brevibacterium. The His degradation pathway was also identified in this group.

Ectoine is produced in significant quantities by both Brevibacterium and Brachybacterium and typically functions as an osmoprotectant for bacterial cells [30]. This process could be crucial for the survival of actinomycetes in the highly saline cheese environment, since it’s not a trend observed between the lactic acid bacteria.

Another amino acid that was produced in great quantity in this study by actinomycetes was ornithine. Ornithine is found to be present in cheese due to the activity of lactic acid bacteria; however, we demonstrate that other bacteria involved in ripening also make a significant contribution [31]. Positive health effects have been shown for the consumption of ornithine [32]. Ornithine can also be transformed into putrescine and other polyamines by bacteria, which could also impact the flavor profile of the cheese [33].

4.2. Role of LAB in Green Cheese Ripening

Lactic acid bacteria in the sample were represented by Lactobacillus delbrueckii, Streptococcus thermophilus, Lactiplantibacillus plantarum, Streptococcus thermophilus and Lactococcus lactis. The bacteria are responsible for metabolizing lactose in milk into lactic acid during the acidification step of cheese production [25]. This process can continue into the ripening phase of the cheese. Taken together, these species represent a significant portion of the microbial community in our study, accounting for 40% of the species abundance. Even though the counts of lactic acid bacteria are high, their contributions to the metabolic functions of the community during the ripening stage appear low in our study, with only Lactobacillus having a high number of pathways for cys and CoA biosynthesis.

4.3. Role of Mold in Green Cheese Ripening

One of the dominant microorganisms on the cheese was the filamentous fungi A. puulaauensis, with 16% of classified reads and 3% of the organism abundance. This mold is shown to grow on a variety of milk and goat cheeses, including Grana and some artisan Italian cheeses [15,16,17]. It’s also shown to produce the mycotoxin sterigmatocystin, but its concentration on cheeses is low, and does not pose a health risk [26,27].

In this study, we found that fatty acid biosynthesis is entirely carried out by Aspergillus puulaauensis, which could contribute to the final flavor profile of the cheese. Additionally, the species plays a significant role in carbohydrate and amino acid metabolism, specifically the degradation of proline.

Interestingly, the dominant organism in this study by sheer abundance was a virus from the family Genomoviridae. This family contains viruses that infect fungi, and it is certainly plausible that one of these viruses infects the species of Aspergillus identified in this sample.

While no typical pathways for amino acid degradation were detected in bacteria, the filamentous fungi Aspergillus puulaauensis was shown to have pathways for the degradation of Pro, GABA, and Met, as well as pathways for the synthesis of putrescin from ornithine. Our findings could mean A. puulaauensis can metabolize compounds produced by the actinomycetes community, rather than directly participating in the proteolysis of casein and other peptides in the cheese.

4.4. Antibiotic Resistance and Virulence Genes in Green Cheese

Upon assembling the sequencing reads, several scaffolds were identified as containing genes that encode resistance to antibiotics, including kanamycin, tetracycline, streptomycin, and others. The resistances to tetracycline and streptomycin, in particular, were found on the same contigs that contain Staphylococcus aureus plasmid origins of replication, strongly implicating them as part of a plasmid. Reads identified as Staphylococcus equorum represented approximately 0.02% of all reads, which is too few to consider the species present in the sample. We speculate this represents either environmental DNA, contamination during handling, or the plasmid being carried by another host. Similarly, Enterococcus faecium was found with only 0.02% of reads.

Antimicrobial resistance has been identified in several species of lactic acid bacteria in cheese, including some of the resistance genes identified in this study [34,35]. E. faecium resistance is also identified in French cheeses [36].

4.5. Comparison Between Taxonomic Profile Results from Shotgun Sequencing and Amplicon Sequencing

A previous study on traditional Bulgarian green cheese employed amplicon metagenomic sequencing using several targeted regions of both the 16S and ITS marker genes [2]. While Brevibacteria was present in most samples, it was only dominant in one of them and could not be classified to the species level using any 16S region. The current shotgun approach elucidates the presence of B. aurantiacum with certainty, and most likely additional species B. limosum, B. spongiae, and B. sandarakinum. We speculate that the genetic distance between the species of this genus is small, which makes the unambiguous assignment of the many sequences shared between the species difficult.

Discrepancies between the taxonomic profiles of the two studies might be attributed to sampling location on the cheese, the length of the maturation period, and the particular 16S region chosen for sequencing, since different regions are known to differ in the outcome of taxonomic annotation [37].

One striking difference between the amplicon sequencing of the cheese is the dominant fungi, which was Debaryomyces hansenii in most samples of the previous study. Considering that the fungi grow on the cheese rind, it is possible that multiple fungi inhabit the cheese surface, and these differences are sample-dependent. This suggestion also lends support to the idea that multiple sub-communities of organisms could form their own distinct clusters within the cheese, as well as on the rind. An important limitation of the current study is the number of samples, which makes it difficult to draw solid conclusions about all microorganisms that contribute to the sensory properties of traditional Bulgarian green cheese.

5. Conclusions

Our study found that the microbial community of traditional Bulgarian green cheese is represented by three main groups of microorganisms – actinomycetes, LAB and mold, all of which most likely contribute to the sensory properties of the food product. However, it’s important to note that the single sample used in this study cannot capture all sources of variability, such as the heterogeneity of the surface microbial community or the variability within cheese batches, which requires a more robust sampling approach. We found the presence of various pathways for the synthesis of useful secondary metabolites, which could point to potential benefits for consumption of the cheese. However, it’s important to note that the presence of these pathways does not assure the presence of their end products, and in this reguard, our study provides a starting point for potential biochemical quantification for the products of these pathways. We were also able to confidently rule out the presence of any harmful microorganisms within the cheese, yet found traces of plasmids carrying microbial resistance genes, the source of which remains speculative. We were also able to provide a potential framework for the interactions of microorganisms groups within the cheese, and how those interactions shape the quality of the end product during ripening. We believe this study helps to elucidate the microbial structure and role of the dominant microorganisms within the surface of traditional Bulgarian green cheese.