Advertisement

The diversity of microbial communities in Chinese milk fan and their effects on volatile organic compound profiles

Open ArchivePublished:December 23, 2020DOI:https://doi.org/10.3168/jds.2020-19053

      ABSTRACT

      Milk fan is a cheese-like fermented milk product produced in Yunnan Province, China. In this study, we characterized the microbial communities of milk fan from 6 distinct geographical origins and investigated their generation of volatile organic compounds (VOC). The microbial communities found in all milk fan samples were dominated by Lactococcus, Lactobacillus, and Raoultella bacteria and Rhodotorula, Torulaspora, and Candida fungi. Samples from the Kunming and Weishan regions had greater bacterial richness, and samples from Xizhou had greater fungal community richness. Sixty prominent VOC (i.e., those having odor activity values ≥1), including esters, acids, alcohols, aldehydes, ketones, and aromatic compounds, were identified by gas chromatography-mass spectrometry analysis of milk fan samples. Pearson correlation analysis revealed that Lactobacillus, Rhodotorula, Lodderomyces, and Debaryomyces had significant correlations with various VOC, revealing a total of 13 compounds that are characteristic of the odor of milk fan. These bacteria and fungi are therefore identified as functional microorganisms that collectively create the complex VOC profile of milk fan. This study provides a comprehensive overview of the microbial community of milk fan and demonstrates its contribution to the unique aroma profile of this fermented milk product

      Key words

      INTRODUCTION

      Milk fan is a cheese-like fermented milk product produced by the Bai people of the Dali Bai Autonomous Prefecture in Yunnan Province, China. It has been widely produced and consumed for more than a thousand years due to its rich nutrition and typical flavor (
      • Liu W.
      • Sun Z.
      • Zhang J.
      • Gao W.
      • Wang W.
      • Wu L.
      • Sun T.
      • Chen W.
      • Liu X.
      • Zhang H.
      Analysis of microbial composition in acid whey for dairy fan making in Yunnan by conventional method and 16s rRNA sequencing.
      ). The prominent aroma features of the milk fan are cheesy, milky, sour, fruit, fatty, and rancid. The overall aroma is typically mild without any obvious pungent odors, such as sulfur and musty (
      • Tian H.
      • Xu X.X.
      • Chen C.
      • Yu H.Y.
      Flavoromics approach to identifying the key aroma compounds in traditional Chinese milk fan.
      ). Milk fan is made through a process that includes preheating, mixing, curdling, agitation, stretching, forming, and drying (Supplemental Figure S1; https://doi.org/10.6084/m9.figshare.13352558.v1), during which it acquires its characteristic paper-fan-like shape (
      • Tian H.
      • Xu X.X.
      • Chen C.
      • Yu H.Y.
      Flavoromics approach to identifying the key aroma compounds in traditional Chinese milk fan.
      ). Similar to Western cheese production techniques that require the addition of rennet and starter cultures, milk fan is produced by coagulation of milk with 10% (vol/vol) sour juice (pH 2.5–3.0), which is a naturally fermented from papaya juice (
      • Liu W.
      • Sun Z.
      • Zhang J.
      • Gao W.
      • Wang W.
      • Wu L.
      • Sun T.
      • Chen W.
      • Liu X.
      • Zhang H.
      Analysis of microbial composition in acid whey for dairy fan making in Yunnan by conventional method and 16s rRNA sequencing.
      ). This traditional processing method and its microbial diversity create the unique combination of sensory attributes and flavor characteristics in milk fan.
      It is well established that microbial community composition influences the flavor, texture, and safety of cheese products (
      • Bertuzzi A.S.
      • Walsh A.M.
      • Sheehan J.J.
      • Cotter P.
      • Crispie D.F.
      • Mcsweeney P.L.H.
      • Kilcawley K.N.
      • Rea M.C.
      Omics-based insights into flavor development and microbial succession within surface-ripened cheese.
      ). Notably, the addition of sour juice to milk fan yields a heterogeneous and locally unique collection of microorganisms (
      • Liu W.
      • Sun Z.
      • Zhang J.
      • Gao W.
      • Wang W.
      • Wu L.
      • Sun T.
      • Chen W.
      • Liu X.
      • Zhang H.
      Analysis of microbial composition in acid whey for dairy fan making in Yunnan by conventional method and 16s rRNA sequencing.
      ). Recently, 16S rRNA gene amplicon sequencing has been used to profile the bacterial communities in milk fan, which revealed that Lactobacillus and Lactococcus were the most abundant genera of bacteria, and that Lactiplantibacillus pentosus, Companilactobacillus crustorum, Levilactobacillus brevis, and Lactobacillus kefiranofaciens were the dominant species (
      • Xue J.
      • Yang Y.
      • Wang Z.
      • Guo Y.
      • Shao Y.
      Bacterial diversity in Chinese Rushan cheese from different geographical origins.
      ). Furthermore, many species of fungi, particularly yeasts, also play a key role in the maturation process of cheese. In some traditional fermented dairy products, Kluyveromyces and Torulaspora have been found to be the dominant yeasts (
      • Gao M.L.
      • Hou H.M.
      • Teng X.X.
      • Zhu Y.L.
      • Hao H.S.
      • Zhang G.L.
      Microbial diversity in raw milk and traditional fermented dairy products (Hurood cheese and Jueke) from Inner Mongolia, China.
      ;
      • Zheng X.
      • Li K.
      • Shi X.
      • Ni Y.
      • Li B.
      • Zhuge B.
      Potential characterization of yeasts isolated from Kazak artisanal cheese to produce flavoring compounds.
      ). However, little research has focused on the fungi present in milk fan; therefore, the key yeast species are, to date, unidentified.
      The aroma of cheese is a crucial product characteristic that influences the sensorial (i.e., organoleptic) experience and consumer preference (
      • Caspia E.L.
      • Coggins P.C.
      • Schilling M.W.
      • Yoon Y.
      • White C.H.
      The relationship between consumer acceptability and descriptive sensory attributes in Cheddar cheese.
      ;
      • Drake S.L.
      • Gerard P.D.
      • Drake M.A.
      Consumer preferences for mild Cheddar cheese flavors.
      ). Microorganisms in cheese convert precursors such as carbohydrates, proteins, and fats into a wide range of volatile compounds responsible for the organoleptic profile of cheese (
      • van Mastrigt O.
      • Gallegos Tejeda D.
      • Kristensen M.N.
      • Abee T.
      • Smid E.J.
      Aroma formation during cheese ripening is best resembled by Lactococcus lactis retentostat cultures.
      ). Many studies of cheese-making processes have documented the relationships between microorganisms and volatile organic compounds (VOC) profiles of different cheeses. For example,
      • Carpino S.
      • Randazzo C.
      • Pino A.
      • Russo N.
      • Rapisarda T.
      • Belvedere G.
      • Caggia C.
      Influence of PDO Ragusano cheese biofilm microbiota on flavour compounds formation.
      characterized the biofilm microbiota of Protected Designation of Origin (PDO) Ragusano cheese and investigated its ability to generate VOC. Some positive correlations were found between certain lactic acid bacteria (LAB) and VOC such as Enterococcus hirae with alcohols; Lactococcus lactis, Lactiplantibacillus plantarum, Lacticaseibacillus casei, and Lactobacillus delbrueckii with aldehydes; and Limosilactobacillus fermentum, Lactobacillus helveticus, and Lentilactobacillus hilgardii with ketones.
      • Duru I.C.
      • Laine P.
      • Andreevskaya M.
      • Paulin L.
      • Kananen S.
      • Tynkkynen S.
      • Auvinen P.
      • Smolander O.-P.
      Metagenomic and metatranscriptomic analysis of the microbial community in Swiss-type Maasdam cheese during ripening.
      found several flavor-forming pathways in Lc. lactis, Lacticaseibacillus rhamnosus, Lb. helveticus, and Propionibacterium freudenreichii, the dominant species in Swiss-type Maasdam cheese. These pathways were predicted to produce methanethiol, acetoin, diacetyl, acetate, ethanol, propionate, and free fatty acids. An investigation of the enzymatic activities of microorganisms isolated from feta cheese brine revealed the proteolytic, peptidolytic, and lipolytic activities of Lacticaseibacillus paracasei, Debaryomyces hansenii, and Saccharomyces cerevisiae (
      • Bintsis T.
      • Vafopoulou-Mastrojiannaki A.
      • Litopoulou-Tzanetaki E.
      • Robinson R.K.
      Protease, peptidase and esterase activities by lactobacilli and yeast isolates from Feta cheese brine.
      ). These studies demonstrate the important role of microorganisms in cheese flavor generation.
      In our previous research, we determined that the key aroma compounds of milk fan are propanoic acid, butanoic acid, octanoic acid, octanal, nonanal, 2-nonanone, and ethyl hexanoate, which make a great contribution to aroma characteristics (
      • Tian H.
      • Xu X.X.
      • Chen C.
      • Yu H.Y.
      Flavoromics approach to identifying the key aroma compounds in traditional Chinese milk fan.
      ). However, research on the relationship between aroma and microbes in milk fan remains scant. In this study, we explored the microbial composition of milk fan and its contributions to the development of aroma profile. First, we used high-throughput sequencing to characterize the composition of bacterial and fungal communities in milk fan collected from different regions of Yunnan Province. Then, we systematically explored the relationships between the core microbes and VOC of milk fan by Pearson correlation analysis. The results of this study can serve as a theoretical basis for the microbial community and its correlation with volatile profiles in milk fan, and can be used to improve the fermentation and aromatic quality of such traditional cheese.

      MATERIALS AND METHODS

      Collection of Milk Fan Samples

      Eighteen milk fan samples were collected from 6 regions (Supplemental Figure S2), including Kunming (K), Yuxi (Y), Dengchuan (D), Xizhou (X), Jianchuan (J), and Weishan (W). Three workshops within a region were chosen and one sample was collected from each site from May 4 to 9, 2019. All samples were made from raw milk produced by local cows and curded using locally fermented sour juice. They were prepared following the same procedures by local residents, and then collected at a similar point in aging (Supplemental Figure S1). Approximately 500 g from each sample was transferred into sealed sterile bags and sent to the Food Microbiology Laboratory of Shanghai Institute of Technology (Shanghai, China). Each sample was divided into 2, with one part designated for microbial community analysis and refrigerated at 4°C for less than 12 h, and the other designated for volatile compounds and chemical composition analysis, and thus frozen at −20°C. The contents of moisture, ash, fat, and protein, and the pH value of milk fan, were determined according to standard methods nos. 926.08, 945.46, 905.02, and 930.29, respectively, described by
      • AOAC International
      Official Methods of Analysis.
      .

      Sample Processing and DNA Extraction

      Before extracting, 1 g of sample was ground in a sterile quartz mortar for 2 min and homogenized with 9 mL of sterile NaCl solution (0.85%, wt/vol) for 10 min. To separate fat and coarse particles, the suspension was centrifuged at 5,000 × g for 10 min at 4°C, and 1 g of the upper sediment was collected into a sterile microcentrifuge tube using a sterile spoon. Total DNA was extracted from the sediment using the Qiagen DNA Mini-Kit (Qiagen, Hilden, Germany) according to the manufacturer's protocol.
      A Nanodrop ND 1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE) and 0.8% agarose gel electrophoresis were used to determine the concentration and purity of the extracted DNA, respectively. Samples were then stored at −20°C.

      PCR Amplification, Quantification, and Sequencing

      For bacteria, the V3–V4 domains of 16S rDNA were amplified by PCR (95°C for 5 min; followed by 25 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 40 s with a final extension of 72°C for 10 min) using primer pairs 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) with specific barcode. For fungi, the ITS1 rDNA regions were amplified by PCR (95°C for 2 min, followed by 30 cycles at 95°C for 30 s, 61°C for 30 s, and 72°C for 45 s with a final extension at 72°C for 10 min) with the primers ITS5F (5′-GGAAGTAAAAGTCGTAACAAGG-3′) and ITS1R (5′-GCTGCGTTCTTCATCGATGC-3′) with specific barcode. The purified amplicons were pooled at equimolar concentrations, and further paired-end sequencing was performed using an Illumina MiSeq instrument (Illumina Inc., San Diego, CA).

      Volatile Profile Analysis of Milk Fan

      Three grams of each milk fan sample was cut into evenly sized cubes and placed in 15-mL headspace bottles with Teflon covers. To each bottle, 100 μL of 2-octanol solution (19 mg/L) was added as an internal standard solution. The volatile compounds were extracted with a solid-phase microextraction (SPME) fiber (Supelco Inc., Bellefonte, PA). After equilibration, a 75-μm carboxen-polydimethylsiloxane fiber was exposed to the headspace of the glass vial for 30 min at 60°C.
      For GC-MS analysis, an Agilent 7890-5973 GC-MS (Agilent Technologies, Santa Clara, CA) and HP-Innowax analytical fused-silica capillary column (60 m × 0.25 mm × 0.25 μm; Agilent) were used. The GC-MS conditions were as follows. First, desorption was performed within 5 min in splitless mode at 250°C. After injection of samples, the oven temperature was held at 40°C for 4 min, increased at 3°C/min to 100°C and held for 2 min, increased to 150°C at 4°C/min, and then increased at 10°C/min to 230°C and held at this temperature for 5 min. The quadrupole mass filter was operated at 150°C. Helium was used as the carrier gas at a constant linear velocity of 1 mL/min. The transfer line temperature was 250°C. The quadrupole mass spectrometer (mass selective detector) was used for chemical identification with electron ionization energy set to 70 eV, the ion source temperature set at 230°C, and the quadrupole in scan mode (i.e., the total ion currents were recorded from m/z of 30 to 450 at a rate of 1 scan per second). Each analysis was performed in triplicate.
      Then, the odor activity values (OAV) were calculated based on the ratio of the concentration of the volatile compound to its detection odor threshold. The thresholds were drawn from the literature (
      • Gemert L.J.V.
      Odour threshold values in water in odour thresholds.
      ). Odor activity values were used to preliminarily evaluate the contribution of each compound to overall milk fan aroma. When the OAV is >1, a compound is considered to contribute to overall aroma.

      Bioinformatics and Statistical Analyses

      Sequences were analyzed using the QIIME-1.9.1 (
      • Caporaso J.G.
      • Kuczynski J.
      • Stombaugh J.
      • Bittinger K.
      • Bushman F.D.
      • Costello E.K.
      • Fierer N.
      • Peña A.G.
      • Goodrich J.K.
      • Gordon J.I.
      • Huttley G.A.
      • Kelley S.T.
      • Knights D.
      • Koenig J.E.
      • Ley R.E.
      • Lozupone C.A.
      • McDonald D.
      • Muegge B.D.
      • Pirrung M.
      • Reeder J.
      • Sevinsky J.R.
      • Turnbaugh P.J.
      • Walters W.A.
      • Widmann J.
      • Yatsunenko T.
      • Zaneveld J.
      • Knight R.
      QIIME allows analysis of high-throughput community sequencing data.
      ) pipeline. First, low-quality sequences were filtered and the remaining sequences were assigned with effective tags to operational taxonomic units (OTU) based on 97% similarity of sequences determined by Uparse software-7.0.1090 (
      • Edgar R.C.
      UPARSE: Highly accurate OTU sequences from microbial amplicon reads.
      ). Representative OTUs with high frequencies of occurrence were selected and annotated for taxonomic information (e.g., phylum, family, and genus levels). The α-diversity indices were calculated to evaluate microbial community richness (according to the Ace and Chao1 estimators) and community diversity (according to the Shannon and Simpson estimators) in milk fan samples of different geographical origins. Principal coordinate analysis (PCoA) based on weighted and unweighted UniFrac distance calculations was used to compare the global microbiota composition in each group at OTU levels. Correlations between VOCs (with OAV >1) present in all samples and bacterial and fungal classes with relative abundance >1% in the milk fan were calculated by Pearson correlation analysis with the Hmisc package in R. Strong correlations (|R| > 0.6, P < 0.05) were further analyzed using Cytoscape software (version 3.5.1; https://cytoscape.org/) with nodes representing VOC and microbial class and line colors representing the positive or negative correlations (
      • Shannon P.
      • Markiel A.
      • Ozier O.
      • Baliga N.S.
      • Wang J.T.
      • Ramage D.
      • Ideker T.
      Cytoscape: a software environment for integrated models of biomolecular interaction networks.
      ). For the microbial classes, nodes were sized to reflect the relative abundance of the microbes, and the thickness of the connection lines shows the strength of the correlation.
      All statistical analyses were executed using R software (version 3.2.5; https://www.R-project.org/). The Kruskal-Wallis test was used to compare the relative abundances of taxa based on the rarefied OTU subset. The Benjamini-Yekutieli method was used to control for multiple tests. The PCoA analyses were performed using the ade4 package (
      • Zapala M.A.
      • Schork N.J.
      Multivariate regression analysis of distance matrices for testing associations between gene expression patterns and related variables.
      ), and the pheatmap package was used for cluster analysis and heatmap construction in R. The data for chemical analysis, α-diversity, and correlation analysis were analyzed using ANOVA with Duncan's multiple range test; P < 0.05 was considered to indicate statistical significance.

      Nucleotide Sequence Accession Number

      The sequence data reported in this paper have been deposited in the National Center for Biotechnology Information database under accession number PRJNA608625 (https://www.ncbi.nlm.nih.gov/nuccore/?term=PRJNA608625).

      RESULTS

      Sample Information

      The results of the chemical composition analysis and sequences data of milk fan samples from 6 different regions are presented in Table 1. Significant differences (P < 0.05) in the contents of ash, protein, fat, and moisture were found in milk fan from different regions, but no difference (P > 0.05) in pH value was observed in these samples. The total numbers of effective sequences of samples ranged from 35,495 to 57,481 for bacteria, and from 62,010 to 72,990 for fungi. All of the sequences were clustered with a 97% sequence identity cut-off. The number of OTU for bacteria ranged from 143 to 336, and for fungi ranged from 37 to 96. Rarefaction curves (Supplemental Figure S3) showed that the number of observed genera was saturated, suggesting that good representation of the microbial diversity had been captured by the sequencing depth used.
      Table 1Chemical composition and sequence data of milk fan samples from 6 locations in China (n = 3 samples per location)
      Values are means ± SD.
      Sampling location (code)Moisture (%)Ash (%)Protein (%)Fat (%)pHNo. of readsNo. of OTU
      Operational taxonomic unit.
      Average length of sequence base
      BacteriaFungiBacteriaFungiBacteriaFungi
      Kunming (K)17.00 ± 1.25
      Values in the same column with different superscripts are significantly different at P < 0.05.
      2.21 ± 0.05
      Values in the same column with different superscripts are significantly different at P < 0.05.
      47.6 ± 0.12
      Values in the same column with different superscripts are significantly different at P < 0.05.
      27.3 ± 0.97
      Values in the same column with different superscripts are significantly different at P < 0.05.
      4.90 ± 0.02
      Values in the same column with different superscripts are significantly different at P < 0.05.
      50,335 ± 4,86572,990 ± 1,670336 ± 9437 ± 6427 ± 0237 ± 0
      Yuxi (Y)13.40 ± 0.89
      Values in the same column with different superscripts are significantly different at P < 0.05.
      2.74 ± 0.07
      Values in the same column with different superscripts are significantly different at P < 0.05.
      37.8 ± 0.59
      Values in the same column with different superscripts are significantly different at P < 0.05.
      36.2 ± 0.21
      Values in the same column with different superscripts are significantly different at P < 0.05.
      4.94 ± 0.03
      Values in the same column with different superscripts are significantly different at P < 0.05.
      34,642 ± 2,38762,010 ± 14,213143 ± 1671 ± 9427 ± 0301 ± 41
      Dengchuan (D)20.60 ± 1.58
      Values in the same column with different superscripts are significantly different at P < 0.05.
      2.05 ± 0.09
      Values in the same column with different superscripts are significantly different at P < 0.05.
      34.1 ± 0.23
      Values in the same column with different superscripts are significantly different at P < 0.05.
      36.3 ± 0.86
      Values in the same column with different superscripts are significantly different at P < 0.05.
      4.95 ± 0.02
      Values in the same column with different superscripts are significantly different at P < 0.05.
      35,495 ± 5,11067,155 ± 5,988174 ± 796 ± 9427 ± 0268 ± 12
      Xizhou (X)19.3 ± 1.14
      Values in the same column with different superscripts are significantly different at P < 0.05.
      2.15 ± 0.11
      Values in the same column with different superscripts are significantly different at P < 0.05.
      45.8 ± 0.09
      Values in the same column with different superscripts are significantly different at P < 0.05.
      24.4 ± 0.33
      Values in the same column with different superscripts are significantly different at P < 0.05.
      4.87 ± 0.03
      Values in the same column with different superscripts are significantly different at P < 0.05.
      57,481 ± 3,16672,403 ± 958270 ± 444 ± 1428 ± 0249 ± 5
      Jianchuan (J)12.6 ± 0.96
      Values in the same column with different superscripts are significantly different at P < 0.05.
      2.72 ± 0.08
      Values in the same column with different superscripts are significantly different at P < 0.05.
      37.1 ± 0.79
      Values in the same column with different superscripts are significantly different at P < 0.05.
      38.7 ± 0.49
      Values in the same column with different superscripts are significantly different at P < 0.05.
      4.94 ± 0.03
      Values in the same column with different superscripts are significantly different at P < 0.05.
      41,233 ± 8,55272,952 ± 2,424336 ± 9466 ± 2426 ± 1222 ± 0
      Weishan (W)16.2 ± 0.61
      Values in the same column with different superscripts are significantly different at P < 0.05.
      2.53 ± 0.04
      Values in the same column with different superscripts are significantly different at P < 0.05.
      34.7 ± 0.41
      Values in the same column with different superscripts are significantly different at P < 0.05.
      39.2 ± 0.11
      Values in the same column with different superscripts are significantly different at P < 0.05.
      4.88 ± 0.03
      Values in the same column with different superscripts are significantly different at P < 0.05.
      46,741 ± 3,10270,786 ± 1,841336 ± 9461 ± 2422 ± 1235 ± 2
      a–f Values in the same column with different superscripts are significantly different at P < 0.05.
      1 Values are means ± SD.
      2 Operational taxonomic unit.

      Microbial Composition and Core Microbiota of Milk Fan

      We first quantified the predominant microbiota in the milk fan samples; taxa with average relative abundance exceeding 1% were identified as dominant (Figure 1A, Figure 1B, Supplemental Tables S1 and S2; https://doi.org/10.6084/m9.figshare.13352546.v1). At the phylum level, the microbial community from all samples consisted mainly of 5 bacterial and 2 fungal phyla. The bacterial community was dominated by Firmicutes, Proteobacteria, Bacteroidetes, Actinobacteria, and Epsilonbacteraeota. The fungal community was dominated by Basidiomycota and Ascomycota. For bacteria (Figure 1A), the most dominant genus was Lactococcus (39.96%), followed by Lactobacillus (16.65%), Raoultella (8.38%), Streptococcus (7.32%), Acetobacter (4.99%), Acinetobacter (3.44%), Klebsiella (2.54%), Leuconostoc (2.44%), Chryseobacterium (2.16%), Enhydrobacter (1.59%), Enterococcus (1.45%), and Pseudomonas (1.28%). For fungi (Figure 1B), the dominant genera were Rhodotorula (36.09%), Torulaspora (17,65%), Candida (16.61%), Lodderomyces (7.82%), Naganishia (5.83%), Cystofilobasidium (2.98%), Trichosporon (2.94%), Pichia (2.85%), Dekkera (2.55%), and Kluyveromyces (2.52%). Furthermore, milk fan samples from different locations had similar compositions of microbiota but different relative abundances. For example, the bacterial genus Lactococcus was dominant in samples K (66.09%), Y (61.07%), X (45.03%), and D (44.29%), but not in sample J (21.38%) or W (1.89%). Rhodotorula was the dominant fungal genus in samples J (72.75%), X (71.26%), and D (48.02%), but Candida dominated in samples K (49.53%) and W (42.82%), and Torulaspora dominated in sample Y (59.48%).
      Figure thumbnail gr1
      Figure 1Relative abundances of (A) bacteria and (B) fungi at the phylum and genus level in milk fan samples. K, Y, D, X, J, and W are milk fan samples collected from Kunming, Yuxi, Dengchuan, Xizhou, Jianchuan, and Weishan, respectively.

      α-Diversity Analysis of Milk Fan Samples from Different Regions

      Analysis of the α-diversity of microbes in milk fan samples was conducted after normalization with the QIIME platform to compare the Ace, Chao1, Shannon, and Simpson indices of the samples (Figure 2, Supplemental Tables S3 and S4). The first 2 indices, Ace and Chao1, represent the richness of the community, whereas Shannon and Simpson indices reflect the species diversity of the community (
      • He W.
      • Chung H.Y.
      Exploring core functional microbiota related with flavor compounds involved in the fermentation of a natural fermented plain sufu (Chinese fermented soybean curd).
      ). For bacteria, the Ace and Chao1 indices of samples K and W were higher than those of the other samples, indicating the greatest bacterial community richness of these samples. For fungi, sample X had the highest fungal community richness among these samples according to highest Ace and Chao1 index values. The Shannon and Simpson indices of all samples were quite scattered and distinct for each sample, revealing no general pattern in either bacterial or fungal community diversity.
      Figure thumbnail gr2
      Figure 2Differences in the microbial α-diversities of milk fan samples: (A, E) Ace, (B, F) Chao1, (C, G) Shannon, and (D, H) Simpson indices, respectively, in bacterial communities (A, B, C, and D) and in fungal communities (E, F, G, and H). K, Y, D, X, J, and W are milk fan samples collected from Kunming, Yuxi, Dengchuan, Xizhou, Jianchuan, and Weishan, respectively. The box represents the middle 50% of the data. The horizontal line inside the box is the median value of the index. The upper and lower ends of the box are the hinges of the distribution of index value. The vertical lines from the ends of box connect the extreme data points to their respective hinges.

      β-Diversity in the Milk Fan Samples

      To analyze the β-diversities of the microbiota in the milk fan samples, PCoA plots were calculated and compared based on weighted and unweighted UniFrac distances of samples collected from different regions (Figure 3). For bacteria, PCoA revealed 5 distinct groups. The weighted UniFrac distance distributions of samples D, X, J, and W were large, which indicated a microbial structural difference among these samples. The weighted UniFrac distance between samples K and Y was relatively small, and the symbols that represent these samples were slightly mixed (Figure 3A). For fungi, PCoA also revealed 5 distinct groups. The weighted UniFrac distance distributions of samples K, Y, J, and W were large. The weighted UniFrac distance between D and X samples was relatively small, and the symbols that represent these samples were slightly overlapped (Figure 3C).
      Figure thumbnail gr3
      Figure 3Microbial β-diversity of milk fan samples in bacteria based on (A) weighted and (B) unweighted UniFrac distance and in fungi based on (C) weighted and (D) unweighted UniFrac distance. An R-value close to 1 indicates that the difference between groups is greater than the difference within groups. K, Y, D, X, J, and W are milk fan samples collected from Kunming, Yuxi, Dengchuan, Xizhou, Jianchuan, and Weishan, respectively.
      To further identify clustering patterns of microbes at the genus level in milk fan samples, the distribution of genera in different samples was analyzed and the results are shown in Figure 4. For bacteria, samples X and K were enriched in Lactococcus, whereas sample J was enriched in Lactobacillus. In contrast, the amounts of Aeromonas, Klebsiella, and Streptococcaceae were significantly lower in samples X, K, and W than in the other samples. For the fungi, Rhodotorula were more abundant in samples X and W than in others, whereas Candida and Lodderomyces were present at higher levels in sample K than in the other samples.
      Figure thumbnail gr4
      Figure 4Heat map of the microbiological diversity at the genus level of (A) bacteria and (B) fungi in milk fan samples. The data used for heat map are the log10 values of the operational taxonomic unit (OTU) number in samples. K, Y, D, X, J, and W are milk fan samples collected from Kunming, Yuxi, Dengchuan, Xizhou, Jianchuan, and Weishan, respectively.

      Volatile Compound Profiling of Milk Fan

      We used GC-MS analysis to determine the VOC profile of milk fan samples. A total of 60 compounds that might contribute to milk fan aroma were identified and quantified (Supplemental Table S5). These compounds were clustered into 6 different groups according to their chemical structures, and comprised 17 esters, 11 carboxylic acids, 8 aromatic compounds, 7 ketones, 7 alcohols, and 5 aldehydes.
      To better illustrate the distinctions of VOC content among various samples, a heatmap of the VOC contents was created (Figure 5). Butyric acid and hexanoic acid were present at high concentrations in all samples, which confirms findings that butyric acid is key to the aroma profile of milk fan (
      • Tian H.
      • Xu X.X.
      • Chen C.
      • Yu H.Y.
      Flavoromics approach to identifying the key aroma compounds in traditional Chinese milk fan.
      ). The concentrations of 2-pentanone, ethyl acetate, and 2-nonanone were greatest in samples K and Y, whereas ethyl hexanoate, ethyl caprylate, and ethyl butyrate concentrations were greatest in samples J and W, and the concentration of propionic acid was greatest in samples D and X. As summarized in Table 2, 13 VOC were determined to be characteristic of milk fan: ethyl butyrate, isobutyl butyrate, ethyl caprylate, ethyl hexanoate, 1-nonanal, propionic acid, butyric acid, valeric acid, hexanoic acid, heptanoic acid, phenol, octanoic acid, and decanoic acid. The OAV of these VOC were >1, indicating their contributions to the aroma of samples. These VOC were present in samples of all regions, showing that these compounds might be widespread in milk fan.
      Figure thumbnail gr5
      Figure 5Heat map of the volatile compound contents of milk fan samples. The data used for heat map are the log10 values of the volatile compound contents in samples. K, Y, D, X, J, and W are milk fan samples collected from Kunming, Yuxi, Dengchuan, Xizhou, Jianchuan, and Weishan, respectively.
      Table 2Odor activity values (OAV
      OAV were calculated based on the ratio of the concentration of the volatile compound to its detection odor threshold.
      ) of 13 characteristic volatile compounds detected in milk fan samples
      No.CompoundSample
      K, Y, D, X, J, and W are milk fan samples collected from Kunming, Yuxi, Dengchuan, Xizhou, Jianchuann, and Weishan regions, respectively.
      T
      Detection threshold from the literature (Gemert, 2011).
      (μg/kg)
      KYDXJW
      1Ethyl butyrate8751676112959
      2Isobutyl butyrate41310489.4
      3Ethyl caprylate1121102102441126
      4Ethyl hexanoate241041181431990
      51-Nonanal1326771163033.1
      6Propionic acid256444810617
      7Butyric acid19614118415741161238.5
      8Valeric acid4097191435758
      9Hexanoic acid82304050103105100
      10Heptanoic acid525681923
      11Phenol52534721
      12Octanoic acid117434213626622415
      13Decanoic acid322661080
      1 OAV were calculated based on the ratio of the concentration of the volatile compound to its detection odor threshold.
      2 K, Y, D, X, J, and W are milk fan samples collected from Kunming, Yuxi, Dengchuan, Xizhou, Jianchuann, and Weishan regions, respectively.
      3 Detection threshold from the literature (
      • Gemert L.J.V.
      Odour threshold values in water in odour thresholds.
      ).

      Correlation Analysis of Core Microbiota and Major Volatile Compounds

      Correlations between detected VOC and microorganisms were explored using Pearson rank correlation coefficient analysis. For simplification, we considered the 20 dominant bacterial taxa, 15 fungal taxa, and 13 characteristic VOC listed above. A total of 38 significant pair-wise correlations (|r| > 0.6, P < 0.05) among these 3 data sets were obtained by Pearson rank correlation test and visualized as a concatenation panel (Figure 6, Supplemental Table S6). Of these, we found 31 positive correlations (red lines) and 7 negative correlations (black lines) between VOCs and microorganisms. In general, bacteria had 16 pair-wise positive associations with VOCs and fungi had 15 pair-wise associations with VOCs, which suggests that both bacteria and fungi are crucial for formation of the volatile compounds of milk fan. As the dominant bacterial genus, Lactobacillus was positively correlated with octanoic acid (r = 0.90, P < 0.01), hexanoic acid (r = 0.79, P = 0.02), ethyl hexanoate (r = 0.63, P = 0.03), butyric acid (r = 0.60, P = 0.04), and ethyl caprylate (r = 0.60, P = 0.04). The dominant fungal genus, Rhodotorula, was positively correlated with heptanoic acid (r = 0.68, P = 0.04), octanoic acid (r = 0.67, P = 0.04), decanoic acid (r = 0.88, P = 0.02), and isobutyl butyrate (r = 0.86, P = 0.03). Notably, some negative correlations were observed, such as Lactococcus with heptanoic acid (r = −0.90, P = 0.01) and decanoic acid (r = −0.83, P = 0.04).
      Figure thumbnail gr6
      Figure 6Correlation analysis between major volatile organic compounds (VOC) and the dominant bacterial and fungal taxa at the genus level in milk fan by Pearson correlation analysis. Yellow diamonds represent the volatile components (found in all samples with odor activity values >1), and green and blue circles represent the most abundant bacterial and fungal taxa. Red (black) lines represent positive (negative) correlations between volatile compounds and the microbiome. In the case of microbial taxa, the size of the node reflects their average relative abundance. Color and thickness of connector lines illustrate the direction and strength of the correlation coefficient.
      Identifying a greater number of connections for certain bacteria and fungi with VOCs suggests the importance of these microorganisms in milk fan. For example, Lactobacillus had 5 pair-wise associations with VOC, and Rhodotorula, Lodderomyces, and Debaryomyces had 4, 4, and 3 pair-wise associations with VOC, respectively. The associations of these 4 microbes could cover all characteristic VOCs in the milk fan samples of the current study, which indicates the fundamental role of these bacteria and fungi.

      DISCUSSION

      Scientific advances in understanding the function of microorganisms have laid the foundation for the industrialization of fermented food manufacturing (
      • Stefanovic E.
      • Kilcawley K.N.
      • Roces C.
      • Rea M.C.
      • O'Sullivan M.
      • Sheehan J.J.
      • McAuliffe O.
      Evaluation of the potential of Lactobacillus paracasei adjuncts for flavor compounds development and diversification in short-aged Cheddar cheese.
      ). Compared with industrialized cheeses, traditional cheeses are still mostly produced and consumed by native peoples using locally available raw food materials based on their inherited knowledge and artisanal techniques (
      • Tamang J.P.
      • Watanabe K.
      • Holzapfel W.H.
      Review: Diversity of microorganisms in global fermented foods and beverages.
      ). Milk fan is one such traditional, cheese-like fermented food. Local production limits the yield and popularity of milk fan, and the unclean production environment means that milk fan is easily contaminated by pathogenic bacteria. Furthermore, milk fan is often unstable and non-uniform in quality between different production batches, and thus inconsistent in taste and flavor. In this study, we explored the core microbes in milk fan at the bacterial and fungal levels, and established relationships between microbial composition and VOCs. These results may contribute to the development and industrialization of this traditional cheese.
      The analysis of bacterial and fungal profiles demonstrated that samples K and W had the greatest bacterial community richness whereas samples X had highest fungal community richness. Bacterial community in samples K and Y and fungal community in samples D and X showed convergent profiles. Most microorganisms in milk fan are introduced by the addition of fermented sour juice. During its fermentation, various microorganisms are unintentionally introduced from the environment (
      • Wang Y.
      • Guo J.
      • Huang A.
      Study on bacterial diversity in traditional sour whey of Yunnan province.
      ). Thus, the local environment, including the climate and position, may affect the microbiological diversity of milk fan samples.
      We confirmed that Lactococcus and Lactobacillus were the most abundant bacteria microbial genera in all samples. These results are consistent with previous studies of the bacterial diversity in milk fan (
      • Liu W.
      • Sun Z.
      • Zhang J.
      • Gao W.
      • Wang W.
      • Wu L.
      • Sun T.
      • Chen W.
      • Liu X.
      • Zhang H.
      Analysis of microbial composition in acid whey for dairy fan making in Yunnan by conventional method and 16s rRNA sequencing.
      ;
      • Xue J.
      • Yang Y.
      • Wang Z.
      • Guo Y.
      • Shao Y.
      Bacterial diversity in Chinese Rushan cheese from different geographical origins.
      ). Another noteworthy feature of the bacterial profile was the dominant role of LAB (e.g., Lactococcus, Lactobacillus, Streptococcus, and Leuconostoc). As stated before, milk fan is produced following the traditional procedures involving sour juice, which is rich in LAB. The addition of sour juice could thus introduce LAB to milk fan (
      • Wang Y.
      • Guo J.
      • Huang A.
      Study on bacterial diversity in traditional sour whey of Yunnan province.
      ). Lactococcus and Lactobacillus are also dominant in other Chinese traditional fermented dairy products, such as Hurood and Jueke cheese in Inner Mongolia, and Kazak artisanal cheese from the Uighur Autonomous Region (
      • Gao M.L.
      • Hou H.M.
      • Teng X.X.
      • Zhu Y.L.
      • Hao H.S.
      • Zhang G.L.
      Microbial diversity in raw milk and traditional fermented dairy products (Hurood cheese and Jueke) from Inner Mongolia, China.
      ;
      • Zheng X.
      • Liu F.
      • Li K.
      • Shi X.
      • Ni Y.
      • Li B.
      • Zhuge B.
      Evaluating the microbial ecology and metabolite profile in Kazak artisanal cheeses from Xinjiang, China.
      ). However, their relative abundance in LAB composition differs. Streptococcus is more abundant in Kazak artisanal cheese, and Leuconostoc was abundant in Hurood and Jueke cheese. In milk fan, these 2 genera were detected at lower abundances in the present study.
      The functional role of fungi in cheese manufacture has been widely reported (
      • Büchl N.R.
      • Seiler H.
      Yeast and molds: Yeasts in milk and dairy products.
      ;
      • Gao M.L.
      • Hou H.M.
      • Teng X.X.
      • Zhu Y.L.
      • Hao H.S.
      • Zhang G.L.
      Microbial diversity in raw milk and traditional fermented dairy products (Hurood cheese and Jueke) from Inner Mongolia, China.
      ), but not for milk fan. Our results revealed that the fungal composition of milk fan samples had a lower degree of diversity than that of bacteria. The most abundant fungal genera were Rhodotorula, Torulaspora, and Candida. Among these, Torulaspora was found in Kazak artisanal cheese and identified as the major flavor-producing taxon (
      • Zheng X.
      • Liu F.
      • Shi X.
      • Wang B.
      • Li K.
      • Li B.
      • Zhuge B.
      Dynamic correlations between microbiota succession and flavor development involved in the ripening of Kazak artisanal cheese.
      ). Candida was identified as dominant yeast genus in hurunge, a traditional starter culture that is used for fermented dairy products by nomadic families in the Inner Mongolia Autonomic Region (
      • Shuangquan
      • Burentegusi
      • Yu B.
      • Miyamoto T.
      Microflora in traditional starter cultures for fermented milk, hurunge, from Inner Mongolia, China.
      ). Rhodotorula has not been reported in Chinese cheese products before but was detected in curd of Spanish blue-veined Cabrales cheese (
      • Álvarez-Martín P.
      • Flórez A.B.
      • López-Díaz T.M.
      • Mayo B.
      Phenotypic and molecular identification of yeast species associated with Spanish blue-veined Cabrales cheese.
      ).
      It is noteworthy that cheese ripening is highly dependent on microbial interactions, and thus cheese flavor-production is due to a mixture of various microorganisms (
      • Carpino S.
      • Randazzo C.
      • Pino A.
      • Russo N.
      • Rapisarda T.
      • Belvedere G.
      • Caggia C.
      Influence of PDO Ragusano cheese biofilm microbiota on flavour compounds formation.
      ). Therefore, the correlations between microorganisms and VOCs were calculated and visualized to further investigate the effect of microorganisms on VOC generation. We observed that Lactobacillus, Rhodotorula, Lodderomyces, and Debaryomyces had more positive correlations with VOC, highlighting their function in the production of volatile compounds. Among them, Lactobacillus was correlated with octanoic acid, hexanoic acid, ethyl hexanoate, butyric acid, and ethyl caprylate. The LAB collectively influence flavor development of cheese through several basic mechanisms (
      • Curtin A.C.
      • Gobbetti M.
      • McSweeney P.L.H.
      Peptidolytic, esterolytic and amino acid catabolic activities of selected bacterial strains from the surface of smear cheese.
      ;
      • Jin H.
      • Mo L.
      • Pan L.
      • Hou Q.
      • Li C.
      • Darima I.
      • Yu J.
      Using PacBio sequencing to investigate the bacterial microbiota of traditional Buryatian cottage cheese and comparison with Italian and Kazakhstan artisanal cheeses.
      ); in particular, Lactobacillus is generally considered to have strong esterase activity (
      • Fox P.F.
      • Wallace J.M.
      Formation of flavor compounds in cheese.
      ;
      • Collins Y.F.
      • McSweeney P.L.H.
      • Wilkinson M.G.
      Lipolysis and free fatty acid catabolism in cheese: A review of current knowledge.
      ). Esterases of Lactobacillus not only hydrolyze the di- and monoglycerides in milk fat to release short-chain fatty acids, but also synthesize esters from glycerides and alcohols via a transferase reaction (
      • Deeth H.C.
      • Touch V.
      Methods for detecting lipase activity in milk and milk products.
      ;
      • Fox P.F.
      • Guinee T.P.
      • Cogan T.M.
      • McSweeney P.L.H.
      Biochemistry of Cheese Ripening.
      ;
      • Holland R.
      • Liu S.Q.
      • Crow V.L.
      • Delabre M.L.
      • Lubbers M.
      • Bennett M.
      • Norris G.
      Esterases of lactic acid bacteria and cheese flavour: Milk fat hydrolysis, alcoholysis and esterification.
      ). Esterase activity could contribute to the formation of these VOC that are correlated with Lactobacillus.
      We also found statistical correlations between fungal genera and VOC. Rhodotorula, Debaryomyces, and Lodderomyces were positively correlated with heptanoic acid, octanoic acid, 1-nonanal, decanoic acid, heptanoic acid, octanoic acid, phenol, propionic acid, valeric acid, and isobutyl butyrate. In cheese, yeasts have a wide range of esterolytic and lipolytic enzymes to generate VOC (
      • Mounier J.
      Microbial interactions in smear-ripened cheeses.
      ). Debaryomyces was found to be a major yeast in most surface-ripened cheeses, and to primarily produce free fatty acids, branched-chain aldehydes, and alcohols (
      • Gori K.
      • Sørensen L.M.
      • Petersen M.A.
      • Jespersen L.
      • Arneborg N.
      Debaryomyces hansenii strains differ in their production of flavor compounds in a cheese-surface model.
      ). Rhodotorula and Lodderomyces are present in many types of cheeses with lipolytic activity (
      • Zheng X.
      • Li K.
      • Shi X.
      • Ni Y.
      • Li B.
      • Zhuge B.
      Potential characterization of yeasts isolated from Kazak artisanal cheese to produce flavoring compounds.
      ). These results agreed with our current observations.
      Lactobacillus, Rhodotorula, Debaryomyces, and Lodderomyces were associated with all VOCs that were characteristic of milk fan odor, suggesting that they can be regarded as the core functional microorganisms that are indispensable in aroma formation in milk fan. However, these correlations were established mainly based on statistical analysis, and it is difficult to determine which volatile compounds are specifically produced by a specific group of microorganisms. Volatile compounds can be produced by a wide variety of microorganisms through highly connected metabolic pathways (
      • Bertuzzi A.S.
      • Walsh A.M.
      • Sheehan J.J.
      • Cotter P.
      • Crispie D.F.
      • Mcsweeney P.L.H.
      • Kilcawley K.N.
      • Rea M.C.
      Omics-based insights into flavor development and microbial succession within surface-ripened cheese.
      ). Furthermore, the number of samples was limited in the present study, although sufficient to support selection of starters for future large-scale production of milk fan. In the future, by including more samples and using particular starter mixtures to include and omit certain groups of microorganisms, we will establish clear microbe–aroma relationships between microorganism groups and particular VOC in dairy products.

      CONCLUSIONS

      We characterized the microbial community diversity of milk fan samples from 6 regions of Yunnan Province and correlated the characteristic VOC with their corresponding producing microbes. The results showed that microbial profiles of milk fan varied among regions. Bacterial taxa including Lactococcus, Lactobacillus, and Raoultella, and fungal taxa, including Rhodotorula, Torulaspora, and Candida were identified as the dominant genera of milk fan. Based on correlations among the microorganisms and VOC, Lactobacillus, Rhodotorula, Lodderomyces, and Debaryomyces are the core functional microorganisms that together facilitate the formation of the complex VOC profile of milk fan. More in-depth studies are required to validate these applications in the future, such as using pure strains of these microorganisms in fermentation to verify their roles in aroma formation.

      ACKNOWLEDGMENTS

      This work was sponsored by the National Natural Science Foundation of China (Award No. 31771943). The samples were sequenced by Shanghai Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China). The authors declare that they have no conflicts of interest.

      REFERENCES

        • Álvarez-Martín P.
        • Flórez A.B.
        • López-Díaz T.M.
        • Mayo B.
        Phenotypic and molecular identification of yeast species associated with Spanish blue-veined Cabrales cheese.
        Int. Dairy J. 2007; 17: 961-967
        • AOAC International
        Official Methods of Analysis.
        17th ed. AOAC International, Washington, DC2000
        • Bertuzzi A.S.
        • Walsh A.M.
        • Sheehan J.J.
        • Cotter P.
        • Crispie D.F.
        • Mcsweeney P.L.H.
        • Kilcawley K.N.
        • Rea M.C.
        Omics-based insights into flavor development and microbial succession within surface-ripened cheese.
        mSystems. 2018; 3 (29404426): e00211-e00217
        • Bintsis T.
        • Vafopoulou-Mastrojiannaki A.
        • Litopoulou-Tzanetaki E.
        • Robinson R.K.
        Protease, peptidase and esterase activities by lactobacilli and yeast isolates from Feta cheese brine.
        J. Appl. Microbiol. 2003; 95 (12807455): 68-77
        • Büchl N.R.
        • Seiler H.
        Yeast and molds: Yeasts in milk and dairy products.
        in: Fuquay J.W. Fox P.F. McSweeney P.L.H. Encyclopedia of Dairy Sciences. Elsevier Ltd., London, UK2011: 744-753
        • Caporaso J.G.
        • Kuczynski J.
        • Stombaugh J.
        • Bittinger K.
        • Bushman F.D.
        • Costello E.K.
        • Fierer N.
        • Peña A.G.
        • Goodrich J.K.
        • Gordon J.I.
        • Huttley G.A.
        • Kelley S.T.
        • Knights D.
        • Koenig J.E.
        • Ley R.E.
        • Lozupone C.A.
        • McDonald D.
        • Muegge B.D.
        • Pirrung M.
        • Reeder J.
        • Sevinsky J.R.
        • Turnbaugh P.J.
        • Walters W.A.
        • Widmann J.
        • Yatsunenko T.
        • Zaneveld J.
        • Knight R.
        QIIME allows analysis of high-throughput community sequencing data.
        Nat. Methods. 2010; 7 (20383131): 335-336
        • Carpino S.
        • Randazzo C.
        • Pino A.
        • Russo N.
        • Rapisarda T.
        • Belvedere G.
        • Caggia C.
        Influence of PDO Ragusano cheese biofilm microbiota on flavour compounds formation.
        Food Microbiol. 2017; 61 (27697162): 126-135
        • Caspia E.L.
        • Coggins P.C.
        • Schilling M.W.
        • Yoon Y.
        • White C.H.
        The relationship between consumer acceptability and descriptive sensory attributes in Cheddar cheese.
        J. Sens. Stud. 2006; 21: 112-127
        • Collins Y.F.
        • McSweeney P.L.H.
        • Wilkinson M.G.
        Lipolysis and free fatty acid catabolism in cheese: A review of current knowledge.
        Int. Dairy J. 2003; 13: 841-866
        • Curtin A.C.
        • Gobbetti M.
        • McSweeney P.L.H.
        Peptidolytic, esterolytic and amino acid catabolic activities of selected bacterial strains from the surface of smear cheese.
        Int. J. Food Microbiol. 2002; 76 (12051480): 231-240
        • Deeth H.C.
        • Touch V.
        Methods for detecting lipase activity in milk and milk products.
        Aust. J. Dairy Technol. 2000; 55: 153-168
        • Drake S.L.
        • Gerard P.D.
        • Drake M.A.
        Consumer preferences for mild Cheddar cheese flavors.
        J. Food Sci. 2008; 73 (19021820): S449-S455
        • Duru I.C.
        • Laine P.
        • Andreevskaya M.
        • Paulin L.
        • Kananen S.
        • Tynkkynen S.
        • Auvinen P.
        • Smolander O.-P.
        Metagenomic and metatranscriptomic analysis of the microbial community in Swiss-type Maasdam cheese during ripening.
        Int. J. Food Microbiol. 2018; 281 (29803134): 10-22
        • Edgar R.C.
        UPARSE: Highly accurate OTU sequences from microbial amplicon reads.
        Nat. Methods. 2013; 10 (23955772): 996-998
        • Fox P.F.
        • Guinee T.P.
        • Cogan T.M.
        • McSweeney P.L.H.
        Biochemistry of Cheese Ripening.
        in: Fox P.F. Guinee T.P. Cogan T.M. McSweeney P.L.H. Fundamentals of Cheese Science. Springer, New York, NY2000: 391-442
        • Fox P.F.
        • Wallace J.M.
        Formation of flavor compounds in cheese.
        Adv. Appl. Microbiol. 1997; 45 (9342826): 17-85
        • Gao M.L.
        • Hou H.M.
        • Teng X.X.
        • Zhu Y.L.
        • Hao H.S.
        • Zhang G.L.
        Microbial diversity in raw milk and traditional fermented dairy products (Hurood cheese and Jueke) from Inner Mongolia, China.
        Genet. Mol. Res. 2017; 16 (28290619): 1-13
        • Gemert L.J.V.
        Odour threshold values in water in odour thresholds.
        in: Compilations of Odour Threshold Values in Air, Water and Other Media. Oliemans Punter & Partners, Utrecht, the Netherlands2011: 157-289
        • Gori K.
        • Sørensen L.M.
        • Petersen M.A.
        • Jespersen L.
        • Arneborg N.
        Debaryomyces hansenii strains differ in their production of flavor compounds in a cheese-surface model.
        MicrobiologyOpen. 2012; 1 (22950022): 161-168
        • He W.
        • Chung H.Y.
        Exploring core functional microbiota related with flavor compounds involved in the fermentation of a natural fermented plain sufu (Chinese fermented soybean curd).
        Food Microbiol. 2020; 90 (32336369)103408
        • Holland R.
        • Liu S.Q.
        • Crow V.L.
        • Delabre M.L.
        • Lubbers M.
        • Bennett M.
        • Norris G.
        Esterases of lactic acid bacteria and cheese flavour: Milk fat hydrolysis, alcoholysis and esterification.
        Int. Dairy J. 2005; 15: 711-718
        • Jin H.
        • Mo L.
        • Pan L.
        • Hou Q.
        • Li C.
        • Darima I.
        • Yu J.
        Using PacBio sequencing to investigate the bacterial microbiota of traditional Buryatian cottage cheese and comparison with Italian and Kazakhstan artisanal cheeses.
        J. Dairy Sci. 2018; 101 (29753477): 6885-6896
        • Liu W.
        • Sun Z.
        • Zhang J.
        • Gao W.
        • Wang W.
        • Wu L.
        • Sun T.
        • Chen W.
        • Liu X.
        • Zhang H.
        Analysis of microbial composition in acid whey for dairy fan making in Yunnan by conventional method and 16s rRNA sequencing.
        Curr. Microbiol. 2009; 59 (19459001): 199-205
        • Mounier J.
        Microbial interactions in smear-ripened cheeses.
        in: Bora N. Dodd C. Desmasures N. Diversity, Dynamics and Functional Role of Actinomycetes on European Smear Ripened Cheeses. Springer, Cham, Switzerland2015: 155-166
        • Shannon P.
        • Markiel A.
        • Ozier O.
        • Baliga N.S.
        • Wang J.T.
        • Ramage D.
        • Ideker T.
        Cytoscape: a software environment for integrated models of biomolecular interaction networks.
        Genome Res. 2003; 13 (14597658): 2498-2504
        • Stefanovic E.
        • Kilcawley K.N.
        • Roces C.
        • Rea M.C.
        • O'Sullivan M.
        • Sheehan J.J.
        • McAuliffe O.
        Evaluation of the potential of Lactobacillus paracasei adjuncts for flavor compounds development and diversification in short-aged Cheddar cheese.
        Front. Microbiol. 2018; 9 (30026739)1506
        • Tamang J.P.
        • Watanabe K.
        • Holzapfel W.H.
        Review: Diversity of microorganisms in global fermented foods and beverages.
        Front. Microbiol. 2016; 7 (27047484): 377
        • Tian H.
        • Xu X.X.
        • Chen C.
        • Yu H.Y.
        Flavoromics approach to identifying the key aroma compounds in traditional Chinese milk fan.
        J. Dairy Sci. 2019; 102 (31477286): 9639-9650
        • van Mastrigt O.
        • Gallegos Tejeda D.
        • Kristensen M.N.
        • Abee T.
        • Smid E.J.
        Aroma formation during cheese ripening is best resembled by Lactococcus lactis retentostat cultures.
        Microb. Cell Fact. 2018; 17 (29973201): 104
        • Wang Y.
        • Guo J.
        • Huang A.
        Study on bacterial diversity in traditional sour whey of Yunnan province.
        J. Food Process. Preserv. 2018; 42e13553
        • Xue J.
        • Yang Y.
        • Wang Z.
        • Guo Y.
        • Shao Y.
        Bacterial diversity in Chinese Rushan cheese from different geographical origins.
        Front. Microbiol. 2018; 9 (30177922)1920
        • Shuangquan
        • Burentegusi
        • Yu B.
        • Miyamoto T.
        Microflora in traditional starter cultures for fermented milk, hurunge, from Inner Mongolia, China.
        Anim. Sci. J. 2006; 77: 235-241
        • Zapala M.A.
        • Schork N.J.
        Multivariate regression analysis of distance matrices for testing associations between gene expression patterns and related variables.
        Proc. Natl. Acad. Sci. USA. 2006; 103 (17146048): 19430-19435
        • Zheng X.
        • Li K.
        • Shi X.
        • Ni Y.
        • Li B.
        • Zhuge B.
        Potential characterization of yeasts isolated from Kazak artisanal cheese to produce flavoring compounds.
        MicrobiologyOpen. 2018; 7 (29277964): 533
        • Zheng X.
        • Liu F.
        • Li K.
        • Shi X.
        • Ni Y.
        • Li B.
        • Zhuge B.
        Evaluating the microbial ecology and metabolite profile in Kazak artisanal cheeses from Xinjiang, China.
        Food Res. Int. 2018; 111 (30007669): 130-136
        • Zheng X.
        • Liu F.
        • Shi X.
        • Wang B.
        • Li K.
        • Li B.
        • Zhuge B.
        Dynamic correlations between microbiota succession and flavor development involved in the ripening of Kazak artisanal cheese.
        Food Res. Int. 2018; 105 (29433268): 733-742