基于高通量測序對建昌板鴨保藏過程中的 細菌群落動態變化的研究

◎ 李家怡，林 巧，刁尚鵬，吳 兵，劉軍華，楊雨欣，任妍菱，陳 瑤，陳亦然，鄧友勝，祝藝嘉

（1.攀西特色作物研究與利用四川省重點實驗室，四川 西昌 615013；2.德昌縣茂源長（童耳朵）食品有限責任公司，四川 德昌 615599）

建昌板鴨由德昌特色建昌鴨經傳統腌制風干工藝處理一個月左右而成。目前尚無科學的鹽水鴨長期貯存方法，易繁殖微生物，導致鴨肉失鮮、失色、失去光澤、發黃甚至腐敗變質[1]。細菌在板鴨產品的貯藏過程中起著關鍵的作用，其決定了建昌板鴨的貨架期長短。以往對建昌板鴨微生物區系的研究大多基于微生物的純化和鑒定，工作量較大，易遺漏一些不能培養的微生物，不能反映建昌板鴨不同時期的細菌組成情況。高通量測序技術有助于更全面、真實地掌握食品保藏過程中微生物的種群結構和變化、更全面地了解建昌板鴨貯藏過程中的菌群變化規律。本實驗以板鴨保藏過程中不同時間段的細菌群落構成為研究對象，采用llumina MiSeq高通量測序技術對11個樣品中的細菌菌群結構及多樣性變化規律進行研究，并對建昌板鴨貨架期和優勢腐敗菌進行了研究[1]。

1 材料與方法

1.1 樣品采集

采用《肉與肉制品 取樣方法》（GB/T 9695.19—2008），對建昌板鴨儲存的第1 d、11 d、13 d、17 d、21 d、25 d、33 d、37 d、41 d、49 d和62 d進行無菌取樣，每次取樣25～30 g，取樣后置于-80 ℃冰箱冷藏。

1.2 儀器與設備

高通量測序儀（百邁客生物科技有限公司）；超凈工作臺（無錫一凈凈化設備有限公司）；HiSeq測序平臺（Illumina公司）。

1.3 高通量測序

委托百邁客生物科技有限公司進行測序。

1.3.1 基因序列擴增及測序

對樣品進行DNA提取后，用PCR儀進行擴增，擴增區域為細菌的16SrRNA基因V3～V4區。擴增完成后，使用2%瓊脂糖凝膠電泳進行檢測，將DNA純化后建庫。最后，利用HiSeq臺進行測序，對測序結果進行處理，對樣品進行多樣性分析、組成分析和樣品差異分析。

1.3.2 序列分析

測序得到的原始數據（Raw Data）中有一定的干擾數據（Dirty Data），將原始數據進行拼接、過濾、去嵌合體后，得到有效數據（Effective Tags），在97%相似度閾值下，用RDP軟件將序列在SILVA和UNITE分類學數據庫中進行對比，再進行物種分類，最后參考各序列之間的距離，對多條序列進行OTU聚類。對OTU的數量、Beta-deversity、Alphadeversity以及菌種在每個分類水平上的具體分布情況進行分析，得到相關的微生物群落結構及組成變化。利用MTHUR軟件對這些微生物群落的多樣性指數和豐度指數進行分析。

1.3.3 數據處理

在Usearch軟件上進行操作分類單元（Operational Taxonomic Unit，OTU）聚類分析。一個物種有一個OTU，因此OTU分析可確定物種數量。根據OTU聚類結果和各分類層次的物種比較數據，利用QIME軟件，根據分類信息確定OTU多樣性的相對指數和序列深度是否足夠，分析各分類層次微生物群落的組成及樣品的組成和多樣性。

2 結果與分析

2.1 16S rDNA V3區測序分析

處理獲得的序列信息，去除低質量片段、檢測出錯誤的序列、連接后引物序列以及無法與數據庫中序列信息比對一致的序列，獲得的序列為有效信息，見表1。有效序列所占比例符合預期[2]。

表1 預處理過程統計表

2.2 測序數據的合理性分析

稀釋性曲線（Rarefaction Curve）從樣本中隨機抽取若干序列，統計這些序列所代表的物種數目，并用序列數與物種數構建曲線，用于驗證測序數據量是否足以反映樣品中的物種多樣性，并間接反映樣本中的物種豐富度[3]。圖1反映了在連續取樣情況下新特征（新物種）出現的速率，在一定范圍內，隨著測序條數的加大，若曲線表現為急劇上升則表示群落中有大量物種被發現；當曲線趨于平緩，則表示此環境中的物種不會隨測序數量的增加而顯著增多。稀釋曲線可作為對各樣本測序量是否充分的判斷，曲線急劇上升表明測序量不足，需增加序列條數；反之，則表明樣品序列充分，可進行數據分析。

為驗證測序量是否能真實反映原始樣品的微生物群落的多樣性，對建昌板鴨第1 d、11 d、13 d、17 d、21 d、25 d、33 d、37 d、41 d、49 d和62 d的微生物組成多樣性利用稀釋曲線進行評估。圖1顯示16S測序的結果中，每個樣品的稀釋曲線在測序量為10 000 bp左右時基本不再變化，表示本次測序的深度已足夠，說明取樣合理，測序深度已基本覆蓋樣品中所有微生物，同時也表明了11個樣本微生物組成具有豐富的多樣性。

圖1 16S測序的樣品稀釋曲線圖

2.3 Venn圖分析

2.3.1 OTU數圖

將超過97%相似的序列劃分為一個OTU，每個OTU對應一個物種。每個樣本的OTU數量通過對樣本進行聚類得到，如圖2所示，樣品中細菌數量在貯藏期間呈先減少后增加的趨勢，同時也能看到樣品中細菌數量在貯藏初始期最高。

圖2 16S測序的OTU數圖

2.3.2 酒醅菌群的花瓣圖

Venn圖也可以用花瓣圖表示，對11個樣本的所有OTU數據進行統計分析，觀察每個樣本共有的OTU和每個樣本獨有的OTU（即OTU只存在于一個樣本中）。如圖3所示，在16S測序的OTU花瓣圖中，11個樣品中有41條共同的OTU數，僅有發酵后期第41 d、49 d、62 d的樣品存在獨有的OTU數，說明各樣品間大部分是相似微生物，即板鴨產品的整個保藏過程中，獨特的細菌種類并不多，但在貯藏期結束時有獨特的細菌參與，出現優勢腐敗菌，其中樣品9、10、11具有特有的OTU數分別為1個、1個、4個，即在板鴨貯藏的末期存在獨特的細菌參與，隨著細菌種類的增多，板鴨的腐敗逐漸加速。

圖3 16S測序的OTU花瓣圖

2.4 物種注釋及分類學分析

2.4.1 OTU聚類結果

先將低豐度的OTU處理后，根據最后得到的OTU計算出每個樣品中各等級能得到的tags數，再計算出在各分類水平下，物種類型在每個樣品中的數目。如表2所示，微生物種類隨時間的推移，基本呈現先減少后增加的趨勢。

表2 16S測序的各等級中的物種數表

2.4.2 菌群多水平上的物種分布圖

（1）門水平。圖4中，色塊的大小說明其物種的相對豐度比例。為得到最好的視圖，僅展示前10個物種的豐富度，后10個物種的總體顯示為Others。11個樣品在門水平上，含細菌微生物最多是厚壁菌門（Firmicutes），其次為變形菌門（Proteobacteria），且隨著時間的變化微生物種類逐漸增多，但在取樣的第62 d的樣品中，放線菌門（Actinobacteria）所占的比例遠高于其他樣品，說明放線菌門是板鴨保藏后期的優勢菌，可能為腐敗菌。

圖4 細菌在門水平上的分布圖

（2）綱水平。由圖5可知，11個樣品在綱水平上，含細菌微生物最多的是芽孢桿菌綱（Bacilli），但在取樣的第62 d的樣品中，放線菌綱（Actinobacteria）所占的比例突增，說明放線菌綱是板鴨保藏后期的優勢菌，可能為腐敗菌。

圖5 細菌在綱水平上的分布圖

（3）目水平。由圖6可知，11個樣品在目水平上，含細菌微生物優勢菌為乳桿菌目（Lactobacillales）和芽孢桿菌目（Bacillales），但在取樣的第62 d的樣品中，棒狀桿菌目（Corynebacteriales）所占的比例突增，說明放線菌綱中的棒狀桿菌目是板鴨腐敗期的優勢菌。

圖6 細菌在目水平上的分布圖

（4）科水平。由圖7可知，11個樣品在科水平上，含細菌微生物優勢菌為葡萄球菌科（Staphylococcaceae）和腸球菌科（Enterococcaceae），第62 d的樣品中，棒狀桿菌科（Corynebacteriaceae）所占的比例突增，說明放線菌綱中的棒狀桿菌目棒狀桿菌科是板鴨腐敗期的優勢菌。

圖7 細菌在科水平上的分布圖

（5）屬水平。由圖8可知，貯藏末期乳桿菌屬（Lactobacillus）逐漸增多，第62 d的樣品中，棒狀桿菌屬（Corynebacterium）所占的比例突增，說明乳桿菌屬、棒狀桿菌屬是板鴨腐敗期的優勢菌。

圖8 細菌在屬水平上的分布圖

（6）菌群豐度聚類熱圖。熱圖是數據矩陣中值的一種圖形表示形式，根據種類或樣本豐度相似性進行聚類。將高豐度和低豐度的物種分塊聚集，并通過顏色梯度和相似度反映群落組成的相似性和差異性。根據每個樣本的種類組成和相對豐度，在門水平上提取物種，用R語言工具繪制，在門水平上進行熱圖的數據聚類分析[4]。

熱圖聚類結果中，縱向聚類表示不同物種在各樣品間豐度的相似情況，兩物種間距離越近，枝長越短，說明這兩個物種在各樣品間的豐度越相似；橫向聚類表示不同樣品的各物種豐度的相似情況，與縱向聚類一樣，兩樣品間距離越近，枝長越短，說明這兩個樣品的各物種豐度越相似。樣品門水平下的物種豐度聚類熱圖見圖9，在門水平上，16S測序中第17 d、62 d的樣品物種豐度最相似。

圖9 16S測序的豐度聚類熱圖

2.5 Alpha多樣性分析

建昌板鴨在儲藏過程中的Alpha多樣性指數變化見表3，chao1指數能正確反應OTU數的大小，且值越大，樣品的多樣性越大。第41 d、49 d、62 d的樣品Shannon指數及Simpson指數逐漸增加，并在最后達到最高值，說明儲存末期的板鴨細菌群落豐富度逐漸增大和均衡化程度逐漸降低。11份樣品的測序覆蓋指數均在0.99以上，說明測序數量已達飽和，結果可以反映微生物群落的多樣性。

2.6 Beta多樣性分析

對板鴨不同貯藏階段的主成分進行主成分分析（PCA），相關矩陣特征值如表3所示，共提取3個特征值不小于1的主成分，其貢獻率分別為74.74%、14.95%、3.43%，主成分分析一般提取包含90%以上信息的主成分，根據主成分特征值及旋轉成分矩陣，計算得到板鴨主成分得分，以PC1為橫坐標、PC2和PC3為縱坐標，將樣本投影到三維坐標系上，得到不同貯藏階段樣品分布情況[5]。由圖10可知，在16S測序中，BY1和BY8、BY3和BY6、BY7和BY5最相似。

表3 16S測序中α多樣性的相關指數表

圖10 建昌板鴨各樣品關鍵風味物質主成分得分圖

3 結論與討論

本實驗通過高通量測序平臺對建昌板鴨貯藏過程中的微生物進行16S rRNA基因分析，探究其多樣性以及群落結構。通過多樣性分析發現，隨著測序量的增加，各多樣性曲線逐漸趨于平行，體現測序達到飽和，測序結果基本能反映不同時期板鴨微生物菌群多樣性組成，說明建昌板鴨貯藏過程中其微生物組成具有豐富的多樣性。通過群落結構分析發現，其細菌主要分布于厚壁菌門、變形菌門、放線菌門和擬桿菌門等11個門，15個綱，92個屬，厚壁菌門為優勢菌門。放線菌門隨著貯存時間的延長也逐漸增多，厚壁桿菌門中芽孢桿菌綱中的乳桿菌屬、放線菌門中棒狀桿菌目棒狀桿菌屬占一定優勢。相比其他研究方法，本方法對貯藏過程中的建昌板鴨細菌群落分析更全面，信息量更大，并且利用高通量測序技術能快速了解貯藏過程中的主要腐敗菌為乳桿菌屬及棒狀桿菌屬，其中乳桿菌屬為主要優勢菌，對于建昌板鴨的貯藏和保鮮和研究優勢腐敗菌的生長變化，深入了解肉品表面微生物群落特征從而延長板鴨的貨架期具有十分重要的意義。