空间异质性对大曲微生物群落的影响

唐慧芳，黄 钧，周荣清,2,*，秦 辉，张宿义，董 异，王 超，王小军，母 雨，潘强林

（1.四川大学轻工科学与工程学院，四川 成都 610065；2.国家固态酿造工程技术研究中心，四川 泸州 646699；3.泸州老窖股份有限公司，四川 泸州 646699）

大曲在我国传统白酒的生产中具有重要作用，一直都是行业关注的热点。大曲是由未灭菌大麦、小麦或豌豆等原料在开放的环境中，籍以源于两者的微生物群落，在驱动因素的胁迫作用下，形成包含功能菌群、粗酶制剂、挥发性物质及前体代谢成分等颇具特色的必需酿酒原料之一[1-5]。这些驱动因子包括环境温度、湿度、曲心温度、水分含量和酸度。后三者可通过适度调节前二者进行调控[6]，因此将前两者视为环境因子，后三者视为内生因子。这些驱动因子调控曲坯微生态体系中种（株）间互作关系、营养网络等成为影响大曲属性的主要因素[5]。

因所使用原料及发酵过程峰值温度的差异，可将大曲分为3 种不同的类型，峰值温度在62～65、60～62 ℃及50 ℃左右的分别是高温、中高温和低温大曲，这些大曲分别用于酱香型、浓香型和清香型白酒的生产[7-8]。高温大曲以Kroppenstedtia、Saccharopolyspora、Bacillus、Weissella和Lactobacillus为主要优势细菌属，真菌则以Thermoascus、Thermomyces和Aspergillus为主[9]。中温大曲优势真菌属主要包括Issatchenkia、Thermoascus、Trichocladium、Aspergillus、Rhizopus和Rhizomucor，优势细菌属为Weissella和Lactobacillus[10]。低温大曲中优势细菌属包括Bacillus和Lactobacillus，优势真菌属包括Rhizopus和Saccharomyces[11]。

近年来，基于功能菌株（群）强化大曲发酵是改善品质的热门研究方向之一，且已取得良好结果。如Wang Peng等[12]接种了Bacillus licheniformis，既改善了大曲风味特征，又提高了淀粉水解酶活力。He Guiqiang等[13]接种B.velezensis和B.subtilis显著提高了淀粉水解力和酯化力，也提高了酯类、吡嗪类和醇类等特征成分含量。宇宙射线诱变是一种独特育种手段，具有效率高、突变范围宽、周期短等特点，应用于酿酒微生物菌种的选育具有广阔前景[14]。

本研究以接种太空诱变曲扩培的强化大曲为对象，探讨在大型曲房中心位置8 层货架式曲架的上、中、下层曲（分别为第1～3、4～5、6～8层）在主发酵期中理化性质与微生物群落结构的变化趋势，研究驱动因子对群落多样性的影响规律。旨在揭示多层发酵房驱动因子与大曲群落的相关性，为大曲生产过程的自动化调控提供理论依据。

1 材料与方法

1.1 材料与试剂

大曲样品取自四川泸州著名酿酒企业制曲中心。

化学试剂均为分析纯，购自成都金山化学试剂有限公司；2×TaqPCR Master Mix、ITS1/ITS4引物、DNase/RNase-free Deionized Water、真菌基因组DNA抽提试剂盒 生工生物工程（上海）股份有限公司。

1.2 仪器与设备

TGL 16M高速冷冻离心机 湘麓离心机仪器公司；S100TMThermal Cycler聚合酶链式反应（polymerase chain reaction，PCR）仪、Gel DocTMXR凝胶成像仪 美国Bio-Rad公司；TU-1901紫外分光光度计 北京普析通用仪器有限责任公司；85-2型数显恒温磁力搅拌器 上海 双捷实验设备有限公司；TES1360A温湿度测试仪 泰仕电子工业股份有限公司；TP3001电子探针温度计 天津凯士达仪器仪表有限公司。

1.3 方法

1.3.1 样品制备与收集

生产流程及设施简图如图1所示。

图1 架子曲生产流程（A）与曲房布局（B）Fig.1 Flow chart of the production process of Daqu (A) and layout of Daqu fermentation room (B)

大曲生产是接种1%（m/m）的母曲，按照DB 5105/T 36—2020《地理标志产品 泸州老窖系列曲药（久香牌）生产技术规范》[15]和Zheng Xiaowei等[16]所述过程及参数控制进行。母曲是搭载在神州11号太空舱并在太空诱变1 个月曲粉逐级扩大的陈曲。翻曲是为了通风散热，调节曲块水分，为曲坯发酵提供适宜的温、湿度环境[17]，架子曲以空气导流调控装置调节室内温、湿度与大曲曲心温度，模拟人工翻曲，实现无翻曲的生产工艺。

取样：分别从曲房中心曲架的上、中、下层随机取2 块大曲，粉碎混匀后堆积成锥型，分别从曲粉堆的四周与中心点5 点共取样300 g，再次混匀后分为两份，其中一份280 g置于4 ℃保藏用于理化检测，另一份20 g则置于-80 ℃保藏用于微生物群落结构分析。样品名称及相应信息如表1所示。

表1 取样点的信息及编号Table 1 Information and codes of sampling sites

1.3.2 微生物群落组成检测

总DNA提取：使用Fast DNA SPIN kit，按照供应商提供操作程序提取后，分别经1%琼脂糖凝胶电泳和Nano Drop ND-1000分光光度计检测其纯度和含量。分别使用引物338F/806R和ITS5/ITS1 PCR扩增细菌的16S rRNA V3～V4区和真菌的ITS1区域。扩增体系25 μL：5×reaction buffer和5×GC buffer各5 μL，0.25 μL DNA聚合酶（5 U/μL，Q5 High-Fidelity），2 μL dNTPs（2.5 mmol/L），正反引物各1 μL（10 μmol/L），2 μL DNA模板，8.75 μL ddH2O。细菌扩增程序：98 ℃预热2 min，98 ℃变性15 s，55 ℃退火30 s，72 ℃延伸30 s，经过25 个循环，最后72 ℃保持5 min。真菌扩增程序：95 ℃预热3 min，95 ℃变性30 s，61 ℃退火30 s，72 ℃延伸45 s，经32 个循环，最后72 ℃保持10 min。PCR产物琼脂糖凝胶电泳进行目的片段纯化回收。纯化样品后委托派森诺生物科技有限公司在Illumina NovaSeq平台完成DNA片段进行双端（Paired-end，2×300）测序。

DADA2序列去噪：首先调用QIIME切除序列的引物片段，弃去未匹配引物的序列；然后调用DADA2进行质控、去噪、拼接、去嵌合体获得扩增子序列变体（amplicon sequence variants，ASV）。接着通过Greengenes/UNITE数据库，将ASV特征序列与数据库中的参考序列相比对，获取每个ASV所对应的分类学信息，去噪后，合并ASV特征序列和ASV表格，并去除singletons ASV。

1.3.3 理化性质的检测

参照QB/T 4257—2011《酿酒大曲通用分析方法》[18]所述方法和步骤测定水分、酸度、液化力、糖化力、发酵力和酯化力。环境温度与湿度使用温湿度测试仪测定，曲心温度使用电子探针温度计测定。所有实验均进行3 次重复。

1.4 数据处理

使用OriginPro 2022对驱动因子以及大曲酶特性进行可视化。使用派森诺基因云平台（https://www.genescloud.cn/home）对微生物群落进行主坐标分析（principal co-ordinates analysis，PCoA）、冗余分析（redundancy analysis，RDA）以及基于average的层次聚类分析（hierarchical clustering analysis，HCA）。使用基迪奥生信云平台（https://www.omicshare.com/tools/home/task/index.html）和图图云平台（https://www.cloudtutu.com/#/login）对驱动因子和优势微生物属（相对丰度＞1%）分别进行普氏分析与Mantel检验。使用IBM SPSS Statistics 25软件对优势微生物属通过Spearman算法（r＞0.70，P＜0.05）构建共现网络，经Cytoscape可视化。

2 结果与分析

2.1 发酵过程空间异质性驱动因子的变化趋势

如图2所示，上、中、下三层大曲的水分含量逐步降低，到10 d均降低至15%以下，达到转房标准[19]。不同层次环境温度、湿度及大曲曲心温度的变化趋势类似。发酵起始2 d后，因菌群繁殖与代谢释放大量热，曲坯的水分大量散发[20]，曲房湿度陡增，热量跟随水蒸气上升，导致下层温度、湿度低于上、中层。随着发酵过程的进行，中、下层大曲的曲心温度略高于上层大曲，可能与环境温、湿度和微生物群落组成有关。发酵初期，曲坯水分、温度适合产酸微生物繁殖，所以大曲酸度陡升[21]，尤其是上层曲坯酸度高达1.7 mmol/10 g，随后又下降直至达到稳定的状态。同期上层曲坯的酸度略高于中、下层曲，可能与上层大曲高温、高湿环境的周期更长，从而产酸更多有关[22]。

图2 驱动因子的变化趋势Fig.2 Changes in driving factors during fermentation

2.2 发酵过程酶活力的变化

由图3A可知，发酵前2 d，各层样品均未检出液化力。第3天时略升，第8天时上层、中层液化力趋于稳定。下层曲的液化力8 d后速增，可能是下层温度偏低，菌群生长缓慢，高湿度环境更适合霉菌生长所致。3 个不同层次大曲的糖化力变化趋势相同，但强度不同（图3B）。糖化酶活力起伏振荡波动后在6 d后逐渐趋于稳定的趋势，类似邢钢等[23]报道的结果。这是原料所含内源酶被激活、后遇高温失活，同时因微生物菌群分泌糖化酶随生物量增高而提高，并伴随高温和酸胁迫影响的综合结果。类似前述的不同层次间温度变化趋势，糖化力也是上层曲的波动程度最大，下层曲的波动最小。且下层曲糖化力明显高于中层和上层曲，可能与下层前期温度低且升温缓、湿度适宜密切相关。

图3 发酵过程中不同层次大曲酶活力的变化趋势Fig.3 Changes in enzymatic activities in different layers of Daqu during fermentation

上、中、下三层曲的发酵力也类似陈晓茹等[24]报道的结果，显著降后陡然升，又降低至趋于平稳的趋势（图3C），与酸度的变化趋势相反（图1）。发酵力表征的是大曲中的酶利用糖转化为CO2与乙醇的能力，与酵母菌的活力密切相关[23]。发酵力的变化趋势类似温度。而酯化力的变化趋势则与液化力类似（图3D）。酯化力同样与霉菌和酵母的繁殖代谢密切相关[25-27]，霉菌和酵母发酵前期生长较为缓慢，随后生长旺盛[27]，酯化力则陡增，随后趋于稳定，且在发酵后期，上、下层样品的酯化酶活力显著高于中层大曲。

2.3 微生物群落结构的差异

样品中检出真菌和细菌有效序列分别是53111～88490和67181～115456，其高质量序列是52943～88244和45465～101691。3 层大曲微生物群落的α多样性指数均呈现波动增加的趋势（图4）。上层、下层大曲细菌的α多样性指数均在8 d达到顶峰，中层则是在6 d，且明显低于上、中层样品。3 层大曲真菌的α多样性指数波动相似，仅上层第8天样品有差异。

图4 不同层次对微生物群落α多样性的影响Fig.4 Effects of different layers on the α-diversity of microbial communities in Daqu

微生物群落主要由Firmicutes和Proteobacteria 2 个细菌门及Ascomycota、Mucoromycota和Basidiomycota 3 个真菌门组成，包括56 个细菌属和59 个真菌属，平均相对丰度大于1%的优势菌属包括Lactobacillus、Weissella、Bacillus、Kosakonia、Staphylococcus、Thermoactinomyces等12 个细菌属（图5A）和Pichia、Thermoascus、Rhizomucor等7 个真菌属（图6A）。不同空间位置环境温度、湿度等参数及其变化显著影响生产过程曲坯的微生物多样性和丰富度[21]。如中层和下层的Firmicutes相对丰度始终在85.16%～99.03%之间，上层曲则是前5 d都是Firmicutes占绝对优势，而从6 d起，Proteobacteria丰度陡增，占据主导地位。

图6 不同层次对真菌群落结构和组成的影响Fig.6 Structure and composition of fungal community in different layers of Daqu

如图5A所示，整个过程中，Weissella和Lactobacillus在发酵过程中的相对丰度交替变化，如0 d 曲坯中Weissella相对丰度高达81.77%，随着发酵过程的进行，其在各层曲的丰度迅速降低至2%后逐渐增高，然后又降至约1%，类似Zhu Min等[28]报道的结果。Lactobacillus经1 d发酵就达到峰值，上、中层曲中相对丰度分别是92.44%和90.08%，随之渐降，8 d降至1%以下，而下层在4 d才达到峰值（约为70%），结束时降至42.89%。这些趋势同发酵过程中酸度增高，水分逐步降低，适合Lactobacillus繁殖密切相关[28]。0 d曲坯中Bacillus相对丰 度低于1%，发酵前期，除1 d下层样品外，中层和上层相对丰度都在1%左右，在发酵后期相对丰度提升，因为Bacillus可产耐热芽孢，可以承受高达60～62 ℃高温[29]，故在发酵8 d中、下层曲中丰度分别达到了57.53%、39.39%，在上层曲中虽不占主导，但也稳定在10%左右。

初始曲坯中，Kosakonia相对丰度仅为8.24%，空间位置显著影响其相对丰度，主发酵结束时，上、中、下层曲相对丰度分别为30%、2.19%和5.38%。Staphylococcus在发酵过程中最高相对丰度分别达5.18%、15.79%、9.75%。类似地，在发酵结束时中层和下层曲中Thermoactinomyces相对丰度分别是27.84%和8.96%，而Pantoea和Kroppenstedtia仅在上层和中层分别约为20%和22%，在其余层次中的相对丰度则小于1%，是非优势属。

细菌群落基于Bray-Curtis距离的PCoA结果表明，大曲微生物群落时空特征差异显著（图5B）。具体地，0 d样品独聚一簇，上层曲1～5 d和6～10 d两个阶段成簇，且相距较远，表明群落结构差异显著。中层和下层曲均聚类成3 簇，前者是1～5、6 d和8～10 d分别成簇，前两簇距离较近；下层则是8 d和10 d分别单独成簇，两者距离较远，1～6 d结构相似且距离紧密。发酵前5 d的3 层曲细菌群落各自聚类在同簇，随着发酵的进行，空间位置则显著影响群落结构，聚类距离趋远。这些结果也被HCA验证，显示出与PCoA相似的群集模式（图5C）。

Pichia、Thermoascus、Rhizomucor等7 个属为优势菌（图6A）。整个主发酵期间，Ascomycota始终占优势，其相对丰度最高可达99.93%。起始时（0 d），Pichia相对丰度高达99.38%，前5 d的3 个层次大曲Pichia相对丰度都维持在96.70%～99.82%之间。6 d后，上层和中层曲中Pichia相对丰度便逐步下降，而下层曲直到8 d，Pichia相对丰度仍高达97.40%，10 d才降至68.88%。后期，中层曲Pichia相对丰度较下层曲低，但仍在60%以上；上层曲则降至10.76%。在上、中层和下层大曲分别经6 d和8 d的发酵后，Thermoascus、Rhizomucor与Rhizopus都成为优势菌。上层曲因中后期温度较高，环境湿度和曲坯水分显著降低，抑制酵母菌和霉菌等热敏真菌繁殖[21,30]，所以Thermoascus[10]和Thermomyces[31]等嗜热真菌[32]变为优势菌，与Xiao Chen等[30]报道的结果相符。

基于真菌多样性的PCoA和HCA都表明（图6B、C），发酵前期曲坯的真菌群落结构均相似，上层和中层发酵0～5 d曲坯聚类成簇、下层则是0～8 d聚集一簇，6～8 d中层曲聚类成另一簇，上、中、下3 层发酵10 d曲坯单独成簇。

2.4 变量对群落结构的驱动作用

使用普氏分析与Mantel检验研究驱动因子对微生物群落演替的贡献（图7）。普氏分析表明驱动因子对细菌群落的贡献显著，而对真菌贡献不显著。Mantel检验显示曲坯水分含量与环境湿度显著影响上层曲坯的微生物组成，而酸度与曲心温度则次之。同时，水分含量与酸度均显著影响中层曲坯的微生物群落，环境湿度对其真菌组成影响较细菌群落大。而在下层大曲中，水分含量对细菌和真菌均有显著影响，而酸度仅对细菌群落影响显著。结合普氏分析与Mantel检验的结果表明，驱动因子对细菌群落的贡献更显著。

图7 微生物群落与环境因子间的普氏分析（A、B）和 Mantel检验（C～E）Fig.7 Procrustes analysis (A,B) and Mantel test (C–E) between microbial communities and driving factors

通过方差分解分析（variance partitioning analysis，VPA）辨析内生因子与环境因子对微生物群落的影响强度，如表2所示。驱动因子对上层大曲的细菌群落结构影响显著，中层次之，而下层大曲细菌结构仅受内生因子的显著影响；仅中层大曲真菌群落受驱动因子显著影响。另外，在上层曲中，内生因子对真菌群落影响强于细菌，而环境因子仅对细菌群落影响显著。在中层大曲中，驱动因子对真菌和细菌群落的贡献度相同，其中内生因子的贡献度略大于环境因子。而在下层大曲中仅内生因子对微生物群落贡献显著，且对细菌和真菌群落的贡献度相同。其结果与普氏分析和Mantel检验一致。

表2 驱动因子与微生物的VPATable 2 VPA between driving factors and microbial communities

使用RDA明确驱动因子与具体优势菌属的相关性。如图8所示，水分和酸度均与Lactobacillus和Pichia呈正相关。而酸度与Bacillus和Rhizopus等呈负相关，不同于Zhu Min等[28]的研究，可能与其优势微生物相对丰度以及生产工艺不同有关。上层曲心温度和室温呈负相关，而在中、下层两者呈强正相关，这可能由上、中、下层空气流动程度不同等条件差异导致。在上层大曲中，曲心温度与Pantoea、Bacillus、Kosakonia、Staphylococcus、Rhizomucor、Thermoascus等呈正相关。中层大曲曲心温度与Thermoactinomyces、Kroppenstedtia、Scopulibacillus、Kosakonia、Pantoea等呈正相关，下层大曲曲心温度Pichia、Thermomyces、Rhizomucor、Thermoascus等呈正相关。这些结果很好证明了Bacillus、Rhizomucor、Thermoascus、Thermomyces等耐高温能力强[30,33]。

图8 微生物群落与驱动因子间的RDAFig.8 RDA between microbial communities and driving factors

2.5 微生物间相关性

如图9所示，在上层大曲中存在81 对互作关系，其中53 对为正相关；在中层大曲与下层大曲中分别存在55 对和43 对互作关系，其中正相关关系分别为35 对和36 对。上、中、下层大曲中共享的仅有10 对关系，包括Lactobacillus与Acinetobacter、Rhizomucor这两对负相关，以及Bacillus与Acinetobacter、Thermoascus、Rhizomucor，Enterobacter与Kosakonia、Pantoea，Pantoea与Kosakonia，Rhizomucor与Thermoascus、Rhizopus这8 对正相关。Lactobacillus与Rhizomucor呈负相关可能是因为乳杆菌生产乳酸能力强，高酸度环境不利于Rhizomucor这类不耐酸菌属生长[34]，在前期Lactobacillus相对丰度高，抑制了Rhizomucor繁殖，后期Lactobacillus相对丰度降低，酸度也下降，提高了Rhizomucor的丰度。Bacillus和Thermoascus呈正相关可能与它们都是耐热菌有关[4,30]。除此之外，8 对相同的关系在上层和中层曲中，上层与下层大曲和中层与下层大曲中分别存在18 对和5 对相同的关系。综上所述，不同层次大曲微生物群落结构差异显著的同时，微生物间的共生共斥关系亦不同，表明不同层次间驱动因素的差异可能也导致微生物间的互作关系发生改变，从而共同促进大曲群落形成。

图9 微生物共现网络在不同层次间的差异Fig.9 Differences in microbial co-occurrence networks among different layers of Daqu

3 结论

以多层货架式工艺生产的主发酵期上、中、下3 层大曲为对象，探讨了其驱动因子和酶活力的变化过程，以及驱动因子对大曲微生物群落的贡献。结果表明，3 层曲的水分均在10 d内逐步降低至15%以下，达到转房标准。不同层间环境温度、湿度及曲心温度的变化程度不同，但波动趋势相同。同期的上层曲坯酸度略高于中层和下层，酶活力则反之。在大曲微生态体系中，微生物群落的时空性特征显著，尤其是发酵中后期，PCoA和HCA也揭示了其规律性。VPA、普氏分析和Mantel检验都表明，驱动因子显著影响细菌群落组成，但不同驱动因子对不同层次大曲的影响程度不同。具体地，驱动因子均显著影响上、中层曲的微生物群落结构，而下层曲仅受内生因子的影响显著。另外，内生因子对上层曲的真菌群落贡献强于细菌，而在中、下层曲中贡献度相同。环境因子仅对上层细菌群落和中层微生物群落具有显著影响。其中，水分和酸度均与Lactobacillus和Pichia呈正相关，但不同层次驱动因子间相关性不同。同时，不同层次大曲中微生物间相互作用也有差异。综上，多层货架式工艺生产大曲过程中，不同层次间大曲酶特性及群落结构的差异与驱动因子显著相关，可通过调控发酵过程中驱动因子从而对大曲质量进行调控，这为大曲的智能制造提供了科学依据。