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.

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 chromatographymass 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 andindependent 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 highthroughput 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 1). 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.

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.

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)  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)  (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 R 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 CaCl 2 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 CaCl 2 . 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 NaH 2 PO 4 . The mass spectra of unknown flavor components were compared with those in the Wiley 275.L (Agilent) database.

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 1A 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 1B). By contrast, the relative abundance of acids increased from days 15 to 30, ranging from 62.70 to 70.34% ( Figure 1B). 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 1B). 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 1C), and more specifically, from −0.36 to 1.77 for ethanol ( Figure 1D). 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 1C) and −0.01 to 2.17 for lactic acid and acetic acid specifically ( Figure 1D).

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 Figure S3 and Tables S1, S2). In addition, the bacterial Shannon diversity index  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 2A). 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 2B). 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.
Based on 16S rRNA and ITS gene amplicons, Unifrac distancebased weighted principal coordinates analysis was conducted to evaluate similarities and differences in the microbiota (prokaryotic and eukaryotic microbiota) among the different samples ( Figure 2C). 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 1). 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 2D) 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 2D). 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 (R 2 = 0.45, P < 0.05) and the moisture content of samples from days 25 to 30 (R 2 = 0.56, P < 0.05) ( Figure 2D). 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 (R 2 = 0.61, P < 0.05) and the pH of samples from days 25 to 30 (R 2 = 0.52, P < 0.01) ( Figure 2E).

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 3 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.

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) 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 4). 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.
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 5). 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 5B,C). 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 5D). 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 5E). The lactic acid content ( Figure 5F) correlated with the functional shift in lactic acid production. PYC and malate dehydrogenase (oxaloacetatedecarboxylating) (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 5G). Conversely, the acetic acid content fluctuated over time ( Figure 5H).

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 2A,B). Compared with previous studies , 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 1). 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 1C,D). There was a significant difference (P < 0.05) between the ethanol and lactic acid contents at different stages of liquor production ( Figure 2D 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 2E). 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 2). 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 4, 5). 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 4, 5D). 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 5C). At this stage, the ACS, PYC, and ALDH genes were not expressed in the genus Lactobacillus (Figure 4). 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 5). 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 4). 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 4, 5E,H). With the genus Lactobacillus being the main contributor under low pH, lactic acid production was the major microbial function in Stage II (Figure 5F).
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 5I), 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.