Response and assembly of abundant and rare taxa in Zaopei under different combination patterns of Daqu and pit mud: from microbial ecology to Baijiu brewing microecosystem

Yu Mu, Jun Hung, Rongqing Zhou,b,, Suyi Zhng, Hui Qin, Hnln Tng, Qinglin Pn, Huifng Tng

a College of Biomass Science and Engineering, Sichuan University, Chengdu 610065, China

b National Engineering Research Center of Solid-State Manufacturing, Luzhou 646000, China

c Luzhou Laojiao Co., Ltd., Luzhou 646699, China

Keywords: Bioaugmented Daqu Artif icial pit mud Abundant and rare taxa Community assembly Metabolic function

ABSTRACT The quality and aroma of strong-flavor Baijiu are mainly dependent on Daqu, pit mud (PM), and the interaction of both. However, little is known about how their combination patterns affect the microbiome and metabolome of Zaopei, especially the metabolic function of rare taxa. Here, an experiment on industrial size was designed to assess the effects of 6 combinations (3 kinds of Daqu × 2 kinds of PM) on the composition and assembly of different taxa, as well as the flavor profile. The results showed that Zaopei’s microbiota was composed of a few abundant taxa and enormous rare taxa, and rare bacterial and abundant fungal subcommunities were signif icantly affected by combination patterns. The assembly processes of abundant/rare taxa and bacterial/fungal communities were distinct, and environmental changes mediated the balance between stochastic and deterministic processes in rare bacteria assembly. Furthermore, specif ic combination patterns improved the f lavor quality of Zaopei by enhancing the interspecies interaction, which was closely related to rare taxa, especially rare bacteria. These f indings highlighted that rare bacteria might be the keystone in involving community interaction and maintaining metabolic function, which provided a scientif ic foundation for better understanding and regulating the brewing microbiota from the viewpoint of microbial ecology.

1. Introduction

Baijiu is one of the most famous traditional alcoholic beverages with a distinctive aroma and bioactive molecules profile in the world, which has a history of more than 2 000 years[1]. Strong-f lavor Baijiu (SFB) has a market share of about 70% in China, its unique skills include using medium-temperatureDaqu, grain steaming and liquor distilling simultaneously, and fermenting in a mud pit, which endows the fresh Baijiu with strong pit f lavor and long aftertaste[2].In the SFB fermentation process,Zaopei(fermentation substrate)is the unique niche in which the succession and assembly of the microbial community and the metabolism of various compounds are completed[3-4].Daqu(starter) provides multiple functional species and enzymes to initiate this process, while pit mud (PM) with various anaerobes is essential for flavor synthesis[3,5-6]. However,Daqufrom different batches often occurs quality f luctuation due to spontaneous fermentation[7], and the evolution of microbial consortia in PM needs more than 20 years[8]. Therefore, the development of bioaugmentedDaquand artificial PM (APM) based on the inoculating of the functional strain or consortia has become one of the effective strategies to regulate the microbiome and metabolome ofZaopei. For example, the inoculation ofBacillus,Saccharomyces,andWickerhamomycesimproved the enzyme activities and flavor characteristics ofDaqu[9-11]. Similarly, the inoculation of hexanoic acid-producing bacteria and/or high-quality PM shortened the maturation period of PM[12-13]. BioaugmentedDaquand APM have been recommended in previous studies due to their beneficial effects such as reducing higher alcohol content, increasing ethanol yield, and improving flavor quality[11-12,14]. In addition, bioaugmentedDaqucan also promote the directional evolution of PM microbiota when used continuously, as the interphase mass transfer affects not only material metabolism but community diversity[8,15]. Therefore, the production of high-quality SFB is highly dependent onDaqu, PM, and their interaction, and the attributes and combinations ofDaquand PM should be considered to achieve better performance.

In recent years, it was comprehensively revealed the microbes involved in SFB fermentation by emerging molecular biology technology.Lactobacilluswas the dominant bacteria inZaopei,whereas the dominant fungi varied withDaququality, pit age,and geographical factors, and their potential sources were also diverse[5,16-18]. Similar to natural ecosystems, the microbiota in brewing micro-ecosystems is skewed, with few abundant taxa coexisting alongside substantial rare taxa[19]. With the development of microbial ecology, there is increasing evidence that rare taxa play a disproportionate role in ecosystems as the “seed bank” of genetic and functional diversity[20]. For example, rare taxa were not only linked to the restoration of soil ecosystem multifunctionality as keystone species but also had a greater impact on crop yield than abundant ones due to their critical role in mediating soil nitrogen cycling[21-22].However, studies on Baijiu fermentation mainly focused on abundant taxa, leaving no information on community interactions and metabolic functions of rare taxa.

Furthermore, the responses of bacteria/fungi communities and abundant/rare taxa to environmental disturbances are often distinct, resulting in their radically different community assembly processes[23-24]. The assembly of species into communities is controlled by deterministic (selection) and stochastic (dispersal and drift)processes[25]. Quantifying the relative importance of deterministic and stochastic processes to community assembly of the specific ecosystem may be beneficial to regulating the functional microbiota.Several studies proved that the community assembly of different microecosystems was governed by distinct ecological processes[23,26],which may be related to the disparate growth habits and dispersal ability of microorganisms, as well as their differential response to environmental changes[27-28]. Additionally, environmental changes mediated the balance between deterministic and stochastic processes in community assembly at the temporal and spatial scales[29]. DifferentDaquand/or PM also altered the physicochemical properties ofZaoepi, such as acidity, reducing sugar (RS) and ethanol content[5,14,30].It is reported that these variations affected the community assembly process ofZaopei, but not considering the unique ecological processes of abundant and rare taxa[31]. This may not be effective in identifying the key driving factors of assembly processes and regulating community and metabolism subsequently.

In this study, the effects of combinations ofDaquand PM with different attributes on the physicochemical properties, volatile metabolites, and microbial communities ofZaopeiwere investigated by polyphase detection approaches. The objectives were to(i) compare the responses of abundant and rare taxa to different combination patterns, (ii) reveal the assembly processes and main driving factors of abundant and rare microbial communities, and(iii) explore the metabolic function of rare taxa.

2. Materials and methods

2.1 Experimental design and sampling

The experiment was conducted in Luzhou Laojiao Co., Ltd.(Luzhou, Sichuan, China), and 6 combination patterns (3 kinds ofDaqu× 2 kinds of PM) were applied to SFB fermentation (n= 48).The 3 kinds ofDaquare spaceDaqu(SD), fortifiedDaqu(FD), and traditionalDaqu(TD). They were produced by the conventional process and stored for 3 months before being used. SD was manufactured via inoculatingMuqupowder at a 1% ratio by weight of wheat, andMuquwas cultured batch by batch via inoculatingDaquthat endured the space flight for a month on the Shenzhou 11 spacecraft. FD was produced by inoculating the mixed suspension ofBacillus velezensiswith a high yield of tetramethylpyrazine andBacillus subtiliswith a high yield of 4-ethylguaiacol at a ratio of 106:106cells/mL. TD was produced from the same batch without inoculation. The differences in their attributes were as follows[32]: SD with the highest liquefaction and esterification powers was dominated byKroppenstedtia,Bacillus, andPichia; FD with the highest saccharification and fermentation powers was characterized byEnterococcus,Bacillus, andThermomyces; while TD with the lowest enzyme activities was mainly composed ofLactobacillus,WeissellaandAspergillus(Table S1).

The two kinds of PM, namely fortified APM (FPM) and conventional APM (CPM), were manufactured in 2019 according to the process described previously[33]. The identical aged PM cultures and differentDaquwere used for cultivation, and FD and TD were inoculated in FPM and CPM, respectively. Their attributes were also disparate significantly,the pH and the abundance of archaea andClostridiumwere higher in the FPM, while the total acid and the abundance of eubacteria andLactobacilluswere higher in the CPM[33](Table S1).

The SFB fermentation process was as follows: (i) theZaopeifrom the last batch of fermentation was mixed with crushed sorghum and rice husk at 20:5:1 for distillation; (ii) the moisture content was adjusted to 53%−58% (m/m) with 80 °C hot water; (iii) 20%Daqupowder was added to the cooled grains and mixed well;(iv) subsequently, the mixture was put into the mud pit for solid-state fermentation. The experiment was performed in December 2021 and samples were collected in March 2022 after 90 days of fermentation.Firstly, 100 g per each point ofZaopeiwas collected from the 4 corners and center points of the upper, middle, and bottom layers, and mixed evenly. Then taken out of 500 g as one sample and divided into three parts for analyzing microbial communities (−80 °C), volatile flavors (4 °C),and physicochemical properties (4 °C), respectively. Each experimental group included eight biological replicates that were collected from different mud pits, and TD + CPM treatment was the control.

2.2 High-throughput sequencing and data processing

Total genomic DNA was extracted using the FastDNA SPIN extraction kit (MP Biomedicals, Santa Ana, CA, USA) according to the manufacturer’s protocols. The quality and quantity of obtained DNA were checked by NanoDrop ND-1000 spectrophotometer(Thermo Scientific, Waltham, MA, USA). The bacterial 16S rRNA V3-V4 region and fungal ITS1 region were respectively amplified by PCR using primers 338F (5’-ACTCCTACGGGAGGCAGCA-3’)/806R (5’-GGACTACHVGGGTWTCTAAT-3’) and ITS5F(5’-GGAAGTAAAAGTCGTAACAAGG-3’)/ITS1R(5’-GCTGCGTTCTTCATCGATGC-3’), and the amplification procedure was described in our previous study[15]. After being separated by 2% agarose gel electrophoresis, the amplicons were purified and quantified, and the final products were sequenced using the Illumina NovaSeq platform (2 × 250 bp).

The raw data were processed by QIIME2 (2019.4). After removing the primers from raw sequences, the remaining sequences were used to generate an amplicon sequence variant (ASV) table by DADA2, which includes quality filter, denoising, merging, and chimera removing[34]. Then, all ASVs were assigned using the naive Bayes taxonomy classifier in the feature-classifier plugin against the Silva (v 132) and UNITE (v 8.0) databases. Raw sequence data for all the samples were uploaded to the NCBI Sequence Read Archive under accession number PRJNA870386.

2.3 Detection of volatile flavors

Trace 1300-TSQ 9000 GC-MS (Thermo Scientific, Waltham,MA, USA) equipped with a VF-WAX-MS capillary column (30.0 m ×0.25 mm, 0.25 μm, Agilent, Santa Clara, USA) was applied to detect the volatile flavors ofZaopei[5]. Briefly, 1.00 g sample and 20 μL methyl octanoate (internal standard, 7.9 mg/100 mL) were accurately added to a 20 mL headspace bottle, and a 50/30 μm divinylbenzene/carboxen/polydimethylsiloxane (DVB/CAR/PDMS)fiber (2 cm, Supelco, Bellefonte, PA, USA) was used to extract volatile flavors at 60 °C for 50 min, and then the extraction head was removed and inserted into the inlet and desorbed for 5 min. GC conditions: inlet temperature, 250 °C; carrier gas, high-purity helium(＞ 99.999%); flow rate, 1 mL/min; spitless. The temperature of the column is as follows: 40 °C for 5 min, 4 °C/min to 100 °C, and 6 °C/min to 230 °C for 10 min. MS conditions: ion source temperature, 250 °C; transmission line temperature, 300 °C; ionization mode, EI (70 eV); scan range, 35–400 amu. After comparing with the standard mass spectrogram from the NIST 2017 library, only compounds with similarity > 80% remained for subsequent analysis.

2.4 Analysis of physicochemical properties

The moisture content was determined using a gravimetric method(105 °C, 4 h). The pH was determined using a pH meter (pHS-3C,INESA, Shanghai, China) after 1:10 (m/V)Zaopei-water mixture extracted for 30 min at room temperature. The acidity was measured by the alkali titration method with a standard NaOH solution(0.1 mol/L). The ethanol content was detected by an alcoholmeter after distillation. The starch and RS content were measured by 3,5-dinitrosalicylic acid (DNS) method.

2.5 Definition of abundant and rare taxa

To analyze the responses of abundant and rare taxa to different combination patterns, all ASVs were classified into 6 categories following recent works[21-22]: (i) always abundant taxa (AAT) with relative abundance ≥ 1% in all samples; (ii) conditionally abundant taxa (CAT) with relative abundance ≥ 1% in some samples and≥ 0.01% in all samples; (iii) always rare taxa (ART) with relative abundance ≤ 0.01% in all samples; (iv) conditionally rare taxa (CRT)with relative abundance ＜ 1% in all samples and ＜ 0.01% in some samples; (v) moderate taxa (MT) with relative abundance between 0.01% and 1% in all samples; (vi) conditionally rare and abundant taxa (CRAT) with relative abundance ranging from rare (＜ 0.01%)to abundant (≥ 1%). Subsequently, AAT, CAT, and CRAT were collectively referred to as abundant taxa, while ART and CAT were as rare taxa for analysis.

2.6 Statistical analysis

All statistical analyses were conducted in R software (v 4.1.3;https://www.r-project.org/) unless otherwise stated. Theα-diversity indices (Shannon and Chao1) and Bray-Curtis distance-based community dissimilarity for abundant and rare taxa were calculated by the “vegan” package. The effects ofDaqu, PM, and their interaction on community diversity and physicochemical properties were determined by two-way ANOVA in the “stats” package. To determine the driving factors of community diversity variation, the Mantel test and the distance-based redundancy analysis (i.e., CAP)were performed by “ggcor” and “vegan” packages, respectively.

The assembly processes of abundant and rare taxa were assessed by beta nearest-taxon index (βNTI) and null-model-based Bray-Curtis-based Raup-Crick (RCbray) indices in the “picante” package[25].A value of |βNTI| > 2 is suggesting that the community assembly is controlled by a deterministic process (βNTI ＜ –2, homogeneous selection; βNTI > 2, heterogeneous selection). When |βNTI| is ＜ 2(stochastic process), RCbray> 0.95 and ＜ –0.95 respectively means the dispersal limitation and homogenizing dispersal, while |RCbray| ＜ 0.95 represents an undominated process mostly comprising weak selection,weak dispersal, diversification, and/or drift[35]. Relationships between βNTI values and the changes in physicochemical properties were evaluated by the linear regression analysis.

The co-occurrence network was constructed based on the molecular ecological network analysis (MENA, http://ieg2.ou.edu/MENA) pipeline[36]. Pearson correlation coefficient among ASVs that were present in at least 3/8 samples was conducted, and the cutoff threshold was determined based on random matrix theory, then the network data was exported to Gephi (v 0.92) for visualization. The flavor profile was compared by principal component analysis (PCA)based on “factoextra” package, and the two-way permutational multivariate analysis of variance (PERMANOVA) test was performed using the function “adonis” in the “vegan” package. The orthogonal partial least squares discrimination analysis (OPLS-DA) model was established by Simca 14.0 (Umetrics AB, Umea, Vasterbotten, Sweden)to identify the differential compounds, and the correlations between differential compounds and the top 20 ASVs of abundant and rare taxa were calculated using “Hmisc” package. The main microbial predictors for the content of butyric acid and hexanoic acid were determined by the random forest model[37]. TheR2andPvalues of models were assessed with 5 000 permutations of the response variable by the “A3” package.The “rfPermute” package was further used to obtain the significance of each predictor on the response variables. The functional compositions of abundant and rare taxa were predicted by PICRUSt2[38].

3. Results

3.1 Distribution and diversity of abundant and rare taxa

After quality control, 3 177 647 and 5 332 748 high-quality sequences were identified for bacteria and fungi in all samples, which were clustered into 2 766 and 2 162 ASVs, respectively (Table S2).Among abundant taxa, 0.94% of the bacterial ASVs and 4.17% of the fungal ASVs accounted for 96.76% and 94.68% of the total sequence number, respectively. By comparison, 99.06% of the bacterial ASVs and 95.83% of the fungal ASVs were identified as rare taxa and their sequence number was only 3.24% and 5.32%, respectively(Fig. 1A). At the genus level, the abundant bacterial subcommunity was only composed ofLactobacillusandAcetobacter, while the rare one contained 540 genera and was also dominated byLactobacillus(Figs. 1B-C). Abundant and rare fungal subcommunities were composed of 36 and 409 genera, respectively, withAspergillus,Pichia, andThermoascusas the dominant (Figs. 1D-E). Sourcetracking analysis showed thatDaquand PM were the primary sources of dominant abundant fungi, but not dominant rare bacteria (Fig. S1).

Fig. 1 Changes in microbial community composition of abundant and rare taxa by different combination patterns. (A) The proportion of each category of microbial taxa to the overall ASV and sequence numbers. (B-E) Relative abundance of abundant (B, D) and rare (C, E) bacterial and fungal genera in different treatments. Each treatment included 8 biological replicates and only the top 10 genera of each subcommunity were shown in the figure. (B, C) Bacteria; (D, E) Fungi.

The Chao1 and Shannon indices of rare taxa were significantly higher than those of abundant ones, and theα-diversity of rare bacterial and abundant fungal subcommunities was remarkably altered by combination patterns (Figs. S2A-B). Two-way ANOVA suggested that rare bacterial subcommunity was affected by the attributes of bothDaquand PM, and SD + FPM treatment exhibited the highest diversity and richness. By contrast, variations in the abundant fungal subcommunity were mainly related toDaquattributes, and itsα-diversity was increased by SD and FD. In addition, the results of community dissimilarity indicated that theβ-diversity of fungal communities also varied with combination patterns, especially the abundant subcommunity (Fig. S2C).

3.2 Environmental responses of abundant and rare taxa

The results of two-way ANOVA showed that the combination patterns ofDaquand PM significantly changed all physicochemical properties ofZaopeiexcept for moisture (Table S3). More specifically, acidity, pH, and starch consumption were closely related toDaquattributes. The former was positively correlated with TD,while the latter two with SD and FD, consistent with the differences in their microbial community and enzyme activity (Table S1)[32]. The RS content was mainly affected by PM attributes and positively correlated with CPM. By contrast,Zaopei’s ethanol content was significantly influenced byDaqu, PM, and their interaction, and SD + FPM and FD + CPM treatments sharply increased the ethanol yield.

We further explored the effects of these environmental changes on the structure and diversity of different subcommunities. Mantel tests based on overall communities indicated that the impacts ofZaopeiproperties on rare bacteria were stronger than on abundant bacteria,while on rare fungi were weaker than on abundant fungi (Table S4).Meanwhile, theα-diversity (Shannon and Chao1) of rare bacterial and abundant fungal subcommunities were significantly influenced byZaopeiproperties, pH/moisture and acidity/starch consumption were associated with the former and latter, respectively (Figs. 2A-B).Moreover, CAP analysis suggested that theβ-diversity of different subcommunities was also affected by environmental changes,especially the abundant fungi that were significantly correlated with pH, acidity, moisture, and starch consumption (Figs. 2C-F). Notably,these environmental changes related to community diversity were primarily induced by theDaquattributes (Table S3).

Fig. 2 Driving factors of microbial diversity of abundant and rare bacterial and fungal subcommunities in Zaopei. (A-B) Physicochemical properties affecting α-diversity were identified by the Mantel test. (C-F) Physicochemical properties influencing β-diversity were determined by constrained analysis of principal coordinates (CAP) based on Bray-Curtis distance at the ASV level. (C) Abundant bacteria; (D) Abundant fungi; (E) Rare bacteria; (F) Rare fungi. ΔStarch, the consumption of starch content during fermentation; RS, reducing sugar. *P ＜ 0.05; **P ＜ 0.01; ***P ＜ 0.001.

3.3 Assembly processes of abundant and rare taxa and their driving factors

The potential roles of deterministic and stochastic processes in the community assembly were quantified using βNTI and RCbray(Fig. S3). Abundant bacterial, abundant and rare fungal subcommunities were governed by the stochastic assembly, with the relative importance of 82.36%, 83.95%, and 70.83%, respectively,whereas rare bacterial subcommunity was mainly driven by the deterministic assembly (74.64%) (Fig. 3A). Different combination patterns altered the community assembly processes, and the contribution ofDaquattributes was higher than PM attributes(Figs. 3B-C). More specifically, SD and FD increased the proportion of heterogeneous selection in abundant fungal subcommunity and decreased in abundant bacterial subcommunity compared with TD,while there was no difference between CPM and FPM. For rare taxa, homogeneous selection (60.71%−92.86%) shaped the bacterial subcommunity, and the relative importance was improved by SD and FD regardless of PM attributes. By contrast, the community assembly of rare fungi was determined by both stochastic (46.43%−85.71%)and deterministic (14.29%−53.57%) processes, and the contribution of the latter was increased by SD and FD, especially in FPM.

Fig. 3 Deterministic and stochastic process in community assembly. (A) Relative importance of deterministic (heterogeneous selection, homogeneous selection)and stochastic (dispersal limitation, homogenizing dispersal, undominated) assembly processes in shaping abundant and rare bacterial and fungal sub-communities.(B) Changes in community assembly processes of abundant and rare taxa by different combination patterns. B1-B4. Abundant bacteria, rare bacteria, abundant fungi, rare fungi, respectively. (C) Impacts of Daqu and PM attributes on the assembly of abundant and rare bacterial and fungal subcommunities. C1-C4. Abundant bacteria, rare bacteria, abundant fungi, rare fungi, respectively.

The linear regression analysis between βNTI and each physicochemical property was conducted to explore the major driving factors of community assembly. The results showed that variations inZaopeiproperties had a greater impact on the assembly process of bacterial communities than that of fungal ones, especially the rare taxa (Fig. 4). The assembly processes of all subcommunities were significantly affected by changes in moisture, which were negatively and positively correlated with the βNTI values of bacterial and fungal communities, respectively. Changes in RS were also negatively associated with βNTI values of abundant bacterial and rare fungal subcommunities. However, environmental changes couldn’t alter the dominance of stochastic process (|βNTI| ＜ 2) in shaping these subcommunities. Conversely, the community assembly of rare bacteria was transitioned from stochastic process to homogeneous selection (βNTI ＜ –2) with increasing pH and moisture and decreasing ethanol. These results indicated that the assembly of rare bacteria inZaopeicould be regulated reasonably.

3.4 Co-occurrence network between abundant and rare taxa

The co-occurrence networks involving abundant and rare taxa were constructed based on the MENA pipeline to reveal the effect of combination patterns on community interactions (Fig. 5). According to Deng et al.[36], networks with similar cutoff thresholds can be compared and further analyzed (Table S5). All networks were mainly composed of rare taxa and positive connections, accounting for 56.93%−71.43% of the total nodes and 77.90%−96.25% of the total edges, respectively. All combination patterns increased the number of nodes and edges compared with the control (TD + CPM),indicating enhanced network complexity and interspecies interactions.Interestingly, TD and SD combined with FPM had a greater network size than those combined with CPM, while FD was the opposite(Fig. 5). The connections in the FD + CPM network were more than twice that of the FD + FPM network despite a similar nodes number, suggesting the better performance of specific combination patterns. Additionally, the proportions of rare-rare, rare-abundant, and abundant-abundant connections in all networks were 46.70%−65.11%,30.27%−44.56%, and 2.97%−12.98%, respectively (Table S5).

3.5 Changes in volatile flavors and its correlation with abundant and rare taxa

We identified 72 volatile flavors fromZaopei, including esters (41),alcohols (12), acids (6), carbonyls (6), phenols (3), and others (4).Among them, esters, acids, and alcohols dominated and accounted for 75.67%−86.40%, 8.18%−19.11%, and 3.33%−4.74% of the total concentration, respectively (Fig. S4A). The results of PCA and two-way PERMANOVA test showed that the flavor profile ofZaopeiwas significantly affected byDaqu, PM, and their interaction(Fig. 6A). For example, SD + FPM and FD + CPM treatments significantly increased the volatiles concentration ofZaopeicompared with the control (TD + CPM), while other treatments had a weak effect(Fig. S4B). Notably, we also observed a significant positive correlation between network size (number of nodes and edges) and total volatile concentration (Figs. S4C-D).

To determine which compounds were closely related to the attributes ofDaquand PM, four OPLS-DA models were constructed and 18 compounds were identified as differential metabolites based on variable importance in projection (VIP) > 1.0 (Fig. S5 and Table S6). Among the 13 compounds related to PM attributes, the contribution of 11 compounds such as hexanoic acid, butyric acid, and phenylethyl alcohol for FPM was higher, while that of ethyl lactate and ethyl palmitate was more important to CPM (Figs. S5A-B). By contrast, 15 compounds were associated withDaquattributes, and only ethyl palmitate dominated in TD (Figs. S5C-F). A total of nine compounds were responsible for the distinction between SD and FD,4 compounds (e.g., hexanoic acid and hexyl hexanoate) showed a greater contribution to the former, while 5 compounds (e.g., acetic acid and ethyl hexanoate) to the latter (Figs. S5G-H). Additionally,the cluster heatmap analysis showed most differential metabolites were enriched in SD + FPM and FD + CPM treatments (Fig. 6B),which further demonstrated the specific combination patterns ofDaquand PM had great potential in flavor formation.

The correlations between the 18 differential metabolites and the top 20 ASVs in each subcommunity were investigated based on Spearman correlation coefficients. As shown in Figs. 6C-F, the bacterial community contributed more to the variation in flavor profile than the fungal community. Moreover, rare bacterial subcommunity exhibited more associations than abundant ones regardless of the positive or negative correlations, while the associations of abundant and rare fungal subcommunities were similar. Therefore, the rare bacterial subcommunity might play a critical role in regulating the flavor profile ofZaopei.

3.6 Functional diversity of abundant and rare taxa

PICRUSt2 was performed to explore the functional diversity of abundant and rare taxa and their contribution to major metabolic pathways in SFB fermentation. The results showed that the functional diversity of abundant and rare fungal subcommunities was similar,which were annotated to 950 and 1 015 enzymes, respectively.Conversely, the functional diversity of rare bacterial subcommunity was significantly higher than that of abundant ones, which were annotated to 2 217 and 783 enzymes, respectively. As shown in Fig. 7, glycolysis/gluconeogenesis, fatty acid biosynthesis, and phenylalanine metabolism pathways were jointly regulated by bacterial and fungal communities, and rare taxa played an essential role in the expression of EC:2.3.1.41 and EC:1.4.3.21. Substrate degradation, ethanol production, and esters synthesis were mainly responsible for the fungal community, and the abundance of four esterases and the cellulose-utilization capacity were higher in rare subcommunity. By contrast, butanoate metabolism and fatty acid elongation pathways were closely related to bacterial communities,in particular, EC:1.1.1.157, EC:1.3.1.44, EC:2.8.3.8, EC:2.3.1.19,EC:2.7.2.7, EC:4.2.1.55, and EC:3.1.2.20 were only detected in rare subcommunity.

Considering the important contribution of butyric acid and hexanoic acid to SFB aroma characteristics, we further identified the main microbial predictors for their content inZaopeiby random forest analysis (Fig. S6). The results suggested that theα-diversity of different subcommunities was potentially more important to hexanoic acid content (R2= 0.281,P＜ 0.001) than butyric acid content (R2= 0.130,P= 0.017). The Shannon and Chao1 indices of rare bacterial subcommunity were found to be the most important variable for predicting the content of butyric acid and hexanoic acid,demonstrating that rare bacteria might play a vital role in the synthesis of butyric acid and hexanoic acid.

4. Discussion

4.1 Effects of combination patterns on the abundant and rare microbial taxa

Generally, the distribution of microorganisms in a natural community is skewed, with few abundant taxa coexisting alongside a large number of rare taxa[20]. Recently, similar observations have been documented in traditional fermented foods, and some rare taxa were identified as core functional microbes[18-19]. In this study, 0.94%of abundant bacteria and 4.17% of abundant fungi accounted for 96.76% and 94.68% of its total sequences, respectively (Fig. 1A),resulting in a significantly higherα-diversity in rare taxa than in abundant ones (Figs. 2A-B), which was consistent with several natural ecosystems[27,39]. The community structures of abundant and rare taxa were mainly composed ofLactobacillus,Acetobacter,Thermoactinomyces,Bacillus,Aspergillus,Pichia, andThermoascus(Figs. 1B-E), these genera were also prevalent in previous studies[4-5,16]. The community structure and diversity of abundant and rare taxa exhibited a distinct response to different combination patterns, and rare bacterial and abundant fungal subcommunities were more sensitive (Figs. 1 and S2). Prior works reported that the response of microbial communities to external disturbances might depend on their environmental breadth (i.e., the resistance to variation in environmental factors)[21,39]. Therefore, this differential response pattern might be attributed to the fact thatZaopeiproperties had a greater impact on rare bacterial and abundant fungal subcommunities(Fig. 2 and Table S4).

Compared with the control (TD + CPM), different combinations decreased the relative abundance of rareLactobacillusand abundantAspergillus, whereas increased that of rareBacillus,Kroppenstedtia, andThermoactinomycesas well as abundantPichiaandThermoascus(Figs. 1C-D). It is reported that excessive amounts and rapid succession ofLactobacillusoften destroyed the balance of community and metabolism during Baijiu fermentation[40-41],while increasedBacillusandPichiawere beneficial for flavor- and ethanol-producing due to their metabolic function[42-43]. Additionally,thermophilic microorganismsKroppenstedtia,Thermoactinomyces,andThermoascusalso promoted the synthesis of multiple volatile flavors[44], andThermoactinomyceswith carboxylesterase genes had previously been identified as the core bacteria in SFB brewing microecosystem[16]. Therefore, an improved flavor quality and ethanol yield might be partly attributed to the enrichment of these microorganisms(Fig. S4B and Table S3). Notably, the source-tracking analysis indicated thatDaquand PM were the important sources of these fungi but not these bacteria (Fig. S1), which was in line with the results reported by Tan et al.[4]and Zhao et al.[18]. Recently, researchers investigated the microbiota of SFB brewing micro-ecosystems comprehensively and the results highlighted that functional bacteria can colonize the environment and contribute to fermentation by transferring intoZaopei[17]. Similarly, Wang et al.[45]also claimed that the brewing environment contributed to most of the bacterial communities inZaopeiduring Baijiu fermentation. Therefore, future research should pay more attention to the exploration and regulation of environmental microorganisms.

4.2 Effects of combination patterns on community assembly of abundant and rare taxa

Elucidating the relative importance of deterministic and stochastic processes to community assembly is a focus for microbial ecology and may contribute to our understanding of community evolution[25].Here, we used βNTI and RCbaryto quantify the ecological processes of community assembly inZaopei(Fig. S3). The assembly pattern of abundant and rare bacterial subcommunities differed significantly,with stochastic and deterministic processes dominating the former and the latter, respectively (Fig. 3A). Similar results have been reported in some natural ecosystems, such as coastal wetlands[24].The distinct assembly patterns between abundant and rare bacterial subcommunities could be attributed to the differential response to external disturbances and differences in niche breadth. Generally, rare taxa have a very narrow environmental niche and are less resistant to environmental changes, whereas abundant taxa occupy a broad niche and are more adaptable to environmental disturbances[20,24].This was confirmed by more significant correlations between the βNTI values of rare bacteria and physicochemical properties(Fig. 4). Specifically, community assembly of rare bacteria was associated with pH, moisture, and ethanol, while that of abundant bacteria with RS and moisture. Moreover, the community diversity of rare bacteria was also influenced by pH and moisture, and that of abundant bacteria by RS (Fig. 2). Changes in pH and RS were closely related to the attributes ofDaquand PM, respectively, whereas variations in moisture might cause by the initial parameter ofZaopei(Table S3). Despite the importance of pH and RS in regulating bacterial community assembly has been reported during Baijiu fermentation[31], our findings highlighted that the community assembly of rare bacteria was more affected by pH, while RS was crucial for the abundant ones. Therefore, different combination patterns ofDaquand PM altered the pH and RS ofZaopeiand further regulated the community assembly of rare and abundant bacterial subcommunities.

Fig. 4 Linear regression analysis between βNTI values of each subcommunity and changes in Zaopei physicochemical properties. Horizontal-dashed lines indicate the βNTI significance thresholds of +2 and −2. (A-H) Abundant bacteria, rare bacteria, rare fungi, abundant bacteria, rare bacteria, rare fungi, abundant fungi, rare bacteria, respectively.

In contrast to the bacterial community, we found that the assembly pattern of abundant and rare fungal subcommunities was similar and mainly driven by stochastic process involving an undominated process and dispersal limitation (Fig. 3A). The stochastic assembly of fungal communities was as expected because variations inZaopeiproperties had a weak effect on the community succession of fungi during SFB fermentation[18,30]. Similarly, prior works suggested that the community assembly of fungi in fermentation micro-ecosystems was mainly determined by stochastic process[31,46]. However, the absolute dominance ofSaccharomyces cerevisiaeresulted in the fungal community assembly governed by homogeneous dispersal during Xiaoqu Baijiu fermentation[31], which was different from the present results. This discrepancy might be related to the process characteristic and microbial composition of the two kinds of Baijiu.The fermentation period of SFB (60−90 days) is much longer than that of Xiaoqu Baijiu (7−15 days), the long-term environmental stress (e.g., no oxygen, high acidity and ethanol) sharply reduced the fungal biomass ofZaopeiat the later stage of fermentation[30], andS. cerevisiaewith tolerance to acid and ethanol was also disappeared at this stage[47]. Therefore, the fungal community assembly in this study was mainly governed by an undominated process related to the random birth and death of species. In addition, all combination patterns except TD + FPM increased the proportion of deterministic process that was more important in the rare fungal subcommunity compared with the control (TD + CPM) (Figs. 3B-C). It has been demonstrated that the accumulation of environmental stress was positively associated with the relative importance of deterministic processes[29]. Therefore, we observed more associations between βNTI values of rare fungi and changes in physicochemical properties(Fig. 4), which further supported that the rare taxa had a narrow environmental niche[20].

4.3 Rare taxa play critical roles in community interaction and flavor formation

More than 60% of the volatile compounds inZaopeiwill be transferred to the fresh Baijiu during distillation, which determines the aroma and quality characteristics of SFB[48]. Two-way PERMANOVA test showed thatDaqu, PM, and their interaction were responsible for these variations, and the concentration of 18 compounds varied with combination patterns, such as hexanoic acid, butyric acid, ethyl hexanoate and ethyl lactate (Figs. 6A-B and Table S6). These volatile compounds have been identified as the aroma-active compounds of SFB[49], which supported the view that high-qualityDaquand PM could improve the flavor profile ofZaopeiand the sensory scores of Baijiu[5,9,12]. Generally, volatile compounds are synthesized under the synergistic effect of multiple species, thereby the flavor profile ofZaopeiwas influenced by community interactions[50]. Moreover,a recent study highlighted that the flavor quality ofZaopeiwas also closely related to the topological properties of co-occurrence networks[51]. In this study, SD + FPM and FD + CPM treatments remarkably enhanced the interspecies interactions in co-occurrence networks (Fig. 5), and a positive correlation between network size and flavor concentration was observed (Figs. S4C-D), which might explain their significantly improved flavor quality (Fig. S4B). These results suggested that optimizing the combination pattern ofDaquand PM might be an effective strategy to improve the flavor quality of SFB.

Fig. 5 Co-occurrence network analysis based on abundant and rare bacteria and fungi in Zaopei. Each node represents an ASV, and its size is proportional to the number of connections. Each edge indicates a significant Pearson correlation, red and green represent positive and negative correlations, respectively.(A) TD + CPM; (B) TD + FPM; (C) FD + CPM; (D) FD + FPM; (E) SD + CPM; (F) SD + FPM.

In addition, rare taxa dominated all co-occurrence networks,and the proportion of rare-rare and rare-abundant connections was obviously larger than that of abundant-abundant connections (Fig. 5 and Table S5). Similar observations have been documented in natural ecosystems and fermentation micro-ecosystems[18,39], suggesting that rare taxa might serve as hubs of interaction networks to maintain the composition and functioning of a microbiome[52]. Considering the significant contribution of rare bacterial subcommunity to variations in flavor profile (Fig. 6), the functional diversity of different subcommunities was predicted by PICRUSt2. We found that the rare bacterial subcommunity had the highest functional diversity, in particular, several key enzymes related to the synthesis of butyric acid and hexanoic acid (e.g., EC:1.1.1.157, EC:2.8.3.8, EC:2.3.1.19,EC:2.7.2.7, and EC:4.2.1.55) were only detected in it (Fig. 7). More importantly, the results of random forest model indicated that itsα-diversity was the most significant variable for predicting the content of butyric acid and hexanoic acid inZaopei(Fig. S6). Similarly, some nitrogen cycling genes in soil were also carried only by rare bacteria,which were more important than abundant ones in predicting crop yield[22]. The synthesis of butyric acid and hexanoic acid was closely related to the PM microbial consortia during SFB fermentation[3],but they could also be generated and accumulated in the absence of PM[40]. The present study indicated that this might be associated with the rare bacterial subcommunity inZaopei. A recent study also provided evidence that rare bacteria were the keystone taxa in fermentation micro-ecosystems, the bioaugmentation of rare species(Komagataeibacter europaeus) improved the community stability and flavor quality[19]. These results highlighted the importance of rare taxa in maintaining ecological functions, enhancing community interactions, and improving metabolic activities, which may be attributed to the fact that rare taxa were the “seed bank” of genetic and functional diversity in ecosystems[20].

Fig. 6 Changes in the volatile flavor of Zaopei and its correlation with abundant and rare bacterial and fungal subcommunity. (A) Effects of Daqu, PM and their interactions on the flavor profile. (B) Cluster heatmap analysis of differential metabolites obtained from OPLS-DA models based on VIP ＞ 1.(C-F) Flavor associations of the top 20 most abundant ASVs in each subcommunity, and only significant associations were shown in the plot (P ＜ 0.05). Red and blue represented positive and negative correlation, respectively, and the shade of the color was proportional to the correlation coefficient.

Fig. 7 Functional diversity of abundant and rare bacterial and fungal subcommunities and their contribution to major metabolic pathways of Zaopei fermentation based on PICRUSt2 prediction.

5. Conclusion

This study is the first to reveal the environmental response,assembly pattern, and metabolic potential of abundant and rare taxa inZaopeifrom the perspective of microbial ecology. The results showed that the composition and diversity of abundant and rare taxa exhibited unique responses to combination patterns ofDaquand PM, and rare bacterial and abundant fungal subcommunities were more sensitive.The community assembly of abundant and rare bacteria was governed by stochastic process and homogeneous selection, respectively, while that of fungi was dominated by stochastic process and less affected by changes in physicochemical properties. Moreover, the co-occurrence network and flavor profile were significantly improved by SD + FPM and FD + CPM treatments, and rare bacteria played a dominant role in maintaining the interspecies interactions and metabolic function of communities, especially the synthesis of hexanoic acid and butyric acid. Our study provides new insights for regulating Baijiu fermentation and improves our understanding of the key roles of rare taxa in brewing micro-ecosystems. However, the current results are partly obtained from functional predictions, for which the specific species at play and their source are unclear. Future work needs to consider the composition and diversity of rare taxa in the SFB brewing environment and further isolation and functional studies on these key rare taxa by multi-omics.

Conflicts of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgement

This research was supported by the Cooperation Project of Luzhou Laojiao Co., Ltd. and Sichuan University (21H0997).

Appendix A. Supplementary data

Supplementary data associated with this article can be found, in the online version, at http://doi.org/10.26599/FSHW.2022.9250121.