吳進菊，李宇昂，王梓杭，鄭婷婷，于 博，余海忠

（湖北文理學(xué)院食品科學(xué)技術(shù)學(xué)院g化學(xué)工程學(xué)院，湖北 襄陽 441053）

湖北襄陽大頭菜歷史悠久，相傳是三國時期諸葛亮所發(fā)明，因此又叫孔明菜、諸葛菜。襄陽大頭菜在2007年獲得“國家地理標志保護產(chǎn)品”，是我國四大名腌菜之一[1]。襄陽大頭菜制作周期長，需經(jīng)歷長達數(shù)月至數(shù)年的發(fā)酵過程才能成熟。在漫長的發(fā)酵過程中，有大量的微生物參與其中，蘊藏了豐富的微生物資源。

近年來，很多研究人員把目光轉(zhuǎn)向發(fā)酵食品的微生物多樣性研究，取得了很大的成果，如Liang Huipeng等[2]研究發(fā)現(xiàn)榨菜發(fā)酵初期乳酸菌生長迅速，然后達到穩(wěn)定的發(fā)酵水平，弧菌、鹽單胞菌、明串珠菌、魏斯氏菌在發(fā)酵過程中數(shù)量減少，片球菌數(shù)量增加。而家庭制作的泡菜中耐酸乳桿菌和短乳桿菌為優(yōu)勢菌屬[3]。青菜泡菜發(fā)酵過程中主要菌屬為乳酸菌屬、假單胞菌屬、弧菌屬和鹽單胞菌屬。大多數(shù)隸屬于變形菌門或擬桿菌門的細菌菌屬在發(fā)酵前期檢測頻率較高，而乳酸菌（隸屬于厚壁菌門）在發(fā)酵后期占主導(dǎo)地位[4]。Zhou Qi等[5]研究發(fā)現(xiàn)東北酸菜發(fā)酵過程中優(yōu)勢微生物為乳桿菌屬、明串珠菌屬、腸桿菌屬、熱袍菌屬、假單胞菌屬等，隨著發(fā)酵的進行，乳酸菌逐漸增多，而腸桿菌逐漸減少。趙慧君[6]、郭壯[7]等對襄陽大頭菜腌制液膜醭中的微生物菌群結(jié)構(gòu)進行了研究，結(jié)果表明，細菌中變形菌門和硬壁菌門為優(yōu)勢菌門，屬水平上鹽單胞菌、色鹽桿菌屬、海桿菌屬、鹽厭氧菌屬、四聯(lián)球菌為優(yōu)勢種屬；真菌中子囊菌門平均含量高達99.99%，屬水平上合囊菌屬、畢赤酵母屬、念珠菌屬為優(yōu)勢菌屬，平均相對含量分別為56.10%、40.96%和2.32%。另外，研究人員對發(fā)酵魚類[8-9]、酒類[10-14]、發(fā)酵乳制品[15-20]、發(fā)酵豆制品[21-24]、發(fā)酵蔬菜[25-30]等發(fā)酵食品中微生物多樣性方面也進行了大量的研究。

目前針對傳統(tǒng)發(fā)酵大頭菜發(fā)酵過程中微生物方面的研究較少，早期研究基本是建立在傳統(tǒng)純培養(yǎng)技術(shù)研究方法上，只能獲得微生物總量中極少的一部分，存在很大的局限性和片面性，從而對發(fā)酵過程中微生物群落結(jié)構(gòu)變化和種群多樣性了解不夠清楚，限制了大頭菜標準化和工業(yè)化生產(chǎn)。本研究采用Illumina MiSeq高通量測序技術(shù)考察了大頭菜發(fā)酵過程中細菌的多樣性和菌群結(jié)構(gòu)的演替規(guī)律，全面揭示了大頭菜發(fā)酵過程中細菌群落結(jié)構(gòu)的動態(tài)變化，為今后大頭菜生產(chǎn)的改進提供了一定的研究基礎(chǔ)。

1 材料與方法

1.1 材料與試劑

24 份大頭菜發(fā)酵液樣本 湖北孔明菜食品有限公司；FastPfu聚合酶 北京全式金生物技術(shù)有限公司；土壤基因組提取試劑盒 美國MP Biomedicals公司；DL2000 DNA Marker 寶生物工程（大連）有限公司。

1.2 儀器與設(shè)備

MiSeq高通量測序平臺 美國Illumina公司；Quanti Fluor?-ST藍色熒光定量系統(tǒng) 美國Promega公司；9700型聚合酶鏈式反應(yīng)（polymerase chain reaction，PCR）儀 美國ABI公司；臺式冷凍離心機 德國Eppendorf公司。

1.3 方法

1.3.1 樣本采集

實驗樣本采自湖北孔明菜食品有限公司，大頭菜采用“三腌、五鹵、六曬”的傳統(tǒng)工藝進行生產(chǎn)，按照加工工序的特點，在每一道腌制和鹵制工序完成后取樣，共采集了8 個不同發(fā)酵時期的大頭菜發(fā)酵液樣本，分別編號為1_1～1_8。每個發(fā)酵時期平行取3 個發(fā)酵液樣本，編號分別為a、b和c。樣本采集后快速運送至實驗室，置于-80 ℃冰箱保存?zhèn)溆谩?/p>

1.3.2 Illumina MiSeq測序和數(shù)據(jù)分析

將24 份大頭菜發(fā)酵液樣本送至上海美吉生物醫(yī)藥科技有限公司進行測序。以微生物總DNA為模板，以338F（5’-ACTCCTACGGGAGGCAGCAG-3’）和806R（5’-GGACTACHVGGGTWTCTAAT-3’）為上下游引物擴增細菌16S rRNA基因的V3-V4可變區(qū)序列。將PCR產(chǎn)物純化后進行熒光定量，進一步構(gòu)建MiSeq文庫后進行雙端測序，測序讀長300 bp。數(shù)據(jù)分析采用美吉I-Sanger云平臺進行在線分析。對優(yōu)化序列提取非重復(fù)序列，按照97%相似性水平對非重復(fù)序列（不含單序列）進行可操作分類單元（operational taxonomic units，OTU）聚類，去除嵌合體，獲得OTU代表序列。利用Mothur軟件計算各樣本的α多樣性，表征細菌群落的豐富度和多樣性。采用群落Bar圖和Heatmap圖反映樣本在各分類學(xué)水平上的物種組成。主成分分析（principal component analysis，PCA），采用降維的方法對數(shù)據(jù)進行簡化分析，利用物種豐度表，基于歐氏距離進行作圖，樣本物種組成越相似，反映在PCA圖中的距離越近。根據(jù)β多樣性距離矩陣進行層級聚類分析，使用非加權(quán)組平均法（unweighted pair-group method with arithmetic mean，UPGMA）算法構(gòu)建樹狀結(jié)構(gòu)，可視化呈現(xiàn)不同環(huán)境樣本中微生物進化的差異程度，距離算法采用Bray-Curtis，樣本層級聚類方式為Average。采用R語言（版本號3.2.5）vegan包進行相似性分析（analysis of similarities，ANOSIM）。ANOSIM又稱為置換多因素方差分析或非參數(shù)多因素方差分析，是一種非參數(shù)檢驗，用來檢驗組間的差異是否顯著大于組內(nèi)差異，從而判斷分組是否有意義。

1.3.3 樣本序列的提交及序列號

將本研究中所有fastq序列文件提交到美國國立生物技術(shù)信息中心SRA數(shù)據(jù)庫，序列號為PRJNA524770。

2 結(jié)果與分析

2.1 測序結(jié)果和α多樣性分析

通過Illumina MiSeq高通量測序，24 份大頭菜發(fā)酵液樣本共獲得1 158 752 條優(yōu)化序列，平均長度為448 bp，各樣本抽平多樣性指數(shù)如表1所示。采用RDP classifier貝葉斯算法對97%相似水平的OTU代表序列進行生物信息統(tǒng)計分析，按最小樣本序列數(shù)進行抽平分析（每個樣本27 407 條有效序列），結(jié)果表明，所有樣本序列歸屬于19 個門、32 個綱、67 個目、133 個科、293 個屬、482 個種、658 個OTU，每個樣本的OTU分類信息見表1。

表1 樣本信息表和α多樣性指數(shù)分析Table 1 Bacterial community composition and alpha diversity estimators of pickle juice samples

由表1可知，每個樣本的Coverage值均為0.996以上，表明測序深度良好，覆蓋率高，可以真實展示樣本中的絕大多數(shù)細菌。通過Shannon指數(shù)和Simpson指數(shù)分析可知，在所有樣本中，1_1的細菌多樣性最高，而1_2細菌多樣性相對較低。這從各樣本門、目、屬、OTU的個數(shù)也能看出，1_1在各分類學(xué)水平的數(shù)目均高于其他樣本，而1_2低于其他樣本。1_1、1_5和1_7中，Chao 1指數(shù)和ACE值較高，表明樣本中群落的豐富度較高。總體來說，發(fā)酵前期和發(fā)酵后期樣本的細菌多樣性沒有呈現(xiàn)出整體的差異性，在8 個不同的發(fā)酵時期，樣本細菌多樣性有的時期高，有的時期低，呈現(xiàn)不規(guī)律變化。

稀釋曲線通常是在97%的相似水平下，從樣本中隨機抽取一定數(shù)量的個體，以測序深度為橫坐標，以各個樣本中的細菌多樣性指數(shù)為縱坐標，在二維水平建立坐標系，將樣本內(nèi)細菌的多樣性隨著測序程度加深的變化情況清楚的展示出來。由圖1可知，24 個樣本隨著抽取的reads條數(shù)的增加，曲線平坦，已接近直線，說明此次取樣的樣本測序數(shù)據(jù)量足夠大，能比較全面地反映樣本中絕大多數(shù)細菌的多樣性信息。

圖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）、互養(yǎng)菌門（Synergistetes）等。其中變形菌門、厚壁菌門、放線菌門、擬桿菌門、螺旋菌門為優(yōu)勢菌門，其他14 個菌門相對豐度較低，在任一發(fā)酵時期中的相對豐度均在1%以下。不同發(fā)酵時期大頭菜樣本細菌群落組成如圖1所示。在整個發(fā)酵過程中，變形菌門和厚壁菌門占據(jù)了絕對的優(yōu)勢，兩者相對豐度之和在1_1中達到75.0%（平均值，下同），在其他發(fā)酵時期均達到了90%以上，這與大頭菜腌制液膜醭中的微生物菌群組成相似[6]。1_1中，變形菌門相對豐度最高，為67.7%，厚壁菌門、放線菌門、擬桿菌門三者相當，分別為7.32%、13.1%和10.4%。而在1_2中，厚壁菌門占據(jù)了絕對的優(yōu)勢，相對豐度達到了66.2%，變形菌門、放線菌門和擬桿菌門的相對豐度均明顯下降。在后續(xù)的發(fā)酵過程中，變形菌門的相對豐度都高于厚壁菌門。

圖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圖根據(jù)物種和樣本間豐度的相似性進行聚類，使高豐度和低豐度的物種分塊聚集，通過顏色變化與相似程度反映不同樣本在分類學(xué)水平上群落組成的相似性和差異性。對屬水平總豐度排在前35的物種進行層級聚類，物種和樣本層級聚類方式均采用平均聚類方式，得到的群落Heatmap如圖3所示。24 個大頭菜發(fā)酵液樣本中共檢測出293 個細菌屬，其中有46 個菌屬相對豐度在1%以上，相對豐度排在前5 位的菌屬有鹽厭氧菌屬（Halanaerobium）、弧菌屬（Vibrio）、鹽單胞菌屬（Halomonas）、乳桿菌屬（Lactobacillus）和色鹽桿菌屬（Chromohalobacter）。在發(fā)酵的前3 個時期，弧菌屬為優(yōu)勢菌屬，相對豐度為27.6%～42.4%，而在發(fā)酵的后5 個時期，弧菌屬相對豐度急劇下降，均在3%以下，而鹽厭氧菌屬相對豐度迅速增加，發(fā)酵前3 個時期相對豐度均在0.1%以下，而發(fā)酵后5 個時期相對豐度增加至31.1%～42.0%，占據(jù)了絕對優(yōu)勢。在1_4～1_8中，除鹽厭氧菌屬外，優(yōu)勢菌屬還有鹽單胞菌屬、色鹽桿菌屬和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%。

關(guān)于乳桿菌屬，在1_1中相對豐度僅為1.4%，而在1_2中提高至55.6%，占據(jù)了絕對優(yōu)勢，隨后在1_3中又下降至4.0%。在1_5～1_8中，乳桿菌屬相對豐度下降至1%以下。乳酸菌是發(fā)酵蔬菜中常見的優(yōu)勢菌屬[1-4]，在蔬菜發(fā)酵過程中起著非常關(guān)鍵的作用，可以抑制硝酸鹽還原細菌的作用，從而降低亞硝酸鹽的形成[3]。在發(fā)酵前期，加入的食鹽量相對較少，所以不耐鹽的弧菌屬、乳桿菌屬等為優(yōu)勢菌屬。而隨著發(fā)酵的進行，大量的食鹽不斷加入，使發(fā)酵液中的食鹽濃度大大提高，在發(fā)酵結(jié)束后發(fā)酵液中的食鹽濃度甚至達到飽和狀態(tài)，因此，發(fā)酵后期不耐鹽的菌屬相對豐度急劇下降，而耐鹽菌屬發(fā)展成為優(yōu)勢菌屬，如鹽厭氧菌屬、鹽單胞菌屬、色鹽桿菌屬等。

由圖3可知，根據(jù)物種間豐度的相似性，總豐度排在前35的細菌菌屬可聚為4 類，鹽厭氧菌屬、norank_f_TBZ33、鹽單胞菌屬和色鹽桿菌屬歸為I類，在發(fā)酵前期相對豐度非常低，而在發(fā)酵后期相對豐度較高，成為優(yōu)勢菌屬。魏斯氏菌屬（Weissella）、弧菌屬和乳桿菌屬歸為II類，在發(fā)酵前期相對豐度較高，為優(yōu)勢菌屬，而在發(fā)酵后期相對豐度有所降低。III類包含了8 個細菌種屬，在發(fā)酵前期相對豐度較低，后期相對豐度有所提高。IV包含了20 個細菌種屬，在發(fā)酵前期相對豐度中等，后期相對豐度進一步降低。在細菌屬水平上根據(jù)樣本間豐度的相似性，24 個發(fā)酵液樣本可以聚為2大類，1_1、1_2和1_3的9 個樣本聚為一類，1_4、1_5、1_6、1_7和1_8的15 個樣本歸為一類，由此可將整個發(fā)酵過程分為2 個階段，發(fā)酵前期和發(fā)酵后期。

2.4 發(fā)酵過程中大頭菜發(fā)酵液樣本細菌群落結(jié)構(gòu)的分析

在對大頭菜發(fā)酵液中門和屬分類學(xué)水平細菌菌群相對豐度進行分析的基礎(chǔ)上，進一步采用PCA和樣本層級聚類分析24 個大頭菜發(fā)酵液樣本細菌群落結(jié)構(gòu)的相似性和差異程度。PCA可以將原有的復(fù)雜數(shù)據(jù)降維，通過分析不同樣本群落組成可以反映樣本間的差異和距離。樣本物種組成越相似，反映在PCA圖中的距離則越近。由圖4可知，在OTU分類水平上，不同樣本細菌群落結(jié)構(gòu)信息主要集中在前2 個主成分，其中PC1的貢獻率為62.9%，PC2的貢獻率為25.4%，累計方差貢獻率為88.3%，因此可以認為PC1和PC2包含了樣本的原始大量細菌組成數(shù)據(jù)信息，可以反映樣本的整體細菌組成信息。24 個大頭菜發(fā)酵液樣本呈現(xiàn)出明顯的聚類趨勢，1_4～1_8的15 個樣本距離較近，處于第2象限，且在PCA圖中沒有完全區(qū)分開，聚為一類。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基礎(chǔ)上，進一步采用樣本層級聚類分析大頭菜發(fā)酵液樣本細菌群落結(jié)構(gòu)的相似性和差異程度。在OTU分類水平上，根據(jù)β多樣性距離矩陣進行樣本層級聚類分析，距離算法采用Bray-Curtis，樹狀結(jié)構(gòu)的構(gòu)建使用UPGMA算法，結(jié)果見圖5。24 個大頭菜發(fā)酵液樣本可聚為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結(jié)果基本一致。由此也可將整個發(fā)酵過程分為2 個階段，發(fā)酵前期和發(fā)酵后期。在發(fā)酵的相同階段，細菌群落組成結(jié)構(gòu)相似。

圖5 OTU水平上樣本細菌群落聚類分析Fig. 5 Cluster analysis of bacterial communities at the OTU level

2.5 不同發(fā)酵階段細菌菌群共有和獨有物種分析

為判斷將所有樣本分為兩組是否有意義，采用ANOSIM進行分析。樣本間的距離采用Bray-Curtis距離算法，置換次數(shù)設(shè)為999 次，結(jié)果如圖6所示。組間距離最大值為276，最小值為136，中位數(shù)為209。發(fā)酵前期組內(nèi)距離最大值為144，最小值為1，中位數(shù)為123.5。發(fā)酵后期組內(nèi)距離最大值為114，最小值為2，中位數(shù)為62。統(tǒng)計值R為0.998 1，說明組間差異遠遠大于組內(nèi)差異。P＝0.001（＜0.05），說明本次檢驗的可信度很高。因此，將整個發(fā)酵過程分為發(fā)酵前期和發(fā)酵后期兩個階段是有意義的、可行的。

圖6 距離箱線圖Fig. 6 Distances box plot

進一步統(tǒng)計分析發(fā)酵前期和發(fā)酵后期樣本中共有的和獨有的細菌菌群數(shù)量，研究不同分類水平樣本的組成相似性和重疊程度。結(jié)果表明，24 個大頭菜發(fā)酵液樣本中共檢出19 個門，發(fā)酵前期獨有菌門有2 個，為芽單胞菌門（Gemmatimonadetes）和酸桿菌門（Acidobacteria），發(fā)酵后期獨有菌門為互養(yǎng)菌門（Synergistetes），且3 個獨有菌門在各樣本中的相對豐度都極低，兩個發(fā)酵階段共有菌門為16 個，占有效序列條數(shù)的99.9%以上。在屬水平上，發(fā)酵前期獨有菌屬有22 個，發(fā)酵后期獨有菌屬有28 個，共有菌屬為243 個。OTU水平上，發(fā)酵前期獨有OTU數(shù)目為67，發(fā)酵后期OTU數(shù)目為80，共有OTU數(shù)目為511。在發(fā)酵的兩個不同階段，各分類水平上的物種在數(shù)目上差異并不大，而且在各個階段獨有的物種相對豐度非常低，因此發(fā)酵兩個階段的細菌菌群結(jié)構(gòu)差異主要表現(xiàn)在共有物種的組成上。

3 結(jié) 論

通過Illumina MiSeq高通量測序，24 份大頭菜發(fā)酵液樣本共獲得1 158 752 條優(yōu)化序列，對97%相似水平的OTU代表序列進行生物信息統(tǒng)計分析發(fā)現(xiàn)，抽平后所有樣本序列歸屬于19 個門、32 個綱、67 個目、133 個科、293 個屬、482 個種、658 個OTU。門水平上，在整個發(fā)酵過程中，變形菌門和厚壁菌門占據(jù)了絕對的優(yōu)勢，兩者相對豐度之和達到75.0%～98.3%。在發(fā)酵前期，弧菌屬為優(yōu)勢菌屬，而隨著發(fā)酵的進行，弧菌屬相對豐度急劇下降，鹽厭氧菌屬成為相對豐度最高的優(yōu)勢菌屬。除此之外，鹽單胞菌屬、色鹽桿菌屬和norank_f_TBZ33豐度也顯著提高，成為發(fā)酵后期的優(yōu)勢菌屬。在發(fā)酵前期和發(fā)酵后期，各分類水平上的物種在數(shù)目上差異并不大，而且在各個階段獨有的物種相對豐度非常低，因此發(fā)酵兩個階段的細菌菌群結(jié)構(gòu)差異主要表現(xiàn)在共有物種的組成上。本研究揭示了襄陽大頭菜發(fā)酵過程中細菌群落結(jié)構(gòu)的動態(tài)演替，為進一步研制大頭菜發(fā)酵劑、改進大頭菜生產(chǎn)工藝提供一定的研究基礎(chǔ)。