遂宁榨菜发酵过程中细菌群落多样性和基因功能预测分析

郝卓莉

（石家庄职业技术学院 食品与药品工程系，河北 石家庄 050081）

榨菜的原料青菜头，学名“茎瘤芥”，分布于我国各地，收藏季节短，不易保存，因此大部分用于制作榨菜[1]。其传统的制作工艺用盐量高，抑制了高品质的乳酸菌活性，同时对产品的风味、营养及安全性等品质产生诸多不良影响，不符合现代社会对于营养、天然和绿色生活方式的追求[2-3]。而高盐的腌制会伴随着大量亚硝酸盐的生产，亚硝酸盐可与蛋白质中的仲胺结合生成亚硝胺，进而促进癌细胞转化和增殖。此外，大量的亚硝酸盐会导致肠原性青紫症[4]。低盐腌制榨菜不仅降低了含盐量，还可以增加酸度及减少亚硝酸盐含量[5]。此外，低盐工艺有利于酵母菌、乳酸菌等有益微生物生长繁殖，可根据理化因素调控进程，以达到改善榨菜品质、安全性及营养价值的目的[6]。随着生活方式向绿色化及健康化转变，榨菜的低盐化是未来的发展趋势[7]。遂宁是我国榨菜的主要产地之一，其通常采用低盐发酵工艺生产榨菜，生产的榨菜已成为当地的特产之一。遂宁榨菜坚持使用传统发酵工艺，采用低盐腌制榨菜，而低盐发酵过程同样也会经历多种微生物的共同演替及复杂的代谢变化。因此，了解遂宁榨菜低盐发酵过程中微生物群落多样性及基因功能的演替显得尤为重要。

近年来，大量分子生物学方法被用来分析食品中的微生物多样性[8]。其中，高通量测序技术是最好的选择之一，因为它能产生大量有用的信息，为研究人员提供了一个相对公正的视角来看待复杂生态系统中的微生物群落组成[9]。LIANG H等[10]通过高通量测序技术发现，乳酸菌是传统榨菜发酵过程中最重要的微生物之一，呈先迅速增加后维持稳定的变化趋势，且在前期以明串珠菌属（Leuconostoc）和魏斯氏菌属（Weissella）为主，中期以乳杆菌属（Lactobacillus）为主，而后期则以乳杆菌属和片球菌属（Pediococcus）为主。尽管发酵过程中微生物复杂多样，但乳酸菌始终是榨菜发酵过程中的优势菌群、主导着榨菜的发酵[9]。目前，很多乳酸菌（如肠膜明串珠菌，短乳杆菌和植物乳杆菌）已被用作发酵剂，以加快榨菜的发酵过程，提高产品品质[11-12]。然而，在遂宁榨菜工业生产中，微生物多样性及其基因功能在发酵过程中的作用仍有待进一步研究。

本研究利用Illumina Miseq高通量测序技术对遂宁榨菜发酵过程中的细菌群落多样性进行解析，同时结合基于标记基因序列来预测功能丰度（phylogenetic investigation of communities by reconstruction of unobserved states，PICRUSt）软件对细菌种群的基因功能的演替进行预测，旨在为工业化生产提供理论依据及筛选特性菌种提供数据参考。

1 材料与方法

1.1 材料与试剂

榨菜发酵液样品：取四川省遂宁市某调味品有限公司两个不同发酵池同一位置发酵1 d、15 d、30 d、45 d和60 d的榨菜发酵液，每个发酵池的不同发酵时间点各取一管，取样后立即装入采样管中封口，用液氮冷冻后带回实验室置于-80 ℃储存，以备后续实验。

1.2 仪器与设备

ABI ProFlex型PCR仪、QuantStudio 5型荧光定量PCR仪、iBright TM FL1500型凝胶成像系统、SL 8型台式离心机、NaneDrop 2000型超微量分光光度计：美国Thermo公司；DYC-Mini1型电泳槽：美国Bio Rad公司；Vortex1型旋涡混匀振荡器：上海沪析实业有限公司；Vetiri梯度基因扩增仪：美国AB公司；Illumina Miseq PE250型测序仪：美国Illumina公司。

1.3 实验方法

1.3.1 pH和总酸的测定

榨菜发酵液中总酸含量采用国标GB 5009.239—2016《食品酸度的测定》的方法测定[13]；使用pH计直接测定榨菜发酵液pH值。

1.3.2 DNA提取及Illumina Miseq高通量测序

1.3.3 测序数据分析

基于Illumina MiseqPE250测序平台，使用FLASH软件[14]通过Overlap对原始数据进行拼接，使用QIIME软件[15]过滤低质量的序列，利用UCHIME软件[14]鉴定并去除嵌合体，得到有效序列。根据UCLUST软件[16]在相似性97%的水平上进行可操作分类单元（operationaltaxonomic unit，OTU）聚类，得到代表序列。使用核糖体数据库项目（ribosomal database project，RDP）classifier[17]软件进行分类学分析，随后使用Silva 132数据库[18]对代表序列比对注释，生成不同分类水平上的物种丰度表，并删除不能注释到门水平的OTU。基于OTU表计算各样品的α-多样性指数与相对丰度，采用加权（Weighted）UniFrac距离算法结合主坐标分析（principal coordinates analysis，PCoA）与层级聚类分析均通过R软件实现。根据代表序列结果，利用PICRUSt软件[19]对榨菜发酵过程中微生物的基因功能进行预测。

2 结果与分析

2.1 发酵过程中pH值及总酸的变化

由图1可知，遂宁榨菜发酵起始pH值为6.5，随着发酵的不断进行，其pH值在发酵1～45 d内迅速下降至pH 3.83；此后，pH值随发酵时间的不断延长而缓慢下降，于发酵60 d稳定在pH 3.62。总酸含量变化趋势则与pH值变化趋势相反，起始总酸含量为0.37 g/L，随着发酵的进行而逐渐上升，于发酵60 d达到最大，其总酸含量为9.43 g/L。

图1 遂宁榨菜发酵过程中pH值和总酸含量的变化Fig.1 Changes of pH value and total acid contents during fermentation process of Suining Zhaicai

2.2 测序数据统计与α-多样性分析

由表1可知，不同发酵时间（10个）的样品测序共获得495 681条有效序列，其覆盖率为98.87%～99.42%，表明本次实验所测得的数据可真实反映样品的细菌群落组成。对各样品测得的有效序列按样品最低序列数进行抽平处理后，可得到318个OTU。Alpha-多样性结果显示，在整个发酵过程中，样品中细菌类群的香农（Shannon）指数和辛普森（Simpson）指数均呈先升后降的变化趋势，且都在发酵15 d时达到最大，分别为3.11和0.91。发现这种现象与遂宁榨菜发酵环境有关，发酵环境的不断酸化是影响原核微生物群落α-多样性的主要因素[10]。

表1 不同发酵时间的样品中检测到的有效序列数、OTU数及α-多样性指数Table 1 Number of valid reads，OTUs and alpha diversity indexes observed in samples with different fermentation time

不同发酵时间的样品中检测到OTU数的韦恩图与柱形图见图2。

图2 不同发酵时间的样品中检测到OTU数的韦恩图（A）与柱形图（B）Fig.2 Venn diagram (A) and bar chart (B) based on OTUs number of samples with different fermentation time

由图2A可知，每个发酵时间点都有特有的OTU，发酵1 d、15 d、30 d、45 d、60 d样品中特有的OTU数分别为36、54、15、21、42个。由图2B可知，不同发酵时间的样品检出的OTU数量各有差异，发酵1 d、15 d、30 d、45 d、60 d样品中检出的OTU数依次为138、147、128、122、104个，在整个发酵阶段一直存在的OTU数为34个，说明这些OTU生长范围广泛，可适应逐渐酸化的发酵环境。上述结果表明，不同发酵时间样品的细菌类群数量存在一定的差异。

2.3 β-多样性分析

基于加权UniFrac距离对所有样品进行主坐标分析（PCoA），结果见图3。由图3可知，PCoA的第一主坐标和第二主坐标分别占总变异的63.72%和25.11%，不同发酵时间样品中细菌菌群结构的组间差异显著（不同发酵时间的样品在其投影平面上的距离较远，大都分散于不同的象限，只有发酵45 d和60 d样品位于同一象限，表明发酵45 d和60 d样品的微生物群落结构较为相似，而与其他发酵时间的样品差异较大，且发酵过程中微生物群落呈抛物线状变化。结果表明，不同发酵时间的样品中细菌类群呈现出明显的时间异质性。LIU Z等[20]研究结果表明，发酵时间会显著影响泡菜发酵过程中微生物群落结构。此外，同一发酵时间样品中细菌群落结构的差异较小（投影平面上的间距较近，且聚集在一起），说明遂宁榨菜发酵具备很强的批次稳定性。

图3 各样品基于加权UniFrac距离的主坐标分析Fig.3 Principal coordinates analysis of each sample based on weighted UniFrac distance

2.4 细菌群落的组成与结构

由图4可知，在门水平上，发酵1 d样品中的优势菌门（平均相对丰度＞1%）是变形菌门（Proteobacteria，76.75%）、厚壁菌门（Firmicutes，16.67%）、类杆菌门（Bacteroidota，4.79%）和弯曲杆菌门（Campilobacterota，1.32%）；而其他发酵时间的样品中的优势菌门却只有变形菌门（Proteobacteria）和厚壁菌门（Firmicutes）。厚壁菌门（Firmicutes）的平均相对丰度随发酵时间的不断延长而持续增长，从16.67%增长至89.06%，而变形菌门（Proteobacteria）的平均相对丰度则与厚壁菌门（Firmicutes）的变化相反，其从76.75%降低至9.26%。这与LIANG H等[10]研究结果相一致。上述结果表明，变形菌门和厚壁菌门主导着遂宁榨菜的发酵，且厚壁菌门（Firmicutes）与变形菌门（Proteobacteria）的比值越大越趋近遂宁榨菜发酵的成熟。

图4 门水平上各样品中细菌类群的相对丰度和结构Fig.4 Relative abundance and structure of bacterial taxa observed from different samples at phylum level

由图5可知，在属水平上，不同发酵时间样品中检测到优势菌属（平均相对丰度＞1%）的结构也具有明显差异。发酵1 d样品以假单胞菌属（Pseudomonas）（24.79%）、Cobetia（24.41%）、弧菌属（Vibrio）（18.28%）、魏斯氏菌属（Weissella）（8.07%）、鞘氨醇杆菌属（Sphingobacterium）（4.63%）、不动杆菌属（Acinetobacter）（3.55%）、明串珠菌属（Leuconostoc）（3.54%）、柠檬酸杆菌属（Citrobacter）（1.95%）、乳球菌属（Lactococcus）（1.64%）、芽孢杆菌属（Bacillus）（1.55%）、Malaciobacter（1.32%）、肠杆菌属（Enterobacter）（1.17%）为优势菌属；发酵15 d样品以明串珠菌属（Leuconostoc）（27.20%）占据首要优势地位，其次是假单胞菌属（17.05%）、乳杆菌属（Lactobacillus）（15.54%）、乳球菌属（Lactococcus）（11.79%）、弧菌属（Vibrio）（18.28%）、不动杆菌属（Acinetobacter）（4.38%）、芽孢杆菌属（Bacillus）（2.44%）、嗜冷杆菌属（Psychrobacter）（2.34%）、盐弧菌属（Salinivibrio）（1.50%）、魏斯氏菌属（Weissella）（1.19%）、盐单胞菌属（Halomonas）（1.16%）、Idiomarina（1.14%）；其他发酵时间样品（30 d，45 d和60 d）则以乳杆菌属（Lactobacillus）（33.91%，64.71%和78.63%）占据首要优势地位，其次乳球菌属（Lactococcus）（32.27%，5.96%和1.19%）、弧菌属（Vibrio）（7.78%，5.17%和1.57%）、乳球菌属（Lactococcus）（6.43%，2.54%和1.16%）、不动杆菌属（Acinetobacter）（2.43%，2.09%和1.29%）、四联球菌属（Tetragenococcus）（1.80%，3.08%和2.95%）、芽孢杆菌属（Bacillus）（1.56%，2.01%和1.18%）、盐单胞菌属（1.51%，2.51%和3.16%）、芽孢杆菌属（Bacillus）（1.56%，2.01%和1.18%）和片球菌属（Pediococcus）（1.22%，1.93%和2.24%）。

图5 属水平上各样品中细菌类群的相对丰度和结构Fig.5 Relative abundance and structure of bacterial taxa observed from different samples at genus level

上述结果表明乳酸菌主导着遂宁榨菜发酵的成熟，且不同种乳酸菌的交替变化共同促进榨菜发酵的进行。乳酸菌通常来自发酵原料和自然环境，其可将原料中的小分子糖类物质转化为酸类化合物，从而赋予榨菜柔和的酸感；部分糖类物质被转化为醇类化合物，这些醇类物质和酸类代谢物反应进一步形成酯类化合物，进而赋予产品具有酯香和醇香等特征[21]。乳球菌属和明串珠菌属的平均相对丰度在整个发酵过程中均呈先增加后降低的变化趋势，而乳杆菌属则持续升高。该结果与张锐等[22]的研究结果相一致，说明明串珠菌属和乳球菌属可为乳杆菌属的繁殖提供合适的生长环境。在低盐的作用下，榨菜表面渗出的汁液非常适合明串珠菌属的生长繁殖，其产生的二氧化碳和酸使pH值迅速下降，从而阻止其他有害微生物的生长繁殖[22]。同时，明串珠菌属可将多余的糖转化为甘露醇和葡聚糖，而这两种物质非常适宜其他乳酸菌按一定的顺序生长，从而起到发酵启动剂的作用[22]。而假单胞菌属在发酵1 d样品中的平均相对丰度可达24.79%，正是由于新鲜蔬菜和发酵器具都可检测到大量的假单胞菌属[23]，随后假单胞菌属的平均相对丰度逐渐降低，直至发酵60 d不足0.01%。此外，随着发酵的不断进行，盐单胞菌属和片球菌属的平均相对丰度均呈逐渐增加的趋势，推测其与生产榨菜“盐渍”工序有关[24]。盐单胞菌属和片球菌属具有耐盐的特性，能分解葡萄糖，产生多种风味物质。而四联球菌属的平均相对丰度则呈先增加后保持稳定的变化趋势，四联球菌属广泛存在与含盐的腌制类发酵食品中，可参与氨基酸的合成及生成醛、醇、酮和酯等挥发性风味物质，从而提升发酵食品的风味和口感[25]。

2.5 细菌群落功能预测

PICRUSt软件分析结果见图6。由图6A可知，遂宁榨菜发酵过程中细菌群落的预测功能可分为4类（平均相对丰度＞1%），其中代谢（metabolism）（43.37%）、环境信息处理（environmental information processing）（18.06%）、遗传信息处理（genetic information processing）（16.06%）和细胞过程（cellularprocesses）（2.39%），说明遂宁榨菜发酵过程中大多数细菌都参与了原料的代谢，也正是多种微生物的生长代谢促进了原料的质地及营养的变化。代谢功能包括碳水化合物代谢（carbohydrate metabolism）（10.62%）、氨基酸代谢（amino acid metabolism）（7.75%）、能量代谢（energy metabolism）（5.65%）、辅助因子和维生素代谢（metabolism of cofactors and vitamins）（3.84%）、核苷酸代谢（nucleotide metabolism）（3.19%）、脂质代谢（lipid metabolism）（2.77%）、外源生物降解和代谢（xenobiotics biodegradation and metabolism）（2.16%）、多聚糖的生物合成和代谢（glycanbiosynthesisandmetabolism）（2.15%）、酶家族（enzyme families）（2.09%）、其他氨基酸代谢（metabolism of other amino acids）（1.62%）和萜类化合物和聚酮化合物代谢（metabolism of terpenoids and polyketides）（1.49%），其中碳水化合物代谢、能量代谢和氨基酸代谢占总代谢功能的55.38%，说明遂宁榨菜发酵过程中细菌群落积极参与了基本的代谢过程，尤其是乳酸菌，其可将原料中糖类物质及蛋白质生成酸类化合物和氨基酸，而发酵环境的不断酸化，使其需要更多的能量来抵御外界环境，确保其能维持正常的生长代谢。此外，环境信息处理功能包括膜转运（membranetransport）（15.70%）和信号转导（signal transduction）（2.36%）；遗传信息处理功能包括复制与修复（replication and repair）（6.91%）、翻译（translation）（3.95%）、转录（transcription）（2.86%）和折叠分选与降解（folding，sorting and degradation）（2.33%）；细胞过程功能只涉及到细胞运动（cell motility）（2.39%）这一条功能通路。同时，不同发酵时间细菌群落结构也会对这些通路的变化造成显著影响。

图6 榨菜发酵过程中细菌群落所预测到KEGG一级通路（A）和KEGG二级通路（B）平均相对丰度和组成结构的变化Fig.6 Changes of average relative abundance and composition in KEGG pathways (level 1) (A) and KEGG pathways (level 2)(B) observed from predicted functional categories of bacterial community during the fermentation process of Zhacai

由图6B可知，碳水化合物代谢和能量代谢等两条二级通路主要富集在发酵1 d样品中，表明该发酵时间的细菌群落主要是为发酵的进行提供能量和降解原料中的大分子物质；酶家族和多聚糖的生物合成和代谢等两条二级通路要富集在发酵15 d 和30 d样品中，揭示这个发酵时期的细菌群落为乳杆菌属的快速生长繁殖提供感应信号及相应底物；氨基酸代谢、脂质代谢、其他氨基酸代谢、萜类化合物和聚酮化合物代谢、核苷酸代谢、脂质代谢、外源生物降解和代谢和辅助因子和维生素代谢等8条二级通路则主要富集在发酵45 d和60 d样品中，表明乳杆菌属的大量繁殖可丰富榨菜产品中的氨基酸和脂质[26]。

3 结论

本研究利用Illumina Miseq高通量测序技术对遂宁榨菜发酵过程中的细菌多样性进行了解析，并结合PICRUSt软件预测菌群功能的变化。研究发现不同发酵时间（10个）样品测序共获得495 681条有效序列，抽平处理后可聚类为318个OTU。样品的微生物群落结构和多样性均存在明显差异，微生物细菌类群的Shannon指数和Simpson指数均呈先增后降的变化趋势，且都在发酵15 d时达到最大，分别为3.11、0.91。门水平上，发酵过程中的优势菌门为变形菌门（Proteobacteria）和厚壁菌门（Firmicutes），其两者相对丰度之和＞93%。属水平上，乳杆菌属（Lactobacillus）成为发酵后期的绝对优势菌属，其相对丰度从0.77%增加至78.63%。此外，PICRUSt预测结果显示，遂宁榨菜发酵过程中细菌群落的基因功能主要与新陈代谢类功能有关，包括氨基酸代谢、碳水化合物代谢、能量代谢等。本研究提高了对遂宁榨菜发酵过程中微生物演替模式的了解，为筛选潜在价值的乳酸菌提供数据参考。