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

1The microbiota in natural environments is complex and often includes thousands of genera from a diverse range of species. In industrial fermentation, the microbiota depends on the conditions and nutrient availability. In natural environments or under artificially controlled conditions, fermentation microbiota have the potential to use different raw materials to produce a variety of metabolites. They are commonly characterized by the presence of specific microorganisms, orderly microbial succession and an unusual functional shift. 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.

2In recent years, the concept of core microbiota has emerged. The core microbiota play important roles in natural ecological systems, human health and industrial processes. However, the core microbiota that influence the quantity and quality of fermented foods in traditional SSF are poorly understood. Traditionally, species were isolated and identified from food fermentation processes by culture-dependent and -independent techniques. Although these studies gained some insight into the microbial populations, our knowledge about the core structural microbiota in fermented foods remain limited. With the development of high-throughput sequencing techniques, researches have focused on the structure of the core microbiota in traditional SSF processes. 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. Similarly, in soy sauce fermentation, 22 frequent genera have been reported to be candidates of the core microbiota.

3Although 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. To address this, research has focused on the function of the core microbiota in traditional food fermentation processes. 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. 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. 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.

4Chinese 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. In addition to the common characteristics of traditional SSF, 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. 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. 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. 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. Based on this information, we explored the correlations between the core microbiota and important metabolites.

5All 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.

6All 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.

7Raw sequences were processed using ‘quantitative insights into microbial ecology’ software (QIIME, version 1.8). 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). 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). 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). Fungal OTUs were manually corrected using the CBS database (Centraalbureau voor Schimmelcultures, Holland) and bacterial OTUs were manually corrected by Ezbiocloud (The Jongsik Chun Lab. and Chunlab, Inc., South Korea). The amplicon databases have been submitted to NCBI SRA and are available under accession numbers SRR5298103 and SRR5298193.

8The 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.

9Raw metatranscriptomic data were processed by removing the rRNA sequences and low quality reads (Q < 0.02). Trinity (version 2.1.1) was employed for de novo transcriptomic assemblies. Using DIAMOND (version 0.8.34) and inputting results into MEGAN (version 5.11.3), we obtained microbial taxonomy information. Then, we compared high quality reads to the non-redundant protein database and the KEGG database to obtain functional annotation information. Heatmap Illustrator (HemI, v.1.0) 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.

10The 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.

11To 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).

12To understand the composition of the microbiota, OTU reads were selected to depict the bacterial and fungal microbiota in liquor production. 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 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 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 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).

13Sequence 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. 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.

14Based 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 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 (R2 = 0.45, P < 0.05) and the moisture content of samples from days 25 to 30 (R2 = 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 (R2 = 0.61, P < 0.05) and the pH of samples from days 25 to 30 (R2 = 0.52, P < 0.01) (Figure 2E).

15To 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 detected 7 positive and 26 negative correlations using Spearman’s correlation coefficient, 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.

16To 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%).

17To 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.

18Furthermore, 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 (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 (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 5G). Conversely, the acetic acid content fluctuated over time (Figure 5H).

19In Chinese soy sauce aroma liquor production, the core microorganisms play critical roles in the characteristics of the flavor metabolites. Compared with previous studies, 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. 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).

20As 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. 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).

21Accordingly, 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. They were able to grow aerobically or anaerobically at low pH and had the ability to limit ethanol and acetic acid production. 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. 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. This genus was able to grow anaerobically at low pH and had the ability to adapt to a large temperature variation (10–45°C). 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).

22It 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. 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. 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. 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. In addition, a variety of heterofermentative Lactobacillus genera had characteristics about rapid growth and effective carbohydrate metabolism. 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. Furthermore, the genus Lactobacillus exerts antifungal effects on food or feed-borne filamentous fungi and yeasts in traditional SSF, such as kimchi fermentation. 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. 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.

23In 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.

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

25The 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.