Moutai liquor production requires raw materials to undergo multi-round of fermentation, and in each round newly prepared Daqu is added and mixed with sorghum (27). In this study, the prepared Daqu samples were collected from Kweichow Moutai Liquor Co., Ltd. (106°22′E, 27°51′N, Guizhou, China). From December 2020 to July 2021, the prepared Daqu in six production cycles were sampled, and 7 biological replicated samples were randomly collected in each production cycle, resulting in 42 samples collected in total (Supplementary Table 1). The collected samples were stored in dry ice and transported back to the laboratory within 24 h for storage at − 80°C.

The 42 prepared Daqu samples were mixed to form a pooled sample. DNA was extracted from the pooled sample using HiPure Bacterial DNA Kits (Magen, Guangzhou, China) following the manufacturer’s instructions. After extraction, the quality of the DNA was analyzed by Qubit and Nanodrop (Thermo Fisher Scientific, MA, United States). Then the extracted DNA was subjected to sequencing using a NovaSeq 6,000 sequencer (Illumina Inc., CA, United States) with pair-end technology (PE 150). The DNA was firstly fragmented to a size of around 350 bp by sonication. Then, the fragments were end-repaired, A-tailed, and adaptor ligated by using the NEBNext ΜLtra DNA Library Prep Kit for Illumina (NEB, MA, United States) following the manufacturer’s instructions. PCR was applied on the DNA fragments with the length of 300–400 bp, and the amplified products were purified using an AMPure XP system (Beckman Coulter, CA, United States). Size distribution of the libraries were analyzed by 2,100 Bioanalyzer (Agilent, CA, United States). The libraries were also quantified using real-time PCR.

Raw data were filtered using FASTP (version 0.18.0) with the following standards: (1) removing reads with ≥ 50% bases having Phred quality scores ≤  20; (2) removing reads with ≥ 10% unidentified nucleotides (N); (3) removing reads aligned to the barcode adapter. The reads were assembled by MEGAHIT (version 1.2.9). MetaGeneMark (version 3.38) was used for gene identification, which output the amino acid sequences of the genes, and the amino acid sequences were used as the protein database for metaproteome analysis. The genes obtained by metagenomic sequencing were aligned by Minimap2 (version 2.17-r941) against the NCBI bacterial reference database and thereby the corresponding proteins were annotated with taxonomic information.

Each sample was weighed at 0.5 g and grinded twice with a mechanical grinder (JXFSTPRP-CL, Shanghai Jingxin Industrial Development Co., Ltd., Shanghai, China). After grinding twice with 5 steel balls (−50°C, 70 Hz, on 120 s, off 120 s), the powder was suspended in 1 mL borax/polyvinylpolypyrrolidone/phenol (BPP) solution (100 mM EDTA, 50 mM borax, 50 mM vitamin C, 30% sucrose, 10 mM Tris-base, 1% Triton-100, 5 mM dithiothreitol, 1% polyvinylpolypyrrolidone) and grinded twice with the mechanical grinder under the same conditions as above described. The samples were then suspended in 1 mL DNA extraction phenol reagent (Beijing Solarbio Science and Technology Co., Ltd., Beijing, China) followed by vortex for 2 min and centrifugation at 12000 g for 20 min at 4°C. The upper phenol phase was taken and mixed with 1 mL BPP solution under the same conditions as above to repeat the lysis step. Then, the upper phenol phase was taken and mixed with 5 mL 0.1 M ammonium acetate in methanol pre-cooled to −20°C, followed by precipitation for 4 h at −20°C. The purified proteins were centrifuged at 12000 g for 20 min at 4°C and washed twice with 1 mL 0.1 M ammonium acetate in methanol pre-cooled to −20°C. After drying at room temperature, the protein precipitates were dissolved in 300 μL of lysate buffer (1% sodium dodecyl sulfate and 8 M urea aqueous solution).

The protein was quantified using Pierce BCA protein assay kit (Thermo Fisher Scientific, MA, United States). For each sample, 200 μg proteins were taken, and 8 M urea solution and 20 μL 1 M triethylammonium bicarbonate buffer (TEAB) solution was added to reach a constant volume of 200 μL. Then, 4 μL 0.5 M tris-(2-carboxyethyl) phosphine hydrochloride (TCEP) solution was added for reduction. The mixture was vortexed at 600 rpm for 1 h at 37°C. After that, 18 μL 0.5 M iodoacetamide solution was added for alkylation for 45 min at 25°C in the dark. Then, 1.2 mL acetone pre-cooled to −20°C was added, followed by precipitation for 4 h at −20°C. The purified proteins were centrifuged at 15000 g for 15 min at 4°C and washed twice with 90% acetone pre-cooled to −20°C. After drying at room temperature, the protein precipitate was dissolved in 200 μL 0.1 M TEAB. Then, 20 μg trypsin (Beijing Wallis Technology Co., Ltd., Beijing, China) per 1 mg protein was added for digestion at 600 rpm for 16 h at 37°C. Peptides from each sample were desalted with MonoSpin C18 column (Shimadzu, Tokyo, Japan) and quantified using the Pierce quantitative colorimetric peptide assay kit (Thermo Fisher Scientific, MA, United States).

For metaproteomic analysis, all samples were analyzed by a nanospray Orbitrap Fusion Lumos Tribrid MS (Thermo Fisher Scientific, MA, United States) with a Nano ACQUITY UPLC system (Waters Corporation, MA, United States). For each sample, 10 μg peptides were redissolved in 30 μL solvent A (0.1% formic acid in water) spiked with 1 × iRT standard peptides (Biognosys AG, Schlieren, Switzerland). Then 1 μg of peptide sample was loaded to an C18 column (Acclaim PepMap, 75 μm × 25 cm, Thermo Fisher Scientific, MA, United States) and separated with a 60-min gradient (Supplementary Table 2), from 2 to 95% solvent B (0.1% formic acid, 20% water and 80% acetonitrile). The flow rate was maintained at 250 nL/min and the column temperature was maintained at 40°C. DIA mode was used to analyze the samples. The parameters were: method duration = 60 min, ion source type = NSI, spraying voltage positive ion = 2,200 V, spraying voltage negative ion = 2,100 V, ion transfer tube temperature = 320°C, cycle time = 3 s, MS desired minimum points across the peak = 9, using wide quad isolation = True, MS orbitrap resolution = 120,000, MS scan range (m/z) = 349.5–1500.5, MS maximum injection time = 50 ms, MS AGC target = 400,000, MS/MS isolation window = 1.6, maximum number of multiplexed ions = 0, CID activation time = 10 ms, activation type = HCD, collision energy = 32%, stepped collision energy = 5%, MS2 PTR reaction time = 50 ms, MS/MS orbitrap resolution = 30,000, MS/MS scan range (m/z) = 200–2000, MS/MS maximum injection time = 72 ms, MS/MS AGC target = 50,000. Sixty variable windows were set for MS/MS acquisition (Supplementary Table 3).

For parallel reaction monitoring (PRM) analysis, results of DIA experiment were used to select proteotypic or protein group specific peptides and to develop PRM assays by SpectroDive 11.6 (Biognosys AG, Switzerland). The inclusion list of the final optimized PRM method was present in Supplementary Data 1. Peptides of every 7 samples from the same production cycle were mixed by equal peptide quantity, resulting in six samples for PRM analysis. The peptides were re-dissolved in solvent A and analyzed by on-line nanospray LC–MS/MS on an Orbitrap Fusion lumos Tribrid MS (Thermo Fisher Scientific, MA, United States) coupled to an EASY-nanoLC 1,000 system (Thermo Fisher Scientific, MA, United States). 2 μL peptide sample was loaded onto a 25 cm analytical column (Acclaim PepMap C18, 75 μm × 25 cm) and separated with a 120 min-gradient starting at 2% solvent B followed by a stepwise increase to 35% B in 95 min, 60% B in 17 min, 100% B in 3 min and stayed there for 5 min. The column flow rate was maintained at 400 nL/min with the column temperature of 55°C. The electrospray voltage was set to 2 kV. PRM settings were as follow: Full MS scans with the mass range from m/z 450 to 1,350 were acquired with a resolution of 60,000, AGC target of 1 × 10⁶ and a maximum injection time of 50 ms. MS2 spectra were acquired with a resolution of 15,000, AGC target of 5 × 10⁴ and a maximum injection time of 80 ms.

The DIA raw data were analyzed by Spectronaut (28) (version 15.4.210913, Biognosys AG, Schlieren, Switzerland). The data were searched by directDIA against a protein sequence database without the need of spectral libraries. The protein sequence database (1,648,851 entries) was built from metagenomic sequencing. Trypsin was used for proteolysis and the maximum number of missed cleavages was 2. Retention time prediction type was set to dynamic iRT. Data extraction was determined by Spectronaut based on the extensive mass calibration. Spectronaut will determine the ideal extraction window dynamically depending on iRT calibration and gradient stability. Q-value (FDR) cutoff on precursor and protein level was applied as 1%. Decoy generation was set to mutated. All selected precursors passing the filters were used for quantification. MS2 interference removed all interfering fragment ions except for the three least interfering ones. The average top 3 filtered peptides which passed the 1% Q-value cutoff were used to calculate the major group quantities. Only the leading protein (with the strongest evidence and ranked first in the result) in each protein group was taken into consideration in all the subsequent analysis.

The PRM raw data were analyzed by SpectroDive 11.6 with the default settings. SpectroDive calculated the ideal mass tolerances for data extraction and scoring based on its extensive mass calibration. Q-value (FDR) cutoff on precursor was applied as 1%. Peptides were manually inspected to verify the matched mass spectra and the peak integration.

Principal component analysis (PCA) and PLS-DA were performed using MetaboAnalyst (16, 29) (version 5.0, https://www.metaboanalyst.ca/) and MetaboAnalystR (version 3.0, https://github.com/xia-lab/MetaboAnalystR). The quantified proteins were annotated with eggNOG (30) (version 4.5.1, http://eggnogdb.embl.de/). COG and EC number annotations were extracted from the eggNOG results. Statistical significance of the differences among the groups of samples was investigated by Kruskal–Wallis test (p-value < 0.05). To explore the patterns of microbial interactions, network analysis was carried out based on Spearman rank correlations. Genera with the relative abundance > 0.1% and appeared in at least 80% of the samples were shown as nodes. A connection stood for significant correlation (p-value < 0.05). Size of each node was proportional to the number of connections. The nodes were colored by genera occupancy. The thickness of edge was proportional to the value of Spearman correlation coefficients.

Data visualization was performed using R (version 4.0.4, https://www.r-project.org/), with the packages ggplot2 (version 3.3.5, https://github.com/tidyverse/ggplot2) and Venn Diagram (version 1.6.20). The interactive platform Gephi (31) (version 0.9.2, https://gephi.org/) was used to generate network diagrams.

The Moutai liquor production process is shown in Figure 1A. In this study, we collected a total of 42 prepared Daqu samples from six production cycles with seven biological replicated samples in each cycle. The prepared Daqu were collected from December 2020 to July 2021 (Supplementary Table 1). As shown in Figure 1B, microbial proteins were extracted from the prepared Daqu using mechanical grinding, and digested into peptides with trypsin. The digested samples were analyzed by label-free DIA metaproteomics. The DIA data were searched by directDIA against a sample-specific protein sequence database (1,648,851 protein entries) constructed by metagenomic sequencing of the same prepared Daqu samples. From the 42 samples, 3,009 proteins and 7,872 peptides were identified and quantified totally, with 2,672 ± 104 proteins and 6,720 ± 469 peptides (mean ± standard deviation) per sample (Supplementary Figure 1; Supplementary Data 2). The numbers of quantified proteins and peptides from the prepared Daqu samples of each production cycle were very close (Figures 2A,B). The consistent protein and peptide numbers indicated the reproducible measurements by the DIA-based metaproteomics in this study.

Based on the quantified proteins, we performed principal component analysis (PCA) to observe how the prepared Daqu samples in the 6 production cycles were clustered by themselves (Figure 2C). In the PCA score plot, the samples in cycles 2, 3 and 4, as well as cycles 5 and 6 were highly overlapped. We also explored the species composition of Daqu microbiota across the 6 production cycles (Figure 2D; Supplementary Data 3). There were in total 82 microbial species identified. Differences of microbial abundances were significant among cycle 1, cycles 2–4, and cycles 5–6. The species in Bacillus were less abundant in cycle 1 compared to cycles 2–4, while those in Kroppenstedtia and Oceanobacillus were more abundant in cycles 5–6 compared to cycle 1 and cycles 2–4.

According to the PCA result and taxonomic compositions, the prepared Daqu samples from 6 production cycles were then divided into three groups: (1) group 1: cycle 1; (2) group 2: cycles 2, 3 and 4; (3) group 3: cycles 5 and 6. Besides PCA, we performed partial least squares–discriminant analysis (PLS-DA) for the classification among the 3 groups (Figures 2E,F). The PLS-DA score plots revealed a good separation of the different groups. The species composition of Daqu microbiota was compared among the 3 groups as well, revealing more pronounced abundance differences than the previous 6-cycle-based comparison (Figure 2G). The quantified proteins and peptides were then counted within each group (Supplementary Figure 2). The 3 groups shared 2,930 proteins and 7,668 peptides, accounting for 97.4 and 97.4%, respectively, of all the quantified proteins and peptides from the 3 groups. All the above results indicated that the metaproteome profiles and taxonomic compositions of the prepared Daqu samples were relatively stable across different production cycles, while can be clustered into different groups revealing the compositional and functional changes of prepared Daqu during the production cycles.

From the 42 prepared Daqu samples, we identified a total of 18 dominant species (belonging to 10 genera) with an average relative abundance greater than 0.5% (Supplementary Data 3). The top 8 dominant species were Kroppenstedtia eburnea (17.83%), Bacillus sonorensis (3.35%), Virgibacillus halodenitrificans (2.81%), Bacillus amyloliquefaciens (2.26%), Bacillus velezensis (1.98%), Oceanobacillus kimchii (1.75%), Bacillus thermoamylovorans (1.69%), and Enterococcus faecium (1.05%). Kruskal–Wallis test was carried out to find the significantly differential species (with relative abundance > 0.15% in at least one group and FDR-corrected p-value < 0.05) among different groups (Figures 3A,B, as well as Supplementary Figure 3; Supplementary Data 4). The result indicates that 9 species in the Bacillus genus and 8 species in the Staphylococcus genus were more abundant in the group 3 compared to the other groups.

We next constructed a microbial community co-occurrence network for the prepared Daqu samples from each production cycle (Figures 3C–E; Supplementary Figure 4; Supplementary Data 5). The co-occurrence networks of microbial communities of the prepared Daqu samples from different production cycles showed different connectivity patterns. In the microbial community of production cycle 1, Bacillus was the top genera with the largest degree of connection (27.78% of the total degree of connection), followed by Staphylococcus (19.44%), Brachybacterium (5.56%), Corynebacterium (5.56%), and Oceanobacillus (2.78%). In production cycle 3, the top 4 genera with the largest degree of connection were Bacillus (26.83%), Staphylococcus (24.39%), Brachybacterium (2.44%), and Corynebacterium (2.44%). In production cycle 6, connectivity among the microbial community was primarily driven by Staphylococcus (31.91%), followed by Bacillus (23.4%) and Brachybacterium (4.26%). We observed higher connectivity of Staphylococcus in group 3, while higher connectivity of Bacillus in group 1 and group 2. In addition, we investigated the alpha diversity of the Daqu microbial communities from different production cycles by computing the Shannon diversity indexes (Figure 3F). The microbiota of group 3 exhibited significantly higher microbial diversity than the other groups.

A total of 306 differential microbial proteins among the three groups of prepared Daqu were screened by PLS-DA (VIP score > 1.5) and submitted for Kyoto Encyclopedia of Genes and Genomes (KEGG) annotation (Figure 4A; Supplementary Data 6). Among the annotated pathways, most changes occurred in the metabolism of carbohydrate, amino acid, energy, cofactors and vitamins, nucleotide, and lipid. The differential proteins also involved in the genetic information processing, including translation, folding, sorting and degradation, as well as transcription.

The abundances of the KEGG pathways were calculated based on the quantification results of all the proteins assigned to each pathway. Then, the differential pathways were screened by Kruskal-Wallis test. There were 93 differential pathways with Kruskal-Wallis p-value < 0.05 among the three groups of samples (Supplementary Figure 5; Supplementary Data 7), mainly focusing on amino acid metabolism (16), carbohydrate metabolism (14) and signal transduction (8). Subsequently, the above differential pathways among the three groups were pairwisely compared by Wilcoxon rank sum test, and the differential pathways with Wilcoxon p-value < 0.05 were extracted (Figure 4B; Supplementary Data 7). The differential pathways among all the samples of the 6 production cycles are also shown in Supplementary Figure 6. The pathways of galactose metabolism, starch and sucrose metabolism, and other glycan degradation in carbohydrate metabolism were significantly decreased in group 3 (cycles 5 and 6). The galactose metabolism pathway was more abundant in group 1 (cycle 1), while the starch and sucrose metabolism pathways were more abundant in group 2 (cycles 2–4). In contrast, amino acid metabolism, organic acid metabolism in carbohydrate metabolism, biosynthesis of secondary metabolites, and metabolism of cofactors and vitamins were significantly elevated in group 3.

Saccharification is the most basic process in the production of Chinese liquor, which involves a variety of enzymes, including amylase, glucoamylase and cellulosase. These enzymes hydrolyze starch and cellulose in the raw materials into low-molecular carbohydrates, which serve as substrates for subsequent fermentations. The quantified proteins in the 42 prepared Daqu samples were annotated with enzymatic function using the eggNOG database (30), and a total of 18 proteins were matched to the 4 key enzymes for starch and cellulose hydrolysis, including α-amylase, 1,4-β-fibrodiglycosidase, glucoamylase, and β-D-glucosidase.

The abundances of the 18 enzymatic proteins in the three groups of prepared Daqu samples are shown in Figure 5 and Supplementary Data 8. The differences of abundances of the 18 enzymatic proteins among the 6 cycles are also shown in Supplementary Figure 7. Most of the proteins were less abundant in group 3 (cycles 5 and 6) compared to group 1 and group 2. PRM analysis was used to verify the enzymes. Fifteen of the 18 proteins were detected by PRM, and their relative quantities among the different cycles and groups by PRM maintained good agreement with the label-free DIA results (Supplementary Figure 8). In addition, we also analyzed the taxonomic source of the proteins. Most proteins (5) of α-amylase were from eukaryote, except 1 from Actinobacteria. Most proteins (7) of 1,4-β-fibrodiglycosidase and all the glucoamylase proteins (2) were from fungi, while 1 protein of 1,4-β-fibrodiglycosidase was from Bacilli. The proteins of β-D-glucosidase were from Actinobacteria (1) and fungi (1). Although the abundance of bacteria was higher than fungi in the prepared Daqu samples, the enzymes related to starch and cellulose hydrolysis were mainly from fungi.