泡菜作为传统发酵蔬菜制品的代表，因其味道酸爽、质地脆生、营养丰富等特点受到广大消费者的青睐[1-2]。目前，泡菜的生产模式逐渐从传统家庭作坊转变为工业化生产，工业泡菜的发酵工艺根据蔬菜的含水量分为两类，一类是含水量较高的蔬菜(如萝卜和青菜)直接用高盐腌制[3]，另一类是用泡菜母水泡制豇豆等含水量较低的蔬菜。与其他发酵蔬菜类似，工业豇豆泡菜也是一个多菌混合发酵体系，复杂的微生物群落对于豇豆泡菜独特的口感、风味和质地起到了极其重要的作用。
近年来，泡菜中微生物群落结构在国内外得到广泛关注，LIANG等[4]利用变性梯度凝胶电泳(polymerase chain reaction-denatured gradient gel electrophoersis，PCR-DGGE)技术探究了工业青菜发酵过程中微生物，结果表明乳杆菌属(Lactobacillus)、假单胞菌属(Pseudomonas)、弧菌属(Vibrio) 和盐单胞菌属(Halomonas)是优势细菌属和德巴利氏酵母属 (Debaryomyces) 是优势真菌属；朱琳等[5]通过基于16S rDNA的高通量测序在萝卜泡菜中检测到肠膜明串珠菌属、乳球菌属、魏斯氏菌属和片球菌属；LEE等[6]运用宏基因组学技术发现明串珠菌属(Leuconostoc)、乳杆菌属(Lactobacillus)和魏斯氏菌属(Weissella)是韩国泡菜发酵过程中的优势菌属。但这些方法仅停留在DNA层面上，已经死亡的微生物和处于非活跃状态的微生物也包括在内，并不能真实反映发酵过程中真正发挥作用的微生物群落组成情况。
宏转录组学技术是以样品中全部微生物RNA信息为分析对象，从转录组水平分析样本微生物的组成及活性基因表达情况，其结果反映的是活性微生物的群落结构[7-8]。近年来，宏转录组学技术已用于传统酸面团[7]和奶酪[9]的研究，本研究以工业豇豆泡菜为研究对象，运用宏转录组学技术对发酵过程中活性微生物群落结构进行探究，以确定其优势活性微生物组成，同时对发酵过程中功能基因的表达情况进行了探究，以期为对工业泡菜发酵机理解析及发酵过程控制提供理论基础。
理化试剂：NaOH、亚铁氰化钾、乙酸锌、四硼酸钠、对氨基苯磺酸、NaNO2、盐酸萘乙二胺(均为分析纯)，成都市科龙化工试剂厂。
测序试剂：RNA PowerSoil® Total RNA Isolation Kit 提取试剂盒，美国 Mobio公司；Ribo-Zero rRNA Removal Kit，美国EPICENTRE公司； TruSeq Stranded mRNA LT Sample Prep Kit，美国Illumina公司； Agilent High Sensitivity DNA Kit，美国安捷伦公司。
PHS-3C酸度计，上海仪电科学仪器股份有限公司；WS202盐度计，上海民仪电子有限公司；LC2030高效液相色谱仪，日本岛津；琼脂糖凝胶电泳仪，中国北京市六一仪器厂；高速冷冻离心机，美国赛默飞公司；Agilent 2100 Bioanalyzer，美国安捷伦公司； PCR仪，美国Bio-Rad公司；Illumina Hiseq 10X测序平台 美国Illumina 公司。
1.3.1 样品采集与处理
分别采集发酵1、15、30、60、90 d的豇豆泡菜液(四川省眉山市某泡菜工厂)，取样后冰袋运回实验室，抽滤过膜富集微生物置于无菌无酶离心管中，液氮速封后转移至-80 ℃冰箱备用。该工厂豇豆泡菜发酵工艺大致为将含盐量约10%(质量分数)的泡菜母水与新鲜豇豆按1∶1的质量比投入发酵池，木棍压实，在相对开放的环境中发酵3个月。
1.3.2 泡菜理化指标测定
pH、盐度分别使用pH计和盐度计进行测定；总酸参照GB/T 12456—2008食品总酸的测定[10]中酸碱滴定法测定，结果以乳酸计；亚硝酸盐含量采用 GB5009.33—2016食品安全国家标准食品中亚硝酸盐与硝酸盐的测定[11]测定。采用高效液相色谱法测定乳酸，乙酸和葡萄糖含量。测定条件为色谱柱Carbomix H-NP，7.8 mm×300 mm，5 μm；流动相0.25 mol/L H2SO4，流速0.6 mL/min, 柱温55 ℃, 注射量10 μL。紫外检测器波长210 nm检测乳酸和乙酸, 示差检测器检测葡萄糖，通过与标准品的保留时间和峰面积的比较，对其进行鉴定和定量。以上实验均做3次平行测定。
1.3.3 泡菜液样品RNA提取及rRNA的去除
使用RNA PowerSoil® Total RNA Isolation Kit试剂盒提取泡菜液样品的总RNA, 用1.5%(质量分数)的琼脂糖凝胶电泳检测RNA的完整性，然后用紫外分光光度计定量。将rRNA 序列特异性的探针与总RNA杂交，用磁珠去除rRNA/探针复合物，再用乙醇沉淀法纯化mRNA。
1.3.4 cDNA文库构建及测序
采用TruSeq Stranded mRNA LT Sample Prep Kit试剂盒构建cDNA文库，用Agilent 2 100 Bioanalyzer对构建的文库质检，再用Promega QuantiFluor对文库定量。将需要上机的文库梯度稀释，并按所需数据量比例混样，混合后的文库变性解螺旋成单链。基于 Illumina Hiseq 10X测序平台，采用全转录组鸟枪法策略，对上述构建的测序文库进行2×150 bp的双端测序，此步骤委托上海派森诺生物公司进行。
1.3.5 测序数据分析
原始数据经过质控，剔除rRNA序列后获得高质量reads, 使用Trinity对所有测序读段进行从头组装[12]，采用CD-HIT对各样本组装拼接后得到的转录本以95%的相似度和90%的覆盖率进行归并去冗余， 并以最长的序列作为该 Unigene 的代表序列。将每个样本的Unigene序列与NCBI-NT数据库进行BLASTN比对(设置期望值E-value<0.001)，利用MEGAN软件采取最低共同祖先算法采用获得物种注释[13]，结合 Unigene 序列在各样本中的表达量数据，统计物种在各个样本中的表达丰度，获得在门和属分类学水平上活性菌群结构组成分布。将Unigene 序列集上传至 KAAS，与KEGG代谢通路数据库比对，从而获得各 Unigene 的功能注释和分类信息[14]。
工业豇豆泡菜发酵过程中理化指标变化如图1所示。发酵第1天pH值为4.78，发酵15 d后迅速下降至4.0以下，随后维持在3.50左右。随着发酵进行，总酸含量明显上升，发酵30 d达到最高10.01 g/kg, 之后略有下降，到发酵结束时总酸含量为8.32 g/kg。与总酸质量浓度变化趋势相似，发酵前30 d，乳酸质量浓度显著增加至9.02 g/L，随后逐渐下降。另外，整体上，乙酸含量逐渐增加而葡萄糖质量浓度从3.76 g/L逐渐降低至0.22 g/L。发酵前期，泡菜母水及豇豆表面存在的乳酸菌等微生物利用新鲜豇豆中的葡萄糖等糖类产生乳酸和乙酸等有机酸使得总酸含量快速增加，pH下降，发酵后期，由于乳酸积累过多，抑制了乳酸菌的生长代谢，同时微生物为了克服乳酸胁迫而消耗乳酸，导致乳酸及总酸含量下降。发酵过程中的盐度从11%逐渐上升直至发酵结束达到13.1%。这可能是该工厂在泡制豇豆时未进行完全封池，随着发酵进行，发酵池中的水分蒸发，盐度上升。亚硝酸盐含量随着发酵的进行逐渐下降，亚硝酸盐主要来自豇豆原料和杂菌代谢产生，随着乳酸菌的生长代谢产生乳酸、 乙酸和酒石酸等有机酸提供酸性发酵环境不仅抑制杂菌生长还能使亚硝酸盐发生氧化还原反应，同时部分微生物产生亚硝酸还原酶和细菌素也对亚硝酸盐的降解起着重要作用[15-16]。本次取样未能检测出亚硝峰，推测可能在1～15 d出现亚硝峰。发酵过程中的亚硝酸盐含量都远低于GB 2762—2017 食品中污染物限量腌渍蔬菜亚硝酸盐限量规定的20 mg/kg[17]。
图1 工业豇豆泡菜发酵过程中理化特性变化
Fig.1 The changes of physicochemical properties of industrial cowpea paocai during the fermentation
测序后的原始数据经过质控，得到高质量序列后由Trinity 软件进行拼接组装，相应统计结果如表1所示。5个豇豆泡菜样本分别获得 183 858、90 693、96 263、48 246和19 955条序列且拼接效果较好。为了比较5个样本的多样性，使用 QIIME 软件对其计算α多样性指数，结果如表2所示。Chao1指数和ACE指数衡量物种丰度，Shannon 指数和Simpson 指数衡量物种多样性和均匀度[18]。由表2可知，发酵第1天的Chao1和ACE指数最大，表明发酵第1天的活性菌群丰富度最高，Shannon指数值和Simpson指数值也比较高，表明发酵第1天活性菌群多样性较高，而发酵90 d所有α多样性指数都最低，表明发酵后期活性菌群多样性和丰富度较低，同时也说明在发酵过程中活性微生物群落演替逐渐趋于稳定。
表1 序列组装结果统计表
Table 1 Statistics of sequence assembly results
发酵时间/d序列总数序列长度/bpN50/bpN90/bp最短长度/bp最长长度/bp1183858629023023452262017139159069336278440425239201829130962632929785429621920133076048246217402055032492019964901995571341193652262013155
注：N50/N90长度越长，说明组装拼接效果越好，得到的宏转录组序列越完整；其中，N50长度尤为重要，是衡量宏转录组组装拼接效果好坏的主要指标之一
表2 α多样性指数分析表
Table 2 Analysis of α diversity index
发酵时间/dSimpsonShannonChao1ACE10.783.331645.451646.511150.803.00827.11825.14300.733.421128.671127.24600.512.49833.00829.46900.401.54458.17460.20
工业豇豆泡菜发酵过程中活性微生物在门水平上群落结构演替如图2所示。发酵第1天，厚壁菌门(Firmicutes)的相对丰度高达88%，此外还检出少量变形菌门(Proteobacteria, 7.9%)。发酵第15天，厚壁菌门的相对丰度下降而子囊菌门(Ascomycota)的相对丰度上升，分别为57.8%和41.7%。发酵30 d后，厚壁菌门以80%的相对丰度再次成为绝对优势菌门。随着发酵进行，厚壁菌门的相对丰度又逐渐降低而子囊菌门的相对丰度显著增加，最终，子囊菌门分别以77%和98% 的相对丰度主导了发酵60 d和发酵90 d。随着发酵时间的延长，原料中的营养物质减少，代谢产物增加，微生物相互竞争，厚壁菌门中部分菌属对泡菜液体系不适应逐渐被淘汰，而子囊菌门中耐酸耐盐的真菌开始在发酵后期发挥作用。
图2 工业豇豆泡菜发酵过程中门水平上活性微生物群落结构
Fig.2 Structure of the active microbial community on phylum level during the industrial cowpea paocai fermentation
从属水平来看(图3)，发酵第1天相对丰度最高的是乳杆菌属(Lactobacillus， 61.8%)，其次是厚壁菌门的片球菌属(Pediococcus)，其相对丰度为19.7%，此外还检测到变形菌门的弧菌属(Vibrio)和厚壁菌门的魏斯氏菌属(Weissella)，其相对丰度分别为7%和4.65%。乳杆菌属(Lactobacillus)、片球菌属(Pediococcus)和魏斯氏菌属(Weissella)主要来自泡菜母水和新鲜豇豆，通过同型或异型发酵产生乳酸和乙酸，使得发酵前期总酸含量迅速上升，对豇豆泡菜酸爽口感的形成有重要作用，这些菌属在腌菜[19]、家庭泡菜[20]和东北酸菜[21]等发酵蔬菜中也都有发现。发酵第15天时，乳杆菌属(Lactobacillus)的相对丰度下降至55.6%,而子囊菌门的Millerozyma的相对丰度上升至33%，Millerozyma是一种罕见的酵母病原体[22]，但也有研究报道Millerozyma的某些菌种具有降解生物胺的作用[23]，Millerozyma在豇豆泡菜发酵过程中的作用有待进一步分析。另外，还发现了少量子囊菌门的Starmerella(4.34%)，而发酵初期相对丰度较高的片球菌属(Pediococcus)、弧菌属(Vibrio)和魏斯氏菌属(Weissella)的相对丰度均低于0.05%，说明这些菌属是工业豇豆泡菜发酵的启动微生物，参与发酵初始阶段，对发酵中后期作用较小，这与之前佟婷婷等[24]发现魏斯氏菌属是泡菜母水作引子发酵蔬菜过程中的启动菌结果保持一致。随着发酵进行，发酵30 d时乳杆菌属(Lactobacillus)再次成为绝对优势菌属，其相对丰度为77.6%，另外，Millerozyma的相对丰度快速下降至7%。发酵第60天，Starmerella的相对丰度迅速增加到69.3%，成为主导菌属，同时乳杆菌属(Lactobacillus)的相对丰度明显下降，但仍以16.6%的相对丰度成为第二优势菌属，在发酵前中期乳杆菌属生长代谢活跃产生大量乳酸，由于乳酸积累过多且发酵后期盐度上升，乳杆菌生长代谢受到抑制。到发酵90 d时，Starmerella以76.6%的相对丰度成为绝对优势菌属，其次是子囊菌门的Sugiyamaella (7.1%)和Millerozyma (6.05%)，仅有1.14%的乳杆菌属(Lactobacillus)在发酵末期被检测出。发酵后期，Starmerella逐渐占据主导地位，有关Starmerella报道较少，其中Starmerella bacillaris产甘油和挥发性酸的能力较强，被用于发酵葡萄酒[25-26]。与本研究结果不同，李恒等[27]发现Kazachstania和假丝酵母属(Candida)是泡菜母水中的主要真菌属，可能是检测方法、泡菜品种、发酵工艺、取样时间及取样地点不同导致。
图3 工业豇豆泡菜发酵过程中属水平上活性微生物群落结构
Fig.3 Structure of the active microbial community on genus level during the industrial cowpea paocai fermentation
各样本在KEGG第一等级上的表达量分布如图4所示。发酵第1天，有关代谢(metabolism) 的基因丰度最高，占比 48.7%，其次是遗传信息处理(genetic information processing) 和环境信息处理 (environmental information processing)，分别占17.3%和11.4%。在发酵15～30 d的转录本中注释到有关代谢的基因最多，分别占 41%和48%，其次是人类疾病(human diseases)，分别占23.9%和13%。随着发酵进行，代谢相关基因的表达量逐渐下降，发酵60 d时，代谢相关基因的相对表达量为32.9%，仍维持较高的表达量。但到发酵结束时，人类疾病相关基因的相对表达量达46%，生物体系统(organismal systems)次之，占比31%，而有关代谢的基因表达量仅占19%。在豇豆泡菜发酵前中期，以乳杆菌属为主的微生物利用新鲜豇豆中的营养物质迅速生长并通过不同代谢途径生成各种影响豇豆泡菜品质的代谢产物，使得此发酵期间代谢相关基因的表达量较高。在发酵后期，人类疾病相关基因的表达量增加，在此阶段真菌丰度也迅速增加，由此推测部分真菌可能与人类疾病密切相关。
图4 工业豇豆泡菜发酵过程中基因表达的 KEGG第一功能等级富集分析
Fig.4 KEGG analysis of gene expression on first level during the industrial cowpea paocai fermentation
对代谢子功能进行具体分析，如图5所示，共注释到11种代谢功能。发酵第1天，有关碳水化合物代谢(carbohydrate metabolism)和能量代谢(energy metabolism)的基因表达量较高，分别为17.9%和16.9%，此外还注释到少量的核酸代谢(nucleotide metabolism)和氨基酸代谢(amino acid metabolism)。发酵第15天，碳水化合物代谢和能量代谢相关基因的相对表达量分别为18.9%和11.8%，另外还注释到3%的氨基酸代谢。发酵第30天，有关碳水化合物代谢和能量代谢的基因相对表达量分别下降至18.2%和7.5%，而氨基酸代谢相关基因的表达量上升至7.4%。发酵60 d，碳水化合物代谢相关基因表达量继续下降至11.1%，能量代谢相关基因表达量开始增加，占9.5%。到发酵90 d时，能量代谢的相关基因表达量占15.9%，而碳水化合物代谢仅注释到1.2%。除发酵第90天外，碳水化合物代谢相关基因表达量在整个发酵过程中都是最高的。发酵第1天，新鲜豇豆含大量以糖类为主的碳水化合物，乳杆菌属，片球菌属和魏斯氏菌属等乳酸菌迅速生长，并通过糖酵解等代谢途径利用葡萄糖产生丙酮酸等代谢产物，使得碳水化合物代谢相关基因表达量较高。发酵15～30 d时，乳杆菌属等微生物继续活跃参与碳水化合物代谢，并进一步将发酵初期积累的丙酮酸及其他中间代谢产物生成乳酸和乙酸等风味成分，豇豆泡菜中乳酸和乙酸含量在此发酵阶段迅速增加，葡萄糖含量显著下降，大量乳酸和乙酸的积累使得发酵前30 d的总酸含量迅速增加，pH快速下降。发酵60 d后，随着葡萄糖等糖类基本耗尽且乳杆菌属丰度明显降低，碳水化合物代谢相关基因表达量逐渐降低，尤其是乳杆菌属含量很低的发酵90 d, 说明豇豆泡菜发酵过程中乳杆菌属是参与碳水化合物代谢的主要菌属。发酵前30 d还有少量氨基酸代谢，这可能是因为微生物将新鲜豇豆中的蛋白质分解成氨基酸，部分氨基酸可能在微生物氨基酸脱羧酶的作用下生成生物胺。此外，在整个发酵过程中能量代谢相关基因表达量较为稳定，这是因为微生物需要维持正常生理活动，能量代谢是细胞进行其他生理活动的基础。
图5 工业豇豆泡菜发酵过程中代谢途径的基因表达富集分析
Fig.5 KEGG analysis of gene expression of metabolism during the industrial cowpea paocai fermentation
本研究利用宏转录组学技术对工业豇豆泡菜活性微生物群落结构进行了解析，结果表明乳杆菌属(Lactobacillus)和片球菌属 (Pediococcus) 是发酵初期主要的活性菌属，发酵中期(15～30 d)的优势活性菌属是乳杆菌属(Lactobacillus)，随着发酵进行，细菌丰度逐渐下降，以Starmerella为主的真菌丰度逐渐增加，最终主导了发酵末期(60～90 d)。同时对其进行了KEGG功能注释，工业豇豆泡菜液中活性微生物主要参与代谢功能，集中在碳水化合物代谢和能量代谢，碳水化合物代谢相关基因表达量占比最高。本研究结果加深了对工业豇豆泡菜微生物群落组成和多样性的认识，发现了发酵过程中起关键作用的菌属，为解析泡菜发酵过程中微生物代谢机理提供参考，后续可进一步分析活菌群落在发酵过程中具体的代谢途径来揭示其对工业豇豆泡菜风味和品质的贡献，以期为发酵蔬菜工业化发展提供借鉴。
[1] 陈功.试论中国泡菜历史与发展[J].食品与发酵科技,2010,46(3):1-5.
[2] 张亚豪,梁会朋,常聪,等.基于实时荧光定量PCR技术监测四川工业泡菜发酵过程中主要细菌的变化[J].中国调味品,2018,43(1):35-38;47.
[3] LIANG Huipeng,CHEN Huiying,JI Chaofan,et al.Dynamic and functional characteristics of predominant species in industrial paocai as revealed by combined DGGE and metagenomic sequencing[J].Frontiers in Microbiology,2018,9:2 416.
[4] LIANG Huipeng,YIN Liguo,Zhang Yahao,et al.Dynamics and diversity of a microbial community during the fermentation of industrialized Qingcai paocai, a traditional Chinese fermented vegetable food, as assessed by Illumina MiSeq sequencing, DGGE and qPCR assay[J].Annals of Microbiology,2015,68(2):111-122.
[5] 朱琳,高凤,曾椿淋,等.萝卜泡菜细菌多样性的高通量测序分析[J].现代食品科技, 2018,34(2):225-231.
[6] LEE M,SONG J H,JUNG M Y,et al.Large-scale targeted metagenomics analysis of bacterial ecological changes in 88 kimchi samples during fermentation[J].Food microbiology,2017,66:173-183.
[7] 张国华,王伟,涂建,等.基于宏转录组学技术解析传统酸面团中微生物代谢机理[J].中国粮油学报,2019,34(11):10-16.
[8] 田露,闵建红,张帝,等.宏组学在发酵食品微生物群落研究中的应用进展[J].生物技术通报,2019,35(4):116-124.
[9] MONNET C,DUGAT-BONY E,SWENNEN D,et al.Investigation of the activity of the microorganisms in a reblochon-style cheese by metatranscriptomic analysis[J].Frontiers in Microbiology,2016,7(337):536.
[10] GB/T 12456—2008,食品中总酸的测定[S].北京:中国标准出版社,2008.
[11] GB 5009.33—2016,食品中亚硝酸盐与硝酸盐的测定[S].北京:中国标准出版社,2016.
[12] GRABHERR M G,HAAS B J,Y ASSOUR M,et al.Trinity: Reconstructing a full-length transcriptome without a genome from RNA-Seq data[J].Nature Biotechnology,2011,29(7):644-652.
[13] HUSON D H,MITRA S,RUSCHEWEYH H J,et al.Integrative analysis of environmental sequences using MEGAN4[J].Genome Research,2011,21(9):1 552-1 560.
[14] KANEHISA M,GOTO S,KAWASHIMA S,et al.The KEGG resource for deciphering the genome[J].Nucleic Acids Res,2004,32: 277-280.
[15] 徐柯, 成林林, 袁美，等.泡豇豆发酵过程中有机酸变化及对亚硝酸盐降解的影响[J].食品与发酵工业,2019,45(17):60-65;72.
[16] 张庆芳,迟乃玉,郑燕,等.乳酸菌降解亚硝酸盐机理的研究[J].食品与发酵工业,2002,28(8):27-31.
[17] GB 2762—2017 食品中污染物限量[S].北京:中国标准出版社,2017
[18] ZANG Jinhong, XU Yanshun, XIA Wenshui,et al. Dynamics and diversity of microbial community succession during fermentation of Suan yu, a Chinese traditional fermented fish, determined by high throughput sequencing[J].Food Res Int,2018,111:565-573.
[19] LIU Zhanggen,LI Junyi,HUANG Tao,et al.Comparison of the bacterial communities in home-made Nanfeng yancai with and without salt[J].Food Res Int,2019,125:108 509.
[20] XIAO Yangsheng,XIONG Tao,PENG Zhen,et al.Correlation between microbiota and flavours in fermentation of Chinese Sichuan Paocai[J].Food Res Int,2019,114:123-132.
[21] HE Zhen,CHEN Huiying,WANG Xinyi,et al.Effects of different temperatures on bacterial diversity and volatile flavor compounds during the fermentation of suancai, a traditional fermented vegetable food from northeastern China[J].LWT - Food Science and Technology,2020,118:108 773.
[22] HONG S I,SUH Y S,KIM H O,et al.Successful treatment of catheter related blood stream infection by millerozyma farinosa with micafungin: A case report[J].Infection & Chemotherapy,2018,50(4):362-366.
[23] 程淑敏,蓝翔,徐莹,等.生物胺降解菌在酱油中的应用[J].食品与发酵工业,2019,45(4):129-134.
[24] 佟婷婷,田丰伟,王刚,等.基于宏基因组分析四川泡菜母水作引子的泡菜发酵过程中细菌多样性变化[J].食品工业科技,2015,36(21):173-177.
[25] VASILEIOS E,KALLIOPI R,FABRIZIO T, et al.Exploitation of the non- Saccharomyces yeast Starmerella bacillaris (synonym Candida zemplinina) in wine fermentation: Physiological and molecular characterizations[J].International Journal of Food Microbiology,2015,199:33-40.
[26] LEMOS J W J F,NADAI C,CREPALDE L T,et al.Potential use of Starmerella bacillaris as fermentation starter for the production of low-alcohol beverages obtained from unripe grapes[J].International journal of food microbiology,2019,303:1-8.
[27] 李恒,陈功,伍亚龙,等.高通量测序方法研究传统四川泡菜母水中微生物群落的动态变化[J].食品科学,2018,39(24):131-138.