风干肉中产脂肪酶瑞士乳杆菌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 个。