Gouda cheese production is a complex interplay of a mixed-strain starter culture, non-starter microorganisms, and time. A typical Gouda cheese starter culture, used in lyophilized or frozen form, contains (undefined) strains of the lactic acid bacteria (LAB) Lactococcus lactis subsp. lactis (Lc. lactis), Lactococcus cremoris subsp. cremoris (Lc. cremoris), Lactococcus lactis subsp. lactis biovar diacetylactis (Lc. lactis biovar diacetylactis), and Leuconostoc mesenteroides and/or Leuconostoc lactis (Düsterhöft et al., 2017). The main role of this mesophilic starter culture is to lower the pH of the fermenting milk by converting lactose into lactate, which is almost completed after 10 h. Besides acidification by the homofermentative Lactococcus, citrate utilization by Lc. lactis biovar diacetylactis and/or leuconostocs occurs, which results into buttery flavor compounds, such as diacetyl, acetoin, and/or 2,3-butanediol (Starrenburg and Hugenholtz, 1991; Smit et al., 2005; Steele et al., 2013; Düsterhöft et al., 2017). In addition, all strains contribute to flavor development through proteolysis, peptidolysis, and amino acid conversions. The heterofermentative leuconostocs produce few eyes because of carbon dioxide formation through heterolactic fermentation and citrate conversion. These subprocesses of fermentation and ripening determine the complex and dynamic flavor profile of Gouda cheese, which is characterized by a sweet and mild taste.

The application of the same LAB starter cultures in a Gouda cheese production environment increases the risk of a bacteriophage infection, which may lead to failed fermentation processes and thus economical losses (Mahony and van Sinderen, 2014). Starter cultures consisting of a complex mixture of strains are in general less phage-sensitive, but a starter culture with only few strains guarantees a more consistent end-product quality (Erkus et al., 2013; Spus et al., 2015). Therefore, most cheese factories apply a starter culture rotation with three to five starter cultures with different phage sensitivity profiles (Høier et al., 2010; Parente et al., 2017; Frantzen et al., 2018). However, little is known about the impact of the use of rotating starter cultures on the final organoleptic properties of Gouda cheese throughout various batch productions.

Although mixed-strain LAB starter cultures are used to steer the Gouda fermentation and ripening process, mesophilic non-starter lactic acid bacteria (NSLAB) originate from the milk or the cheese production equipment or environment (e.g., brine bath and ripening room) and are able to grow to counts as high as 10⁸ CFU/g during the ripening phase (Fitzsimons et al., 2001; Van Hoorde et al., 2008; Fox et al., 2017; Gobbetti et al., 2018). In general, Lacticaseibacillus paracasei and Lactiplantibacillus plantarum are the main NSLAB across various cheese types ripened for 4 months or longer, including Gouda cheese (Gobbetti et al., 2015). Positive contributions of NSLAB to cheese flavor are well known, as they hydrolyze peptides and convert the resulting free amino acids into flavor compounds, which however depends on the production batch (Montel et al., 2014; Gobbetti et al., 2015). Moreover, in some cases, they cause off-flavors or undesired gas production (Martley and Crow, 1993; Oberg et al., 2016). To take benefit from the NSLAB regarding flavor formation, but to avoid batch variability, selected NSLAB are added as adjunct cultures in some cheese productions (Antonsson et al., 2003; Van Hoorde et al., 2010b; O’Brien et al., 2017; Randazzo et al., 2021). Additionally, Enterococcus durans, Enterococcus faecalis, Enterococcus faecium, Lacticaseibacillus rhamnosus, Latilactobacillus curvatus, Lentilactobacillus parabuchneri, and Pediococcus pentosaceus may be part of the Gouda cheese NSLAB communities (Van Hoorde et al., 2008, 2010a). Also, culture-dependent methods detect a higher relative abundance of NSLAB in aged Gouda cheese compared to culture-independent approaches, because the starter culture strains occur in a viable but not culturable (VBNC) state during ripening that can only be detected through DNA sequencing (Ganesan et al., 2007; Ruggirello et al., 2014). Finally, given the reclassification of the Lactobacillus genus (Zheng et al., 2020), an update is needed to assess the full diversity of lactobacilli in Gouda cheese.

Whereas culture-dependent methods are more laborious and only isolate microorganisms that can grow on the agar media and under the incubation conditions applied, resulting into an inherent selection bias, culture-independent methods are more suited to investigate the taxonomical diversity of a microbial ecosystem (Ercolini, 2013; De Filippis et al., 2017; Levante et al., 2021). Most metagenetic high-throughput amplicon-based sequencing studies on cheese target only a part of the 16S rRNA gene, as has been performed, for instance, for Irish artisan cheeses (Quigley et al., 2012), soft Herve cheese (Delcenserie et al., 2014), Portuguese Pico cheese (Riquelme et al., 2015), traditional Brazilian cheeses (Kamimura et al., 2019), Italian mozzarella (Marino et al., 2019), Gouda cheeses (Park et al., 2019), and Grana Padano cheeses (Zago et al., 2021). Alternatively, the microbiota of some cheeses has been investigated by high-throughput sequencing of the whole 16S rRNA gene, using PacBio technology, which allows a more accurate species-level identification (Jin et al., 2018; Yang et al., 2020). The latter study has also attempted to link the most prevailing species with sensory profiles generated with an electric nose (Yang et al., 2020). However, since the metagenetic sequence reads were clustered into operational taxonomic units (OTUs) instead of amplicon sequence variants (ASVs), taxonomic assignment could not be performed below species level (Callahan et al., 2017). Indeed, ASVs allow a single-nucleotide resolution that can resolve taxonomical identification at species and sometimes even strain level (Callahan et al., 2019; Johnson et al., 2019). In that perspective, a meta-analysis of cheese microbiomes has indicated a significantly different volatile profile for strains of the same species, encompassing Brevibacterium linens, Lc. lactis, and Streptococcus thermophilus involved in artisan Irish cheese production (Walsh et al., 2020), which emphasizes the relevance of obtaining a reliable taxonomical resolution below species level.

Besides an accurate taxonomic assignment of the microbial dynamics during Gouda cheese ripening, a thorough metabolomic analysis as a function of the ripening time is needed, in particular to be able to make the link between the microbial taxa and flavor compounds produced unequivocally. For example, research on smear-ripened cheddar cheeses has correlated both Debaryomyces hansenii and Glutamicibacter arilaitensis with alcohols and carboxylic acids, and the latter species with ketones as well, whereas B. linens and Geotrichum candidum have been correlated with sulfur compounds (Bertuzzi et al., 2018). An in-depth study of Gouda cheese has investigated the sensory and chemical properties of 36 productions from five countries and indicated the age of the cheeses as most determining factor to explain differences in concentrations of flavor compounds formed during ripening (Jo et al., 2018). Lactic acid has been identified as the main organic acid in Gouda cheeses. Acetic acid, diacetyl, 2-and 3-methylbutanal, and ethyl butyrate are among the most important aroma-active compounds. This is in line with other studies that have indicated acetoin (buttery notes), lactones (fruity/sweet/buttery notes), methylketones (floral/nutty notes), and short-chain fatty acids (cheesy/rancid notes) as additional major constituents of Gouda cheeses (Dirinck and De Winne, 1999; Van Leuven et al., 2008; Van Hoorde et al., 2010b; Jung et al., 2013; Ruyssen et al., 2013). Sulfur compounds, such as dimethyl sulfides and methional, are present at relatively low concentrations, but do have an impact on the Gouda cheese flavor (Van Leuven et al., 2008; Jo et al., 2018). In addition to their contribution to the basic cheese taste (Toelstede and Hofmann, 2008), amino acids are precursors for volatile organic compounds that contribute to the cheese aroma, such as higher alcohols, higher aldehydes, and esters (Smit et al., 2005). However, they can also be the precursors of biogenic amines, which have a toxicological potential, especially histamine and tyramine (Komprda et al., 2007; Linares et al., 2012). All these compounds can be produced by a wide range of LAB species, which is considered to be a strain-dependent trait (Barbieri et al., 2019; Walsh et al., 2020). However, few studies have focused on potential differences in microbial and metabolite compositions of the core and rind of the cheeses investigated (Calasso et al., 2016; Dugat-Bony et al., 2016; Salazar et al., 2018).

Whereas several studies on Gouda cheese have focused on its microbial composition with limited focus on the metabolites produced, other studies have investigated a restricted part of the metabolome with limited attention for the microbial composition. Therefore, the present study aimed at a systematic assessment of the overall microbiome (below species level) and metabolome (including endogenous and microbial substrates and metabolites, both volatile and non-volatile ones) as well as their link with the application of rotated starter cultures and/or cheese cores and rinds of Gouda cheeses from 23 different production batches, throughout their whole ripening time and made in the same dairy factory.

The present study deals with the investigation of batch-to-batch variations of a Gouda-type cheese produced on large scale in a European dairy factory, with special attention for the rotated starter culture mixtures applied, the ripening time, and the microbial and metabolite compositions of the cores and rinds of the cheese wheels. Although the commercial supplier of the mixed-strain starter cultures guaranteed that the starter culture mixtures A, B, and C used in the present study lead to similar end-products, significant differences were found between the Gouda A, B, and C cheeses produced, not only regarding their microbiome but also their metabolome. The main acidifiers, Lc. lactis and Lc. cremoris, did occur in all cheeses, albeit not in the same relative abundances, but species of the heterofermenter Leuconostoc did not, at least not to the same extent or according to the starter culture mixture used. The data of the present study further indicated that the dairy factory involved harbored a house microbiota that was mainly composed of Lacc. paracasei, T. halophilus, and Loil. rennini.

For the bacterial metagenetic analysis, a high-throughput, full-length 16S rRNA gene sequencing and ASV analysis was applied, which allowed accurate identifications at both species and strain level (Callahan et al., 2019). The big advantage of such a metagenetics approach lays in the single-nucleotide resolution, which even resulted in the detection of highly correlated, non-identical ASVs, probably belonging to one strain or cluster of strains. However, this deeper resolution also led to new challenges. For example, clusters with shared ASVs made further analysis less straightforward. Additionally, regardless of the precision, the use of only one gene might give a biased view on the bacterial diversity regarding subspecies level.

Most differences in microbiome and metabolome among the Gouda cheeses produced with the three different mixed-strain starter cultures applied were found after 36 weeks of ripening, as illustrated by the bacterial beta diversity and the scaled PCA of metabolites that grouped these cheese production batches according to the starter culture mixture used. At later ripening times, differences between the cheese production batches made with the same starter culture mixture increased and became more important than differences between the three groups of Gouda cheeses produced. This corroborated with a higher relative abundance of NSLAB that had a low diversity on species level across batches, but a striking diversity on ASV level. A similar situation has been described for a cheddar cheese factory, where cheese batches produced in different vats during the same day, or in the same vat during different days, have less than half of the NSLAB strains encountered as identical (Williams et al., 2002). The existence of a house microbiota and its role in the propagation of NSLAB in the factory and cheeses produced has been described frequently (Mounier et al., 2006; Bokulich and Mills, 2013a; Delcenserie et al., 2014; Calasso et al., 2016). Indeed, cheeses produced later during the day tend to have a higher bacterial diversity, suggesting a gradual accumulation of NSLAB in the production facility (O’Sullivan et al., 2015). This house microbiota was probably a main source of NSLAB contamination and associated variability for the 23 Gouda cheeses investigated.

The acidifying lactococci, Lc. cremoris and Lc. lactis, were the most abundant bacterial species in the Gouda cheeses produced, the former LAB species being present in higher relative abundances in the starter culture mixtures applied. The typical differences regarding salt tolerance and arginine conversion between these homofermentative LAB species were not always confirmed. For example, all Gouda cheese production batches had the same ASV profiles in the cores and rinds after 36 weeks of ripening, although it could be expected that the more salty rind would contain more ASVs of the more salt-tolerant Lc. lactis. Additionally, high arginine concentrations were not correlated with Lc. cremoris, although this species can not metabolize arginine (Schleifer et al., 1985). Yet, the higher arginine concentrations of the Gouda B cheeses could be linked with their slightly higher relative abundances of Lc. cremoris. Alternatively, in most cheese production batches, the relative abundance of Lc. cremoris decreased in favor of Lc. lactis, in particular between 45 and 100 weeks of ripening, which has to be ascribed to the faster speed of lysis of Lc. cremoris than Lc. lactis (Erkus et al., 2013; Fox et al., 2017).

Two ASVs of Lc. cremoris, Lc. cremoris_01 and Lc. cremoris_06, showed a high correlation with acetoin. However, acetoin is normally produced from citrate by Lc. lactis biovar diacetylactis, likely suggesting that the acetoin-producing strains lysed before 26 weeks of ripening. Alternatively, citrate-metabolizing and non-metabolizing Lc. lactis strains could share the same ASVs, which made it impossible to find meaningful correlations. Possibly, a substantial fraction of the Lactococcus reads was from cells in a VBNC state, as serine was one of the few amino acids positively correlated with Lactococcus and serine and methionine are particularly produced when Lactococcus enters a VBNC state (Ganesan et al., 2007).

In general, the different Lactococcus ASVs identified in the present study allowed to link each cheese production batch with its starter culture mixture applied, which can be of use for food authentication, although the differences in ASVs could not explain metabolic differences between the Gouda cheese production batches. However, it is well known that the Lactococcus genotype is often not in line with its phenotype (Tailliez et al., 1998; Fernández et al., 2011; Cavanagh et al., 2015), which has been reflected in the different reclassifications of Lactococcus species over time (Orla-Jensen, 1919; Garvie and Farrow, 1982; Schleifer et al., 1985; Li et al., 2021).

Considering the non-lactococci in the starter cultures used, a striking difference was the very low relative abundance of leuconostocs in the Gouda B cheeses, although they were present in the original starter culture mixtures applied. The Gouda A cheeses had a higher and approximately the same level of leuconostocs as the Gouda C cheeses. Yet, the presence of the Leuc. pseudomesenteroides_01 and Leuc. pseudomesenteroides_02 ASVs in the starter culture mixtures was lower than 0.008% (based on the total number of sequence reads). In general, undefined Gouda cheese starter cultures contain a low share of leuconostocs (Smid et al., 2014; Porcellato and Skeie, 2016; Frantzen et al., 2017). Overall, Leuc. mesenteroides was more abundant in the starter culture mixtures used, but it was marginally present in the Gouda cheeses produced compared to Leuc. pseudomesenteroides. In older literature and product information sheets, Leuc. mesenteroides subsp. cremoris and Leuconostoc lactis are considered as the main (sub)species in Gouda cheeses and mesophilic Gouda cheese starter cultures, but during the last 15 years, isolation of Leuc. mesenteroides subsp. mesenteroides, Leuc. mesenteroides subsp. dextranicum, and Leuc. pseudomesenteroides is more common (Frantzen et al., 2017). Of course, some Leuc. pseudomesenteroides isolates have also been wrongly identified as Leuc. mesenteroides, due to the use of solely partial 16S rRNA gene sequencing and/or phenotypical tests (Frantzen et al., 2017). Apart from a different effect on cheese metabolites, Leuc. pseudomesenteroides grows better and has a larger genome size compared to Leuc. mesenteroides subsp. cremoris, suggesting a higher survival potential during cheese ripening (Pedersen et al., 2016; Frantzen et al., 2017).

The three different ASV clusters of Leuc. pseudomesenteroides found in the Gouda cheeses of the present study correlated with either acidification and the lack of an efficient citrate metabolism (ASV clusters 1 and 3) or amino acid expulsion and flavor formation (ASV cluster 2). However, the lactococci associated with each Leuconostoc ASV cluster could be responsible for these differences as well. Indeed, leuconostocs are not believed to directly contribute to lactose fermentation (Ardö and Varming, 2010). Also, they lack cell envelope proteinases, making them dependent on peptide release from lactococci, and although they have only limited aminopeptidolytic activity, they can increase the concentrations of some amino acids extracellularly (Pedersen et al., 2016). Further, whereas the concentrations of most amino acids were lower in the Gouda B cheeses compared to the other ones, reflecting the lower relative abundance of leuconostocs, these cheeses contained the highest acetoin and 2,3-butanedione concentrations. This could be ascribed to the citrate-metabolizing Lc. lactis strains abundant or present in those cheeses, questioning the need of Leuconostoc for the production of these aroma compounds. Indeed, Lc. lactis biovar diacetylactis can even produce more acetoin compared to Leuconostoc (Starrenburg and Hugenholtz, 1991; Hugenholtz, 1993). However, the Gouda B cheeses with a low relative abundance of Leuconostoc contained more NSLAB, suggesting that inhibition of the growth of undesired bacteria seems to be an interesting contribution of Leuconostoc, an effect that has been described before (Martley and Crow, 1993). This effect was stronger for strains of the ASV cluster 1 in the Gouda C cheeses, compared to those of ASV cluster 2 in the Gouda A cheeses.

The increase of the relative abundances of the NSLAB species from 5 to 20% between 26 and 45 weeks of Gouda cheese ripening was in line with other culture-independent Gouda cheese studies (Porcellato and Skeie, 2016; Park et al., 2019). Lacticaseibacillus paracasei was the most abundant NSLAB species during the first 45 weeks of ripening, confirming previous reports (Van Hoorde et al., 2008; Kołakowski et al., 2012). It was represented by 10 of ASVs, differing between the Gouda cheese production batches and ripening times, but often occurring in a similar ratio, suggesting the presence of non-identical copies of the 16S rRNA gene within one strain. Possibly, the cheese manufacturing environment contained a high diversity of Lacc. paracasei and each production batch got inoculated with some strains that prevailed during ripening and contributed to the release of amino acids and the production of volatile organic compounds. In support of these beneficial properties, the use of Lacc. paracasei as adjunct culture for cheese production is well-known (Antonsson et al., 2003; Van Hoorde et al., 2010b; Stefanovic et al., 2018; Asahina et al., 2020; Bancalari et al., 2020; Randazzo et al., 2021).

Tetragenococcus halophilus was the most abundant NSLAB species in multiple Gouda cheeses after 75 and 100 weeks of ripening. This LAB species has been found culture-independently in Polish Oscypek cheese (Alegría et al., 2012), Spanish Cabrales and Manchego-type cheeses (Diaz et al., 2016; Botello-Morte et al., 2022), the rind of Gouda cheese (Salazar et al., 2018), and Irish hard cheese (only genus level; Quigley et al., 2012), but always at low relative abundances. Culture-dependently, it has been isolated from Mexican Cotija cheese (Morales et al., 2011), Brie de Meaux cheese (Unno et al., 2020), and Spanish blue-veined cheeses (Rodríguez et al., 2022). In these culture-dependent studies, T. halophilus has only been picked from media enriched with 2.5, 5.0, or 7.0% (m/v) NaCl. Besides a halophilic character, this LAB species displays an alkaliphilic behavior, as growth is better at pH 8.0 than at pH 7.0, whereas no growth occurs at pH 5.0 (Justé et al., 2012; Unno et al., 2020). Except for its salt dependency in agar media, its rare isolation could be ascribed to its prevalence in cheese rinds or long-ripened cheeses, which give lower counts with a culture-dependent analysis and fewer cells for proper DNA extraction with a culture-independent analysis.

The presence of T. halophilus in cheeses probably originates from the brine bath, as its occurrence in the latter was shown both culture-dependently and culture-independently (data not shown). Indeed, it has been found as a prevalent species in Belgian, Italian, and Danish brines (Marino et al., 2017; Haastrup et al., 2018; Vermote et al., 2018). One Italian brine shows a relative abundance of more than 50% (on sequence read basis) of Tetragenococcus and has a pH of 4.6, whereas the pH of the other brines varies between 4.8 and 5.7, suggesting that either growth is still possible at a pH lower than 5.0 or that no growth but survival occurs thanks to its salt tolerance. Since T. halophilus was already abundant in the Gouda cheese rinds of the present study after 36 weeks of ripening, it could likely proliferate and migrate from the brine to the cheese rind and finally to the core upon further cheese ripening. Furthermore, strong correlations were found between T. halophilus and the concentrations of free amino acids and short-chain fatty acids, which was in line with the proteolytic and lipolytic properties of isolates from Mexican cheeses (Morales et al., 2011). These properties could make T. halophilus interesting as an adjunct starter culture for (Gouda) cheese production. This LAB species is already used as primary starter culture for the production of oriental fermented salted food products, such as soy sauce, fish sauce, and shrimp paste (Justé et al., 2014). However, some T. halophilus strains harbor histidine and/or tyrosine decarboxylase genes, allowing histamine and/or tyramine production, respectively (Diaz et al., 2016; Chun et al., 2019). Alternatively, well-chosen strains would even lower biogenic amine formation, probably indirectly through the repression of other tetragenococci that produce biogenic amines (Udomsil et al., 2011; Kim et al., 2019), which supports the need to perform a genome sequencing and functional analysis approach of potential adjunct culture strains before their application.

Although T. halophilus was represented by more than 100 different ASVs in the Gouda cheeses examined, only two main ASV clusters occurred after 75 weeks of ripening, which confirmed its strain-dependent properties. Whereas the T. halophilus ASV cluster 1 was mainly found in cheeses with higher pH (which may have been caused by biogenic amine formation) and higher concentrations of short-chain fatty acids, the T. halophilus ASV cluster 2 was mainly found in cheeses with higher dry mass, reflecting the acid-sensitive and lipolytic character of strains of the former cluster and the osmotolerance of those of the latter one. In addition, the T. halophilus ASV cluster 1 strongly correlated with Loil. rennini and high and low concentrations of several biogenic amines and their precursor amino acids, respectively, after 75 weeks of ripening. The T. halophilus ASV cluster 2 was negatively correlated with histidine and positively correlated with histamine, suggesting the presence of a histidine decarboxylase gene that was expressed after 75 and 100 weeks of ripening. Indeed, histidine decarboxylase genes can be present in T. halophilus genomes, whereas other decarboxylase genes have not been found in the pangenome composed of 15 T. halophilus strains (Satomi et al., 2011; Chun et al., 2019).

Putrescine and cadaverine negatively impacted the sensory scores of the Gouda cheeses examined, whereas histamine did not. However, consumption of foods with a high histamine concentration can result in histamine poisoning (Özogul and Özogul, 2019). In the Gouda cheeses of the present study, the histamine concentrations were low up to 45 weeks of ripening, but they were around the Codex Alimentarius limit of 200 mg/kg (wet mass) established for fish and fish products after 75 and 100 weeks of ripening. Hence, because of a correlation of both ASV clusters 1 and 2 with the presence of several biogenic amines, a careful selection should be made when strains of these clusters are considered to be used as adjunct cultures. Indeed, in the case of the T. halophilus ASV cluster 1, the presence of biogenic amines could be related to the presence of Loil. rennini, and in the case of the T. halophilus ASV cluster 2, histamine production was most likely related to T. halophilus. However, the use of T. halophilus as adjunct culture could be made possible after curing of the plasmid encoding histamine production (Satomi et al., 2011; Botello-Morte et al., 2022).

The salt-tolerant Loil. rennini, which was prevalent in the Gouda cheese rinds after 36 weeks of ripening and in the cores of some cheese production batches after 75 and 100 weeks of ripening, was another prevalent NSLAB species. It was first isolated from rennet and held responsible for cracks and off-flavor defects in Dutch cheeses (Chenoll et al., 2006). It has also been found in Greek Kopanisti and Mana cheeses and in cheese brine (Asteri et al., 2009; Vermote et al., 2018). This LAB species, represented by a limited number of ASVs in the Gouda cheeses of the present study, correlated positively with different biogenic amines. Indeed, our investigation of the genome of a strain isolated from a 2-year Greek Kopanisti cheese (Kazou et al., 2017) revealed the presence of different decarboxylase genes, such as those coding for glutamate decarboxylase and ornithine decarboxylase, which yield 4-aminobutyric acid and putrescine, respectively. Moreover, a high correlation was found between Loil. rennini and a high pH after 75 weeks of ripening, likely to be ascribed to these decarboxylations.

The fungi found in the Gouda cheeses of the present study were in accordance with other reports on Debaryomyces as the most prevailing yeast in Gouda cheeses (Welthagen and Viljoen, 1998; Banjara et al., 2015). Whereas the latter studies did not differentiate between cheese cores and rinds, the present study showed that Debaryomyces was the most abundant yeast in the rind and Saccharomyces in the core of the Gouda cheeses produced. Debaryomyces has been found as a prevalent yeast in cheese brines and cheese production environments before (Mounier et al., 2006; Bokulich and Mills, 2013a; Vermote et al., 2018). Its higher salt tolerance and oxygen dependency than that of Saccharomyces (Sánchez et al., 2006) was reflected by the higher relative abundance of Debaryomyces in the Gouda cheese rinds and of Saccharomyces in the anaerobic cheese cores. The occurrence of yeasts in the Gouda cheeses of the present study may be linked with the production of different types of cheeses in the same dairy factory and their presence in cheese brines (Vermote et al., 2018). However, whereas some studies have mentioned yeast counts in Gouda cheeses up to 10⁵ CFU/g (Welthagen and Viljoen, 1998; Kołakowski et al., 2012), plating of some cheese samples of the present study on yeast-peptone-dextrose agar did not result in growth (data not shown), indicating that yeasts were present in very low numbers.

Batch-to-batch variation in mature Gouda cheeses was significant on microbial and metabolite level. Even when the rotated starter culture mixtures (A, B, and C) were considered to be similar, differences in microbial compositions and metabolite concentrations were found. A metagenetic approach based on high-throughput full-length 16S rRNA gene sequencing, accompanied with an ASV analysis, allowed differentiation of the bacterial communities on intra-species level. Starter culture mixture B did not result in prevailing leuconostocs in the corresponding cheeses, which led to a higher relative abundance of NSLAB. Hence, the presence of leuconostocs seemed to inhibit the growth of NSLAB, an effect that especially occurred in the Gouda C cheeses. Batch-to-batch differences generally increased with ripening time and became then less related to the starter culture mixture applied. Yet, metabolite differences between cores and rinds of the Gouda cheeses were found at all ripening times. The bacterial composition of the rinds at 36 weeks of ripening resembled those of the cores at 75 and 100 weeks of ripening. Among the NSLAB, which originate from the house microbiota that might vary, but will most probably be independent of the season given the constant temperature under which the factory is operating, Lacc. paracasei was most abundant until 45 weeks of ripening. This was taken over by T. halophilus and to a lesser extent Loil. rennini, which both could be responsible for biogenic amine formation and a bad flavor. Strong correlations of T. halophilus with desirable metabolites regarding aroma formation may make it a candidate adjunct culture. Finally, to increase Gouda cheese quality standardization, the use of only one starter culture mixture is proposed, together with a potential adjunct culture. Starter culture mixture rotation could still be applied, which each time should lead to another cheese type; for instance, starter culture mixture A should only be used for young Gouda cheeses, B for matured Gouda cheeses, and C for extra-matured Gouda cheeses. Also, the metagenetic ASV approach applied in the present study will allow differentiation of mixed-strain starter cultures across Gouda cheese producers in particular and contribute to food authentication in general.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.ebi.ac.uk/ena, PRJEB58546.

HD carried out the research and drafted the manuscript. SW and LDV supervised the research and revised and edited the manuscript. LDV was responsible for funding. All authors contributed to the article and approved the submitted version.

This work was supported by the Research Council of the Vrije Universiteit Brussel (SRP7 and IOF3017 projects) and the Agency Flanders Innovation and Entrepreneurship (VLAIO).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.