Omics-Based Insights into Flavor Development and Microbial Succession within Surface-Ripened Cheese

Two smear-culture mixes, D4 and S5, were used for the cheese trials and contained, as outlined in the supplier specification sheet, Brevibacterium linens, Debaryomyces hansenii, Cyberlindnera jadinii, and Brevibacterium casei (for D4) or Staphylococcus xylosus, B. linens, D. hansenii, Geotrichum candidum, and Glutamicibacter arilaitensis (previously classified as Arthrobacter arilaitensis) (for D5). Using metagenomic analysis, performed with Kaiju (11), the relative abundances of the individual species within the mixes were determined (Fig. 1). Overall, Kaiju was able to assign 81.7% ± 1.5% of reads from the starter mix samples to the species level. The proportion of assigned reads for each starter mixture sample is presented in Fig. S1 in the supplemental material. B. casei (60.83%) and C. jadinii (15%) were the most abundant bacterial and yeast species in D4, while B. linens and D. hansenii were minor components in the smear-culture mix, with relative abundances of 5.25% and 1.92%, respectively (Fig. 1; Table S1). In the S5 mix, G. arilaitensis (64.25%), D. hansenii (14.56%), and G. candidum (11.83%) were the most abundant bacteria and yeasts; S. xylosus (0.59%) and B. linens (3.52%) were present at lower relative abundances. Other species, not specified by the suppliers, were identified at low relative abundances in the smear-culture mixes D4 and S5 and are reported in Table S1.

Two test cheeses, D4 and S5, were prepared by smearing young Cheddar cheese curd with the two aforementioned commercial smear-culture mixes and ripened for 30 days at 15°C. Kaiju was used to determine the bacterial and yeast compositions of the cheese surfaces at days 0, 18, 24, and 30 for both the control cheese (unsmeared and ripened under vacuum) and the two test cheeses (11). Overall, Kaiju was able to assign 57.5% ± 8.3% of reads from the cheese samples to the species level. The proportions of assigned reads for each cheese sample are presented in Fig. S2. Compositional data of the cheese surface were analyzed by a one-way analysis of variance (ANOVA), designed with SAS 9.3, to determine the significant differences in the proportions of the individual species present over time (12). The metagenomic sequences of the bacteria used as starter cultures in the Cheddar cheese curd (Lactococcus lactis and Streptococcus thermophilus) and as smearing cultures (B. linens, S. xylosus, and G. arilaitensis) were compared at the strain level, using PanPhlAn, to determine the presence/absence of the inoculated bacterial strains on the cheeses throughout ripening (13).

As expected, lactic acid bacteria dominated the surfaces of all samples at day 0, and their relative abundances on the surface of the control cheese did not significantly change throughout the 30 days of ripening (Fig. 2). L. lactis and S. thermophilus were identified in all samples analyzed (D4, S5, and control) (Fig. 3). L. lactis was the dominant species in the control, constituting 75.85% of the initial population at day 0 and decreasing to 65.99% at day 30. S. thermophilus increased from 19.65% at day 0 to 28.21% at day 30, while the relative abundance of Lactobacillus (Lb.) helveticus was low throughout the ripening period (2.12% at day 0 and 2.72% at day 30) (Table S2). However, over the course of 30 days of ripening, the smearing processes clearly influenced the microbial populations of the cheese surfaces of both test cheeses, D4 and S5, causing a significant reduction in the relative abundances of Lb. helveticus (P < 0.03) and L. lactis (P < 0.0001). From day 0 to day 18, the population on the surface of D4 changed from predominately lactic acid bacteria to Debaryomyces hansenii and Glutamicibacter arilaitensis (Fig. 2). Subsequently, over the course of ripening, the relative abundance of D. hansenii significantly decreased (P < 0.0001) from 34.12% at day 18 to 4.14% at day 30 (Table S2). In parallel, the relative abundance of G. arilaitensis significantly increased (P < 0.0001) from 30.9% at day 18 to become the dominant population on the cheese surface (73.75%) at day 30 (Table S2). Using PanPhlAn, it was determined that the strain of G. arilaitensis detected on the cheese surface of D4 was different from the G. arilaitensis strain used in the smear-culture mix inoculated onto the surface of S5, confirming that the growth of this strain on D4 did not result from cross contamination of the two cheeses during inoculation or ripening (Fig. S3 and S4). However, the G. arilaitensis strain detected on D4 did appear to be more closely related to the strain on the control cheeses (Fig. S4). The secondary microbial population (individually between 1% and 3% of the population) of the D4 surface was composed of species not included in the initial smear-culture mix and included Arthrobacter sp., Corynebacterium variabile, Debaryomyces fabryi, G. candidum, Staphylococcus equorum, and Staphylococcus saprophyticus (Table S2). In addition, some species present in the initial smear-culture mix (C. jadinii and B. casei) were not detected during ripening, while the inoculated B. linens strain was detected at only at a very low relative abundance on the cheese surface of D4 throughout ripening (Table S2).

By comparison, the microbiota was more diverse in cheese S5 (Fig. 2; Table S2). On the cheese surface of S5, the relative abundances of the lactobacilli decreased, while that of B. linens increased significantly (P < 0.004) from day 18 to day 24, reaching 37.05% before decreasing, but not significantly, to 22.84% at day 30 (Table S2). The strain detected was confirmed by PanPhlAn to be that inoculated within the S5 mix (Fig. S3 and S4). The yeasts D. hansenii and G. candidum (components of the S5 mix) were the most abundant populations on the cheese surface at day 18, comprising 21.2% and 37.54% of the microbiota, respectively, but their relative abundances significantly decreased (P < 0.04) by day 24 to 9.57% and 17.6%, respectively, without showing further significant reductions at day 30 (Table S2). S. xylosus, did not correspond to the strain present in the S5 mix (Fig. S3 and S4) and was detected at 9.08% at day 18 but did not change significantly throughout the ripening period (Table S2). In addition, a secondary microbial population, comprising D. fabryi (detected in the S5 mix) (Table S1) and Psychrobacter sp. (not detected in the S5 mix) (Table S1), developed at low relative abundance (1 to 2%) on the surface of cheese S5 (Table S2) over the course of the ripening period. However, some inoculated species were either not detected (S. equorum) at any stage throughout ripening or detected at a very low relative abundance (G. arilaitensis, ~0.44%) on the cheese surface during ripening (Table S2).

Headspace solid-phase microextraction (HS-SPME) gas chromatography-mass spectrometry (GC-MS) was used to analyze the development of volatile compounds at days 0, 18, 24, and 30 of ripening for both the control and test cheeses. In total, 53 volatile compounds that potentially contributed to the flavor development were detected on the cheese surfaces. These compounds are predicted to arise from a variety of substrates and consisted of 8 alcohols, 6 aldehydes, 10 carboxylic acids, 10 esters, 13 ketones, 2 S-thioesters, and 4 sulfur compounds (i.e., a total of 53 compounds) (Fig. 3). As expected, given the microbial diversity on the surface, there was a greater variety and intensity of volatile compounds detected than on the control cheese, on which only 23 of the aforementioned 53 compounds were detected (Fig. 3). In all cheeses, the levels of all volatile compounds detected increased throughout the ripening period, apart from those of 2,3-butanediol, hexanal, heptanal, octanal, nonanal, 2,3-butanedione, and dimethylsulfone (Fig. 3).

The correlation analysis between the relative abundances of microbial species and the abundances of volatile compounds detected on the cheese surface was performed using the Spearman correlation test, as described previously by Walsh et al. (1). From the results of the metagenomic analysis (performed with Kaiju) and the volatile analysis, it was possible to associate both yeasts and bacteria, at the species level, with specific volatile compounds. Figure 4 demonstrates the degrees of correlation between the volatile compounds and the organisms detected.

There was a strong correlation between B. linens and G. candidum and sulfur compounds and 2-methyl-1-butanol. S. xylosus was correlated with sulfur compounds, 2-methyl-1-butanol, and some ethyl esters; Corynebacterium variablile was correlated with ketones. D. hansenii was correlated with acids and alcohols, G. arilaitensis was correlated with ketones, alcohols, and acids, and S. saprophyticus was correlated with ketones, esters, acids, and alcohols (Fig. 4; Table 1).

Using SUPER-FOCUS (https://edwards.sdsu.edu/SUPERFOCUS) (14), whole-metagenome shotgun sequencing was used to characterize the functional potential of the whole microbial community on the cheese surfaces at different stages of ripening. Overall, SUPER-FOCUS was able to assign 62.5% ± 10.9% of reads from the cheese samples to a function. The proportions of assigned reads for each cheese sample are presented in Fig. S2. The functional clusters analyzed were initially organized into three different levels, in relation to the specificity of the metabolic pathways. Pathway data were analyzed to determine the significant differences of the individual metabolic clusters by ANOVA, using SAS 9.3 (12), with the selection of 16 specific functional clusters with relative abundances significantly higher (P < 0.05) on the cheese surfaces of S5 and D4 than on that of the control (Fig. 5).

pH and color analyses were performed on the three cheese types, and the resultant data were examined using a split-plot test, designed with SAS 9.3 (12). A significant interactive effect (P < 0.0001) between smear treatments and ripening time was observed for pH. At days 18, 24, and 30, the pH was significantly higher (P < 0.0001) on the surfaces of S5 and D4 than on that of the control. In addition, the pH was significantly higher (P < 0.0001) on the surface of S5 than on that of D4 from day 18 onwards (Fig. S5).

A significant interactive effect (P < 0.0001) between time and smear treatments was observed for L*, a*, and b* values. The L* value measures the visual lightness (as values increase from 0 to 100), the a* value measures the redness to greenness (positive to negative values, respectively), and the b* value measures the yellowness to blueness (positive to negative values, respectively). At days 18, 24, and 30, the a* value was significantly higher (P < 0.0001) for the surfaces of S5 and D4 than for the surface of the control. At day 30, the a* value was also significantly higher (P < 0.02) on the surface of D4 than on that of S5 (Fig. S6).

Free amino acid (FAA) and free fatty acid (FFA) analyses were performed on the three cheese types, and the experimental results were examined by one-way ANOVA, using SAS 9.3 (12). The concentrations of total FAAs on the surfaces of S5 (15,158 ± 1,683 µg ⋅ mg⁻¹) and D4 (11,914 ± 1,769 µg ⋅ mg⁻¹) were significantly higher (P < 0.05) than those on the control surface (6,605 ± 819 µg ⋅ mg⁻¹). In addition, the concentrations of some individual FAAs, such as tyrosine, proline, and histidine, were significantly higher (P < 0.05) on the surface of S5 than on the surfaces of D4 and the control (Fig. S7).

The concentrations of total FFAs on the surfaces of S5 (22,069 ± 3,875 µg ⋅ mg⁻¹) and D4 (26,562 ± 2,606 µg ⋅ mg⁻¹) were significantly higher (P < 0.05) than on the surface of the control (1,336 ± 70 µg ⋅ mg⁻¹). The concentrations of some individual FFAs, such as C_4:0, C_8:0, C_10:0, C_12:0, C_14:0, and C_18:0, were significantly higher (P < 0.05) on the surface of D4 than on that of S5 or the control (Fig. S7).

A block of commercial Cheddar cheese <24 h after manufacture was aseptically cut into smaller blocks (~8 by 6.5 by 30 cm) and washed with smearing solutions, as described in our previous study (15). Two commercial smear-culture mixes comprising G. candidum, D. hansenii, B. linens, G. arilaitensis, and S. xylosus (S5 mix) (Sacco, Cadorago, Italy) and D. hansenii, C. jadinii, B. casei, and B. linens (D4 mix) (DuPont Danisco, Beaminster, Dorset, United Kingdom) were used to inoculate the surfaces of the cheese curds. The blocks of cheese were washed with the smearing solutions and placed in sterile racks inside sterile plastic bags (Südpack Verpackungen, Ochsenhausen, Germany), as previously described (15). The cheese was ripened for 30 days at 15°C, with a relative humidity of ~97%. At days 7, 10, and 15 of ripening, the cheese blocks were brushed with a sterile sponge that had been soaked in a sterile brine solution (5% NaCl) to uniformly spread the smear microbiota on the cheese surface. As a control, unsmeared cheese blocks were vacuum packed in sterile bags and incubated at 15°C, as with the test cheeses.

Three replicate cheese trials were performed at different times during Cheddar cheese making season. All data presented are the results of the analysis performed on samples taken from the cheese surface (at a depth of ~0.5 cm). All analyses were performed in triplicate.

The pH level was measured on days 0, 18, 24, and 30 using a standard pH meter (MP220; Mettler-Toledo, Schwerzenbach, Switzerland) (52). The data were analyzed by one-way analysis of variance (ANOVA) using SAS 9.3 (12).

At days 0, 18, 24, and 30 of ripening, the color was measured on the cheese surface at room temperature, using a Minolta CR-300 colorimeter (Minolta Camera, Osaka, Japan). The instrument was calibrated on white tile, and the color of the cheese surface was measured using L*, a*, and b* values. The L* value measures the visual lightness (as values increase from 0 to 100), the a* value measures the redness to greenness (positive to negative values, respectively), and the b* value measures the yellowness to blueness (positive to negative values, respectively).

The total DNA was extracted from the smear culture mixes and the cheese samples using the PowerSoil DNA isolation kit as described in the manufacturer’s protocol (Cambio, Cambridge, United Kingdom). For the DNA extraction from the cheese surface, at days 0, 18, 24, and 30, a pretreatment step was included as follows. Samples were removed from different parts of the cheese block and pooled to give a representative sample of 5 g. The cheese was placed in a stomacher bag with 50 ml of 2% trisodium citrate and homogenized using a masticator mixer (IUL SA, Barcelona, Spain) for 5 min.

Fifteen milliliters of the smear-culture mix, or the cheese solution, was placed into sterile Falcon tubes and centrifuged for 30 min at 4,500 × g. After centrifugation, the supernatant was discarded and the pellet was placed in a 2-ml Eppendorf tube. The pellet was washed several times with sterile phosphate-buffered saline (PBS) by centrifuging it at 14,500 × g for 1 min, until the supernatant was completely clear. The pellet was then added to PowerBead tubes (Cambio, Cambridge, United Kingdom) provided with the kit as described in the protocol and homogenized by shaking on the TissueLyser II (Qiagen, West Sussex, United Kingdom) at 20 Hz for 10 min. The DNA was then purified according to the protocol of the standard PowerSoil DNA isolation kit (Cambio, Cambridge, United Kingdom).

Total DNA was initially qualified and quantified by gel electrophoresis and the NanoDrop 1000 (Bio-Sciences, Dublin, Ireland) before more-accurate quantification with the Qubit high-sensitivity DNA assay (Bio-Sciences, Dublin, Ireland).

Whole-metagenome shotgun libraries were prepared in accordance with the Nextera XT DNA library preparation guide from Illumina (53). Libraries for the starter mixture samples were sequenced on the Illumina MiSeq with a 2× 300-bp cycle v3 kit. Libraries for the cheese samples were sequenced on the Illumina NextSeq 500 with a v2 NextSeq 500/550 high-output reagent kit (300 cycles). All sequencing was done in the Teagasc sequencing facility in accordance with standard Illumina sequencing protocols.

Raw whole-metagenome shotgun sequencing reads were processed on the basis of quality and quantity using a combination of Picard tools (https://github.com/broadinstitute/picard) and SAMtools (54). Processing of raw sequence data produced a total of 3,214,480 ± 841,719 filtered reads for samples sequenced on the MiSeq and 19,210,475 ± 12,478,696 filtered reads for samples sequenced on the NextSeq. The metagenomic binning tool Kaiju (11) was used to determine the species-level microbial compositions of samples. The NCBI nonredundant protein database (55) was used with Kaiju. PanPhlAn (13) was used for strain-level analysis of species of interest. PanPhlAn works by aligning sequencing reads against a species pangenome database, built from reference genomes, to identify the gene families present in strains from metagenomic samples. The reference genomes included for each pangenome database are outlined in Table S3. SUPER-FOCUS (14) was used to characterize the microbial metabolic potential of samples. SUPER-FOCUS measures the abundances of subsystems, or groups of proteins with shared functionality, by aligning sequencing reads against a reduced SEED (56) database.

FAA analysis was performed at the end of the ripening (day 30) on the soluble N extracts using a JEOL JLC-500V AA analyzer fitted with a JEOL Na⁺ high-performance cation-exchange column (JEOL, Garden City, Herts, United Kingdom) (57). The chromatographic analyses were conducted at pH 2.2. Results are expressed as micrograms per milligram of cheese.

FFA extractions were performed at the end of the ripening (day 30) according to the method outlined by De Jong and Badings (58). The FFA extracts were derivatized as methyl esters as described by Mannion et al. (59). Fatty acid methyl ester extracts were analyzed using a Varian CP3800 gas chromatograph (Aquilant, Dublin, Ireland) with a CP84000 autosampler and flame ionization detector and a Varian, Inc., 1079 injector (Aquilant, Dublin, Ireland). Results are expressed as micrograms per milligram of cheese.

The volatile compounds were analyzed at days 0, 18, 24, and 30. The surface of the cheese was removed, wrapped in foil, and stored vacuum packed at −20°C until analysis. Before analysis, the samples were defrosted and grated, and 4 g of the cheese surface was used. Analysis was carried out as outlined by Bertuzzi et al. (15).

Statistical analysis was done with SAS 9.3 (12) and R-3.2.2 (60). The R packages ggplot2 and pheatmap were used for data visualization. The vegan package was used to calculate the Bray-Curtis dissimilarity between samples, while the Hmisc package was used for correlation analysis.

Sequencing reads have been deposited in the European Nucleotide Archive under the project accession number PRJEB15423.