Unraveling Core Functional Microbiota in Traditional Solid-State Fermentation by High-Throughput Amplicons and Metatranscriptomics Sequencing.

Fermentation microbiota is specific microorganisms that generate different types of metabolites in many productions. In traditional solid-state fermentation, the structural composition and functional capacity of the core microbiota determine the quality and quantity of products. As a typical example of food fermentation, Chinese Maotai-flavor liquor production involves a complex of various microorganisms and a wide variety of metabolites. However, the microbial succession and functional shift of the core microbiota in this traditional food fermentation remain unclear. Here, high-throughput amplicons (16S rRNA gene amplicon sequencing and internal transcribed space amplicon sequencing) and metatranscriptomics sequencing technologies were combined to reveal the structure and function of the core microbiota in Chinese soy sauce aroma type liquor production. In addition, ultra-performance liquid chromatography and headspace-solid phase microextraction-gas chromatography-mass spectrometry were employed to provide qualitative and quantitative analysis of the major flavor metabolites. A total of 10 fungal and 11 bacterial genera were identified as the core microbiota. In addition, metatranscriptomic analysis revealed pyruvate metabolism in yeasts (genera Pichia, Schizosaccharomyces, Saccharomyces, and Zygosaccharomyces) and lactic acid bacteria (genus Lactobacillus) classified into two stages in the production of flavor components. Stage I involved high-level alcohol (ethanol) production, with the genus Schizosaccharomyces serving as the core functional microorganism. Stage II involved high-level acid (lactic acid and acetic acid) production, with the genus Lactobacillus serving as the core functional microorganism. The functional shift from the genus Schizosaccharomyces to the genus Lactobacillus drives flavor component conversion from alcohol (ethanol) to acid (lactic acid and acetic acid) in Chinese Maotai-flavor liquor production. Our findings provide insight into the effects of the core functional microbiota in soy sauce aroma type liquor production and the characteristics of the fermentation microbiota under different environmental conditions.

Introduction

The microbiota in natural environments is complex and often includes thousands of genera from a diverse range of species (Astudillo-Garcia et al., 2017). In industrial fermentation, the microbiota depends on the conditions and nutrient availability (Sekwati-Monang et al., 2012). In natural environments or under artificially controlled conditions, fermentation microbiota have the potential to use different raw materials to produce a variety of metabolites (Xu and Ji, 2012; Lu et al., 2016; Wang et al., 2016). They are commonly characterized by the presence of specific microorganisms, orderly microbial succession and an unusual functional shift (Azwar et al., 2014; Di Cagno et al., 2014; Kong et al., 2014; Qi et al., 2014). As an example of a typical process involving fermentation microbiota, traditional solid-state fermentation (SSF) generates fermented foods that have a high nutritional value and confer health benefits due to the involvement of multiple microorganisms and biochemical reactions (Giraffa, 2004; Kong et al., 2014; Tamang et al., 2016).

In recent years, the concept of core microbiota has emerged (Astudillo-Garcia et al., 2017). The core microbiota play important roles in natural ecological systems, human health and industrial processes (Schmitt et al., 2012; Zhang et al., 2015; Kable et al., 2016). However, the core microbiota that influence the quantity and quality of fermented foods in traditional SSF are poorly understood (Rantsiou and Cocolin, 2006; Cocolin et al., 2013). Traditionally, species were isolated and identified from food fermentation processes by culture-dependent and -independent techniques (Papalexandratou et al., 2013; Chen et al., 2014; Ricciardi et al., 2015). Although these studies gained some insight into the microbial populations, our knowledge about the core structural microbiota in fermented foods remain limited (Wu et al., 2013). With the development of high-throughput sequencing techniques, researches have focused on the structure of the core microbiota in traditional SSF processes (Ercolini, 2013; Aldrete-Tapia et al., 2014). For example, in cheese production, 10 fungal genera and 14 bacterial genera (taxa > 1%) are considered the core microbiota by 16S rRNA gene and internal transcribed spacer (ITS) amplicon sequencing (Wolfe et al., 2014). Similarly, in soy sauce fermentation, 22 frequent genera have been reported to be candidates of the core microbiota (Sulaiman et al., 2014).

Although the above-mentioned studies illustrate the structure of the core microbiota in traditional fermented foods, the correlation between the microbiota and metabolites during different processes remains unclear (Papagianni, 2014). To address this, research has focused on the function of the core microbiota in traditional food fermentation processes (Garrote et al., 2015). Recently, the high-throughput meta-omics sequencing technique has been demonstrated to be informative in discovering microbial succession patterns and the function shift that are characteristic of the fermentation microbiota and are crucial for their function (Bokulich and Mills, 2012; Solieri et al., 2013). For example, in cheese production, metatranscriptomic sequencing demonstrated that some lactic acid bacteria (LAB) are considered the main functional contributors in the processes of proteolysis, lipolysis, and amino acid/lipid catabolism (De Filippis et al., 2016). Similarly, metagenomic sequencing revealed that in pickle production, the genera Novosphingobium and Sphingomonas are starter microorganisms, which are subsequently replaced by the genera Leuconostoc, Lactobacillus, and Weissella. In addition, metatranscriptomic sequencing analysis revealed that the genes related to carbohydrate transport, hydrolysis, and pyruvate metabolism are actively expressed and highly improve the metabolic capacity of lactic acid production (Jung et al., 2013).

Chinese soy sauce aroma type (Maotai-flavor) liquor production (Supplementary Figure S1) is a spontaneous and repeated batch process thought to originate from ancient fermented beverages (McGovern et al., 2004). In addition to the common characteristics of traditional SSF (Thomas et al., 2013), soy sauce aroma liquor production also displays some specific characteristic features influenced by endogenous factors (such as temperature, pH, moisture, ethanol, lactic acid, and acetic acid) in the pit fermentation process (Table ). The fermentation mechanism required for Chinese soy sauce aroma liquor production therefore involves complex microbiota and metabolites (Xu et al., 2010). Moreover, the fermentation microbiotas use cereals (major sorghum and wheat) as raw materials to form a variety of flavor components through a series of biochemical reactions, including various carbon, nitrogen, and sulfur metabolites (Zhu et al., 2007; Xu et al., 2010; Wu et al., 2012). However, the functional correlation between the core microbiota and important metabolites remains to be established in soy sauce aroma liquor production. More specifically, microbial succession and the functional shift in the core microbiota have not been clarified to date. Therefore, in this study, we employed several high-throughput sequencing technologies (16S rRNA gene amplicon sequencing, ITS amplicon sequencing and metatranscriptomics sequencing) to explore the structure and function of core microorganisms in fermentation microbiota (Schoch et al., 2012; Fierer et al., 2013; Glass et al., 2015). In addition, we combined ultra-performance liquid chromatography (UPLC) and headspace-solid phase microextraction-gas chromatography-mass spectrometer (HS-SPME-GC-MS) to explore the fluctuations in major flavor components (Cordovez et al., 2015; Zothanpuia et al., 2017). Based on this information, we explored the correlations between the core microbiota and important metabolites.

Materials and Methods

Sample Collection

All samples were collected in Xishui County (28.14 N, 106.18 E), Guizhou Province, China. The fermentation process was carried out in four distinct phases: starter (Daqu) making, stacking fermentation, pit fermentation, and distillation. Samples were taken during the pit fermentation phase of liquor production (Supplementary Figure S1). Raw materials (15 tons) were mixed with sorghum and wheat in a pit (3 m × 2.5 m × 4 m). In addition, they were also selected from different locations (points A, B, and C) in the pit. Samples were collected at 5-day intervals in duplicate, one for genomic DNA extraction and total RNA extraction, the other for flavor component analysis. Samples from the three points (1, 2, and 3) were mixed in the same layer to form one sample to reduce the volatility of samples before extraction and analysis.

DNA Extraction, 16S rRNA Gene Amplicon Sequencing, and ITS Amplicon Sequencing

All samples were treated with sterile phosphate-buffered saline (PBS, 0.1 mol/L), then centrifuged at 30 × g for 10 min. The supernatant was then centrifuged again at 10,000 rpm for 10 min. The precipitated cells were milled with liquid nitrogen and genomic DNA was extracted using sodium laurate buffer (sodium laurate 10 g/L, Tris-HCl 0.1 mol/L, NaCl 0.1 mol/L, ethylenediaminetetraacetic acid (EDTA) 0.02 mol/L) containing phenol:chloroform:isoamyl alcohol (25:24:1). Bacterial V3-V4 and the fungal ITS1 of the rRNA were amplified using forward primers (5′-GTACTCCTACGGGAGGCAGCA-3′, 5′-CTTGGTCATTTAGAGGAAGTAA-3′) and the reverse primer (5′-GTGGACTACHVGGGTWTCTAAT-3′, 5′-TGCGTTCTTCATCGATGC-3′), respectively. Amplicons were then sequenced using an Illumina MiSeq platform (Illumina, San Diego, CA, United States) at AuwiGene Technology Co., Ltd., Beijing, China.

16S rRNA Gene and ITS Amplicon Data Analyses

Raw sequences were processed using ‘quantitative insights into microbial ecology’ software (QIIME, version 1.8) (Caporaso et al., 2010). Briefly, sequences with ambiguous bases (‘N’) were removed, and reads were quality-filtered to remove bases with Q < 0.02 by Trimmomatic (version 0.32) (Bolger et al., 2014). Then, overlapping reads merged by Fastq-join (version 1.8) and primer sequences were removed and complete assembled reads were used for further analysis. The sequences were clustered into operational taxonomic units (OTUs) at 97% sequence similarity and chimeras were removed by USEARCH (version 8.1.1861) (Angly et al., 2016). In addition, OTUs of 16S rRNA gene amplicons were mapped to Silva (version 119) (Zhu et al., 2015) and OTUs of ITS amplicons were mapped to Unite (version 7.0) (Kõljalg et al., 2013). Fungal OTUs were manually corrected using the CBS database (Centraalbureau voor Schimmelcultures, Holland) (Hallin and Ussery, 2004) and bacterial OTUs were manually corrected by Ezbiocloud (The Jongsik Chun Lab. and Chunlab, Inc., South Korea) (Gulati et al., 2015). The amplicon databases have been submitted to NCBI SRA and are available under accession numbers SRR5298103 and SRR5298193.

Total RNA Extraction and Metatranscriptomic Sequencing

The initial steps of total RNA extraction were the same as those for DNA extraction. After centrifugation, precipitated cells were milled with liquid nitrogen and total RNA was extracted with sodium laurate buffer (sodium laurate 10 g/L, Tris-HCl 0.1 mol/L, NaCl 0.1 mol/L, EDTA 0.02 mol/L) containing TRIzol (Sigma-Aldrich, St. Louis, MO, United States). Ribo-ZeroTM rRNA Removal Kits (Bacteria) and Ribo-ZeroTM Magnetic Gold Kit (Yeast) (Epicentre, San Diego, CA, United States) were used to remove rRNA from the total RNA. Metatranscriptomic libraries were constructed according to the NEBNext® UltraTM RNA Library Prep Kit (Illumina, New England Biolabs), and were then sequenced on an Illumina Hiseq 2500 platform at the Biomarker Technologies Co., Ltd., Beijing, China (Ziesemer et al., 2015).

Metatranscriptomic Data Analyses

Raw metatranscriptomic data were processed by removing the rRNA sequences and low quality reads (Q < 0.02). Trinity (version 2.1.1) (Haas et al., 2013) was employed for de novo transcriptomic assemblies. Using DIAMOND (version 0.8.34) (Buchfink et al., 2015) and inputting results into MEGAN (version 5.11.3) (Huson et al., 2007), we obtained microbial taxonomy information. Then, we compared high quality reads to the non-redundant protein database (Liu et al., 2013) and the KEGG database (Kanehisa et al., 2004) to obtain functional annotation information. Heatmap Illustrator (HemI, v.1.0) (Deng et al., 2014) was used to illustrate different gene expression patterns in the main metabolic pathways. SPSS (version 19.0, Web Atlas, Paris, France) and HemI were used to establish correlations generated by the metatranscriptomic database analysis (P < 0.05). Other statistical analyses and graphics were executed using Microsoft® Excel and R software (version 3.3.2). The metatranscriptomic database has been submitted to NCBI SRA and is accessible under SRR5306242.

Flavor Components Analyses

The flavor components in samples were analyzed by HS-SPME-GC-MS and UPLC. For HS-SPME-GC-MS analysis, samples were treated as follows: 4 g of each sample was selected and placed in 50 ml centrifuge tubes. Firstly, 40 ml of sterile distilled water and 0.4 g of CaCl2 were added to the sample and mixed. Secondly, tubes were placed in ice water for 30 min. Finally, the supernatant was collected by centrifuging at 8000 rpm for 10 min. The automatic headspace sampling system (Multi-Purpose Sample MPS 2 with a SPME adapter) (GERSTEL Inc., Baltimore, MD, United States) with a 50/30 μm DVB/CAR/PDMS fiber (Supelco Inc., Bellefonte, PA, United States) was used in SPME. Samples were preheated for 5 min and extracted for 45 min at 50°C. The Agilent 6890N GC coupled with an Agilent 5975 Mass Selective Detector was used in GC-MS. The GC-MS conditions were: a starting temperature of 50°C (held for 2 min), then increased to 230°C at a rate of 4°C/min and held at 230°C for 15 min. For UPLC analysis, we used the same method as previously described without 0.4 g CaCl2. The UPLC system consisted of an Acquity UPLC (Waters, Milford, MA, United States) coupled with an Acquity UPLC Tunable UV Detector (Waters). The column was an Acquity UPLC HSS T3 column (2.1 mm × 100 mm, 1.8 μm) and the eluent was NaH2PO4. The mass spectra of unknown flavor components were compared with those in the Wiley 275.L (Agilent) database.

Results

Major Alcohols and Acids in the Flavor Components of Soy Sauce Aroma Liquor Production

To examine the flavor components of liquor production, important non-volatile and volatile components were qualified and quantified by UPLC and HS-SPME-GC-MS across samples. A total of 61 major components were identified, including 10 alcohols (mainly ethanol, 1-hexanol and 3-methyl-1-butanol), 12 acids (mainly lactic acid, acetic acid and succinic acid) and 28 esters (mainly ethyl acetate, isobutyl hexanoate, and diethyl succinate) (Figure and Supplementary Figure S2). It was apparent that ethanol was the most abundant metabolite (18.87 ± 14.98 – 43.43 ± 8.21 g/kg fermented grain), followed by lactic acid (28.05 ± 4.75 – 36.20 ± 6.20 g/kg fermented grain) and acetic acid (7.40 ± 2.05 – 23.98 ± 7.18 g/kg fermented grain) (Supplementary Figure S2). The relative abundance of alcohols increased from days 0 to 15, with the average ranging from 24.82 to 32.35%. Among the alcohols, the relative abundance of ethanol increased from days 5 to 15, with the average ranging from 18.02 to 28.18% (Figure ). By contrast, the relative abundance of acids increased from days 15 to 30, ranging from 62.70 to 70.34% (Figure ). Among the acids, the relative abundance of lactic acid and acetic acid increased from days 15 to 30, with the average ranging between 34.81 and 42.62% (Figure ). The production rate of the major flavor components is shown in Supplementary Figure S3. The production rate of alcohols decreased from days 0 to 10, with the average ranging from 0.85 to 1.54 for alcohols (Figure ), and more specifically, from -0.36 to 1.77 for ethanol (Figure ). By contrast, the production rate of acids increased from days 20 to 30, with the average ranging from -0.34 to 4.91 for acids generally (Figure ) and -0.01 to 2.17 for lactic acid and acetic acid specifically (Figure ).

The average content of flavor components across all the samples. Samples were sorted based on the fermentation time. (A) Average content of alcohols and acids across all the samples (n = 18, each bar n = 3). (B) Average relative abundance of major alcohols and acids across all the samples (n = 18, each bar n = 3). (C) The production rate of major alcohols and acidsacross all the samples (n = 18). (D) The production rate of ethanol and lactic acidacross all the samples (n = 18).

Endogenous Factors Distinctly Reflected Microbial Succession in the Fermentation Microbiota

To understand the composition of the microbiota, OTU reads were selected to depict the bacterial and fungal microbiota in liquor production (Stegen et al., 2016). A total of 712,610 high quality reads were generated from bacterial 16S rDNA V3–V4 sequences in 18 samples, ranging from 15,272 to 106,976 reads. We also obtained 912,794 high quality reads from fugal ITS1 sequences in 18 samples, ranging from 27,240 to 77,221 reads. In addition, 16S rRNA gene amplicons sequences clustered into 522 OTUs and ITS amplicons sequences clustered into 552 OTUs at a 97% similarity level. Bacterial sequence numbers were normalized to 12,426 reads and fungal sequence numbers were normalized to 27,470 reads for describing the composition characteristics (Supplementary Tables S1, S2). The Good’s coverage estimator (Long et al., 2016a) revealed that 99.55–99.74% of bacterial and fungal OTUs were obtained in all samples, suggesting that major bacterial and fungal OTUs had been captured. The bacterial Chao1 richness estimator (Long et al., 2016b) decreased with a range from 296.22 ± 23.22 to 160.07 ± 3.03 at different sampling times, whereas the fungal Chao1 richness estimator increased with a range from 153.91 ± 37.31 to 262.16 ± 4.60 (Supplementary Figure S3 and Tables S1, S2). In addition, the bacterial Shannon diversity index (Chen et al., 2015) decreased with a range from 4.54 ± 0.69 to 0.26 ± 0.04 at different sampling times, whereas the fungal Shannon diversity index increased with a range between 1.77 ± 0.12 and 3.09 ± 0.55 (Supplementary Figure S3 and Tables S1, S2).

Sequence similarity thresholds (80 to 100%) were used to define the OTUs. The range of similarity cut-offs, 80, 90, 97, and 99%, nominally estimated phylum, class, genus, and species, respectively (Whiteson et al., 2014). Firmicutes and Ascomycota were dominant at the phylum level, comprising about 90% of bacteria and 95% of fungi (Figure ). The largest groups of bacterial genera were Lactobacillus, Kroppenstedtia, Acidithiobacillus, Acetobacter, and Pediococcus, whereas the largest groups of fungal genera were Pichia, Aspergillus, Saccharomycopsis, Penicillium, and Zygosaccharomyces (Figure ). The genus Lactobacillus (on average accounting for 44.86 to 98.83% of bacteria) was the most abundant bacterial genus. The genus Pichia (on average accounting for 91.34 to 64.59% of fungi) was the most abundant fungal genus in each sample. The overall abundances of bacterial genera showed that the average relative abundance of the genus Lactobacillus increased but the abundance of the genera Kroppenstedtia (12.92–0.12%) and Acidithiobacillus (7.55–0.12%) decreased at different sampling times. Simultaneously, the overall relative abundances of fungal genera also changed. The relative abundance of the genus Pichia decreased but the genera Aspergillus (0.50–8.42%) and Saccharomycopsis (0.35–10.61%) increased at different sampling times.

Microbiota and major flavor components analysis across all the samples. Average bacterial (A) and fungal (B) distribution at the genus-level in microbiota based on 16S rRNA gene and ITS amplicons (n = 18, each bar n = 3). (C) Amplicons analysis represented the similarities of microbial compositions based on principal component analysis (PCA). (D) Ethanol, lactic acid, and acetic acid production had different stages. Bars represented mean (±SE). Asterisk indicates significant differences by Mann–Whitney U test (P < 0.05). (E) Plots of PC1 versus four factors showed that ethanol and lactic acid had significant correlation with microbiota in different stages by Spearman correlation coefficient (P < 0.05).

Based on 16S rRNA and ITS gene amplicons, Unifrac distance-based weighted principal coordinates analysis was conducted to evaluate similarities and differences in the microbiota (prokaryotic and eukaryotic microbiota) among the different samples (Figure ). Although there was a slight difference in the microbiota among samples, similar clustering patterns were observed. Samples from days 15, 20, and 25 formed clusters. Samples from days 15 and 20 were distributed similarly to samples from day 25, whereas samples from days 5 and 30 were dispersed within other samples. The differences between samples on days 5 and 30 were larger than the differences between other samples. To demonstrate the major factors associated with microbial succession at the different stages of liquor production, multiple endogenous factors (temperature, pH, moisture, ethanol, lactic acid, and acetic acid) were selected and analyzed to reveal their connection with the microbiota (Table ). We found that the temperature and acetic acid content of samples showed no significant difference (P > 0.05) as determined by a Mann–Whitney U test (Supplementary Figure S4). However, the ethanol content of samples from days 0 to 10 differed significantly (P < 0.05) compared with samples from days 10 to 30. The moisture levels of samples also differed significantly for days 0 to 15 (P < 0.05) compared with samples from days 15 to 30 (Figure ) as determined by a Mann–Whitney U test. The Mann–Whitney U test also revealed a significant difference in the lactic acid content of samples from days 0 to 20 (P < 0.05) compared with samples from days 20 to 30 and pH of samples from days 0 to 15 was significantly different (P < 0.01) compared with that of samples from days 15 to 30 (Figure ). Ethanol and moisture were the best indicators in the early stage of liquor production, with principal coordinate 1 (PC1) being significantly associated with the ethanol content of samples from days 5 to 10 (R2 = 0.45, P < 0.05) and the moisture content of samples from days 25 to 30 (R2 = 0.56, P < 0.05) (Figure ). Lactic acid and pH were the best indicators in the late stage of liquor production, with principal coordinate 1 (PC1) being significantly correlated with the lactic acid content of samples from days 25 to 30 (R2 = 0.61, P < 0.05) and the pH of samples from days 25 to 30 (R2 = 0.52, P < 0.01) (Figure ).

Core Microbiota Showed a Significant Correlation with Major Endogenous Factors at Different Stages of Production

To further illustrate the correlation between the core microbiota and major endogenous factors, amplicons and the major endogenous factor database were combined to explore positive and negative correlations between the microbiota and major endogenous factors using Gephi (version 0.9.1) (Figure and Supplementary Tables S3, S4). In total, 33 pairs of significant and robust relationships (edges) were identified from 33 genera (nodes) (P < 0.05). Gephi network analysis (Jacomy et al., 2009) detected 7 positive and 26 negative correlations using Spearman’s correlation coefficient (Borkowf, 2002), revealing that the genera Sulfobacillus, Clostridium, and Aspergillus negatively correlated with ethanol (P < 0.05) in the early stage of liquor production (Stage I). The genera Proteiniphilum and Schizosaccharomyces showed a positive correlation with ethanol (P < 0.05) in Stage I. Bacterial genera (Lactobacillus and Oscillibacter) and fungal genera (Saccharomyces, Schizosaccharomyces, Purpureocillium, Pseudeurotium, Candida, Fusarium, Rhizomucor, Discosia, and Acremonium) displayed a negative correlation with moisture (P < 0.05) in Stage I. By contrast, in the late stage of liquor production (Stage II), the genus Corynebacterium showed a negative correlation with lactic acid, whereas the genus Lactobacillus showed a positive correlation with lactic acid (P < 0.05). The bacterial genus Lactobacillus and various fungal genera (Penicillium, Zygosaccharomyces, Pleurotheciella, Naumovozyma, Kazachstania, Trichosporon, Pseudeurotium, Veronaea, Discosia, and Pseudogymnoascus) displayed a negative correlation with pH (P < 0.05) in Stage II. The bacterial genera (Kroppenstedtia and Corynebacterium) and the fungal genera (Saccharomycopsis and Schizosaccharomyces) showed a positive correlation with pH (P < 0.05) in Stage II.

Correlation network between microbiota and endogenous factors. Spearman correlation coefficient depicted significant negative and positive correlations (P < 0.05) between microbiota and endogenous factors in Stage I (A) and Stage II (B) (relative abundance > 0.01%). Edge thickness represented the value. Edge color represented the positive (red) or negative (black) correlation. Edge length and point sizes had no meaning for the present study.

Variations in Pyruvate Metabolism Reflected a Functional Shift in the Fermentation Microbiota

To explore microbial function in Chinese soy sauce aroma liquor production, the KEGG database was used to explore the overall metabolic capability of microbiota in situ by HemI (Supplementary Figure S5). Among 64,349 unigenes (i.e., transcripts that originate from the same transcription locus) related to major KEGG metabolic pathways, 9,003 unigenes (13.99%) were attributed to metabolic pathways involved in maintaining the basic survival of microorganisms (ribosomes, cell cycle, and RNA transport), whereas the others (86.01%) were attributed to different types of flavor components. These unigenes were attributed to carbon metabolic pathways (27,788, 43.18%), nitrogen metabolic pathways (20,667, 32.12%) and other metabolic pathways (6,891, 10.71%). The most abundant type of metabolism in carbon metabolic pathways was glycolysis/gluconeogenesis (1,777, 2.76%), followed by pyruvate metabolism (1,480, 2.30%) and the tricarboxylic acid cycle (TCA cycle) (1,433, 2.23%).

To explain the correlations between core microorganisms and major endogenous factors in the microbiota, the KEGG and non-redundant protein databases were combined to illustrate pyruvate metabolism related to ethanol, lactic acid and acetic acid production. In our study, the genera Pichia, Saccharomyces, Schizosaccharomyces, Zygosaccharomyces, and Lactobacillus were the dominant functional contributors to the fermentation microbiota (Supplementary Figure S6). Functional genes related to the genera Pichia, Saccharomyces, Schizosaccharomyces, Zygosaccharomyces, and Lactobacillus were expressed at high levels in pyruvate metabolism (Figure ). Lactate dehydrogenase (cytochrome) (LDHC, K00101, and K00102) and pyruvate dehydrogenase E1 component (PDH, K00161, K00162 and K00163) peaked in expression in the genus Pichia during pyruvate metabolism. Pyruvate kinase (PYK, K00873), aldehyde dehydrogenase (e.g., ADH, K00128, K00129) and alcohol dehydrogenase (e.g., ADH, K00001, and K00002) showed the highest levels of gene expression in the genera Schizosaccharomyces and Lactobacillus during carbohydrate metabolism. Although the genera Pichia, Schizosaccharomyces, Saccharomyces, Zygosaccharomyces, and Lactobacillus shared the same pathway, high level gene expression occurred in different samples. Moreover, to maintain the high-level expression of pyruvate metabolism, PYK became the predominant KEGG gene in the genera Pichia, Schizosaccharomyces, Saccharomyces, Zygosaccharomyces, and Lactobacillus. Consistent with the need to use lactic acid and ensure the conversion of pyruvate to acetaldehyde and acetyl-CoA, our results showed that LDH, PDH and pyruvate decarboxylase (PDC, K01568) became the major KEGG genes in the genus Pichia. KEGG expression genes related to the genus Schizosaccharomyces showed high metabolic capacity via PDH, acetyl-CoA synthetase (ACS, K01895, K01067, etc.), ALDH and ADH, which converted pyruvate to acetaldehyde through acetyl-CoA and acetate, and then converted acetaldehyde to ethanol in the fermentation process. In addition, KEGG genes related to genus Lactobacillus showed high metabolic capacity via LDH, xylulose-5-phosphate/fructose-6-phosphate phosphoketolase (XFP, K01621), phosphate acetyltransferase (PTA, K00625) and ADH in liquor production.

Relative abundance of KEGG expression genes in pyruvate metabolism related to genera Pichia, Schizosaccharomyces, Saccharomyces, Zygosaccharomyces, and Lactobacillus in different stages (Stage I: samples from days 5 to 15 and Stage II: samples from days 20 to 30 (n = 6). Color depth represented the proportion of total reads mapping to KEGG metabolic pathway in all KEGG annotated reads.

Furthermore, as a major metabolic pathway, pyruvate metabolism became increasingly active during the temporal course, as determined by the FPKM (Fragments per Kilobase of transcript per Million mapped reads) method (Guo et al., 2013) (Figure ). ADH in molds (genera Aspergillus and Penicillium), yeasts (genera Pichia, Schizosaccharomyces, Saccharomyces, and Zygosaccharomyces) and LAB (genus Lactobacillus) displayed different metabolic capacities to form ethanol in different samples, which explained the variation in ethanol content between samples during production (Figures ). LDHC could convert lactic acid to pyruvate and its gene expression levels were obviously downregulated from 105.05 to 1.20 in mold (genus Byssochlamys) and yeasts (genera Pichia, Zygosaccharomyces, and Saccharomyces) from days 10 to 30 (Figure ). By contrast, LDH could convert pyruvate to lactic acid and its gene expression levels were upregulated from 23.86 to 830.18 in yeast (genus Schizosaccharomyces) and LAB (genus Lactobacillus) from days 5 to 30 (Figure ). The lactic acid content (Figure ) correlated with the functional shift in lactic acid production. PYC and malate dehydrogenase (oxaloacetate-decarboxylating) (e.g., MAE, K00027, K00028) could convert pyruvate in the TCA cycle and their gene expression levels fluctuated in different samples. XFP could convert glucose to acetyl phosphate without pyruvate metabolism and expression of this gene was upregulated from 556.43 to 13,668.3 in genus Lactobacillus from days 5 to 30 (Figure ). Conversely, the acetic acid content fluctuated over time (Figure ).

KEGG expression genes in pyruvate metabolism associated with the measured content of ethanol, lactic acid, and acetic acid in different stages (n = 6). (A) FPKM (Fragments per Kilobase of transcript per Million mapped reads) of KEGG expression genes related to pyruvate. (B) FPKM of KEGG expression genes related to metabolic pathway from acetaldehyde to ethanol. (C) Ethanol content from Stage I to Stage II. (D) FPKM of KEGG expression genes related to metabolic pathway from lactic acid to pyruvate. (E) FPKM of KEGG expression genes from pyruvate to lactic acid. (F) Lactic acid content from Stage I to Stage II. (G) FPKM of KEGG expression genes from pyruvate to acetyl-phosphate. (H) Acetic acid content from Stage I to Stage II.

Discussion

In Chinese soy sauce aroma liquor production, the core microorganisms play critical roles in the characteristics of the flavor metabolites. Compared with previous studies (Wu et al., 2013; Chen et al., 2014), we found that 11 bacterial genera and 10 fungal genera were the core microbiota in the microbial succession of liquor production, which affected the production of flavor components (Figures ). Compared with previous studies (Wang X. et al., 2014), we found that samples with a high alcohol content differed from those with a high acid content in the production process. The alcohol content increased quickly in samples from days 0 to 10, whereas the acid content increased rapidly from days 20 to 30 (Figure ). In addition, the production rate of alcohol content decreased quickly in samples from days 0 to 15, whereas the acid content increased quickly at samples from days 15 to 30 (Figures ). There was a significant difference (P < 0.05) between the ethanol and lactic acid contents at different stages of liquor production (Figure and Supplementary Figure S4). Ethanol and moisture level significantly correlated with the fermentation microbiota in the early stage of production, whereas lactic acid and pH significantly correlated with the fermentation microbiota in the late stage of production (Figure ). In a previous study, to illustrate the variety of structures and functions within the fermentation microbiota, the fermentation process was separated into several stages (Tao et al., 2014). According to this report, liquor production could be separated into two distinct stages. Stage I was a period of high-level alcohol production in the samples from days 0 to 15. During this stage, the ethanol content increased quickly from 18.87 ± 14.98 g/kg to 37.03 ± 16.36 g/kg fermented grain. The species richness and diversity within the bacteria microbiota decreased, whereas diversity within the fungal microbiota increased (Supplementary Tables S1, S2). Stage II was a period of high-level acid production in the 30.59 ± 6.37 g/kg to 36.20 ± 6.21 g/kg fermented grain. The species richness and diversity of the bacterial microbiota decreased, whereas that of the fungal microbiota remained unchanged during the second stage (Supplementary Tables S1, S2).

As a result of the complex endogenous factors present in a closed system (i.e., the pit fermentation process), it is a challenge to link the core microbiota to major flavor metabolites (such as alcohols and acids) in liquor production (Supplementary Figure S1). Many studies had described the distribution of species in the microbiota and possible correlations between the microbiota and products in different environmental conditions, but these studies have not identified the critical metabolic pathways involved (Graham et al., 2014; Wang Z. et al., 2014). Therefore, in our study, amplicons and metatranscriptomic data were combined to investigate potential correlations between microorganisms, functions and products. In our study, the core bacteria (genera Lactobacillus and Kroppenstedtia) and core fungi (genera Saccharomyces, Schizosaccharomyces, Zygosaccharomyces, and Saccharomycopsis) were significantly correlated with major endogenous factors at different stages of liquor production (P < 0.05) (Figure ). Moreover, metatranscriptomic analysis indicated the relative abundance of unigene expression related to glycolysis/gluconeogenesis, the TCA cycle and pyruvate metabolism increased from 6.79 to 8.24% in the microbiota (Supplementary Figure S5). In addition, KEGG data illustrated that the bacterial genus (Lactobacillus) and fungal genera (Schizosaccharomyces, Pichia, Zygosaccharomyces, and Saccharomyces) were the core functional microbiota in liquor production (Supplementary Figure S6). These results showed that yeasts (genera Schizosaccharomyces, Pichia, Zygosaccharomyces, and Saccharomyces) and LAB (genus Lactobacillus) were dominant functional contributors to major alcohol and acid production (ethanol, lactic acid, and acetic acid production).

Accordingly, pyruvate metabolism related to ethanol, lactic acid, and acetic acid production in yeasts (genera Pichia, Saccharomyces, Schizosaccharomyces, and Zygosaccharomyces) and LAB (genus Lactobacillus) was selected to illustrate the functional shift (functional conversion from ethanol to lactic acid production) in different stages of liquor production (Figures ). Yeasts (genera Pichia, Saccharomyces, Schizosaccharomyces, and Zygosaccharomyces) in Stage I expressed high levels of genes related to pyruvate metabolism. The genera Pichia, Saccharomyces, and Zygosaccharomyces were facultative aerobic microorganisms that were separated from food or wines (Pitt and Hocking, 2009). They were able to grow aerobically or anaerobically at low pH and had the ability to limit ethanol and acetic acid production (Kawahata et al., 2006; Yuangsaard et al., 2013; Lindahl et al., 2016). During this stage, the genera Pichia, Saccharomyces, and Zygosaccharomyces converted lactic acid to pyruvate. Moreover, they even converted pyruvate to acetyl-CoA, acetaldehyde in liquor production (Figures ). Acetyl-CoA is an important intermediate metabolite from the conversion of pyruvate to acetic acid under aerobic or anaerobic conditions (Edirisinghe et al., 2016). The genes PDH E1, ACS, ALDH, and ADH encode a series of consecutively used enzymes that convert pyruvate to ethanol via acetyl-CoA, acetate, and acetaldehyde. With increasing levels of the genus Schizosaccharomyces as the main contributor under aerobic conditions and high moisture, ethanol production was the major microbial function in Stage I (Figure ). At this stage, the ACS, PYC, and ALDH genes were not expressed in the genus Lactobacillus (Figure ). Although the genes related to pyruvate oxidase (POXL) might be expressed at high levels, but it had not activity at anaerobic conditions in liquor production. Generally, enzymes related to pyruvate metabolism were highly expressed in yeasts, but were expressed at low-level in the genus Lactobacillus (Figure ). In addition, the genus Lactobacillus in Stage II expressed high levels of the genes related to pyruvate metabolism. The genus Lactobacillus comprises facultative anaerobic microorganisms identified from cheeses or wines (Lahtinen et al., 2012; Su et al., 2014). This genus was able to grow anaerobically at low pH and had the ability to adapt to a large temperature variation (10–45°C) (Fröhlich-Wyder et al., 2015; Johanningsmeier and McFeeters, 2015; Wang et al., 2015). During this stage, the genus Lactobacillus increased lactic acid productivity and maintained ethanol productivity simultaneously (Figure ). LDH in the genus Lactobacillus converted pyruvate to lactic acid, while at the same time, XFP in the genus Lactobacillus converted glucose to acetyl phosphate directly in heterofermentative fermentation (Figures ). With the genus Lactobacillus being the main contributor under low pH, lactic acid production was the major microbial function in Stage II (Figure ).

It was intriguing to understand why the genera Schizosaccharomyces and Lactobacillus were the major functional contributors among the core microbiota. With the initiation of ethanol production in Stage I, the resulting high ethanol content may inhibit some core bacteria (genera Bacillus, Acetobacter, and Pediococcus) in the microbiota (Araque et al., 2013; Kim et al., 2013; Sowards et al., 2014). Similarly, with the increase in lactic acid content in Stage II, the high content of lactic acid may inhibit some core yeasts (genera Pichia, Schizosaccharomyces, and Zygosaccharomyces) suppressing ethanol production by these genera in the microbiota (Crowley et al., 2013; Su et al., 2014). Our study found that the genes of four key enzymes in heterolactic fermentation and three key enzymes in homolactic fermentation were expressed (Supplementary Figure S7). Besides, the gene expressions of four key enzymes in heterolactic fermentation were much higher than three key enzymes in homolactic fermentation. This phenomenon showed that the most Lactobacillus species had ability to produce lactic acid, ethanol and acetic acid. Previous research showed that L. acidophilus, L. delbrueckii, L. helveticus, and L. salivarius were homofermentative lactobacilli. Except these, most Lactobacillus genera were facultative or obligatory heterofermentative lactobacilli (Lahtinen et al., 2012; Gänzle, 2015). In Supplementary Figures S8, S9, homofermentative and heterofermentative metabolism in the genus Lactobacillus showed that many lactobacilli played their significant function in two-type metabolism and most lactobacilli were facultative or obligatory heterofermentative lactobacilli in Chinese soy sauce aroma liquor production. Therefore, the genus Lactobacillus had developed several survival systems in various acidic environments to prevent cell damage resulting from acid stress by the acid product during the fermentation process (Liu et al., 2015). In addition, a variety of heterofermentative Lactobacillus genera had characteristics about rapid growth (Tamarit et al., 2015) and effective carbohydrate metabolism (Gänzle and Ripari, 2016). Moreover, some homofermentative Lactobacillus genera could produce antibacterial peptides and bacteriocins that enable these organisms to compete against other bacteria in the microbiota, which confers a competitive advantage over other microorganisms in many food industries by reducing the risks of food poisoning and protecting human health against food spoilage (Messaoudi et al., 2013; Garrote et al., 2015). Furthermore, the genus Lactobacillus exerts antifungal effects on food or feed-borne filamentous fungi and yeasts in traditional SSF, such as kimchi fermentation (Ryu et al., 2014). For the reasons stated above, the genus Lactobacillus became the predominant microorganism in liquor production and the major lactic acid and acetic acid producer under strict anaerobic conditions during liquor production. However, the acetic acid content decreased between Stages I and II (Figure ), for which there are several possible reasons. Firstly, low expression of the acetate kinase (ACK, K00925) and ACS gene might indicate that it encodes the key regulatory enzyme in acetic acid production. Secondly, ethyl acetate was the most abundant ester among the flavor components, the synthesis of which required large amounts of ethanol and acetic acid, thus decreasing the acetic acid content (Lin et al., 2012). Therefore, acetic acid was still responsible for acid production in Stage II and the genus Lactobacillus was still the microbial indicator of the major functional shift from Stage I to Stage II in liquor production. As suggested above, a variety of core functional microbiota can be attributed to the different stages of liquor production. Although 11 bacterial and 10 fungal genera were recognized as the core microbiota, only the genera Schizosaccharomyces and Lactobacillus possessed the metabolic capacity to form ethanol, lactic acid and acetic acid in the flavor metabolites.

In our study, we unveiled the structure and function of the core microbiota in Chinese soy sauce aroma type liquor production by high-throughput sequencing technologies and identified the type and quantity of flavor components by UPLC and HS-SPME-GC-MS. Amplicon analysis revealed fundamental information about microbial succession and the positive correlations between the core microbiota structure and major endogenous factors. Furthermore, metatranscriptomic database analysis demonstrated a large amount of information regarding the variety of core functional contributors in the fermentation microbiota, thus making it possible to correlate the core functional microbiota with important flavor metabolites. Finally, we identified many microorganisms as members of the core structural microbiota that had the metabolic capacity to perform specific functional roles. Only a few microorganisms, constituting the core functional microbiota, played a critical role in driving the variety of microbiota. According to the best of our knowledge, this is the first report to analyze the core microbiota in Chinese soy sauce aroma type liquor production systemically and to define the core functional microbiota in different fermentation processes. These findings enrich our knowledge of the core functional microbiota in the fermentation process under artificial conditions or natural environments.

Author Contributions

ZS, HD, YZ, and YX designed this research. ZS, HD, YZ, and YX executed the experiments and analyzed the data. ZS wrote the paper.

Conflict of Interest Statement

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.

Table 1

Correlation between microbiota and endogenous factors.

Sample		Temperature (°C)		pH		Moisture (%)		Ethanol (g/kg fermented grain)		Lactic acid (g/kg fermented grain)		Acetic acid (g/kg fermented grain)
		Range	Pa	Range	Pa	Range	Pa	Range	Pa	Range	Pa	Range	Pa
Time	Day 5 (n = 3)	36.3–42.8	0.800	3.65–3.70	0.008	55.65–55.69	0.001	10.30–34.49	0.381	27.07–28.80	0.046	15.75–23.92	0.235
	Day 10 (n = 3)	38.2–41.7		3.64–3.69		55.62–55.65		17.60–49.99		26.94–29.23		15.74–18.24
	Day 15 (n = 3)	32.7–41.2		3.64–3.70		55.44–55.51		26.87–55.91		24.28–30.47		19.33–24.13
	Day 20 (n = 3)	30.5–40.7		3.61–3.66		54.71–57.73		28.02–49.25		26.71–32.72		16.85–21.08
	Day 25 (n = 3)	29.8–40.5		3.61–3.64		55.69–55.64		34.55–50.74		27.74–30.30		16.30–18.17
	Day 30 (n = 3)	28.7–39.8		3.55–3.59		54.68–54.72		24.16–50.24		30.85–43.00		19.39–32.25