基于宏轉錄組學解析不同貯藏時期冷鮮灘羊肉微生物代謝過程

趙珺怡，楊 波，羅瑞明,*，張赫宇，蘇春霞，胡倩倩，趙曉策

（1.寧夏大學農學院，寧夏 銀川 750021；2.寧夏大學生命科學學院，寧夏 銀川 750021）

冷鮮灘羊肉是生鮮灘羊肉的一種重要產品，在常規0～4 ℃冷藏過程中，冷鮮灘羊肉會被外源性腐敗微生物[1-2]污染，從而導致肉的腐敗變質，甚至爆發食源性疾病乃至食物中毒[3]。有研究表明，微生物的生長繁殖不僅取決于貯藏環境中氧氣含量、溫度高低，也與生長基質的pH值有關[4]。不同屠宰、分割工序[5]與不同包裝形式[6]所含有的微生物有一定的差異。目前，針對冷鮮肉微生物的研究熱點多集中于采用選擇性培養基、聚合酶鏈式反應（polymerase chain reaction，PCR）技術、酶聯免疫法等分析冷鮮肉的生物多樣性，通過致腐微生物生長值與模型參數擬合，結合肉品質感官指標與理化指標，建立貨架期預測模型等[7-11]。但對于冷鮮肉微生物與微環境之間的分子生態關系與群落結構動態演替的分子機制卻鮮有報道。

隨著后基因組時代的來臨，組學技術在食品微生物領域也得到廣泛關注。Dugat-Bony等[12]應用宏基因組和宏轉錄組學技術研究了由9 種微生物組成的奶酪表皮在4 周成熟期中的變化。通過對其轉錄本表達量的分析，發現氨基酸代謝與脂肪代謝途徑是最顯著的代謝途徑，相關差異基因源于假絲地霉酵母、干酪棒桿菌和蜂房哈夫尼亞菌。Jeong等[13]運用宏轉錄組學技術研究了韓國泡菜發酵水中的微生物群落演替，發現發酵泡菜水中初期微生物群落中含有明串珠菌、乳酸菌、假單胞菌屬和泛菌屬，而發酵第3天直至發酵結束明串珠菌占主導地位。在冷鮮肉微生物領域，張赫宇等[14]采用高通量測序技術對不同冷藏階段灘羊肉微生物的16S rRNA基因V3區進行測序，研究發現假單胞菌、腸桿菌、熱死環絲菌和乳酸菌為冷藏過程中的優勢細菌，而優勢真菌主要為酵母菌。張鐵華等[15]研究了冷鮮牛肉微生物菌相變化，結果表明在低溫貯藏過程中，G－小桿菌系是導致無包裝鮮牛肉腐敗變質的主要表面菌系，而G+球菌系是導致真空包裝鮮牛肉腐敗變質的主要菌系。

宏轉錄組學興起于宏基因組之后，以微生物群落的總RNA為研究對象[16]，探索特定時期下群體全部基因組轉錄情況及轉錄規律。宏轉錄組學技術能從轉錄水平研究復雜微生物群落變化，能更好的挖掘潛在的新基因和差異表達基因功能，探究各微生物成員對整體環境的貢獻以及微生物功能與功能之間相互影響調控的作用機制。

本研究以冷鮮灘羊肉表面微生物為研究對象，采用宏轉錄組學技術對0、4、8 d微生物基因組差異轉錄、代謝途徑與3 個貯藏時期微生物群落演替的關系進行探究，通過比較3 個貯藏時期的變化推測微生物優勢菌優勢生長的原因，為冷鮮肉保鮮分子生態學理論構建打下基礎，為冷鮮灘羊肉的貯藏加工提供理論指導。

1 材料與方法

1.1 材料與試劑

灘羊后腿肉（當天屠宰），取自寧夏鹽池縣大夏牧場食品有限公司（灘羊飼養期間不做任何生物性防疫工作，且定期進行清糞、消毒處理），托盤密封包裝。

TRIzol試劑盒 美國ThermoFisher Scientific公司；GeneReadTM試劑盒 德國Qiagen公司；三磷酸脫氧核糖核苷酸（dNTPs）、fragmentation buffer、核糖核酸酶H、DNA聚合酶I 美國TaKaRa Bio公司；AMPure XP純化磁珠 英國Beckman Coultier公司。

1.2 儀器與設備

SW-CJ-1G超凈工作臺 上海博訊實業公司；Sorvall Legend Micro 17常溫離心機、Forma 994 -80 ℃超低溫冰箱 美國Thermo公司；HiSeq高通量測序儀 中國Novogene公司；Tanon-2500凝膠成像系統 中國天能公司；vortex-genie2 G-560E旋渦混合儀 美國SI公司；Qubit 2.0熒光計 上海Life Technologies公司；ScanSpeed MiniVac Alpha冷凍離心機 美國Scan Speed公司；GeneAmp PCR system 9700普通PCR儀 美國Life公司；Implen GmbH納米光度計 德國Schatzbogen公司。

1.3 方法

1.3.1 樣品采集

選用寧夏鹽池縣大夏牧區的灘羊作為實驗材料，在0～4 ℃條件下貯藏的冷鮮灘羊肉上采集樣品。采集樣品前，用體積分數70%乙醇溶液對冷凍收集管、鑷子、刀等取樣器進行消毒，干燥后密封于無菌密封袋。用消毒刀分別在貯藏時間為0、4 d和8 d共9 個灘羊腿（每個貯藏時間段取3 個羊腿）上進行刮擦，得到刮屑。采集的樣本包括肌肉、結締組織、血液和脂肪等。

在第0天時，快速刮取3 個羊腿上的同一位置獲得3 個平行性樣本（≥150 mg），分別放入3 個冷凍管并標記為C0_1、C0_2和C0_3；然后將其表面刮擦物含菌落、血液、肉汁液等放置在-80 ℃的液氮中冷凍1 h。冷鮮貯藏4 d的3 個羊腿與0 d取樣相同位置以同樣方式獲得樣本，分別標記為C4_1、C4_2、C4_3；冷鮮貯藏8 d的3 個羊腿也以同樣的方式取樣分別標記為C8_1、C8_2和C8_3。將所有樣品放入冰箱（0～4 ℃）中冷藏。每組樣本分別進行3 個生物學重復實驗。

1.3.2 RNA的提取

將上述樣本封裝于滅菌容器，在無菌條件玻璃球振蕩分散，制備蒸餾水懸浮液，質量濃度1 g/100 mL。取1.0×109（1 mL菌液OD600nm為1～1.5）的細菌培養物放入2 mL離心管中，室溫12 000×g離心30 s，除去上清液。使用TRIzol試劑盒，按照試劑盒說明分別提取樣本的總RNA，通過瓊脂糖凝膠電泳檢測RNA降解情況（圖1）。確保RNA的純度和濃度符合后續實驗要求，即RNA完整值（RNA integrity number，RIN）大于等于8.0。

圖1 瓊脂糖凝膠電泳的總RNA質量分析Fig.1 Agarose gel electrophoresis for quality analysis of total bacterial RNA

1.3.3 cDNA文庫的構建及Illumina HiSeq測序

使用評估污染物試劑盒對1.3.2節所得RNA進行評估，確定樣品未發生污染后，使用GeneReadTM試劑盒除去rRNA。并加入碎片緩沖液使mRNA分解成小片段。

使用隨機六聚體引物以上述mRNA為模板合成單鏈cDNA，使用TaKaRa試劑盒以單鏈cDNA為模板合成雙鏈cDNA。使用AMPure XP純化磁珠純化雙鏈cDNA后，加入A尾，修復末端。通過PCR擴增修飾的雙鏈cDNA，用AMPure XP磁珠濃縮回收PCR產物，構建基因文庫。

使用Qubit2.0熒光計初步定量每個文庫。將每個文庫質量濃度稀釋至2 ng/μL。使用Q-PCR方法確定文庫的有效濃度（確定為＞2 nmol/L）。

使用Agilent Bioanalyzer 2100系統檢測并確定DNA插入物的大小是否符合預期。在Illumina HiSeq平臺進行高通量測序。

1.3.4 數據預處理

測序產生的原始數據存在一定比例低質量數據，為了保證后續分析的結果準確可靠，首先對原始的測序數據進行預處理，獲取用于后續分析的有效數據。預處理方法參見方法。處理步驟如下：1）去除質量值≤5的堿基數達到一定比例的reads（默認reads長度的40%，設置為40）；2）去除含N的堿基數目達到一定比例的reads（默認reads長度的10%，設置為10）；3）去除Adapter污染（默認Adapter序列與reads序列有15 bp的overlap，設置為15）；4）在有宿主污染可能性的前提下，需與宿主數據庫進行比對，過濾掉可能是宿主污染的reads（默認設置比對一致性大于等于90%的reads為宿主污染）。

1.3.5 De novo組裝

針對每個樣品經預處理得到的clean reads，先使用NCBI的rRNA、tRNA以及SILVA數據庫進行比對分離出來宏基因組中rRNA序列，剩下的mRNA序列則使用拼接軟件Trinity（version: r20140413p1）分別進行從頭組裝，然后對所有樣品的序列整合并使用CD-HIT-EST去冗余（設定序列一致性閾值為0.95），得到unigene集合。

1.3.6 生物信息學分析

將得到的unigene集合進行一系列差異表達基因分析：差異基因的基因本體（gene ontology，GO）注釋、差異基因京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）注釋、GO顯著性富集分析及KEGG顯著性富集分析。

2 結果與分析

2.1 數據預處理結果

表1 數據預處理統計Table 1 Preprocessing results of metatranscriptomic data

限定序列質量后，共獲得有效數據57.02 GB。由表1可知，每個樣品的Q30堿基的占比不小于91.48%，此結果表明序列數據可靠，可用于下一步分析。此外，整個基因組的GC含量（鳥嘌呤和胞嘧啶的比例）在48.9%到55.15%之間，這表明在基因組測序過程中，覆蓋的整個區域沒有偏差，保證了后續拼接和注釋的可靠性。

2.2 De novo組裝結果

單基因長度區間分布在301～500 bp之間，單基因總數達77 709 個，占所有單基因長度區間的60.25%。統計結果如圖2所示，4 746 個單基因長度大于2 000 bp，占所有單基因長度區間的0.04%，平均長度643 bp，連續N50長度為672 bp。164 836 個轉錄序列來自序列重疊組，平均長度733 bp，連續N50長度為871 bp。單基因長度主要分布在300～2 000 bp之間，占所有單基因長度區間的99.94%，8 829 個基因長度大于2 000 bp，占所有單基因長度區間的0.06%。總之，單基因數目豐富，轉錄長度分布合理，可以滿足后續實驗的需要。

圖2 單基因序列和轉錄序列長度分布圖Fig.2 Unigene and transcript length distribution

2.3 物種注釋結果

通過與Nr（NCBI non-redundant protein sequences）庫進行BLAST比對（e-value≤10-5），由于每一條序列可能會有多個比對結果，得到多個不同的分類級別，為了保證其生物意義，采取LCA算法（應用MEGAN軟件的系統分類）[17]，將出現第一個分支前的分類級別，作為該序列的物種注釋信息。從門、屬和種水平上的相對豐度表出發，選取出在各樣品中的最大相對豐度排名前10的門類，并將其余的物種設置為其他，繪制出各樣品對應的物種注釋結果在其水平的統計圖，如圖3所示。

圖3 物種注釋在門（a）、屬（b）和種（c）水平的結果Fig.3 Results of species annotation at the levels of phylum (a),genus (b) and species (c)

根據所有樣品在各水平的物種注釋及豐度信息，選取豐度排名前35 名的屬及其在每個樣品中的豐度信息繪制熱圖，并從分類信息和樣品間差異兩個層面進行聚類，找出研究樣品中聚集較多的物種或樣品，結果展示見圖4。

由圖3與圖4所知，Proteus hauseri（變形桿菌）、Macrococcus caseolyticus（巨大球菌屬）、Pseudomonas（假單胞菌）、Aspergillus flavus（黃曲霉）、Terrabacter（腫大地桿菌屬）、Aeromonas hydrophila（嗜水氣單胞菌）在C4_1樣品中富集，C4_3樣品中的優勢菌株是Trichodesmium erythraeum（紅色毛癬菌）。C8_1樣品中的優勢菌是Lactococcus lactis（乳酸菌），而在C8_2樣品中占主導地位是Gammaproteobacteria bacterium（伽馬蛋白細菌）和L. lactis。在4 d的樣品中發現了3 組在不同樣品中具有相似豐度變化的生物，第1組是P. hauseri和M. caseolyticus；第2組是A. flavus；第3組是Terrabactersp.28和A. hydrophila（嗜水氣單胞菌）。由此可以看出，貯藏至第8天時，冷鮮灘灘羊肉的主要優勢微生物為變形菌門、乳酸菌屬等。其中，變形菌門中的假單胞菌和腸桿菌是造成冷鮮肉腐敗變質的主要微生物，這與已有研究結果一致[18-19]。

灘羊肉冷藏第8天的優勢菌株為G. bacterium（變形菌）和L. lactis（乳酸菌）。在灘羊肉冷藏初期，革蘭氏陰性細菌如P. hauseri、N. spumigena、M. caseolyticus、G. prolifera、A. flavus、A. phagocytophilum、Terrabactersp. 28、A. hydrophila、G. bacterium占據主導地位。

圖4 物種豐度聚類圖在門（a）、屬（b）和種（c）水平的結果Fig.4 Clustering maps showing species abundance at the levels of phylum (a), genus (b) and species (c)

假單胞菌是好氧微生物，生長幾乎不受其他微生物影響[20]。在灘羊肉的冷藏過程中，假單胞菌優先利用貯藏環境中的氧氣，當氧氣消耗殆盡時，灘羊肉成熟過程中產生的葡萄糖成為其生長所需的能量來源[21]。假單胞菌在繁殖和新陳代謝的過程中會產生具有特殊氣味的硫化物，酯、酸和其他變質代謝產物[22]會導致灘羊肉的腐敗變質。與假單胞菌相比，腸桿菌的增長速度較緩慢。腸桿菌為兼性厭氧菌，可以發酵葡萄糖并產酸產氣。在冷藏環境下，灘羊肉中的葡萄糖與中間代謝產物6-磷酸葡萄糖作為碳源被腸桿菌利用，產生胺類代謝物質導致肉的腐敗[23-24]。隨著貯藏時間的延長，灘羊肉細胞逐漸衰老死亡，可提供的能量和氧氣逐漸變少，革蘭氏陽性菌特別是乳酸菌（L. lactis）逐漸成為優勢微生物。乳酸菌分泌一種名為nisin的抗菌肽，該物質能與脂II分子結合在革蘭氏陰性菌細胞膜上形成孔洞，致其細胞死亡[25]，從而確保自己的優勢地位。

2.4 差異基因的篩選結果

通過對0～4、4～8、0～8 d的冷鮮灘羊肉中微生物群落的基因表達情況檢測，采用基于負二項分布的DESeq[26]進行分析。統計顯著性（signature，S）以q-value小于0.005和|log2（fold change）|大于1為篩選條件判斷基因是否為差異基因[27]，TRUE：為差異基因，FALSE意為不是差異基因。將篩選出來的差異表達基因用火山圖展示。橫坐標代表基因在不同實驗組中表達倍數的變化；縱坐標代表基因表達量變化的統計學顯著程度，q-value越小，-lg（q-value）越大，則差異越顯著，結果如圖5所示。

圖5 差異表達基因的火山圖Fig.5 Volcano plots of differentially expressed genes

篩選結果顯示，從貯藏時間為4 d與0 d、8 d與4 d、8 d與0 d的冷鮮灘羊肉中測得的差異表達基因數分別為429、15、1 529 個。在C4與C0中，上調基因為247 個，下調基因為182 個；在C8與C4中，上調基因為2 個，下調基因為13 個；在C8與C0中，上調基因為727 個，下調基因為802 個。上調基因分別占基因總數的42.42%、13.33%和47.54%。下調基因分別占基因總數的57.58%、86.67%和52.46%。

2.5 差異基因的GO富集分析

圖6 差異基因GO顯著富集圖Fig.6 Significantly enriched GO terms of differentially expressed genes

采用非中央超幾何分布的GOseq[28]方法將P＜0.05作為閾值篩選GO條目，如圖6所示，具體地顯示了擁有較多基因富集的GO富集分析的條目。其中，在C0與C4對比組內的60 個分子功能、74 個生物過程和27 個細胞組成分別于329、838、277 個unigenes有關；在C4與C8對比組內8 個分子功能、35 個生物過程和17 個細胞組成分別與10、55、18 個unigenes有關；在C0與C8對比組內的70 個分子功能、121 個生物過程和57 個細胞組成分別與981、1 811、2 860 個unigenes相關。除去細胞結構的差異基因，C0與C4對比組，富集差異基因較多的分子功能是受體、蛋白質、GTP和核苷酸的結合，肽酶的活性以及結構分子的活性；富集差異基因較多的生物過程有信號轉導、蛋白質聚合與水解和光合作用等；C4與C8對比組，富集差異基因較多的生物過程有離子轉運和線粒體轉運；C0與C8對比組中，富集差異基因較多的分子功能是不同氨基酸酶的活性，富集差異基因較多的生物過程主要是核苷酸代謝、蛋白質水解與組裝、蛋白質代謝、磷酸化作用和離子轉運等。

2.6 差異表達基因的EggNOG/COG注釋

通過EggNOG分析將25 類unigenes簇進行分類，如圖7所示。其中包含最多基因的功能基因簇分別為R—一般功能預測（8 150 個匹配基因）、S—功能未知（7 900 個基因）和T—信號轉導機制（6 000 個基因）。此外，同源基因簇內的基因具有相似的結構和生物學功能，被認為來自同一祖先。而E—氨基酸轉運和代謝、F—核苷酸轉運和代謝、G—碳水化合物的運輸和代謝、H—輔酶運輸和代謝與I—脂質運輸和代謝也占有一定份額，這與GO富集分析得到的結果較為一致。

圖7 差異基因單基因EggNOG 注釋結果的統計圖Fig.7 EggNOG annotation results of differentially expressed unigenes

2.7 差異基因的KEGG富集分析

在生物體內，不同基因相互協調行使其生物學功能，通過Pathway顯著性富集能確定差異表達基因參與的最主要生化代謝途徑和信號轉導途徑。KEGG是有關Pathway的主要公共數據庫[29]。將生物代謝通路劃分為6 類，分別為：細胞過程、環境信息處理、遺傳信息處理、人類疾病、新陳代謝、生物體系統，將3 個不同貯藏時期的微生物所富集的差異基因進行KEGG注釋，如圖8所示。

由圖8可知，富集基因較多的代謝通路有核苷酸代謝、脂質代謝、能量代謝、碳水化合物代謝以及氨基酸代謝等。作為微生物的基礎代謝，微生物普遍存在的核苷酸水解酶對核苷酸代謝分解具有促進作用。核苷酸水解酶先將核苷酸水解，生成核苷，再通過酶的催化作用下生成堿基和磷酸核糖[30]，磷酸核糖參與磷酸戊糖途徑生成6-磷酸果糖并且參與糖酵解，最終產物丙酮酸進入三羧酸循環途徑進行反應，將能量代謝與物質代謝聯系在一起，增加了微生物可利用的碳源，增強了微生物的增殖能力。

圖8 不同貯藏時期差異基因表達的KEGG Pathway富集結果Fig.8 KEGG pathway enrichment analysis of differentially expressed genes at different chilled storage period

氨基酸代謝在微生物演替中的作用主要包括2 個方面，一方面主要用于合成自身獨特的蛋白質、多肽和其他含氮物質，以產生相應的氮源為微生物利用；另一方面可通過脫氨、轉氨、聯合脫氨的方式分解為α-酮酸、胺類及二氧化碳[31]。其中，α-酮酸可通過胺化反應轉化為糖、脂質或非必需氨基酸，也可參與三羧酸循環完全氧化成二氧化碳和水，并釋放出大量能量[32]。乳酸菌以游離氨基酸作為底物，通過氨基酸的脫羧反應產生相應的堿性生物胺，維持自身的生長繁殖[33]。且氨基酸代謝產生的大量二氧化碳會降低貯藏環境的氧氣濃度，正好滿足了兼性厭氧的乳酸菌的生長條件，使得乳酸菌成為冷鮮灘羊肉貯藏期內的優勢微生物之一。

能量代謝屬于新陳代謝的一種，碳水化合物代謝也是能量代謝的過程。碳水化合物代謝過程中會形成磷酸鹽化合物，在連接磷酸鹽基團和其他分子的化學鍵中存儲大量能量[34]。隨著能量代謝與碳水化合物代謝的顯著富集表明冷鮮灘羊肉表面微生物的代謝活躍，貯藏后期的微生物種類和數量相比于原始菌相產生了較大的增長，與能量代謝有關的差異基因表達富集較為顯著。經過不同時間的冷鮮貯藏，灘羊肉表面微生物也呈現出一定的差異，但是就總體3 個貯藏時間段而言，因為原始菌相存在的基礎使微生物群落結構的相似性更大，差異性則相對較小，而差異性主要是由于冷鮮肉微生物的生長代謝導致的。

在0～4 ℃的貯藏條件下，微生物的生命活動顯著降低，一些嗜冷菌可以正常增殖，假單胞菌和腸桿菌都可以在這個溫度條件下進行生長繁殖。宰前活體灘羊肉的pH值為中性，但在宰后成熟階段，灘羊肉肌肉組織中的糖原在糖酵解酶作用下生成乳酸。同時，呼吸作用并沒有停止，產生了大量能量和磷酸，在乳酸和磷酸的共同作用下，灘羊肉的pH值下降，造成胴體僵硬。此時鈣激活酶作用于僵硬的胴體使其自溶，pH值降低至5.6～6.0[35]，這種酸性環境最適合乳酸菌的生長。在貯藏過程中，冷鮮灘羊肉中的脂肪成分被附著在表面的優勢微生物脂肪酶催化水解成脂肪酸和甘油[36]，由甘油氧化形成的三碳醇酸是絲氨酸降解的中間產物，經磷酸化后生成3-磷酸甘油酯，可進一步參與糖異生或糖酵解反應。1-磷酸葡萄糖是糖異生途徑產物，在儲藏前期，假單胞菌和腸桿菌先利用氧氣進行生長繁殖，在貯藏后期時，假單胞菌和腸桿菌則利用之前反應產生的葡萄糖為底物繼續繁殖生長，成為主要優勢微生物。而在貯藏中期豐度較高的不動桿菌屬為專性需氧菌，在氧氣含量較少的貯藏后期大幅減少；微小桿菌屬和放線菌屬為非嗜冷菌，在低溫貯藏環境下的生長繁殖相對收到抑制，且乳酸菌在生長代謝中產生的乳酸也會抑制這些微生物的生長繁殖，使它們隨著貯藏時間的延長，逐漸變為劣勢微生物。

3 結 論

為系統揭示冷鮮肉微生物群落的變化，本研究利用宏轉錄組學技術對0、4、8 d三個冷鮮貯藏時間段灘羊肉表面微生物群落的結構進行了解析，共獲得57.02 GB的轉錄本有效數據。研究表明在3 個貯藏時期中，C4與C0組的上調基因為247 個，下調基因為182 個；C8與C4組的上調基因為2 個，下調基因為13 個；而在C8與C0組中，上調基因和下調基因分別為727 個和802 個。這些差異表達基因主要富集在核苷酸代謝、脂質代謝、能量代謝、碳水化合物代謝以及氨基酸代謝等代謝途徑。通過利用這些代謝途徑所產生的代謝物和能量作為微生物生長繁殖的底物，假單胞菌和腸桿菌成為貯藏期間冷鮮灘羊肉表面的優勢微生物。在貯藏后期，由于貯藏環境中殘存的少量氧氣與灘羊肉較低的pH值，乳酸菌豐度漲幅較大，也成為優勢微生物之一。本研究對貯藏過程中的冷鮮灘羊肉微生物進行分析，可為冷鮮灘羊肉貯藏期間微生物及肉品質量預報、控制技術研發提供理論參考。