魔芋软腐病病原菌果胶杆菌ZX67菌株全基因组测序及分析

夏 涛，李梦飞，张 雯，陈雪英，史银连，王红莹，李 刚，祝 伟

（国家民委重点实验室/中南民族大学生命科学学院，武汉 430074）

0 引言

魔芋是中国中西部十余省市山区人民的重要经济作物，是山区人民的主要经济来源之一，魔芋粉是难得的优良保健食品，据估计中国魔芋种植业的价值在2000亿元以上，而近几年来，受到魔芋软腐病的危害，使得曾经大规模种植的魔芋种植面积锐减，山区人民的收入受到了严重影响。虽然自20世纪80年代以来全国众多的高校、科研院所投入了大量的人力和财力开展了软腐病的防治研究，但至今仍然没有取得有效的技术方法来遏制软腐病的发生和泛滥。

软腐病是魔芋病害中影响最为严重的一种细菌性病害，公认是由Pectobacterium引起的[1-2]，被认为是阻碍魔芋加工业蓬勃发展的主要原因。软腐病菌不仅在田间、而且在储存期间都能对魔芋球茎侵染导致腐烂、造成重大损失[1-2]。魔芋软腐病的病原菌包括Pectobacterium属的胡萝卜果胶杆菌(Pectobacterium carotovorum)和菊果胶杆菌(P.chrysanthemi)[3-5]。早在1945年果胶杆菌属被分类为软腐欧文氏菌属的腐烂病原菌[3-4]，现在包含的菌株涵盖了多种植物病原菌[4]，会引起这些植物的软腐、茎腐烂、枯萎和黑腿病[4-5]。该病原菌是一类营腐生性寄生菌，可从寄主的自然孔口和外部伤口侵入，不过大多数病原菌最主要的侵染途径是从机械伤口或虫伤侵入。虽然全球范围内已投入了大量人力物力进行研究，但到目前为止，在防治魔芋软腐病方面始终没有突破性进展，还没有找到有效的防治技术方法。病原菌ZX67菌株分离自魔芋之乡竹溪县感染软腐病的魔芋球茎深处，在此前研究中，依据细胞形态、生理生化、16S rDNA比对分析的结果将其鉴定为胡萝卜果胶杆菌胡萝卜亚种(P.carotovorumsubsp.carotovorum)，为了能更好地研发魔芋软腐病的有效防治方法，进一步了解软腐病病原菌的特性，在用酚-氯仿法提取了全基因组DNA后对其进行全基因组测序和分析，为深入研究软腐病的防治方法提供更丰富的数据。

1 材料与方法

1.1 实验时间、地点

室内实验于2019年9月—2021年11月在中南民族大学生命科学学院实验室进行。

1.2 实验材料

1.2.1 菌株 从腐烂魔芋球茎深处分离出的病原菌胡萝卜果胶杆菌胡萝卜亚种ZX67菌株。

1.2.2 培养基 LB培养基。蛋白胨2.0 g，酵母提取物1.0 g，NaCl 1.0 g，H2O 200 mL，固体平板加1.5%琼脂。121℃，灭菌25 min。

1.2.3 培养条件 固体平板和液体三角瓶均置于30℃培养24 h。

1.3 病原菌基因组DNA的提取方法

将菌株ZX67接种到50 mL液体培养基中，放置于30℃摇床200 r/min过夜培养后，取菌液于离心管中在4℃、11000 g离心7 min，移出上清液，细胞沉淀加缓冲液(10 mmol/L Tris-HCl pH 7.5，1 mmol/L EDTA)重悬，按比例加入蛋白酶K(1 mg/mL)和SDS(10%w/v)于50℃保温1 h破细胞，以碱裂解法抽提基因组DNA[6]，纯化基因组DNA是使用酚-氯仿-异戊醇抽提、70%乙醇沉淀的方法，用琼脂糖凝胶电泳和NanodropOne(Thermo Fisher Scientific,USA)测量光吸收值来评估纯化的DNA纯度、质量和浓度。

1.4 全基因组测序及数据处理与分析

将提取的基因组DNA样品利用第三代测序技术进行测序，测序仪器为第三代测序仪PacBio RSII(Pacific Biosciences,USA)[7]。经过质量检测合格后，使用建库试剂SMRTbellTM Template prep kit 1.0完成对样本的文库构建工作，用测序试剂盒DNA Binding kit P6 V2和DNA Sequencing Kit 4.0 V2进行测序[7]。用de novo（基因组从头测序）即在不需要任何参考序列的情况下，对基因组进行测序，然后使用生物信息学软件和分析方法对测序得到的数据进行拼接和组装。将蛋白序列基于BLAST与各数据库进行比对(blastp,evalue≤1e-5)后得到对应的功能注释信息。所有的注释均使用BLAST软件结合各个数据库的特点完成，目前已知基因功能的数据库有KEGG(Kyoto Encyclopedia of Genes and Genomes)[8]，COG(Cluster of Orthologous Groups of proteins)[9]，GO(Gene Ontology) (http://geneontology.org/)[10]，NR(Non-Redundant Protein Database)[11]，Swissprot[12]，Trembl[13]，Pfam[14]。为了鉴定RNA基因，用tRNAscan-SE软件预测得到tRNA[15]；rRNA分别由通过与近缘参考序列的rRNA库比对，及使用rRNAmmer软件预测得到[16]；sRNA首先基于Rfam软件与Rfam数据库进行比对注释[17]，然后通过cmsearch程序进行最终确定[18]。基于其序列组成，由IslandPath-DIOMB软件预测得到基因岛[19]，该软件主要通过检测序列中核苷酸偏向性和移动性基因（如转座酶或整合酶）来预测基因岛及潜在的水平基因转移。实验中的转座子主要由TransposonPSI软件进行预测，该软件是基于PSI-Blast的同源比对来对基因组或者蛋白序列进行转座子预测的，适用于细菌和真菌。本实验通过PHAST软件预测前噬菌体[20]。同时，分别使用RepeatMasker软件[21]和TRF(Tandem repeats finder)软件[22]对散在重复序列和串联重复序列进行预测。并采用CRISPRFinder对样品基因组进行CRISPR预测[23]。同时将测序得到的序列对其数据进行整理分析，并查阅相关文献后确认基因组上存在的与降解细胞壁相关的酶基因和蛋白分泌系统，并查找与致病性相关的毒力因子[24]。

2 结果与分析

2.1 病原菌基因组DNA的提取和质量分析

采用碱裂解法提取了菌株ZX67的基因组DNA，经过反复纯化后得到200 μL的样品，取2 μL DNA样品用于1.0%的琼脂糖凝胶电泳检测，在电泳后的凝胶上显示所提取DNA条带清晰完整（图1），取1 μL该样品在Nanodrop One上检测，测得其DNA浓度约为107.6 ng/μL，质量为21.5 µg，260:230为2.3，260:280为1.87，样品质量、浓度、纯度均符合测序要求。

图1 基因组DNA电泳图谱

2.2 全基因组组装结果统计

从感染软腐病的腐烂魔芋球茎中分离到的菌株ZX67，其16S rDNA序列与胡萝卜果胶杆菌(P.carotovora)具有99%以上的同源性，P.carotovora在魔芋、土豆等许多植物中都能引起软腐病[25]。利用pacbio rsii技术（单分子实时测序系统）[26]对该菌株的全基因组进行测序，用de novo进行了组装。最后得出的基因组参考序列总长度为4,909,724 bp，GC含量为51.27%（表1）。组装后，将原始数据与组装后的参考基因组序列进行比较，并计算参考基因组序列每个位置的覆盖深度来评估基因组的完整性、装配和排序均匀性，从覆盖分布来看，每个区域的覆盖深度都很高，平均覆盖深度为107 x，一般测序深度在30 x以上，基因组覆盖度和测序错误率的控制得以保证，该结果表明所获得的基因组被高度覆盖，装配结果具有良好的完整性。针对测序样品的组装基因组序列，结合编码基因的预测结果，使用Circos软件[27]对样品基因组进行展示，全基因组图谱如图2，从外向内分别展示了碱基数目对应位置、正链编码基因的对应位置、负链编码基因对应位置、基因组GC含量分布、基因组GC skew（偏移）值分布。菌株ZX67全基因组序列在基因库(GenBank)中的登录号为：CP034211。

表1 组装结果统计表

图2 全基因组圈图

2.3 基因组注释结果和特性分析

P.carotovorumZX67菌株基因组序列包含了4,909,724 bp，与数据库中的已有P.carotovorum序列相比较，碱基数目较高，测序的覆盖率较高。其基因组的GC含量为51.27%，一般的GC%越高，则DNA结构就趋于高稳定性，而较高的AT%易于解链、可能与较高的基因突变率相关，这也可能与该类菌种在20～30℃左右的生存环境温度相适应，在一定程度上也暗示了作为病原菌需要适应宿主防卫行为的变化尤其是增强时适者生存的碱基高自由能状态向高突变率状态转变的可跨越阈值。采用Prodigal v2.6.3[19,26]软件对该基因组分析，总共计算和预测了4,977个编码基因，编码基因的平均长度为1,041 bp，其序列长度占总基因组的85.25%，预测的非编码RNA(ncRNA)主要包括tRNA（74个），rRNA（22个）及sRNA（100个）。基因岛有10个，每个基因岛的平均长度为36,027 bp。预测到转座子6个，每个转座子的平均长度为363 bp。预测的散在重复序列有261个以及串联重复序列103个（表2）。预测存在2个前噬菌体，每个前噬菌体的平均长度为26,258 bp，有5个CRISPR序列，每个序列的平均长度为621 bp（表2、表3）。

表2 基因组特征数据统计表

该菌株的全基因组测序后产生了大量数据，为了解析、注释其包含的编码基因数目、功能等，笔者采用了目前普遍采用比对的方法对预测出来的编码基因进行功能注释，但基因不同于蛋白，不能通过结构来预测功能，所以分别通过与已知的基因功能数据库NR、Swiss-Prot、GO、GOG、KEGG、Trembl和Pfam的比对来推测了该基因组的编码基因数目、功能等；而与每一个数据库比对所得到的编码基因数目、种类等数据结果存在较大的差异。虽然在得到的结果中，一般都是采用氨基酸序列与数据库进行相似性比对，但其中通过与GO和KEGG数据库比对的结果对后续的实验研究具有更重要的作用，因为它们分别在基因功能和代谢通路研究中占据重要地位，这对于后续的改造、寻找防控的靶点等研究很重要。

在已经预测的10个基因岛内，通过数据查找分析发现在果胶杆菌ZX67菌株的基因组中存在有5种与致病性相关的蛋白质分泌系统，其中III型分泌系统位于预测的第4个基因岛内，而IV型分泌系统位于预测的第5个基因岛内，通过在抗性基因数据库(https://card.mcmaster.ca/analyze)进行比对，预测到基因岛中还可能含有可抑制细菌蛋白质合成的黄色霉素(kirromycin ARO:3003368)的抗性基因（相似性为92%），表明基因岛可能与细菌的适应性和致病性相关[28]。

而在基因组中预测到的两个前噬菌体区域不是完整的前噬菌体，第一个前噬菌体序列的长度为38,859 bp，GC含量为48.79%，包含48个CDS编码区域，可编码外壳蛋白和整合酶。另一个前噬菌体长度为13,658 bp，GC含量为53.51%，包含18个CDS编码区域，能编码噬菌体主要的尾管蛋白。此两个噬菌体序列均为不完整的噬菌体，不具有典型噬菌体基因组的完整特性，将这两个前噬菌体序列与NCBI数据库中公布的71种果胶杆菌和2种欧文氏菌的噬菌体基因组序列进行比对，结果显示预测的前噬菌体与其他噬菌体的同源性都不高，均为30%左右，一方面可能是由于序列不完整，另一方面也可能证实了噬菌体基因具有很高的遗传多样性。

另外，此次分离的病原菌基因组上含有5个长度不一的CRISPR区，最长的CRISPR序列长度是1,529 bp，含有25条间隔序列，最短的CRISPR序列长度是98 bp，只含有1条间隔序列（表3）。这预测到的5种CRISPR序列，在数据库中的Pectobacterium属的多个种、Acinetobactersp.NEB 394、Clostridioides difficilestrain CD10010、Proteus terraesubsp.cibarius strain HNCF44W等菌种的全基因组序列上都可以搜索到高度相似的序列。

在全球范围内由果胶杆菌侵染导致的植物软腐病均有发生[12]，为探究菌株ZX67与其他胡萝卜果胶杆菌之间的差异，笔者将菌株ZX67与NCBI数据库中已完成测序的果胶杆菌属的部分菌株信息进行比较（表4），所选菌株分别分离自中国、日本、白俄罗斯、波兰、挪威、巴基斯坦以及南韩等不同国家，它们的寄主也有所不同，除了ZX67是从魔芋中分离出来的，其他菌株则是以山葵、黄瓜、马铃薯和萝卜为寄主，其中来源于马铃薯的居多，如P.versatile14A、P.parmentieriIFB5604、P.parmentieriIFB5432、P.polarisNIBIO1006和P.punjabenseSS95，这表明马铃薯软腐病是由这几种果胶杆菌单独或混合侵染而发生软腐病。菌株ZX67在基因组大小与GC含量上与这些菌株的差异并不明显，但CDS和CRISPR的数量比其他菌株都要多，表明其基因组上能编码的蛋白数量和种类相对要多，这些蛋白可能与细胞壁的降解及致病性相关。

表4 菌株ZX67与NCBI数据库中果胶杆菌属菌株的比较

2.4 基于基因组序列的重要代谢途径分析

果胶杆菌感染的特征在于分泌大量降解植物细胞壁的酶引起的诱发腐烂症状[29]。果胶杆菌可编码多种酶协同作用以降解植物细胞壁。果胶酶是致病过程中最重要的酶，受许多调节因素的控制，以多种形式存在并利用中间薄层和果胶质，使细胞瓦解、内含物外流。在预测的全基因组编码的基因组成中，发现胡萝卜果胶杆菌ZX67菌株的基因组中包含有18种果胶酸裂解酶基因（pelABCLWX、pel123、KKH3_37720、KKH3_37730、KKH3_37740、KKH3_21450、KKH3_21500、KKH3_28390、KKH3_17430、KKH3_08410和PCC21_038530），2种果胶裂解酶(pnl、KKH3_12920)，1种多聚半乳糖醛酸酶(KKH3_42200)，2种鼠李糖半乳糖醛酸裂解酶(rhiE、KKH3_05330)，1种果胶甲基酯酶(pemA)，2 种果胶乙酰化酶(KKH3_29960、KKH3_18740)，3 种纤维素酶(celA1、KKH3_20670、KKH3_14460)，5种葡聚糖酶(bcsZ、cel5、celS、KKH3_22720、KKH3_41050)。此外，该基因组还含有5个基因编码的纤维素酶prt1、prtCDEF，其中一些酶可能在植物细胞壁降解中起作用。这些纤维素酶主要表现内切、外切葡聚糖酶和β-Glucosidase活性，首先降解纤维素释放寡糖、纤维二糖和葡萄糖，然后降解植物的细胞壁（表5）。有人研究和比较了在菊果胶杆菌和胡萝卜果胶杆菌中一些调控基因参与生产这些细胞壁降解酶，如菊果胶杆菌上的pir、pecS-pecM、fur、hns、crp、pecT、expI-expR、kdgR，胡萝卜果胶杆菌中的aepA、rexZ、hor、pehR-pehS、expM、rpoS、rsmA-rsmB-rsmC、hexA、kdgR和expI-expR等基因[29]；在ZX67的全基因组上则发现了kdgR、crp、pir、pecS和rpoS调控基因的存在，这5个调控基因中，既有均出现在2种病原菌中的kdgR，也有分别出现在菊果胶杆菌上的crp、pir、pecS和胡萝卜果胶杆菌中的rpoS，具有一定的复杂性。在果胶酶表达分泌的过程中，会受到一些调节因素的调控，主要的调节蛋白有kdgR、crp、pir和pecS。除果胶酶外，其他与细胞壁降解相关的酶如果胶裂解酶(pel)、多聚半乳糖醛酸酶(peh)、纤维素酶(pel)和蛋白酶(prt)等也被认为是致病的主要因素，这些酶大量产生并分泌到外部环境中[29]。

表5 ZX67基因组中与植物细胞壁降解相关的酶

根据对基因组编码的基因组成分析，在果胶杆菌ZX67菌株的基因组中存在5种类型的蛋白质分泌系统（表6）。蛋白质分泌在软腐病的发病机制中起着重要作用[30]。菌株ZX67基因组上的这5种类型蛋白质分泌系统分述如下：Ⅰ型分泌系统包括6个基因，分别为 PC_1299、PC_3006、PC_3169、PC_3686、PC_3688、PC_3689，编码permease/ATPase和外膜蛋白。Ⅱ型分泌系统包括 15个基因，PC_783、PC_784、PC_1484～1493、PC_1495～1497，是果胶杆菌最重要的分泌系统，又称out系统(outcdefghijklmn)。Ⅱ型分泌系统的15个基因编码转运蛋白HofC同源物和Ⅱ型分泌系统蛋白各组分：C、D、E、F、G、H、I、J、K、M和N。II型分泌系统可将各种蛋白质从周质空间运输到细胞外环境。它负责分泌大多数植物细胞壁降解酶，如纤维素酶、果胶酶和其他一些毒力因子，从而破坏宿主细胞，引起组织坏死和疾病[30]。Ⅲ型分泌系统包括19个基因，编码发病机制、蛋白质分泌、膜和质膜的组成部分。在菌株ZX67的基因组中，还发现存在IV型分泌系统，它包括16个基因，编码ATPase、转糖基化酶、P-type DNA transfer ATPase VirB11、Rhs蛋白、DotU蛋白等基因。与菌株ZX67不同的是，在胡萝卜果胶杆菌SCC1菌株的基因组中没有发现IV型分泌系统(T4SS)[30]。一些病原菌中的T4SS成分通过细菌与宿主细胞的粘附起到毒力和感染的作用[31]。VI型分泌系统编码包括hcp1、vgrGA蛋白、icmF等17个基因。预测结果显示，VI型分泌系统基因在果胶杆菌ZX67基因组中呈零星分布，但在胡萝卜果胶杆菌SCC1菌株的基因组中，VI型分泌系统的基因是成簇分布的[30]。据报道，VI型分泌系统对某些果胶杆菌的毒力几乎没有影响[30]。但在P.carotovorumZX67基因组中，VI蛋白分泌系统成分hcp1是分泌性细胞毒素hcp。

表6 菌株ZX67的基因组中蛋白质分泌系统

2.5 基于全基因组序列的菌株ZX67分类地位分析

病原菌ZX67全基因组测序完成后，在NCBI数据库中找到大约20条已公布测序结果的P.carotovorum菌株的基因组信息(https://www.ncbi.nlm.nih.gov/genome/genomes/1799)，建树观察ZX67基因组与各胡萝卜果胶杆菌菌株之间的亲缘关系（图3）。在胡萝卜果胶杆菌种中，虽都是同一个种的菌株，但其亲缘关系的远近都存在差异，菌株ZX67与P.carotovorumsubsp.carotovorumPC1[25]的亲缘关系紧密，而与菌株DSM 30168[25]的关系较远，由于测序的年代和技术差异，部分测序结果未能覆盖全基因组，不同环境分离出来的菌株其基因组也可能存在差异，依据这些差异，把P.carotovorum又划分为不同的亚种，其中P.carotovorumsubsp.carotovorum是最早提出来的亚种。不同来源的亚种基因组上的差异可能是适应寄主植物的环境变化选择的结果。

图3 菌株ZX67在P.carotovorum种的基因组中的分布位置

3 讨论与结论

胡萝卜果胶杆菌是马铃薯、胡萝卜、白菜等经济作物的常见病原菌，曾经甚至目前还在继续造成经济损失，早已引起世界各国的重视。本实验是对分离自竹溪县高山上发生软腐病的魔芋球茎中的病原菌菌株ZX67进行了全基因组测序，测序的结果覆盖全面，组装了基因组全长序列，用多个基因分析软件对基因组进行了分析，与已发布的同一个种的基因组序列及其文献、数据相比较，更加深入、全面、先进。目前，在NCBI数据库中，已经有大约60个胡萝卜果胶杆菌的全基因组测序结果，其中大约20条P.carotovorum的全基因组测序结果比较完整，另外的40余条由于测序的年代较早和技术差异，测序结果未能覆盖全基因组。在20条全序列比对中，菌株ZX67与菌株DSM 30168的关系较远，而与P.carotovorumsubsp.carotovorumPC1的亲缘关系紧密。这些菌株由于来源不同又存在着部分的差异，依据这些差异，已把P.carotovorum划分为不同的亚种，其中P.carotovorumsubsp.carotovorum是最早提出来的亚种。不同来源的亚种基因组上的差异可能与适应不同的环境变化的结果。

通过对病原菌全基因组序列分析，发现在菌株ZX67基因组上存在大量与降解细胞壁相关的酶基因，并且与其他的菌株有一定的差异，例如胡萝卜果胶杆菌PCC1基因组上有1种果胶裂解酶，4种多聚半乳糖醛酸酶，1种鼠李糖半乳糖醛酸裂解酶，2种果胶甲基酯酶，2种果胶乙酰化酶，3种纤维素酶和11种果胶酸裂解酶，各酶系在数量上与ZX67基因组上的有所不同（表5），ZX67基因组上有8种纤维素酶，有18种果胶酸裂解酶等。其次，二者的分泌蛋白系统也有所差异，除了菌株PCC1含有的分泌系统外，在菌株ZX67基因组上还发现了IV型分泌系统的多个基因；另外VI型分泌系统在分布上也有所不同，有些菌株是成簇分布的，还有一些菌株是分散的基因。在测序方面，本次测序重点突出了对CRISPR、重复序列、前噬菌体等与噬菌体的侵染、防御相关的基因序列的测序与分析，为后续利用噬菌体防治魔芋软腐病的技术改造提供了一个良好的基础。目前，在尚无有效方法能够攻克魔芋软腐病防治难题的困境下，加强CRISPR系统和前噬菌体的相关研究将可能是研发新的防治方法的有效切入点。因细菌基因组中包含的CRISPR序列是病原菌识别并抵御噬菌体侵染的重要元件，通常与CRISPR相关蛋白一起构成自身的获得性免疫系统，CRISPR序列由短重复序列和间区序列组成，重复序列大多高度保守[23,25]；而间区序列是从外源DNA序列中获得的，具有很高的多态性，当含有相同序列的外源DNA入侵时，就会被机体识别，对其进行剪切使之沉默而消除了噬菌体侵染的能力。CRISPR系统能够为细胞提供抵抗噬菌体侵染的免疫能力，相应地，噬菌体以其基因组的高频突变方式去躲避被识别被裂解、以anti-CRISPRs（抗CRISPR蛋白）反攻击宿主细胞，双方在互杀互搏的生存斗争中不断进化，通过激烈争夺生存压力制高点的殊死搏击恰恰使双方达到了互动平衡，并推动了二者基因组的快速进化。因此，若能够利用这些序列作为改造的靶点，将有可能规避宿主细胞的防卫系统而达到防治的目标。

在这项研究中，笔者提出了魔芋软腐病病原菌P.carotovorumZX67菌株的注释基因组序列，其由一个总长度为4,909,724 bp、GC含量为51.27%的片段组成，是从发病魔芋球茎深处分离出来的，根据病原菌生活方式，发现菌株ZX67的基因组具有与其他已测序的Pectobacteriumgenomes类似的大量植物细胞壁降解酶，有与侵染相关的毒力和蛋白分泌系统，并重点分析了CRISPR系统和前噬菌体等有助于建立防治新方法的序列。对胡萝卜果胶杆菌ZX67菌株全基因组的测序及分析的丰富结果，为进一步了解软腐病菌的特性提供了机会，也为研究新的防治措施提供了新的出发点和思路。