Spontaneous SSF of Daqu starters at an industrial scale was performed in the fermentation room of a traditional vinegar production factory in Shanxi province, China. Initially, cereal materials of approximately 4200 kg barley and 1800 kg wheat were ground and mixed. Then the mixtures were stirred with the addition of 36–37% water, and were shaped into bricks (28 cm × 18 cm × 5 cm) and layer-by-layer piled in the fermentation rooms. The stacked layers of Daqu blocks were incubated for 24 days with strict temperature control. According to traditional SSF techniques, the variation of room temperature and core temperature of Daqu SSF process was controlled by forced ventilation, and the piles of Daqu blocks were manual turned every 2 days at the thermophilic and cooling stages of SSF process to allow adequate aeration and to decrease the inoculation temperature. To investigate the influence of environmental factors on microbiota dynamics during Daqu SSF process, a MTSSF process and a LTSSF process of Daqu starters were conducted with the highest inoculation temperature reached 48–53°C and 53–55°C (Figure 1, Stage B), respectively. The MTSSF and LTSSF processes of Daqu starters were both conducted through two independent experiments. Daqu samples were separately collected at days 1, 2, 5, 10, 14, and 24 according to temperature evolution during the SSF process (Figure 1). To obtain adequate information and representation, Daqu blocks from each stage were randomly selected from the upper, middle, and lower locations of two individual processes in triplicates, which were then ground, mixed, and pooled into sterile Stomacher bags (Stomacher Lab System, London, United Kingdom) to provide an experimental Daqu powder sample (approximately 500 g). All of the samples were stored at -20°C for further analysis.

The room temperature was monitored every day by a thermo-hygrometer (GEMlead TH339, Fuzhou, China) stuck on the wall of fermentation room, and the core temperature of Daqu (pile temperature) was determined every day using a thermograph (Neo Thermo TVS-700, Nippon Avionics, Tokyo, Japan) inserted into Daqu blocks with a depth of 1 cm. The moisture of Daqu was determined by dry/wet weight measurement method at 105°C. The pH was measured with a pH meter (Sartorius PB-10, Germany). The total titratable acidity was determined by titration with 0.02 M NaOH exhibiting a titration endpoint of pH 8.2. The reducing sugar content was determined by DNS method (Miller, 1959). The amino acid nitrogen content and protease activity were determined according to the national professional standard methods (SB/T 10317-1999, 1999; QB/T 4257-2011, 2011). One unit of protease activity was defined as the amount of amino acid nitrogen liberated per hour by 1 g Daqu under the assay conditions. Amylase and glucoamylase activities were determined as our previously described (Li et al., 2015). One unit of amylase activity was defined as the amount of starch liquefied per hour by 1 g Daqu in sodium acetate buffer (50 mM, pH 4.6) at 35°C. One unit of glucoamylase activity was defined as the amount of glucose liberated per hour by 1 g Daqu in sodium acetate buffer (50 mM, pH 4.6) at 40°C.

DNA extraction from Daqu samples was performed using the Soil DNA Kit (Omega Bio-Tek, Norcross, GA, USA) according to the manufacturer’s instructions. qPCR analysis was performed in quadruplicate using the commercial kit (SYBR^® Premix Ex TaqTM II, Takara, Dalian, China) with an ABI 7500 Real Time PCR System (Applied Biosystems). Primers pairs P1/P2, Lac1/Lac2, B1/B2, and Y1/Y2 were used to quantify the specific gene abundance of total bacteria, LAB, Bacillus, and fungi (Muyzer et al., 1993; Xu et al., 2011), respectively. Each reaction was performed in a 25 μL volume containing 12.5 μL SYBR Premix Ex Taq (Takara, Dalian, China), 0.5 μL of each primer (10 mM) and 2 μL of 10-fold (Bacillus) or 100-fold (bacteria, LAB and fungi) dilution DNA template. The qPCR thermocycling steps were as follows: 95°C for 30 s, 40 cycles of 95°C for 5 s, 55°C for 34 s, 72°C for 30 s. The fluorescent products were detected at the annealing step of each cycle. Melting curve analysis and agarose gel electrophoresis were performed to confirm the specificity of the amplification. The amplification efficiency and correlation coefficient (R²) for the amplification of specific gene of bacteria, Bacillus, LAB and fungi were 99.3% and 0.992, 97.1% and 0.990, 94.4% and 0.991, 98.1% and 0.991, respectively.

The V4 regions of bacterial 16S rRNA gene and ITS1 regions of fungal rRNA genes were amplified used the specific primers 515f/806r (Peiffer et al., 2013) and ITS5-1737F/ITS2-2043R (Huang et al., 2016) with the barcodes, respectively. All PCR reactions were carried out in triplicate 30 μL reactions with 15 μL of Phusion^® High-Fidelity PCR Master Mix (New England BioLabs), 0.2 μM of each primer and 10 ng DNA templates. Thermal cycling consisted of initial activation at 98°C for 1 min, followed by 30 cycles of denaturation at 98°C for 10 s, annealing at 50°C for 30 s, and elongation at 72°C for 60 s and finally elongation at 72°C for 5 min. Negative control were treated similarly with the exclusion of template DNA and failed to produce visible PCR products. PCR products were mixed in equimolar ratios and mixture PCR products were purified with QIAquick Gel Extraction Kit (QIAGEN, Dusseldorf, Germany). Sequencing libraries were generated using TruSeq^® DNA PCR-Free Sample Preparation Kit (Illumina, USA) following manufacturer’s recommendations and index adaptors were added. The library quality was assessed on the Qubit@ 2.0 Fluorometer (Thermo Scientific) and Agilent Bioanalyzer 2100 system. Finally, the library was sequenced on an Illumina HiSeq2500 platform, at Novogene, Beijing, China.

Paired-end reads from the original DNA fragments were merged using FLASH (Zhang et al., 2014) and assigned to each sample with the unique barcodes. UPARSE software package (Uparse v7.0.1001) with the UPARSE-OTU and UPARSE-OTUref algorithms was used to pick operational taxonomic units (OTUs) at the 97% similarity (Edgar, 2013). Representative sequences were picked for each OTU, and RDP classifier (Version 2.2) was used to annotate the taxonomic information for each representative sequence. Alpha diversity indices Chao1, Shannon, Simpson, ACE and Goods coverage were performed in QIIME (Version 1.7.0) to reflect the diversity and richness of microbial community in different samples (Caporaso et al., 2010). For beta diversity, QIIME calculated both the unweighted and weighted UniFrac distances (Lozupone et al., 2011). Sequencing results are available through the GenBank sequence read archive database under accession number PRJNA316566.

The statistical significance (P ≤ 0.05) of the difference among different batches were identified using a one-way analysis of variance (ANOVA). Pearson’s test was performed to reveal the correlations between environmental variables and abundant classes using SPSS Statistics 19.0. Principal coordinate analysis (PCoA) was performed with the weighted UniFrac distance. Paired t-test and wilcoxon tests within the stats R package (Version 2.15.3) were performed to test whether there was a significant difference in the alpha and beta diversity indices among the two batches. In addition, Analysis of similarities (ANOSIM) (Clarke, 1993) and multi-response permutation procedure (MRPP) (He et al., 2010) analyses were further employed to examine the community difference among the two batches. Canonical correspondence analysis (CCA) between Daqu microbial community and measured variables was performed with Canoco 5.0 software.

Temperature was a universal indicator used to monitor and control the Daqu SSF process. Dynamics of this parameter and other physicochemical characteristics throughout the SSF process are shown in Figure 1. The thermal profile was followed the typical evolution of Daqu SSF process, and allowed to distinguish four stages during the MTSSF and LTSSF processes: the mesophilic stage (days 1 to 2), the thermophilic stage (days 2 to 10), the cooling stage (days 10 to 14), and the maturation stage (days 14 to 24) (Figure 1A). However, the temperature of MTSSF process was generally 3–7°C and 5–12°C higher than that in LTSSF process at thermophilic stage and cooling stage (Figure 1A), respectively. The temperature increased from 48 to 53°C and 44 to 50°C (Phase A of thermophilic stage), and maintained around 55°C and 50°C (Phase B of thermophilic stage) at thermophilic stage of MTSSF and LTSSF processes, respectively.

Generally, the evolution profiles of other physicochemical parameters were similar during the MTSSF and LTSSF processes. The moisture slightly declined, and then markedly decreased after 2 days and finally fell to 12.18 and 11.38% (Figure 1B), respectively. The reducing sugar contents shown a substantially decreased tendency during the MTSSF process with the exception of cooling stage, however, a noteworthy increase of reducing sugar contents (32.94 mg g^-1 dry starter) could be observed at phase A and followed by a quick depletion (dropped to 15.33 mg g^-1 dry starter) at phase B of the thermophilic stage during the LTSSF process (Figure 1B). Despite the variation in initial pH and titratable acidity, the pH firstly decreased to about 4.22 and the titratable acidity increased to approximately 0.20 mmol g^-1 at mesophilic stage, but afterward the pH obviously increased and the titratable acidity strongly decreased at the thermophilic stage, and finally the pH and titratable acidity both maintained near 6.40 and 0.06 mmol g^-1 until the end of MTSSF and LTSSF processes (Figure 1C). Nevertheless, significantly (P < 0.05) higher in pH and lower in titratable acidity were displayed at the mesophilic and thermophilic stages of the MTSSF process compare with the LTSSF process (Figure 1C). The pH and titratable acidity on day 5 were 5.28 and 4.87, 0.11, and 0.17 mmol g^-1 during MTSSF and LTSSF processes, respectively.

Dynamics of enzymatic activities throughout the SSF process are shown in Figure 2. Overall, the evolution profiles of enzymatic activities were similar during the MTSSF and LTSSF processes. The protease activities increased from 0.09 to 0.82 mg amino acid nitrogen g^-1 h^-1 on days 1 to 14 and fell to 0.66 mg amino acid nitrogen g^-1 h^-1 at the end of MTSS process, however, the protease activities shown a gradual increased tendency from initial 0.11 to 0.72 mg amino acid nitrogen g^-1 h^-1 at the end of the LTSS process (Figure 2A). The contents of amino acid nitrogen slightly increased within the first 2 days, and significantly (P < 0.05) increased from days 2 to 5 and fluctuated around a peak on day 5. Moreover, significantly (P < 0.05) higher in protease activity and amino acid nitrogen content were displayed on day 14 during MTSSF process compare with the LTSSF process (Figure 2A). The amylase activities slightly declined within the first 2 days, and then increased rapidly from days 2 to 5, and maintained near 1.37 g liquefied starch g^-1 h^-1 until the end of MTSSF and LTSSF processes (Figure 2B). However, the glucoamylase activity slightly decreased within the first 2 days, and then obviously declined throughout the MTSSF and LTSSF processes (Figure 2B). The glucoamylase activities in MTSSF process dropped steadily from 902.71 to 394.33 mg glucose g^-1 h^-1 on days 5 to 24, in LTSSF process from 512.25 to 268.69 mg glucose g^-1 h^-1. Significantly (P < 0.05) higher in glucoamylase activities were noted after the mesophilic stage of the MTSSF process compare with the LTSSF process (Figure 2B).

Dynamics of biomass of total bacteria, LAB, Bacillus and fungi by qPCR throughout the SSF process are shown in Figure 3. Despite the variation in initial quantities, the levels of total bacteria, LAB, Bacillus and fungi firstly decreased significantly (P < 0.05) at phase A of thermophilic stages during the SSF process (Figure 3). Afterward, their quantities increased significantly (P < 0.05) at phase B of thermophilic stages. Notably, higher levels of total bacteria, LAB, Bacillus and fungi were displayed at the thermophilic stages of the MTSSF process compare with the LTSSF process. Afterward, the microbial quantity displayed decreased and increased trends at cooling stages during MTSSF and LTSSF processes, respectively. The levels of total bacteria, LAB, Bacillus and fungi during the MTSSF and LTSSF processes ranged from 9.40 to 11.24 Log copies g^-1, 8.6 to 10.30 Log copies g^-1, 7.58 to 9.90 Log copies g^-1 and 8.69 to 10.10 Log copies g^-1, respectively. Overall, microbial evolution throughout SSF process is a much more dynamic process, and the quantities of total bacteria, LAB, Bacillus and fungi exhibited similar trends during the MTSSF and LTSSF processes.

Illumina HiSeq sequencing was used to investigate the microbial composition and dynamics during the MTSSF and LTSSF processes. After filtering the low-quality reads and chimeras, for bacterial communities, a total of 633,751 effective tags with an average length of 248 bp from all samples were obtained, and each sample contained 34,386 to 66,238 effective tags with different phylogenetic OTUs ranging from 175 to 485 were generated (Supplementary Table S1). For fungal communities, a total of 714,806 effective tags with an average length of 270 bp from all samples were obtained, and each sample contained 26,404 to 73,805 effective tags with different phylogenetic OTUs ranging from 32 to 134 were generated (Supplementary Table S2), via 97% sequence identity cutoff. Moreover, the rarefaction curves approached the saturation plateau (Supplementary Figure S1), which indicated that almost all bacterial and fungal communities could be well represented.

In general, highly similar and successional dynamics of microbial communities were exhibited during the MTSSF and LTSSF processes at order level (Figures 4A,B) and genus level (Supplementary Figure S2). Groups of Rickettsiales and Streptophyta only dominated the initial 2 days but were retrieved at low frequencies (<1% of total sequences) until the end of SSF process (Figure 4A). Meanwhile, thermophilic Eurotiales first became prominent order on day 10 and predominated until the end of MTSSF and LTSSF processes. However, phylotypes of Enterobacteriales, Lactobacillales, Bacillales, Saccharomycetales, and Mucorales both prevailed the whole MTSSF and LTSSF processes (Figures 4A,B). Among these prevailing orders, despite some variation at some moment (Figures 4A,B), (i) the relative abundance of Lactobacillales and Saccharomycetales first increased from days 1 to 5, then decreased from days 5 to 10, and finally increased from days 10 to 24; (iii) the relative abundance of Enterobacteriales and Bacillales first increased from days 1 to 5, then decreased from days 5 to 14, and finally increased from days 14 to 24; (iii) the relative abundance of Mucorales first increased from days 1 to 10, then decreased from days 10 to 14, and finally increased from days 14 to 24. Nevertheless, several facts should be highlighted (Figures 4A,B): (i) the relative abundance of Enterobacteriales in MTSSF process constituted from 10.30% to 71.73%, while obviously lower in LTSSF process and the percentages ranged from 3.16 to 41.06%; (ii) the relative abundance of Lactobacillales in MTSSF process constituted from 2.34 to 16.68%, while obviously higher in LTSSF process and the percentages ranged from 8.43 to 57.39%; (iii) the relative abundance of Eurotiales in MTSSF process decreased steadily from 36.10 to 28.63% on days 10 to 24, while obviously higher in LTSSF process and the percentages gradual increased from 52.00 to 72.97% on days 10 to 24.

Variations in some small proportions of the top 10 bacterial groups of Xanthomonadales (0.10 to 3.94%), Pseudomonadales (0.17 to 1.85%), Clostridiales (0.02 to 0.98%), Actinomycetales (0.03 to 0.68%), and Burkholderiales (0.03 to 0.57%), and fungal groups of Hypocreales (0.00 to 0.42%), Pleosporales (0.00 to 0.41%), Sordariales (0.00 to 0.25%), and Ascomycota sp. (0.01 to 0.10%) were observed during the MTSSF and LTSSF processes. Other bacterial orders of Rhizobiales (0.01 to 0.14%), Rhodobacterales (0.00 to 0.29%), Bacteroidales (0.03 to 0.22%), Flavobacteriales (0.00 to 0.17%), iii1–15 (0.00 to 0.16%), Desulfobacterales (0.00 to 0.16%), Rhodospirillales (0.00 to 0.14%), and Campylobacterales (0.00 to 0.14%), and other fungal orders of Capnodiales (0.00 to 0.17%), Pezizales (0.00 to 0.07%), and Mortierellales (0.00 to 0.05%) were also observed during the MTSSF and LTSSF processes.

Alpha diversity indexes were conducted to evaluate the microbial richness and diversity varied during the MTSSF and LTSSF processes (Supplementary Tables S1 and S2). For bacterial richness and diversity (Supplementary Table S1), the observed OTUs indexes during the MTSSF process generally declined from 261 to 148, in contrast, the observed OTUs indexes during the LTSSF process markedly increased from 147 to 407 during the LTSSF process. Opposite evolution patterns of Chao 1 and ACE values were also observed. However, for fungal richness and diversity (Supplementary Table S2), similar evolution patterns of alpha diversity indexes were observed during the MTSSF and LTSSF processes, for example, the observed OTUs indexes decreased during the MTSSF and LTSSF processes from days 1 to 10 and days 1 to 5, and afterward increased. Notably, obviously lower in bacterial alpha diversity indexes of Shannon, observed OTUs, Chao 1, and ACE indexes during the MTSSF process compare with LTSSF process were observed after days 2, 5, and 10, respectively.

In addition, the box and whisker plots showed that obviously lower (P > 0.05) and significantly lower (P < 0.05) in bacterial alpha diversity indexes of observed OTUs and Shannon indexes were observed during the MTSSF process compare with the LTSSF process, respectively (Figures 5A,B). However, commonly higher in fungal alpha diversity indexes were shown during the MTSSF process compare with LTSSF process. Furthermore, the box and whisker plots showed that dramatically higher (P > 0.05) in fungal alpha diversity indexes of observed OTUs and Shannon indexes were observed during the MTSSF process compare with the LTSSF process (Figures 5A,B).

Principal coordinate analysis was conducted to evaluate the similarities in microbial communities using the weighted UniFrac distance during the MTSSF and LTSSF processes. For bacterial communities, samples from MTSSF process and LTSSF process tended to form clusters (Figure 4C), the weighted UniFrac distance in MTSSF process were significantly lower (P < 0.01) than that in LTSSF process (Figure 5C). However, for fungal communities, clusters separately tended to form days 1 and 2, days 5, and days 10, 14, 24 both during the MTSSF and LTSSF processes (Figure 4D), no significant difference in the weighted UniFrac distance in the MTSSF and LTSSF processes (Figure 5F). Moreover, both ANOSIM and MRPP tests indicated that obviously differences (R > 0, P > 0.05) in bacterial communities but highly similarities (R < 0, P > 0.05) in fungal communities were exhibited in the MTSSF and LTSSF processes (Table 1).

The effect of environmental factors on the bacterial (Figure 6A) and fungal (Figure 6B) communities distribution was evaluated by CCA. Overall, the two axes explained 90.24 and 77.43% of the variation in bacterial and fungal community differentiation, respectively, suggesting the remarkable correlation between microbial structure and environmental factors. Moisture and acidity showed strongly positively correlation with bacterial and fungal composition at the mesophilic stages, but negative correlation with bacterial and fungal composition at cooling and maturation stages (Figure 6). However, core temperature (pile temperature) and pH showed negative correlation with bacterial and fungal composition at the mesophilic stages, but positively correlation with bacterial and fungal composition at cooling and maturation stages (Figure 6). Furthermore, Daqu samples from days 5 were distributed between days 1 to 2 and days 10 to 24 (Figure 6). Pearson’s correlation analysis revealed that moisture, acidity and glucose were mainly positively correlated with the relative abundances of Enterobacteriales, Lactobacillales, Streptophyta, Rickettsiales, Bacillales, and Saccharomycetales but negatively correlated with Eurotiales and Mucorales (Supplementary Table S3). However, pH showed contrasting correlations (Supplementary Table S3). Temperature were mainly positively correlated with the relative abundances of Enterobacteriales, Lactobacillales, Bacillales, Saccharomycetales, and Mucorales but negatively correlated with Streptophyta, Rickettsiales, and Eurotiales (Supplementary Table S3).