襄阳大头菜发酵过程中细菌的多样性

吴进菊，李宇昂，王梓杭，郑婷婷，于 博，余海忠

（湖北文理学院食品科学技术学院g化学工程学院，湖北 襄阳 441053）

湖北襄阳大头菜历史悠久，相传是三国时期诸葛亮所发明，因此又叫孔明菜、诸葛菜。襄阳大头菜在2007年获得“国家地理标志保护产品”，是我国四大名腌菜之一[1]。襄阳大头菜制作周期长，需经历长达数月至数年的发酵过程才能成熟。在漫长的发酵过程中，有大量的微生物参与其中，蕴藏了丰富的微生物资源。

近年来，很多研究人员把目光转向发酵食品的微生物多样性研究，取得了很大的成果，如Liang Huipeng等[2]研究发现榨菜发酵初期乳酸菌生长迅速，然后达到稳定的发酵水平，弧菌、盐单胞菌、明串珠菌、魏斯氏菌在发酵过程中数量减少，片球菌数量增加。而家庭制作的泡菜中耐酸乳杆菌和短乳杆菌为优势菌属[3]。青菜泡菜发酵过程中主要菌属为乳酸菌属、假单胞菌属、弧菌属和盐单胞菌属。大多数隶属于变形菌门或拟杆菌门的细菌菌属在发酵前期检测频率较高，而乳酸菌（隶属于厚壁菌门）在发酵后期占主导地位[4]。Zhou Qi等[5]研究发现东北酸菜发酵过程中优势微生物为乳杆菌属、明串珠菌属、肠杆菌属、热袍菌属、假单胞菌属等，随着发酵的进行，乳酸菌逐渐增多，而肠杆菌逐渐减少。赵慧君[6]、郭壮[7]等对襄阳大头菜腌制液膜醭中的微生物菌群结构进行了研究，结果表明，细菌中变形菌门和硬壁菌门为优势菌门，属水平上盐单胞菌、色盐杆菌属、海杆菌属、盐厌氧菌属、四联球菌为优势种属；真菌中子囊菌门平均含量高达99.99%，属水平上合囊菌属、毕赤酵母属、念珠菌属为优势菌属，平均相对含量分别为56.10%、40.96%和2.32%。另外，研究人员对发酵鱼类[8-9]、酒类[10-14]、发酵乳制品[15-20]、发酵豆制品[21-24]、发酵蔬菜[25-30]等发酵食品中微生物多样性方面也进行了大量的研究。

目前针对传统发酵大头菜发酵过程中微生物方面的研究较少，早期研究基本是建立在传统纯培养技术研究方法上，只能获得微生物总量中极少的一部分，存在很大的局限性和片面性，从而对发酵过程中微生物群落结构变化和种群多样性了解不够清楚，限制了大头菜标准化和工业化生产。本研究采用Illumina MiSeq高通量测序技术考察了大头菜发酵过程中细菌的多样性和菌群结构的演替规律，全面揭示了大头菜发酵过程中细菌群落结构的动态变化，为今后大头菜生产的改进提供了一定的研究基础。

1 材料与方法

1.1 材料与试剂

24 份大头菜发酵液样本 湖北孔明菜食品有限公司；FastPfu聚合酶 北京全式金生物技术有限公司；土壤基因组提取试剂盒 美国MP Biomedicals公司；DL2000 DNA Marker 宝生物工程（大连）有限公司。

1.2 仪器与设备

MiSeq高通量测序平台 美国Illumina公司；Quanti Fluor™-ST蓝色荧光定量系统 美国Promega公司；9700型聚合酶链式反应（polymerase chain reaction，PCR）仪 美国ABI公司；台式冷冻离心机 德国Eppendorf公司。

1.3 方法

1.3.1 样本采集

实验样本采自湖北孔明菜食品有限公司，大头菜采用“三腌、五卤、六晒”的传统工艺进行生产，按照加工工序的特点，在每一道腌制和卤制工序完成后取样，共采集了8 个不同发酵时期的大头菜发酵液样本，分别编号为1_1～1_8。每个发酵时期平行取3 个发酵液样本，编号分别为a、b和c。样本采集后快速运送至实验室，置于-80 ℃冰箱保存备用。

1.3.2 Illumina MiSeq测序和数据分析

将24 份大头菜发酵液样本送至上海美吉生物医药科技有限公司进行测序。以微生物总DNA为模板，以338F（5’-ACTCCTACGGGAGGCAGCAG-3’）和806R（5’-GGACTACHVGGGTWTCTAAT-3’）为上下游引物扩增细菌16S rRNA基因的V3-V4可变区序列。将PCR产物纯化后进行荧光定量，进一步构建MiSeq文库后进行双端测序，测序读长300 bp。数据分析采用美吉I-Sanger云平台进行在线分析。对优化序列提取非重复序列，按照97%相似性水平对非重复序列（不含单序列）进行可操作分类单元（operational taxonomic units，OTU）聚类，去除嵌合体，获得OTU代表序列。利用Mothur软件计算各样本的α多样性，表征细菌群落的丰富度和多样性。采用群落Bar图和Heatmap图反映样本在各分类学水平上的物种组成。主成分分析（principal component analysis，PCA），采用降维的方法对数据进行简化分析，利用物种丰度表，基于欧氏距离进行作图，样本物种组成越相似，反映在PCA图中的距离越近。根据β多样性距离矩阵进行层级聚类分析，使用非加权组平均法（unweighted pair-group method with arithmetic mean，UPGMA）算法构建树状结构，可视化呈现不同环境样本中微生物进化的差异程度，距离算法采用Bray-Curtis，样本层级聚类方式为Average。采用R语言（版本号3.2.5）vegan包进行相似性分析（analysis of similarities，ANOSIM）。ANOSIM又称为置换多因素方差分析或非参数多因素方差分析，是一种非参数检验，用来检验组间的差异是否显著大于组内差异，从而判断分组是否有意义。

1.3.3 样本序列的提交及序列号

将本研究中所有fastq序列文件提交到美国国立生物技术信息中心SRA数据库，序列号为PRJNA524770。

2 结果与分析

2.1 测序结果和α多样性分析

通过Illumina MiSeq高通量测序，24 份大头菜发酵液样本共获得1 158 752 条优化序列，平均长度为448 bp，各样本抽平多样性指数如表1所示。采用RDP classifier贝叶斯算法对97%相似水平的OTU代表序列进行生物信息统计分析，按最小样本序列数进行抽平分析（每个样本27 407 条有效序列），结果表明，所有样本序列归属于19 个门、32 个纲、67 个目、133 个科、293 个属、482 个种、658 个OTU，每个样本的OTU分类信息见表1。

表1 样本信息表和α多样性指数分析Table 1 Bacterial community composition and alpha diversity estimators of pickle juice samples

由表1可知，每个样本的Coverage值均为0.996以上，表明测序深度良好，覆盖率高，可以真实展示样本中的绝大多数细菌。通过Shannon指数和Simpson指数分析可知，在所有样本中，1_1的细菌多样性最高，而1_2细菌多样性相对较低。这从各样本门、目、属、OTU的个数也能看出，1_1在各分类学水平的数目均高于其他样本，而1_2低于其他样本。1_1、1_5和1_7中，Chao 1指数和ACE值较高，表明样本中群落的丰富度较高。总体来说，发酵前期和发酵后期样本的细菌多样性没有呈现出整体的差异性，在8 个不同的发酵时期，样本细菌多样性有的时期高，有的时期低，呈现不规律变化。

稀释曲线通常是在97%的相似水平下，从样本中随机抽取一定数量的个体，以测序深度为横坐标，以各个样本中的细菌多样性指数为纵坐标，在二维水平建立坐标系，将样本内细菌的多样性随着测序程度加深的变化情况清楚的展示出来。由图1可知，24 个样本随着抽取的reads条数的增加，曲线平坦，已接近直线，说明此次取样的样本测序数据量足够大，能比较全面地反映样本中绝大多数细菌的多样性信息。

图1 稀释曲线Fig. 1 Rarefaction curves

2.2 门水平细菌群落组成分析

24 个样本中共检测出19 个细菌门，包括变形菌门（Proteobacteria）、放线菌门（Actinobacteria）、拟杆菌门（Bacteroidetes）、厚壁菌门（Firmicutes）、螺旋菌门（Spirochaetae）、螺旋体菌门（Saccharibacteria）、绿弯菌门（Chloroflexi）、梭杆菌门（Fusobacteria）、软壁菌门（Tenericutes）、栖热菌门（Deinococcus-Thermus）、疣微菌门（Verrucomicrobia）、酸杆菌门（Acidobacteria）、热袍菌门（Thermotogae）、FBP、浮霉菌门（Planctomycetes）、芽单胞菌门（Gemmatimonadetes）、互养菌门（Synergistetes）等。其中变形菌门、厚壁菌门、放线菌门、拟杆菌门、螺旋菌门为优势菌门，其他14 个菌门相对丰度较低，在任一发酵时期中的相对丰度均在1%以下。不同发酵时期大头菜样本细菌群落组成如图1所示。在整个发酵过程中，变形菌门和厚壁菌门占据了绝对的优势，两者相对丰度之和在1_1中达到75.0%（平均值，下同），在其他发酵时期均达到了90%以上，这与大头菜腌制液膜醭中的微生物菌群组成相似[6]。1_1中，变形菌门相对丰度最高，为67.7%，厚壁菌门、放线菌门、拟杆菌门三者相当，分别为7.32%、13.1%和10.4%。而在1_2中，厚壁菌门占据了绝对的优势，相对丰度达到了66.2%，变形菌门、放线菌门和拟杆菌门的相对丰度均明显下降。在后续的发酵过程中，变形菌门的相对丰度都高于厚壁菌门。

图2 门水平上细菌群落组成分析Fig. 2 Bacterial community composition at the phylum level

2.3 属水平细菌群落组成分析

图3 属水平上细菌群落组成Heatmap图Fig. 3 Heatmap of bacterial community composition at the genus level

Heatmap图根据物种和样本间丰度的相似性进行聚类，使高丰度和低丰度的物种分块聚集，通过颜色变化与相似程度反映不同样本在分类学水平上群落组成的相似性和差异性。对属水平总丰度排在前35的物种进行层级聚类，物种和样本层级聚类方式均采用平均聚类方式，得到的群落Heatmap如图3所示。24 个大头菜发酵液样本中共检测出293 个细菌属，其中有46 个菌属相对丰度在1%以上，相对丰度排在前5 位的菌属有盐厌氧菌属（Halanaerobium）、弧菌属（Vibrio）、盐单胞菌属（Halomonas）、乳杆菌属（Lactobacillus）和色盐杆菌属（Chromohalobacter）。在发酵的前3 个时期，弧菌属为优势菌属，相对丰度为27.6%～42.4%，而在发酵的后5 个时期，弧菌属相对丰度急剧下降，均在3%以下，而盐厌氧菌属相对丰度迅速增加，发酵前3 个时期相对丰度均在0.1%以下，而发酵后5 个时期相对丰度增加至31.1%～42.0%，占据了绝对优势。在1_4～1_8中，除盐厌氧菌属外，优势菌属还有盐单胞菌属、色盐杆菌属和norank_f_TBZ33，而这3 个菌属在1_1～1_3中相对丰度都非常低，甚至没有检测出。嗜冷杆菌属（Psychrobacter）在1_1和1_2中相对丰度较低，分别为0.61%和0.62%，在1_3～1_8中相对丰度有所增加，达到了1.2%～5.3%。

关于乳杆菌属，在1_1中相对丰度仅为1.4%，而在1_2中提高至55.6%，占据了绝对优势，随后在1_3中又下降至4.0%。在1_5～1_8中，乳杆菌属相对丰度下降至1%以下。乳酸菌是发酵蔬菜中常见的优势菌属[1-4]，在蔬菜发酵过程中起着非常关键的作用，可以抑制硝酸盐还原细菌的作用，从而降低亚硝酸盐的形成[3]。在发酵前期，加入的食盐量相对较少，所以不耐盐的弧菌属、乳杆菌属等为优势菌属。而随着发酵的进行，大量的食盐不断加入，使发酵液中的食盐浓度大大提高，在发酵结束后发酵液中的食盐浓度甚至达到饱和状态，因此，发酵后期不耐盐的菌属相对丰度急剧下降，而耐盐菌属发展成为优势菌属，如盐厌氧菌属、盐单胞菌属、色盐杆菌属等。

由图3可知，根据物种间丰度的相似性，总丰度排在前35的细菌菌属可聚为4 类，盐厌氧菌属、norank_f_TBZ33、盐单胞菌属和色盐杆菌属归为I类，在发酵前期相对丰度非常低，而在发酵后期相对丰度较高，成为优势菌属。魏斯氏菌属（Weissella）、弧菌属和乳杆菌属归为II类，在发酵前期相对丰度较高，为优势菌属，而在发酵后期相对丰度有所降低。III类包含了8 个细菌种属，在发酵前期相对丰度较低，后期相对丰度有所提高。IV包含了20 个细菌种属，在发酵前期相对丰度中等，后期相对丰度进一步降低。在细菌属水平上根据样本间丰度的相似性，24 个发酵液样本可以聚为2大类，1_1、1_2和1_3的9 个样本聚为一类，1_4、1_5、1_6、1_7和1_8的15 个样本归为一类，由此可将整个发酵过程分为2 个阶段，发酵前期和发酵后期。

2.4 发酵过程中大头菜发酵液样本细菌群落结构的分析

在对大头菜发酵液中门和属分类学水平细菌菌群相对丰度进行分析的基础上，进一步采用PCA和样本层级聚类分析24 个大头菜发酵液样本细菌群落结构的相似性和差异程度。PCA可以将原有的复杂数据降维，通过分析不同样本群落组成可以反映样本间的差异和距离。样本物种组成越相似，反映在PCA图中的距离则越近。由图4可知，在OTU分类水平上，不同样本细菌群落结构信息主要集中在前2 个主成分，其中PC1的贡献率为62.9%，PC2的贡献率为25.4%，累计方差贡献率为88.3%，因此可以认为PC1和PC2包含了样本的原始大量细菌组成数据信息，可以反映样本的整体细菌组成信息。24 个大头菜发酵液样本呈现出明显的聚类趋势，1_4～1_8的15 个样本距离较近，处于第2象限，且在PCA图中没有完全区分开，聚为一类。1_1和1_3的6 个样本距离较近，处于第3象限，聚为一类。1_2的3 个样本距离较近，处于第1象限，单独聚为一类。

图4 OTU水平上样本细菌群落PCAFig. 4 Principal component analysis of bacterial communities at the OTU level

在PCA基础上，进一步采用样本层级聚类分析大头菜发酵液样本细菌群落结构的相似性和差异程度。在OTU分类水平上，根据β多样性距离矩阵进行样本层级聚类分析，距离算法采用Bray-Curtis，树状结构的构建使用UPGMA算法，结果见图5。24 个大头菜发酵液样本可聚为2 个大类，其中1_4～1_8的15 个样本聚为A大类，1_1～1_3的9 个样本聚为B大类，而B大类又可分为两个小类，其中1_2的3 个样本聚为B1类，1_1和1_3的6 个样本聚为B2类。这与Heatmap图和PCA结果基本一致。由此也可将整个发酵过程分为2 个阶段，发酵前期和发酵后期。在发酵的相同阶段，细菌群落组成结构相似。

图5 OTU水平上样本细菌群落聚类分析Fig. 5 Cluster analysis of bacterial communities at the OTU level

2.5 不同发酵阶段细菌菌群共有和独有物种分析

为判断将所有样本分为两组是否有意义，采用ANOSIM进行分析。样本间的距离采用Bray-Curtis距离算法，置换次数设为999 次，结果如图6所示。组间距离最大值为276，最小值为136，中位数为209。发酵前期组内距离最大值为144，最小值为1，中位数为123.5。发酵后期组内距离最大值为114，最小值为2，中位数为62。统计值R为0.998 1，说明组间差异远远大于组内差异。P＝0.001（＜0.05），说明本次检验的可信度很高。因此，将整个发酵过程分为发酵前期和发酵后期两个阶段是有意义的、可行的。

图6 距离箱线图Fig. 6 Distances box plot

进一步统计分析发酵前期和发酵后期样本中共有的和独有的细菌菌群数量，研究不同分类水平样本的组成相似性和重叠程度。结果表明，24 个大头菜发酵液样本中共检出19 个门，发酵前期独有菌门有2 个，为芽单胞菌门（Gemmatimonadetes）和酸杆菌门（Acidobacteria），发酵后期独有菌门为互养菌门（Synergistetes），且3 个独有菌门在各样本中的相对丰度都极低，两个发酵阶段共有菌门为16 个，占有效序列条数的99.9%以上。在属水平上，发酵前期独有菌属有22 个，发酵后期独有菌属有28 个，共有菌属为243 个。OTU水平上，发酵前期独有OTU数目为67，发酵后期OTU数目为80，共有OTU数目为511。在发酵的两个不同阶段，各分类水平上的物种在数目上差异并不大，而且在各个阶段独有的物种相对丰度非常低，因此发酵两个阶段的细菌菌群结构差异主要表现在共有物种的组成上。

3 结 论

通过Illumina MiSeq高通量测序，24 份大头菜发酵液样本共获得1 158 752 条优化序列，对97%相似水平的OTU代表序列进行生物信息统计分析发现，抽平后所有样本序列归属于19 个门、32 个纲、67 个目、133 个科、293 个属、482 个种、658 个OTU。门水平上，在整个发酵过程中，变形菌门和厚壁菌门占据了绝对的优势，两者相对丰度之和达到75.0%～98.3%。在发酵前期，弧菌属为优势菌属，而随着发酵的进行，弧菌属相对丰度急剧下降，盐厌氧菌属成为相对丰度最高的优势菌属。除此之外，盐单胞菌属、色盐杆菌属和norank_f_TBZ33丰度也显著提高，成为发酵后期的优势菌属。在发酵前期和发酵后期，各分类水平上的物种在数目上差异并不大，而且在各个阶段独有的物种相对丰度非常低，因此发酵两个阶段的细菌菌群结构差异主要表现在共有物种的组成上。本研究揭示了襄阳大头菜发酵过程中细菌群落结构的动态演替，为进一步研制大头菜发酵剂、改进大头菜生产工艺提供一定的研究基础。