風干肉中產脂肪酶瑞士乳桿菌TR13全基因組測序及序列分析

田建軍，張開屏，趙艷紅，李權威，馬牧然，馬俊杰，曹凱慧，靳 燁,*

（1.內蒙古農業大學食品科學與工程學院，內蒙古 呼和浩特 010018；2.內蒙古商貿職業學院食品工程系，內蒙古 呼和浩特 010070）

乳酸菌是一種在自然界廣泛存在的無芽孢、厭氧或兼性厭氧、主要用于發酵各種生食和飼料、代謝產物主要以乳酸為主的革蘭氏陽性細菌[1-3]。由于其在發酵食品生產中的作用和作為益生菌的用途，是食品微生物學研究中最重要的細菌之一[4-5]。瑞士乳桿菌（Lactobacillus helveticus）是具有良好益生特性的一種乳酸菌，在發酵牛奶的過程中可產生抗高血壓的肽[6]，已被廣泛用于發酵和功能性食品當中[7-10]。

傳統風干肉是以牛、羊、馬等反芻動物新鮮肉為主要原料，在適宜溫度和通風良好的條件下，將牛羊肉割成小條，掛在陰涼處，自然晾曬、風干，去除肉中部分水分而形成的風味獨特的肉制品，它屬于自然條件下，依靠內源酶和有益微生物的作用，發生一系列生化變化而制成的一類發酵肉制品[11]。發酵肉制品中含有豐富的微生物資源，且具有高度的微生物群落多樣性[12]。乳酸菌是最早從發酵肉制品中分離出來的微生物，是傳統發酵肉制品中重要的優勢菌群，與產品品質密切相關[13-14]。近幾年也有一些從發酵肉制品中篩選肉用乳酸菌發酵劑的相關報道[15-17]。脂肪酶能夠使發酵肉制品產生獨特的風味。Fernández等[18]發現將脂肪酶添加到中式香腸中能夠有效加速香腸中脂肪降解，促進香腸中脂質來源的揮發性風味物質的生成、減少成熟時間。Zalacain等[19]在發酵香腸中添加外源脂酶，發現脂 酶可促進飽和脂肪酸釋放，改變游離脂肪酸組成且不飽和脂肪酸含量顯著增加。Horchani等[20]將1 株產脂肪酶的葡萄球菌用作香腸發酵劑以此來改善風味品質。Lorenzo等[21]把產脂肪酶的4 種不同商業發酵劑添加到發酵香腸中，對其游離脂肪酸含量和種類進行分，發現加工過程中脂肪酸有顯著的逐漸釋放的現象，進一步研究發現，產脂肪酶微生物菌株不僅僅可改善風味成分的變化，對其肉中油脂的特性也產生了一定的影響。

全基因組測序技術可以將細菌DNA（完整基因組序列）全部檢測出來，便于揭示生物體的完整DNA組成，通過序列分析及其功能注釋，能夠更好地理解物種內部和物種之間的變異。嗜血流感桿菌（Haemophilus influenzae）是世界上第1個測序完成全基因組序列的細菌[22]。Callanan等[23]完成了第1株L. helveticus DPC 4571的全基因組測序。Schmid等[24]對3 株L. helveticus進行了全基因組測序分析。盡管完整的基因組序列在準確描述核心基因組、鑒定特異基因以及作為功能基因組學研究方面具有特殊的價值，但它們在公共數據庫中的代表性仍然很低。NCBI（https://www.ncbi.nlm.nih.gov/genome/browse#!/eukaryotes/）數據庫顯示，截止2019年9月，乳桿菌屬（Lactobacillus）全基因組測序菌株為2 531 株，其中L. helveticus測序菌株僅有57 株，含有質粒的L. helveticus有14 株。到目前為止，就L. helveticus而言，國內外的研究主要集中在其發酵產生的生物活性物質對機體健康的促進作用以及對腸道菌群的影響[25-26]，有關L. helveticus全基因組測序分析為主題的國內外文獻及研究報道較少。

本課題組從自然風干羊肉中篩選到1 株L. helveticus TR13，前期研究表明，該菌株產脂肪酶能力較強且發酵特性良好，接種量2%、36 ℃發酵培養48 h，酶活性為10.89 U/mL，所產脂肪酶對中短鏈及長鏈底物均有降解作用，且對硝基苯酚棕櫚酸酯（C16）降解效果較好。為深入研究TR13的代謝機制，采用PacBio Illumina HiSeq測序平臺對細菌TR13進行全基因組測序分析和GO（Gene Ontology）、COG（Cluster of Orthologous Groups）、KEGG（Kyoto Encyclopedia of Genes and Genomes）數據庫基因組基本功能注釋，以期為深入了解TR13在羊肉發酵香腸生成過程中代謝機理的研究提供生物信息學基礎。

1 材料與方法

1.1 材料與試劑

L. helveticus TR13，內蒙古農業大學食品科學與工程學院肉品科學技術團隊前期從自然風干肉制品中分離、篩選獲得。

DNA抽提試劑盒（細菌）Wizard?Genomic DNA Purification Kit 美國Promega公司；DNA純化試劑盒DNA Clean Beads 中國諾唯贊公司；AP151-03 TransStart TopTaq DNA Polymerase 中國TransGene公司；二代建庫試劑盒NEXTflex? Rapid DNA-Seq Kit美國Biooscientific公司；三代建庫試劑盒SMRTbell Sample and Template Preparation Reagent Kits 美國Pacbio公司。

1.2 儀器與設備

MISEQ測序儀 美國Illumina公司；Pacbio RSII測序儀 美國Pacific Biosciences公司；GeneAmp?9700聚合酶鏈式反應（polymerase chain reaction，PCR）儀美國ABI公司；5424R高速臺式冷凍離心機 德國Eppendorf公司；TBS-380熒光儀 中國Sunnyvale公司；Covaris M220超聲 波破碎儀 中國Gene Company Limited。

1.3 方法

1.3.1L. helveticusTR13基因組DNA的提取

用M R S固體培養基（含1 5%瓊脂粉）活化L. helveticusTR13，挑取單菌落接種于10 mL滅菌MRS液體培養基中，37 ℃培養24 h，按2%接種量接種于500 mL MRS液體培養基中進行擴大培養，24 h后8 000 r/min離心8 min收集菌體，根據Bacterial Genomic DNA Isolation Kit試劑盒說明書，提取菌體基因組，并用0.8%瓊脂糖凝膠電泳檢測基因組質量，保證進行后續實驗的DNA質量足夠高。純化的基因組DNA采用TBS-380熒光儀進行定量。高質量的DNA（OD260nm/OD280nm=1.8～2.0，DNA總量≥15 μg，質量濃度≥50 ng/μL）用于后面的建庫測序。

1.3.2 文庫構建

基因組測序使用PacBio RS II單分子實時測序（SMRT）和Illumina測序平臺組合，Illumina的數據被用來評估基因組雜合度、基因組重復度、基因組大小、是否存在質粒及污染，同時保證更完整更精確的組裝。

Illumina文庫構建：取至少1 μg基因組DNA，利用超聲波破碎儀進行基因組DNA片段化，將DNA樣本剪切成小于400 bp的片段。使用NEXT flex?Rapid DNA-Seq試劑盒進行文庫制備。具體步驟如下：連接A&B接頭；篩選去除接頭自連片段；瓊脂糖凝膠電泳進行片段篩選，保留一端是A接頭、一端是B接頭的片段；氫氧化鈉變性，產生單鏈DNA片段；橋式PCR擴增。

單分子文庫建庫：取至少15 μg基因組DNA，利用超聲波破碎儀將基因組DNA處理成0～10 kb的片段，然后根據PacBio說明書進行片段純化，末端補平，兩端分別連接SMRT bell測序接頭。

1.3.3 Illumina預測

制備的文庫在Illumina HiSeq X Ten儀器上進行雙端測序（2×150 bp），利用PacBio RS II和Illumina平臺生成的數據進行生物信息學分析。經過一系列的質量剪切之后得到的高質量reads，稱之為clean data。利用canu及HGAP軟件進行PacBio數據組裝，將reads組裝成contigs，然后人工判斷成環，得到完整的染色體和質粒基因組。最后利用Illumina測序數據對組裝結果進行校正。本次測序委托上海美吉生物醫藥科技有限公司完成。

1.3.4 基因預測與注釋

利用Glimmer對基因組中的編碼序列（coding sequence，CDS）進行預測，質粒基因采用GeneMarkS軟件預測，tRNAscan-SE進行tRNA預測，Barrnap進行rRNA預測。利用BLAST、Diamond、HMMER等序列比對工具，將獲得的基因提交COG[27]、GO[28]、KEGG[29-30]數據庫進行比對，獲得功能注釋信息。

2 結果與分析

2.1 基因組組裝與預測

采用PacBio和Illumina HiSeq測序平臺，對細菌TR13進行完成圖測序，測序及預測結果如表1所示。TR13染色體基因組為環狀分子，染色體總數為1，基因組序列總長度為2 172 224 bp，GeneMarkS軟件預測沒有發現質粒基因。平均GC含量為36.82%。對細菌的編碼基因進行預測，獲得編碼基因（CDS數）2 356 個，編碼基因總長度為1 883 532 bp，平均長度為799.46 bp，基因密度（即基因組中每1 kb的堿基中所對應的基因數目）為1.08 個/kb，占基因組百分比為86.71%，編碼基因GC含量為37.47%，基因間區長度為288 692 bp，基因間區GC含量為32.57%，基因間區占基因百分比為13.29%。

表1 基因組統計及預測Table 1 Genome statistics and prediction

表2 編碼基因長度分布Table 2 Length distribution of protein-coding genes

如表2所示，序列長度分布在1 000 bp以上的基因數最多，數值為678，占全部編碼基因的28.78%。其中最小基因序列長度為114 bp，對應蛋白質序列長度為37 aa；最大基因序列長度為5 184 bp，對應蛋白質序列長度為1 727 aa。對應蛋白質平均序列長度為265 aa。TR13基因組中含串聯重復序列36 個，重復序列占基因組的總長度5 249 bp，占基因組百分比為0.28%。重復序列長度最小為40 bp，最大為702 bp，36 個重復序列平均長度為145.8 bp。

2 356 個編碼基因中，有1 891 個基因的起始密碼子為ATG，其中有1 286 個基因的終止密碼子為TAA、385 個為TAG、220 個為TGA；154 個基因的起始密碼子為GTG，其中有97 個基因的終止密碼子為TAA、42 個為TAG、15 個為TGA；311 個基因的起始密碼子為TTG，其中有175 個基因的終止密碼子為TAA、84 個為TAG、52 個為TGA。

TR13基因組中含有15 個rRNA操縱子基因，分別由5 個5S rRNA、5 個16S rRNA、5 個23S rRNA組成，其中3 個分布于負義鏈，5S rRNA、16S rRNA、23S rRNA各1 個。TR13基因組中含有63 個tRNA基因，tRNA類型及其數量如表3所示。

表3 tRNA預測結果統計Table 3 Statistics of tRNA prediction

由表3可知，tRNA基因分別轉運Arg、Gly、Leu、Ser等20 種不同的氨基酸，其中轉運Arg、Gly、Leu、Ser的基因數最多，均為5 個，轉 運Cys和Trp的基因最少僅有1 個。

2.2 基因組圈圖分析

基因組圈圖可以全面展示基因組的特征，如基因在正、反義鏈上的分布情況、基因的COG功能分類情況、GC含量、基因組島、同源基因等。通過CGView軟件繪制基因組圈圖，如圖1所示。

CGView圈圖最外面一圈為基因組大小的標識，L. helveticusTR13基因組為1 套環狀無質粒分子，基因總長度為2 172 224 bp。第2圈和第5圈為正鏈、負鏈上的CDS，不同的顏色表示不同的COG功能分類，編碼基因的數量為2 356 個，其中COG注釋的基因為1 942 個，占編碼基因的82.43%，分別從新陳代謝、信息存儲與處理、細胞作用與信號傳遞和無明顯特征4 大分類，復制、重組和修復，翻譯、核糖體結構和生物發生，碳水化合物運輸和代謝等18 個類型上得到了注釋。第3圈和第4圈分別為正鏈負鏈上的CDS、tRNA、rRNA，TR13細菌基因組中含有15 個rRNA操縱子，分別由5 個5S rRNA、5 個16S rRNA、5 個23S rRNA組成，含有63 個tRNA基因，分別轉運Arg、Gly、Leu、Ser等20 種不同的氨基酸。第6圈為GC含量，向外的部分表示該區域GC含量高于全基因組平均GC含量，峰值越高表示與平均GC含量差值越大，向內的部分表示該區域GC含量低于全基因組平均GC含量，峰值越高表示與平均GC含量差值越大，基因組平均GC含量為36.82%。第7圈為GC-skew值，具體算法為（G－C）/（G+C），在生物意義上該值為正值時正鏈更傾向于轉錄CDS，為負值時負鏈更傾向于轉錄CDS；最內一圈為基因組大小標識。基因組圈圖可以使研究者對菌株基因組的特征有更全面、更直觀的認識。

圖1 基因組CGViieeww圈圖Fig. 1 Graphical representation of the genome of strain TR13

2.3 基因組的功能注釋

對預測得到的編碼基因進行基礎的功能注釋，通過與COG、GO和KEGG數據庫比對進行功能注釋。TR13基因組中所有基因能夠注釋到GO信息的基因數目為1 972 個，能夠注釋到COG信息的基因數目為1 942 個，與KEGG數據庫比對能夠定位到具體Pathway的基因數目有1 143 個。

2.3.1 GO注 釋結果

根據對比結果，采用BLAST2go對比軟件，在數據庫中對TR13進行GO注釋，如 圖2所示。TR13基因組中含細胞組成、生物過程、分子功能3大類型，在GO分類中共有1 792 個基因得到了功能注釋，占所有編碼基因（2 356 個）的百分比為76.06%。其中：細胞組成相關的基因數量878 個；生物過程相關的基因數量1 295 個；分子功能相關的基因數量1 478 個。1 792 個基因共注釋了40 種功能特性。

在生物過程這一層面分別得到代謝過程、細胞過程和單一生物過程等18 種功能注釋，其中和代謝有關基因數量最多為1 019 個，與細胞過程相關基因855 個，與單一生物過程相關基因609 個。

在細胞組成層面得到10 種功能注釋，其中細胞膜和細胞膜組分功能相關的基因最多。在分子功能層面分共得到12 種功能注釋，其中與催化活性有關的基因最多為1 083 個，其次為細胞附著相關基因861 個。功能基因中與抗氧化活性相關的基因有硫氧還原蛋白-二硫鍵還原酶、谷胱甘肽還原酶、NAD（FAD）依賴性脫氫酶等7 個基因，信息如表4所示。

2.3.2 COG注釋結果

采用Diamond軟件，參數設置為E-value 10－5，在數據庫eggNOG（http://www.ncbi.nlm.nih.gov/COG/）中對TR13進行COG注釋，如表5所示。

表5 COG注釋結果詳情Table 5 COG annotation results

TR13編碼基因在COG中的分類數量為4，COG的類型數量為18，編碼基因組中共有1 942 個基因在COG數據庫中得到了注釋，注釋基因占比為82.43%。

COG分類分別為新陳代謝、信息存儲與處理、細胞作用與信號傳遞和無明顯特征4 種類型，其中和新陳代謝有關的基因數為615 個，分為8 種類型。信息存儲與處理相關基因565 個，分為3 種類型。細胞作用與信號傳遞相關基因261 個，分為6 種類型，特征較弱的未知功能基因有512 個，歸為1 種類型。

表6 脂質的運輸與代謝相關基因Table 6 Lipid transport- and metabolism-related genes

新陳代謝COG聚類主要集中于碳水化合物的運輸與代謝（相關基因146 個）、氨基酸的運輸與代謝（相關基因143 個）、核苷酸的運輸與代謝（相關基因90 個）、無機離子的運輸與代謝（相關基因82 個）、脂質的運輸與代謝（相關基因44 個）。其中脂質的運輸與代謝如表6所示，gene0183、gene0428、gene2011分別為3 種不同酯酶的編碼基因，gene0893作為三酰甘油脂肪酶基因具有分解脂肪的作用。gene1050（mvaK2）、gene1052（mvaK1）、gene1915（dagK）為3 種不同的激酶基因，是一類可以從高能供體分子（如ATP）轉移磷酸基團到特定靶分子（底物） 的酶，此過程為磷酸化。gene0547（bcrC）為磷酸酶基因，與激酶的作用相反。gene0301（acpS）為holo-酰基載體蛋白合成酶基因，gene0931（acpP）和gene2286（acpP）為酰基載體蛋白質基因，具有合成脂肪酸的作用。

2.3.3 KEGG注釋結果

KEGG（http://www.genome.jp/kegg/）是基因組研究方面的公共數據庫。該數據庫是系統分析基因功能，聯系基因組信息和功能信息的大型知識庫。采用Diamond對比軟件，參數設置為E-value 10－5，在KEGG數據庫中對TR13進行KEGG注釋。

L. helveticusTR13在KEGG數據庫中共有1 144 個基因分別在代謝、遺傳信息處理、環境信息處理、細胞過程、生物體系統、人類疾病6大功能35 個通路上得到功能注釋，結果如圖3所示。樣本中涉及的代謝通路總數以及參與每個代謝通路的基因組數目，分析發現在代謝途徑層面特別是碳水化合物代謝、核苷酸代謝和氨基酸代謝得到較多的基因功能注釋。其中634 個基因在代謝通路上得到注釋，占KEGG中全部注釋基因的55.42%，12 個代謝通路中，碳水化合物代謝相關的基因為132 個，占代謝通路注釋基因（634 個）的20.82%，與脂質代謝通路有關的基因有48 個。110 個基因在環境信息處理層面得到注釋，其中與膜運輸有關的基因為86 個，與信號傳導有關的基因24 個。

圖3 L. helveticus TR13的KEGG功能分類圖Fig. 3 KEGG functional classification of L. helveticus TR13

表7 脂質代謝通路及其相關基因Table 7 Lipid metabolism pathways and related genes

如表7所示，TR13基因組中主要包含了脂肪酸合成、脂肪酸降解、不飽和脂肪酸合成、初級膽汁酸生物合成、次生膽汁酸生物合成、酮體的合成與降解、甘油脂類代謝、甘油磷脂代謝、鞘脂類代謝共9 條通路信息。其中ko00061通路為脂肪酸合成通路，TR13基因組序列分析共獲得gene0428、gene0931、gene1746等15 個相關基因，ko00071為脂肪酸生物降解通路，在TR13基因組中包含gene1800和gene1984相關基因2 個。

表8 脂肪酸生物合成與降解通路信息Table 8 Information about fatty acid biosynthesis and degradation pathways

如表8所示，ko0061為脂肪酸合成代謝通路，包含gene0428、gene0931、gene1746等15 個基因，基因大小分布于243～1 380 bp之間，其中gene0931和gene2286為酰基載體蛋白質acpP基因，酰基載體蛋白是脂肪酸合成中的關鍵蛋白質，作為脂酰基的載體將脂酰基從一個酶反應轉移到另一個酶反應。Ko00071為脂肪酸分解代謝通路，TR13基因組中包含了與其相關的基因有gene1800乙酰CoA酰基轉移酶和乙醛脫氫酶/乙醇脫氫酶gene1984。ko01040為不飽和脂肪酸的生物合成代謝通路，TR13基因組中包含了gene1746和gene2284，均為酰基載體蛋白還原酶基因（fabG）。

2.3.4 毒力基因預測結果

來源于微生物并對微生物自身侵染和引起特定宿主疾病具有促進作用的物質稱為毒力因子，主要包括細菌毒素、調節細菌黏附作用的細胞表面蛋白、對細菌本身具有保護作用的蛋白、細胞表面碳水化合物以及具有細菌致病性的水解酶等。毒力因子數據庫VFDB，是一個全面一體化的病原菌毒力因子信息管理網絡數據庫。利用軟件Diamond比對VFDB核心數據庫數據庫，獲得毒力因子基因注釋概況并進行統計結果如表9所示。L. helveticusTR13有對病毒相關基因的調控、防守毒力因子、進攻毒力因子和非特異性毒性因子4 個類型的毒力基因，其中鐵吸收系統毒力因子基因數最多為24 個，其次為黏附毒力基因17 個，共預測到可能的毒力基因共75 個。

表9 毒力基因預測分類統計Table 9 Classi fication statistics of the predicted virulence genes

2.3.5 耐藥基因預測結果

利用軟件Diamond，通過CARD數據庫注釋，獲得每個基因組中包含的耐藥基因的信息。L. helveticusTR13共預測到150 個耐藥基因，耐藥基因預測分類統計如圖4所示。主要有抗生素耐藥基因、變異或突變體耐藥基因、糖肽耐藥基因、通過分子旁路作用產生抗生素耐藥性的基因、調節抗生素外排的基因、四環素抗性基因、氨基糖苷類耐藥基因、內酰胺耐藥基因等。

圖4 耐藥基因預測分類統計Fig. 4 Classi fication statistics of the predicted drug resistance genes

3 結 論

L. helveticus TR13基因組為1 套環狀無質粒分子，總長度為2 172 224 bp。GC含量為36.82%，預測到2 536個編碼基因，其中編碼基因總長度為1 883 532 bp，平均長度為799.46 bp，平均密度為1.08 個/kb，對應蛋白質平均序列長度為265 aa，占基因組百分比為86.71%，基因中包含15 個rRNA操縱子和63 個tRNA。

TR13基因組中能夠注釋到GO信息的基因數目為1 972 個，包含了40 種功能特性，占所有編碼基因（2 356 個）的76.1%。能夠注釋到COG信息的基因數目為1 942 個，注釋基因占比為82.43%。在KEGG數據庫中共有1 144 個基因分別在代謝、遺傳信息處理、環境信息處理、細胞過程、生物體系統、人類疾病6 大功能35 個通路上得到功能注釋。預測到可能的毒力基因75 個，耐藥基因150 個。

脂質代謝的通路信息分析顯示，TR13基因組中主要包含了脂肪酸合成、脂肪酸降解、不飽和脂肪酸合成、初級膽汁酸生物合成、次生膽汁酸生物合成、酮體的合成與降解、甘油脂類代謝、甘油磷脂代謝、鞘脂類代謝共9 條通路信息。TR13基因組中包含了ko00061脂肪酸合成通路相關的15 個基因，ko00071脂肪酸生物降解通路相關的gene1800和gene1984基因2 個。