基于Illumina MiSeq測(cè)序技術(shù)分析大鯢肉冷藏過程中微生物菌群演替規(guī)律

趙 萍，金文剛,3, ，蘭阿峰，劉俊霞，裴金金，陳德經(jīng),3

（1.陜西理工大學(xué)生物科學(xué)與工程學(xué)院，秦巴生物資源與生態(tài)環(huán)境省部共建培育國家重點(diǎn)實(shí)驗(yàn)室，陜西 漢中 723001；2.陜西理工大學(xué) 陜西省資源生物重點(diǎn)實(shí)驗(yàn)室，陜西 漢中 723001；3.陜西理工大學(xué) 陜南秦巴山區(qū)生物資源綜合開發(fā)協(xié)同創(chuàng)新中心，陜西 漢中 723001）

大鯢（）是我國人工養(yǎng)殖成熟的重點(diǎn)水產(chǎn)動(dòng)物品種之一，已在多個(gè)省、市實(shí)現(xiàn)了規(guī)模化人工養(yǎng)殖，其食用和藥用價(jià)值極高，含有豐富的氨基酸、人體必需脂肪酸、礦物質(zhì)、微量元素及維生素。隨著大鯢養(yǎng)殖規(guī)模的不斷擴(kuò)大，商品鯢價(jià)格不斷下跌，給養(yǎng)殖企業(yè)、養(yǎng)殖戶造成了較大經(jīng)濟(jì)損失，大鯢精深加工與開發(fā)利用已成為制約產(chǎn)業(yè)發(fā)展的瓶頸。迄今，研究人員已在大鯢營養(yǎng)評(píng)價(jià)、分割加工與貯藏保鮮、生物活性成分、熱加工方便食品等方面進(jìn)行了較多嘗試，促進(jìn)了一些高附加值產(chǎn)品的上市以及產(chǎn)業(yè)鏈的延伸。目前，大鯢低溫肉制品（冷鮮和速凍肉）已成為最重要的初加工產(chǎn)品形態(tài)，但冷鮮肉存在易腐敗變質(zhì)、貨架期短的缺陷。肉類腐敗變質(zhì)主要是由微生物的生長繁殖導(dǎo)致，而引起腐敗變質(zhì)的菌群主要為少數(shù)優(yōu)勢(shì)特定腐敗菌。查明并抑制引起大鯢肉腐敗變質(zhì)的優(yōu)勢(shì)菌群，是冷藏過程中實(shí)現(xiàn)品質(zhì)控制和延長貨架期的根本。

Illumina MiSeq高通量測(cè)序技術(shù)是研究微生物多樣性的一種分子生物學(xué)手段，該技術(shù)以細(xì)菌基因序列信息為基礎(chǔ)，通過細(xì)菌的基因序列研究水產(chǎn)品貯藏過程中微生物（尤其是優(yōu)勢(shì)腐敗菌）的組成及變化。傳統(tǒng)的微生物群落多樣性研究，主要是通過生理生化實(shí)驗(yàn)及表型進(jìn)行鑒定，耗費(fèi)時(shí)間且不能對(duì)其進(jìn)行精確鑒定，難以獲得微生物群落多樣性變化的真實(shí)情況。與傳統(tǒng)研究方法相比，Illumina MiSeq高通量測(cè)序技術(shù)具有無需分離純化、測(cè)定速度快、結(jié)果精準(zhǔn)、利用微量樣品就能實(shí)現(xiàn)所有微生物檢測(cè)、安全性和自動(dòng)化程度高、分析全面等優(yōu)點(diǎn)。其中16S rDNA被認(rèn)為是最適于細(xì)菌系統(tǒng)發(fā)育和分類鑒定的指標(biāo)，通常作為揭示生物物種的特征核苷酸序列，16S rDNA擴(kuò)增子測(cè)序是對(duì)高變區(qū)進(jìn)行測(cè)序分析和菌種鑒定，既可以提高效率還可以節(jié)約成本。目前，高通量測(cè)序技術(shù)已成功應(yīng)用于深水玫瑰蝦、鯛魚、貽貝、鯧魚等水產(chǎn)品低溫貯藏過程中微生物多樣性分析。

隨著大鯢分割加工的產(chǎn)業(yè)化推進(jìn)，有關(guān)大鯢分割肉氣調(diào)包裝、冰藏和微凍貯藏理化品質(zhì)指標(biāo)的研究已有報(bào)道。冷鮮肉被譽(yù)為21世紀(jì)鮮肉消費(fèi)的主流，課題組前期對(duì)大鯢肉冷藏過程中主要理化指標(biāo)和揮發(fā)性成分進(jìn)行了初步探討，然而大鯢肉冷藏過程中微生物菌群組成和特定腐敗菌的研究，還鮮見報(bào)道。為此，本研究采用高通量測(cè)序技術(shù)探究大鯢肉冷藏期間微生物菌群演替規(guī)律，旨在為后期靶向抑菌及冷藏保鮮奠定基礎(chǔ)。

1 材料與方法

1.1 材料與試劑

鮮活健康子二代大鯢3 尾（（2.51±0.36）kg），購自陜西漢中龍頭山水產(chǎn)養(yǎng)殖開發(fā)有限公司大鯢生態(tài)養(yǎng)殖基地。

三氯乙酸 上海阿拉丁生化科技股份有限公司；甲基紅、亞甲基藍(lán) 北京鼎國生物技術(shù)有限責(zé)任公司；氧化鎂天津市百世化工有限公司；硼酸 天津市福晨化學(xué)試劑廠；鹽酸 杭州匯普化工儀器有限公司；瓊脂 北京奧博星生物技術(shù)有限責(zé)任公司；Qubit dsDNA Assay Kit試劑盒美國Life Technologies公司；Magpure Soil DNA LQ Kit試劑盒 廣州Magen公司；Tks Gflex DNA ploymerase試劑盒 寶日醫(yī)生物技術(shù)（北京）有限公司。

1.2 儀器與設(shè)備

AD18型分散均質(zhì)機(jī) 上海昂尼儀器儀表有限公司；SW-CJ-IFD型潔凈工作臺(tái) 上海博迅實(shí)業(yè)公司；LQ-C1003型電子天平 深圳市飛亞衡器有限公司；FA3204B型電子天平 上海精科天美科學(xué)儀器有限公司；COHS-250型生化培養(yǎng)箱 常州金壇良友儀器有限公司；Centrifuge 5418型臺(tái)式高速離心機(jī) 德國Eppendorf公司；580BR10905型聚合酶鏈?zhǔn)椒磻?yīng)（polymerase chain reaction，PCR）儀 美國Bio-Rad公司；SN 002358型QIAxtractor高通量核酸提取儀 德國QIAGEN公司；2100型生物分析儀 美國Aglient公司；2000型NanoDrop微量分光光度計(jì) 美國Thermo Fisher公司；2500型凝膠成像儀 北京Tanon公司。

1.3 方法

1.3.1 大鯢肉樣品的制備

鮮活大鯢經(jīng)放血、熱燙、刮黏液、去內(nèi)臟和清洗后，用聚乙烯袋20 min內(nèi)運(yùn)回實(shí)驗(yàn)室，去頭、皮、四肢、尾，并切成約5.0 cm×2.0 cm×0.5 cm的肉塊，放入黑色托盤中并用保鮮膜密封4 ℃冷藏，于宰后第0、2、4、6、8天取樣，每個(gè)時(shí)間點(diǎn)取6 個(gè)平行，分別記為T0（T0-1、T0-2、T0-3、T0-4、T0-5、T0-6），T2、T4、T6、T8的標(biāo)記與T0相似，定時(shí)將樣品取出，置于超凈工作臺(tái)上，去除保鮮膜，用事先滅菌好的剪刀將肉樣剪碎，混勻，用鑷子裝在15 mL離心管中，置于-80 ℃冰箱中貯藏，用于細(xì)菌菌落總數(shù)、揮發(fā)性鹽基氮（total volatile basic nitrogen，TVB-N）測(cè)定和高通量測(cè)序。

1.3.2 細(xì)菌菌落總數(shù)和TVB-N測(cè)定

按照GB 4789.2—2016《食品微生物學(xué)檢驗(yàn) 菌落總數(shù)測(cè)定》和GB 5009.228—2016《食品中揮發(fā)性鹽基氮的測(cè)定》方法，對(duì)不同冷藏時(shí)間大鯢肉中菌落總數(shù)和TVB-N含量進(jìn)行測(cè)定。

1.3.3 高通量測(cè)序

將制取的不同冷藏時(shí)間大鯢肉樣本用干冰保存運(yùn)送至上海歐易生物醫(yī)學(xué)科技有限公司，基于Illumina MiSeq平臺(tái)進(jìn)行高通量測(cè)序。

1.3.3.1 基因組DNA提取

按照Magpure Soil DNA LQ Kit試劑盒說明書提取樣本的基因組DNA，在2 mL的離心管中，加入約500 mg的玻璃珠，再加入0.25～0.5 g大鯢肉樣品和0.6～0.8 mL Buffer SOL，在渦旋儀上以最高速度渦旋5～10 min，再加入60 μL Buffer SDS至樣品中，渦旋混勻30 s，詳細(xì)提取方法按照試劑盒說明書進(jìn)行。之后利用1%瓊脂糖凝膠電泳和NanoDrop 2000微量分光光度計(jì)檢測(cè)DNA濃度。

1.3.3.2 PCR擴(kuò)增與建庫

以基因組DNA為模板，根據(jù)測(cè)序區(qū)域的選擇，使用帶Barcode的特異引物，TaKaRa公司的Tks Gflex DNA polymerase進(jìn)行1輪PCR，確保擴(kuò)增效率和準(zhǔn)確性。第1輪PCR產(chǎn)物使用電泳檢測(cè)，檢測(cè)后使用磁珠純化，純化后作為2輪PCR模板，并進(jìn)行2輪PCR擴(kuò)增，并再次使用電泳檢測(cè)，檢測(cè)后使用磁珠純化，純化后對(duì)PCR產(chǎn)物進(jìn)行瓊脂糖凝膠電泳檢測(cè)和Qubit定量。根據(jù)PCR產(chǎn)物濃度進(jìn)行等量混樣，并利用Illumina Novaseq 6000測(cè)序系統(tǒng)上機(jī)測(cè)序。細(xì)菌多樣性鑒定對(duì)應(yīng)區(qū)域16S V3-V4區(qū)。引物序列：343F 5′-TACGGRAGGCAGCAG-3′；798R 5′-AGGGTATCTAATCCT-3′。

1輪PCR體系：15 μL 2×Gflex PCR Buffer、1 μL 5 pmol/μL Primer F、1 μL 5 pmol/μL Primer R、≥1 μL（50 ng）Template DNA、0.6 μL TKS Gflex DNA Polymerase（1.25 U/μL），加水補(bǔ)足至30 μL。PCR程序：94 ℃預(yù)變性5 min，循環(huán)1 次；94 ℃變性30 s；56 ℃退火30 s，循環(huán)26 次；72 ℃延伸20 s；72 ℃延伸5 min，循環(huán)1 次；4 ℃保藏。

2輪PCR體系：15 μL 2×Gflex PCR Buffer、0.6 μL TKS Gflex DNA Polymerase（1.25 U/μL）、1 μL Adapter I5、1 μL Adapter I7、50 ng第1輪產(chǎn)物，加水補(bǔ)足至30 μL。PCR程序：94 ℃預(yù)變性5 min，循環(huán)1 次；94 ℃變性30 s；56 ℃退火30 s，循環(huán)7 次；72 ℃延伸20 s；72 ℃延伸5 min，循環(huán)1 次；4 ℃保藏。

1.4 數(shù)據(jù)統(tǒng)計(jì)

測(cè)序完成后，使用Trimmomatic軟件對(duì)原始雙端系列進(jìn)行去雜，去雜后的雙端系列利用FLASH軟件進(jìn)行拼接。為保證結(jié)果的準(zhǔn)確性，可進(jìn)行精準(zhǔn)去雜，去除含有模糊堿基、單堿基高重復(fù)區(qū)的序列與及長度過短的序列，同時(shí)，利用UCHIME軟件檢測(cè)并去除序列中的嵌合體。處理后形成的有效序列，采用Vsearch軟件，根據(jù)序列的相似性，將相似度≥97%的序列歸為多個(gè)操作分類單元（operational taxonomic unit，OTU）。使用QIIME軟件包挑選出各個(gè)OTU的代表序列，并將所有代表序列與數(shù)據(jù)庫進(jìn)行對(duì)比注釋。16S使用Silva（version132）數(shù)據(jù)庫對(duì)比，物種對(duì)比注釋使用RDP classifier軟件，保留置信區(qū)間＞0.7的注釋結(jié)果，物種比對(duì)注釋使用BLAST軟件?；诜诸愋畔W(xué)，探討大鯢肉在冷藏期間微生物群落結(jié)構(gòu)變化。使用Excel對(duì)原始數(shù)據(jù)進(jìn)行處理、匯總，采用SPSS 25進(jìn)行顯著性分析，其中多重比較采用Duncan法，＜0.05，差異顯著，利用Origin 2021進(jìn)行繪圖。

2 結(jié)果與分析

2.1 大鯢肉冷藏過程中菌落總數(shù)和TVB-N值的變化

通常水產(chǎn)魚類菌落總數(shù)≥6（lg（CFU/g））時(shí)達(dá)到腐敗，由圖1A可知，隨著冷藏時(shí)間的延長，菌落總數(shù)總體呈上升的趨勢(shì)，在0～4 d增加緩慢，4～8 d增加較快，在第6天超過水產(chǎn)品可接受的上限，第8天達(dá)到最大值8.41（lg（CFU/g）），已經(jīng)超出了可食用水產(chǎn)品的微生物閾值。在冷藏初期，菌落總數(shù)增長速率較慢，主要是低溫、pH值的降低與及初始腐敗菌落數(shù)較低。隨著冷藏時(shí)間的延長，菌落總數(shù)增長速率加快，主要是由于蛋白質(zhì)的自溶及代謝物的積累有利于腐敗菌的生長繁殖。TVB-N是由具有揮發(fā)性的氨、伯胺、仲胺及叔胺等低級(jí)堿性含氮化合物組成，主要是由于食品在酶和微生物的作用下，使蛋白質(zhì)及非蛋白的含氮化合物降解產(chǎn)生氨及胺類等揮發(fā)性堿性含氮化合物。根據(jù)G B/T 18108—2008 規(guī)定，魚肉中TVB-N 值達(dá)到30 mg/100 g，是消費(fèi)者可接受的上限，由圖1B可知，大鯢肉在冷藏期間TVB-N值呈增加的趨勢(shì)，在前4 d增加緩慢，TVB-N值由6.61 mg/100 g增加到12 mg/100 g，之后快速增加并在第6天接近可接受上限，在第8天超出可接受上限達(dá)到最大值47.8 mg/100 g。大鯢肉在發(fā)生自溶生化變化前期，蛋白質(zhì)降解速度較慢，引起蛋白質(zhì)降解的主要因素為內(nèi)源酶，在自溶階段，蛋白質(zhì)的分解為微生物的生長繁殖提供了營養(yǎng)物質(zhì)，使得微生物迅速增長，然后微生物產(chǎn)生各種蛋白酶作用于蛋白質(zhì)，又使得蛋白質(zhì)分解，所以在自溶階段后TVB-N值迅速增加。

圖1 大鯢肉冷藏過程中菌落總數(shù)（A）和TVB-N值（B）的變化Fig.1 Changes in total bacterial count (A) and TVB-N content (B) of giant salamander meat during cold storage

2.2 PCR擴(kuò)增結(jié)果

如圖2所示，所有樣品的PCR擴(kuò)增產(chǎn)物條帶介于750～500 bp之間，特異性和亮度較好，無副帶和拖帶，說明各樣品純度較高，PCR擴(kuò)增條件合適，可用于后續(xù)的高通量測(cè)序分析。

圖2 PCR擴(kuò)增產(chǎn)物電泳圖Fig.2 Electrophoretograms of PCR-amplified products

2.3 測(cè)序樣本數(shù)據(jù)分析

高通量測(cè)序獲得不同冷藏期大鯢肉30 個(gè)樣品原始DNA序列，經(jīng)質(zhì)控和優(yōu)化后得到有效序列，相關(guān)信息見表1。可知，不同冷藏期樣品的有效序列數(shù)均在40 000以上，有效序列百分比達(dá)73%以上，平均長度范圍為412.55～425.01 bp，表明樣本的有效序列滿足后續(xù)微生物多樣性分析要求。

表1 樣品測(cè)序信息Table 1 Results of sample sequencing

將有效序列按照97%的相似度進(jìn)行OTU分類，并選取每個(gè)OTU中豐度最大的序列為該OTU的代表序列，繪制Venn圖，如圖3所示。0、2、4、6、8 d樣品中共有的OTU數(shù)分別為3 016、3 040、1 979、1 456、1 576，獨(dú)有的OTU數(shù)分別為715、721、119、24、37；而0 d-2 d、2 d-4 d、4 d-6 d、6 d-8 d中共有的OTU數(shù)分別為1 985、1 629、1 165、1 056。隨著冷藏時(shí)間的延長樣本中總OTU數(shù)呈降低的趨勢(shì)，且后三個(gè)時(shí)間段顯著降低，說明冷藏前期大鯢肉的微生物多樣性可能比冷藏后期更豐富，一些微生物只能在富氧的環(huán)境中生存。同時(shí)樣本中獨(dú)有的OTU與總OTU變化趨勢(shì)相似，但前后時(shí)間段之間共有OTU變化幅度較小，說明大鯢肉在冷藏過程中其菌群結(jié)構(gòu)變化很大，但主要菌群變化較小。根據(jù)OTU的特點(diǎn)能將5 個(gè)時(shí)間段劃分為3 個(gè)區(qū)間，即冷藏前期（0、2 d）、冷藏中期（4 d）、冷藏后期（6、8 d），各區(qū)之間變化顯著。

圖3 樣品中微生物群落共享和獨(dú)有的OTU Venn分析Fig.3 Venn analysis of shared and unique OTUs of microbial communities between samples

2.4 α多樣性分析

多樣性是根據(jù)97%相似度水平下的OTU信息，通過Chao1指數(shù)、Shannon指數(shù)、Simpson指數(shù)等對(duì)樣品微生物物種豐富度和多樣性進(jìn)行評(píng)估。Chao1指數(shù)反映樣品的菌群豐度，估計(jì)群落中含OTU的數(shù)目；Shannon指數(shù)和Simpson指數(shù)反映樣品中微生物群落的多樣性，Shannon指數(shù)越大，說明群落多樣性越高，Simpson指數(shù)則與之相反，值越大，說明群落多樣性越低。本研究得到的大鯢肉不同冷藏時(shí)間樣品多樣性相關(guān)指數(shù)如表2所示。不同冷藏時(shí)間樣本中微生物的覆蓋率均大于99%，表明測(cè)序結(jié)果可以反映實(shí)際情況。在0 d時(shí)Chao1指數(shù)最大，表明此時(shí)菌群豐度最高，隨著冷藏時(shí)間的延長，Chao1指數(shù)逐漸降低，說明微生物的豐富度逐漸降低，且冷藏中、后期較前期降低顯著（＜0.05）。隨著冷藏時(shí)間的延長，Shannon指數(shù)逐漸降低，Simpson指數(shù)則與之相反，表明大鯢肉的微生物菌群多樣性逐漸降低。冷藏前期氧含量充足、微生物生長繁殖所需的營養(yǎng)物質(zhì)較少、同時(shí)各種微生物存在生存適應(yīng)階段，因此，在冷藏前期微生物豐度和多樣性較高；隨著冷藏時(shí)間的延長，環(huán)境中氧含量降低、蛋白質(zhì)等大分子物質(zhì)在內(nèi)源酶和外源酶的作用下降解，為微生物的生長繁殖提供充足的營養(yǎng)物質(zhì)、優(yōu)勢(shì)物種大量繁殖并產(chǎn)生代謝物和毒素等，引起部分物種死亡導(dǎo)致冷藏后期微生物多樣性和豐富度降低，與Huang Wenbo等的研究結(jié)果相似。

表2 樣品微生物多樣性指數(shù)Table 2 Microbial diversity indexes of samples

圖4A為樣本物種累積曲線，用于描述隨著抽樣量的增大物種增加狀況，在生物多樣性和群落調(diào)查中，被廣泛用于抽樣量充分性的判斷以及物種豐富度估計(jì)，隨著抽樣量的增加，曲線趨于平緩，表明此環(huán)境中的物種并不會(huì)隨著樣本量的增加而顯著增多，抽樣充分。圖4B為樣本Rank Abundance曲線，用于同時(shí)解釋樣品多樣性的兩個(gè)方面，即樣品所含物種的豐富度（指一個(gè)群落或生境中物種數(shù)目的多寡）和均勻程度（指一個(gè)群落或生境中全部物種個(gè)體數(shù)目的分布均勻程度），物種豐富程度由曲線在橫坐標(biāo)上的長度反映，曲線越寬（橫軸的跨度越大），表示物種的組成越豐富；物種組成的均勻程度由曲線的形狀反應(yīng)，曲線越平坦（縱軸的跨度越大），表示物種組成的均勻程度越高。由圖4B可知，隨著OTU數(shù)量的增多，曲線逐漸平坦，表示樣本中物種組成均勻度逐漸提高，即冷藏前期較冷藏中、后期物種豐富度和多樣性高。

圖4 樣本物種累積曲線（A）和Rank Abundance曲線（B）Fig.4 Species accumulation curves (A) and Rank Abundance curves (B)

2.5 群落多樣性分析

2.5.1 樣品在門水平上的多樣性

如圖5所示，在門分類水平上，0、2 d中主要菌門為擬桿菌門（Bacteroidetes）、厚壁菌門（Firmicutes），相對(duì)豐度均在3 7%～4 2%之間；其次為變形菌門（Proteobacteria），相對(duì)豐度分別為14.2%、11.2%；Fibrobacteres、Actinobacteria、Spirochaetes相對(duì)豐度在3%～1%之間，其他菌門相對(duì)豐度均小于1%。4、6、8 d較0、2 d變化顯著，絕對(duì)優(yōu)勢(shì)菌門為變形菌門，相對(duì)豐度在84.7%以上，顯著高于冷藏前期，而擬桿菌門、厚壁菌門，相對(duì)豐度在7%～2%之間，顯著低于冷藏前期，其他菌門均小于1%。變形菌門中包含水產(chǎn)品常見的腐敗菌，如希瓦氏菌屬、假單胞菌屬?？偠灾谡麄€(gè)冷藏過程中，擬桿菌門、厚壁菌門、變形菌門為優(yōu)勢(shì)菌門，說明這3 種微生物對(duì)大鯢肉的品質(zhì)起著重要作用，從門水平上講，細(xì)菌的種類大致相同，而豐度差別比較大。

圖5 樣本中微生物在門水平相對(duì)豐度分布Fig.5 Relative abundance and distribution of microorganisms in samples at the phylum level

2.5.2 樣品中細(xì)菌在屬水平上的多樣性

為細(xì)化研究不同樣品中微生物多樣性的差異，進(jìn)一步在屬水平上分析大鯢肉冷藏過程中細(xì)菌群落多樣性。如圖6所示，0、2 d中主要菌屬為擬桿菌屬（），相對(duì)豐度分別為17.6%、19.1%，隨著時(shí)間的延長，相對(duì)豐度明顯下降，在4、6、8 d相對(duì)豐度分別下降至2.8%、0.7%、1.3%；其次為，相對(duì)豐度分別為8.5%、9.1%，在冷藏中后期豐度均低于1.5%，顯著下降；其余菌屬相對(duì)豐度均小于1%。而在4、6、8 d中假單胞菌屬（）為絕對(duì)優(yōu)勢(shì)菌屬，相對(duì)豐度分別為43.5%、56.4%、47.8%，較冷藏前期大幅增加，說明冷藏后期的條件更有利于其生長；其次為氣單胞菌屬（）、相對(duì)豐度在10%左右，沙雷氏菌屬（）相對(duì)豐度在5%左右，在4、6、8 d相對(duì)豐度分別為6.4%、1.3%、1.6%，這些菌屬較冷藏前期明顯增加；其余均屬相對(duì)豐度均小于1%。在屬水平上，不同冷藏期所含細(xì)菌的種類多樣性和豐富度存在較大的差異。此外，大鯢肉在冷藏初期，菌群種類極為豐富，隨著冷藏時(shí)間的延長，物種的種類逐漸減少。

圖6 樣本中微生物在屬水平相對(duì)豐度分布Fig.6 Relative abundance and distribution of microorganisms in samples at the genus level

2.6 樣本差異分析

基于由物種組成（某些距離矩陣算法也會(huì)考慮到物種進(jìn)化關(guān)系）計(jì)算得到的距離矩陣進(jìn)行主坐標(biāo)分析（principal coordinate analysis，PCoA）。如圖7A所示，PC1為54.19%，PC2為26.73%，PC1和PC2疊加解釋度達(dá)80.92%。在所有樣品中0、2、8 d組內(nèi)樣本之間距離較小，表明組內(nèi)微生物群落結(jié)構(gòu)差異較小，且0、2 d組內(nèi)微生物群落結(jié)構(gòu)相似度較8 d高；而4 d和6 d組內(nèi)樣品之間距離較大，表明組內(nèi)微生物結(jié)構(gòu)差異較大，主要是由于菌群在大鯢肉表面分布不均導(dǎo)致。整體看，0、2 d樣品之間相互重疊且集中分布在第2象限，8 d樣品集中分布在第4象限，4、6 d樣品分散在第1、4象限且6 d在第4象限的樣品更多，綜上表明，0-2、4、6、8 d之間微生物結(jié)構(gòu)差異較大。

圖7 樣本中微生物群落PCoA（A）和LEfSe分析差異物種注釋分支圖（B）Fig.7 PCoA analysis (A) and LEfSe analysis of differential species annotated clade (B) of microbial communities in samples

LEfSe分析可以分析組間菌群差異，找出組間差異的微生物種類。如圖7B所示，不同顏色表示不同分組，紅色節(jié)點(diǎn)表示在紅色組別中豐度相對(duì)較高的差異顯著物種，綠色節(jié)點(diǎn)表示在綠色組別中豐度相對(duì)較高的差異顯著物種，黃色節(jié)點(diǎn)表示在兩組比較中并無顯著差異的物種，節(jié)點(diǎn)直徑大小與相對(duì)豐度大小呈正比，每層節(jié)點(diǎn)由內(nèi)向外分別表示門、綱、目、科、屬。在微生物門水平上，5 組中差異顯著的種類包括5 個(gè)，分別是變形菌門、放線菌門（Actinobacteria）、擬桿菌門、纖維桿菌門（Fibrobacteres）、厚壁菌門。在微生物屬水平上，5 組樣本中顯著差異的菌屬共鑒定出16 個(gè)，冷藏前期假黃單胞菌屬（）、Rikenellaceae_RC9_gut_group、uncultured_bacterium、9、擬桿菌屬、纖維桿菌屬（）、毛螺菌屬（）、、顯著上調(diào)。冷藏中期沙雷氏菌屬（）、不動(dòng)桿菌屬（）顯著上調(diào)。冷藏后期氣單胞菌屬、假單胞菌屬、ambiguous_taxa、、1顯著上調(diào)。

2.7 累加統(tǒng)計(jì)分析

如圖8所示，A、D兩組變化趨勢(shì)相同，均是冷藏前期豐度較高，中、后期明顯下降，A組在5 個(gè)冷藏節(jié)點(diǎn)的值分別為46.32%、46.17%、7.22%、2.36%、2.99%，D組為51.60%、51.17%、12.72%、7.57%、10.54%。其中，擬桿菌屬、9、毛螺菌屬為專性厭氧菌，Rikenellaceae_RC9_gut_group為瘤胃微生物主要生活在厭氧環(huán)境中，該類菌屬，在冷藏初期，貯藏環(huán)境中氧含量較高，不利于生長繁殖，在冷藏中后期貯藏環(huán)境中氧含量較低，能進(jìn)行生長繁殖，但不具有生長優(yōu)勢(shì)，生長繁殖受到抑制，導(dǎo)致相對(duì)豐度明顯下降。B組在整個(gè)冷藏期間波動(dòng)相對(duì)較小，但在冷藏中期豐度較高，在5 個(gè)冷藏節(jié)點(diǎn)的值分別為0.46%、0.66%、13.49%、5.42%、6.46%。其中沙雷氏菌屬為革蘭氏陰性桿狀兼性厭氧菌，菌株能產(chǎn)生粉色、紅色或深紅色色素，對(duì)大鯢肉顏色有顯著影響，不動(dòng)桿菌屬為革蘭氏陰性好氧菌，屬于條件致病菌，是機(jī)體抵抗力降低時(shí)易感染的革蘭氏陰性菌，在冷藏中期相對(duì)豐度顯著增加，表明大鯢肌肉組織自身抗菌活性顯著下降，在冷藏后期快速下降，主要是由于假單胞菌屬的快速增殖競(jìng)爭營養(yǎng)成分和產(chǎn)生大量代謝物和毒素，冷藏環(huán)境中氧含量降低，對(duì)其繼續(xù)增殖具有抑制作用。

圖8 不同菌群累加占比分析Fig.8 Change in relative abundance of different bacterial groups during storage

C組在冷藏中、后期菌群豐度大幅增加，可以判定為主要致腐菌屬，在5 個(gè)冷藏節(jié)點(diǎn)的值分別為1.62%、2.00%、66.57%、84.65%、80.02%。其中假單胞菌屬為革蘭氏陰性需氧桿狀細(xì)菌，分解蛋白質(zhì)的能力極強(qiáng)，容易導(dǎo)致肉類變黏，代謝產(chǎn)物硫化氫、過氧化氫會(huì)使肉類變色，生長的最適pH值為7.0～8.5，能利用H和CO作為能源且能在冰箱中生長繁殖，不能在pH值6或6以下生長，是低溫貯藏肉制品中常見的腐敗菌，其致腐性受群體感應(yīng)系統(tǒng)的調(diào)控。結(jié)合圖6可知，該菌在冷藏中后期相對(duì)豐度顯著增加至50%左右，為主要致腐菌屬。在大鯢肉發(fā)生自溶生化變化前期，蛋白質(zhì)等大分子降解較少，pH值較低，微生物種類較多，該菌的生長繁殖速度相對(duì)較慢；在自溶階段及之后，蛋白質(zhì)等大分子在酶的作用下分解為小分子，貯藏環(huán)境中氧含量降低，部分菌種因不適應(yīng)環(huán)境而死亡，pH值逐漸升高，有利于該菌的生長繁殖，導(dǎo)致冷藏中后期大幅增加；但在第8天相對(duì)豐度較第6天略微下降，主要是由于微生物大量生長繁殖，導(dǎo)致代謝物和毒素積累，對(duì)自身生長具有抑制作用。氣單胞菌屬能代謝產(chǎn)生氧化三甲胺和硫化氫，與假單胞菌屬類似，是大鯢肉在冷藏過程中變色、發(fā)黏并產(chǎn)生不良?xì)馕兜脑蛑唬?、ambiguous_taxa、1在冷藏中后期相對(duì)豐度增加，其中相對(duì)豐度值在中后期達(dá)到10%左右，但其致腐性有待研究。

2.8 進(jìn)化分析

系統(tǒng)發(fā)育學(xué)也稱系統(tǒng)發(fā)生學(xué)，通過某一分類水平上序列間堿基的差異構(gòu)建進(jìn)化樹能夠揭示出有關(guān)生物進(jìn)化過程的順序，了解生物進(jìn)化歷史和機(jī)制，是研究物種形成和進(jìn)化關(guān)系的有力工具。使用FastTree軟件繪制組合圖，選擇豐度在前50的OTU，根據(jù)最大似然法構(gòu)建進(jìn)化樹，并以熱圖展示OTU在不同樣本的豐度。由圖9可知，進(jìn)化樹可以分為5 個(gè)大的分支，第1個(gè)分支為變形菌門中的部分屬，主要為假單胞菌屬，絕大部分菌屬在冷藏中、后期豐度較高，是影響大鯢肉品質(zhì)的關(guān)鍵菌群。第2個(gè)分支主要包括了厚壁菌門中的部分菌屬，第3個(gè)分支主要包括梭桿菌門中的梭菌屬，第4個(gè)分支主要包括了纖維桿菌門中的纖維桿菌屬，第5個(gè)分支包括擬桿菌門中的部分屬，主要為擬桿菌屬。后4個(gè)分支除1，其余菌屬均在冷藏前期豐度較高。通過進(jìn)化分析可知，菌落演替與冷藏時(shí)間具有良好的相關(guān)性。

圖9 TOP 50物種進(jìn)化樹及OTU豐度組合圖Fig.9 Phylogenetic tree and OTU abundance combinations of TOP 50 species

2.9 基于16S的KEGG功能預(yù)測(cè)

京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）是有關(guān)生物系統(tǒng)比較完善的數(shù)據(jù)庫，關(guān)聯(lián)基因組信息和功能信息的知識(shí)庫。其由基因蛋白序列（KEGG Genes）、具有內(nèi)源性和外源性的化學(xué)物質(zhì)（KEGG Ligand）、分子相互作用和代謝通路圖（KEGG Pathway）、各種生物之間的層次關(guān)系（KEGG Brite）構(gòu)成。KEGG在L2水平共41 種代謝通路，其中相對(duì)豐度值占比居前15的代謝通路如表3所示?？梢钥闯?，膜運(yùn)輸、特征差、細(xì)胞過程和信號(hào)傳導(dǎo)、脂代謝、新陳代謝在冷藏中后期表現(xiàn)為增加的趨勢(shì)。其中膜運(yùn)輸?shù)南鄬?duì)豐度最大，達(dá)到13%左右，膜運(yùn)輸包括ABC轉(zhuǎn)運(yùn)器、磷酸轉(zhuǎn)移酶系統(tǒng)、細(xì)菌分泌系統(tǒng)3 條代謝通路，其中，磷酸轉(zhuǎn)移酶系統(tǒng)是細(xì)菌吸收碳水化合物的主要機(jī)制；ABC轉(zhuǎn)運(yùn)器廣泛存在于細(xì)菌、古細(xì)菌和真核生物中，主要是將ATP水解與多種底物主動(dòng)轉(zhuǎn)運(yùn)結(jié)合起來；細(xì)菌分泌系統(tǒng)在革蘭氏陰性細(xì)菌中具有高度的保守性。由3 條代謝通路的特性可以看出，膜運(yùn)輸在冷藏中后期豐度的升高，與革蘭氏陰性菌的增加密切相關(guān)。

表3 KEGG功能預(yù)測(cè)Table 3 KEGG function prediction

碳水化合物代謝、復(fù)制和修復(fù)、能量代謝、翻譯、輔助因子和維生素代謝、核苷酸代謝、轉(zhuǎn)錄、遺傳信息處理、折疊、分類和降解隨著冷藏時(shí)間的延長均表現(xiàn)為降低的趨勢(shì)。其中碳水化合物代謝、能量代謝、輔助因子和維生素代謝等能為細(xì)菌的生長繁殖提供必需的營養(yǎng)物質(zhì)；復(fù)制與修復(fù)、轉(zhuǎn)錄、遺傳信息處理、折疊等是微生物自身合成遺傳物質(zhì)的主要代謝通路；該類代謝通路相對(duì)豐度的變化趨勢(shì)與圖8中主要致腐菌的變化趨勢(shì)相反，表明致腐菌的大量增殖和微生物群落結(jié)構(gòu)的變化對(duì)代謝通路環(huán)境有明顯的影響，弓湃等研究指出假單胞菌屬中CbrAB-CrcZ-Crc參與碳代謝抑制調(diào)控，能為假單胞菌適應(yīng)環(huán)境變化和提高在自然界中的競(jìng)爭力做出代謝優(yōu)化調(diào)控，賦予假單胞菌相比同一環(huán)境中其他微生物而言更多的競(jìng)爭優(yōu)勢(shì)，更有利于假單胞菌在環(huán)境中生存。

氨基酸代謝通路在整個(gè)冷藏過程中波動(dòng)較小，氨基酸代謝在微生物演替中的作用主要包括兩個(gè)方面，一方面是氨基酸的合成代謝，用于合成自身獨(dú)特的蛋白質(zhì)、多肽和其他含氮物質(zhì)，以產(chǎn)生相應(yīng)的氮源為微生物利用；另一方面是氨基酸的分解代謝，通過脫氨基、轉(zhuǎn)氨基、聯(lián)合脫氨基的方式分解為-酮酸、胺類和二氧化碳。其中-酮酸可合成非必需氨基酸、轉(zhuǎn)變?yōu)樘羌磅ヮ?、氧化產(chǎn)生能量；二氧化碳大部分直接排到細(xì)胞外。同時(shí)能量代謝中糖酵解和三羧酸循環(huán)產(chǎn)生二氧化碳都會(huì)導(dǎo)致環(huán)境中二氧化碳濃度的升高，有利于兼性厭氧菌屬的生長繁殖，是冷藏中后期乳球菌屬（）、沙雷氏菌屬、拉恩氏菌屬（）、檸檬酸桿菌屬、希瓦氏菌屬等豐度明顯增高的主要原因。

3 結(jié)論

通過Illumina MiSeq測(cè)序技術(shù)對(duì)托盤包裝大鯢肉冷藏過程中的5 個(gè)時(shí)期30 個(gè)樣品進(jìn)行分析，初步明確了大鯢肉冷藏期間菌相組成及變化規(guī)律，主要結(jié)論如下：1）大鯢肉在冷藏過程中微生物豐度隨著冷藏時(shí)間的延長逐漸降低；門水平分析發(fā)現(xiàn)，細(xì)菌優(yōu)勢(shì)菌門從冷藏前期（0、2 d）的擬桿菌門、厚壁菌門逐漸轉(zhuǎn)變?yōu)橹校? d）、后期（6、8 d）的變形菌門；屬水平分析表明，細(xì)菌優(yōu)勢(shì)菌屬從冷藏前期的擬桿菌屬、逐漸轉(zhuǎn)變?yōu)橹?、后期的假單胞菌屬、氣單胞菌屬、、沙雷氏菌屬?）PCoA顯示0-2、4、6、8 d之間微生物結(jié)構(gòu)差異較大，兩個(gè)PC疊加解釋度達(dá)80.92%；LEfSe分析表明，引起大鯢肉各冷藏期差異顯著的菌門主要為變形菌門、放線菌門、擬桿菌門、纖維桿菌門、厚壁菌門，菌屬主要為假黃單胞菌屬、、沙雷氏菌屬、不動(dòng)桿菌屬、氣單胞菌屬、假單胞菌屬、ambiguous_taxa、、Rikenellaceae_RC9_gut_group、_9、擬桿菌屬、纖維桿菌屬、毛螺菌屬、、_1。3）進(jìn)化分析表明大鯢肉冷藏過程中微生物菌群的演替與冷藏時(shí)間具有較強(qiáng)相關(guān)性。綜合分析，引起冷藏過程中托盤包裝大鯢肉腐敗變質(zhì)的微生物可能為假單胞菌屬、氣單胞菌屬、和沙雷氏菌屬。該研究為今后大鯢肉冷藏過程中靶向抑菌保鮮及貨架期延長提供參考。