iMetaOmics | 江苏省农科院植物细菌团队-解析中国梨火疫菌特征及溯源分析

全基因组学揭示中国梨火疫菌株遗传变异、质粒多样性及溯源分析

iMeta主页：http://www.imeta.science

●  对中国的13株梨火疫病菌以及哈萨克斯坦和吉尔吉斯坦的各1株梨火疫病菌进行了全基因组测序和注释；

●  通过对本研究中15菌株与全球其他区域梨火疫菌株进行全基因组进化树分析，推测梨火疫病菌传播途径；

●  与模式菌株 CFBP1430 基因组比较分析，进一步诠释了该15株梨火疫病菌 基因组的共线性、序列变异（包括SNPs以及DIPs）以及质粒分布等问题;

●  通过各菌株接种梨幼果探讨了其致病性的差异。

引  言

梨火疫病菌（Erwinia amylovora）是世界十大最重要的植物病原细菌之一，常导致苹果、梨等蔷薇科寄主植物的大面积死亡，带来巨大的经济损失。该病害于 1780年首次发生于美国纽约哈德逊河流域，这也是科学史上首个被证实能引起植物病害的细菌。目前已传播到全球60多个国家和地区，包括美洲、欧洲的大部分地区以及非洲、亚洲的部分地区。在过去的100余年中，对E. amylovora的研究工作已从多维度展开。早期工作侧重于该病原菌形态学、传播途径及化学防控等方面，而近二十年来，研究的焦点逐渐转向其致病性的分子机制。该病原菌侵染寄主主要依赖于III型分泌系统（T3SS）以及胞外多糖。其中T3SS主要由Hrp类蛋白（Hypersensitive response and pathogenicity）和Hrc类蛋白（Hrp-conserved）组成，对应的编码基因均聚集在PAI毒力岛上（Pathogenicity islands）。该PAI编码了约60个毒力相关基因，分为四个区域，包括hrp/hrc区域、HAE区域（Hrp-associated enzyme）、HEE区域（Hrp effectors and elicitors）以及IT区域（Island transfer）。E. amylovora通过T3SS分泌多种蛋白，包括效应蛋白DspA/E、Eop1、Eop2、Eop3和Eop4（AvrRpt2EA）、伴侣蛋白DspB/F以及Hrp蛋白（hrpN、hrpW等）。其中，DspA/E和HrpN是E. amylovora重要的致病蛋白。胞外多糖EPS（Exopolysaccharides, EPS）是一种在分泌到细菌表面的长链多糖。已知与梨火疫菌致病性相关的EPS包括梨火疫菌素（amylovoran）和果聚糖（levan）。EPS具有多重功能，包括保护E. amylovora免受脱水和营养匮乏的影响、促进菌体渗透到寄主植物中、阻断寄主维管系统的运输，并保护菌体不被寄主防御系统所识别等。

在中国，2015年首次在新疆维吾尔自治区伊犁哈萨克自治州苹果树上发现梨火疫病，并在2016年严重损害苹果产量。2017年，梨火疫病蔓延至库尔勒市，对‘库尔勒香梨’产业造成严重破坏，2020年该病害蔓延至甘肃。为深入了解我们中国梨火疫病菌遗传信息，我们首先对中国新疆、甘肃不同地区分离的13株（KL20-28, KL17-17, ZY20-3-2，AKS17-22, KL18-6, KL19-1, KL21-15, KS22-2, YL16-01, YL16-15, ZYGT-2, 31070-3, 1506-4）以及邻国哈萨克斯坦KAAL06和吉尔吉斯斯坦KyPL01菌株进行了全基因组测序（图1A；表S1）。每个菌株的染色体大小约为3.8 Mb，且其中质粒有不同分布模式。将所有测序菌株基因组与模式菌株CFBP1430进行比较。分析了单核苷酸多态性（SNPs）、短缺失、插入和其他多态性（DIPs）。此外，通过系统发育分析，阐明了全球不同地理区域梨火疫菌株之间的进化关系。这些结果表明，所有中国菌株都有一个共同的祖先，并与哈萨克斯坦菌株和吉尔吉斯斯坦菌株密切相关。因此，我们认为梨火疫病向中国传播的主要途径可能遵循从中东到中亚的轨迹，最终至新疆及甘肃。

结果与讨论

测序菌株全基因组（图1A; S1−S2）变化范围为3,798,299~3,805,877 bp，基因组GC含量均为53.6%（表S1−S2）。利用Glimmer 3.0对每个菌株全基因组内的蛋白编码序列个数进行预测，范围为3376~3457个（图S2; 表S2）。采用平均核苷酸同一性（ANI）分析来评估基因组的核苷酸序列相似性，结果显示所有15测序菌株的ANI值与CFBP1430的相似性均超过99%（表S3），均被归类属于E. amylovora（图1B）。以KL20-28的基因组为例，对其特征进行了全面描述（图1C；表S5）。基因组测序比较结果显示，KL20-28和CFBP1430中匹配序列的相似性为99.87%。

测序梨火疫菌株基因注释和质粒鉴定

利用COG、KEGG、GO、Swiss-Prot、NR等5个数据库对每个菌株的预测基因进行了注释。例如在KL20-28中，共获得了3403个基因的功能注释，占总基因的99.91%（表S4）。基于COG数据库，将蛋白质分为25类（图S4）。基于GO系统（图S5；表S6）显示，细胞类别含有2133个基因，在细胞成分中占据基因数目最多。在23个KEGG聚类分析中，碳水化合物代谢含有203个基因，各类别中占据基因数最多（图S5；表S7）。此外，我们对致病性相关因子进行了初步分析，包括效应因子DspA/E、Eop1、Eop3、AvrRpt2和HopPtoC ；harpin蛋白HrpN；伴侣蛋白DspB/F；梨火疫菌素合成簇Ams（图S11）。结果表明，上述所有因子在测序15菌株中均高度保守。除了在3’端Eop1有微小的差异外，其余序列的相似性均在99.9%以上（图S11）。

pEA29（~28kb）作为保守质粒在菌株KL20-28, KL17-17, KL18-6, KL19-1, ZYGT-2, 1506-4, 31070-3, KyPL01和 KAAL06均存在，但ZY20-3-2, AKS17-22, KL21-15, KS22-2, YL16-01和YL16-15均未检测到该质粒（图S6；表S8-S9）。此外，在菌株KL17-17中鉴定到一个新的质粒pEA3-like（表S8，S10）。进一步利用VFDB（致病菌毒力因子数据库）对pEA29质粒中的毒力因子基因进行了鉴定和注释。在pEA29质粒共注释到4个相同的毒力相关基因，即甲基受体趋化蛋白cheD、tar/cheM甲基受体趋化蛋白II、pilW转座酶和Ⅳ型纤毛生物合成蛋白PilW（表S15）。然而，通过VFDB分析在质粒pEA3中未发现毒力相关基因。

与CFBP1430相比，所有已测序的菌株均存在染色体倒置现象（图S7；表S11），范围从0.4M到3.7M。此外，15菌株中共存在5087个SNPs和1996个DIPs（表S12）。在所有菌株中共鉴定出172个共享SNPs位点（共2580个）（图1D），约占总SNPs的50.7%。在DIPs方面，共确定了28个共享位点（共420个）（图1E），约占总DIPs的21%。说明DIPs在每个菌株的基因组中表现出相对较高水平的特异性。

SNPs在不同菌株中分布相对均匀，表现出一致的密度模式（图2A；表S12）。DIPs的分布没有表现出明显的模式，且DIPs的数量在不同菌株之间存在较大差异（图2B；表S12）。以CFBP1430基因组为参考，进一步分析了与中国梨火疫菌株致病性相关基因的SNP变异（表S16）。发现这些变异所在基因分别属于细菌分泌系统、鞭毛组装、细菌趋化性、植物-病原体相互作用和次生代谢产物的生物合成等方面，其中32个SNPs在15个已测序的菌株中的变异具有多样性，包括2个同义突变和30个非同义突变（表S17）。

与参考菌株CFBP1430相比，15个已测序菌株的全基因组序列相似性超过99.98%。大多数染色体经历了倒置事件的分析结果，与之前学者报道一致。虽然全基因组结构的E. amylovora菌株在中国目前几乎相同（存在SNPs或DIPs变化），但不排除在未来会因为染色体重组产生更多的菌株类型。

图1. 15梨火疫菌株全基因组测序和基因组多样性分析

（A）标注区域代表分离菌株所在地区。采样点使用https://sgis.kostat.go.kr网络工具地图绘制；（B）本研究15测序菌株与CFBP1430D的ANI值。ANI采用Jspecies软件MUMer算法计算；（C） KL20-28菌株全基因组环状图。从外圈到内圈依次代表：KL20-28基因组序列位置坐标、正链基因、负链基因、ncRNA（黑色表示tRNA，红色表示rRNA）、GC含量（红色表示平均值以上，蓝色表示平均值以下）、GC歪斜值（用于测量圆形染色体中G和C的相对含量），GC歪斜值=（G-C）/（G+C），紫色表示0以上，橙色表示0以下；（D） UpSetR图描述了每个梨火疫菌株特异和共享SNPs的数量。（E）UpSetR图描述了每个梨火疫菌株特异和共享DIPs的数量。

中国梨火疫菌株溯源分析

将本研究测序的15个菌株的信息结合NCBI中132个梨火疫菌株和10个亚洲梨火疫菌株的基因组参考信息（表S13-S14），采用diamond法对所有参与分析的菌株进行比对。进一步使用OrthoMCL法进行相似性聚类，构建进化树。结果显示，所有梨火疫菌株大致可以分为两个分支（图2C；表S13-S14）。其中第一类属于绿色表示的北美洲分支的梨火疫菌株；另一部分则是由不同地理来源组成相对复杂的分支，红色表示已测序的13个中国菌株就归在此分支中。分析表明13个国内测序进化树聚类在同一亚分支，说明菌株来源一致。并且，中国菌株与哈萨克斯坦菌株KAAL06 吉尔吉斯斯坦菌株KyPL01的聚类关系最近。其次与中东黎巴嫩菌株之间的进化关系较近，再次与欧洲各菌株的进化关系较近（图2D）。值得注意的是，韩国的所有菌株均与美国菌株UT5P4聚集在一起，与中国菌株和亚洲其他菌株进化距离相对较远。此外，本研究进一步对已测序15菌株与全球其他255菌株之间进行系统发育分析，结果与上述一致（图S8）。推测梨火疫传播到中国的途径为欧洲→中东→西亚→新疆。

自21世纪初以来，与中国相邻的国家，包括哈萨克斯坦（2008年）、吉尔吉斯斯坦（2008年）和韩国（2015年），都报道了梨火疫病的爆发。2015年，中国新疆伊犁的苹果上首次发生梨火疫病。本研究ANI分析表明与KAAL06和KyPL01相似性最高的3株中国菌株均处于在边境的伊犁哈萨克自治州（YL16-01、1506-4、31070-3）。综合以上信息可以说明：1. 虽然这15个菌株是从中国新疆、甘肃不同地区以及中亚邻国中分离出来的，但所有这些菌株都存在一个共同祖先。2. 截止目前，中国境内的梨火疫菌株是由单一源头引入的，而非多个。结合本课题组之前的报道（Sun, 2023）。中方伊犁口岸对应哈萨克斯坦阿拉木图州的“霍尔果斯”口岸，距中方口岸仅 15 公里。该口岸属公路口岸，具有国际联运地位。梨火疫病的风险在这些地区发生并扩散的可能很高，中国梨火疫的来源也很可能来自此处。

图2测序梨火疫菌株SNP/DIP染色体密度分析、溯源分析以及致病性差异分析

（A）比例尺表示全基因组大小（单位Mb）SNPs出现的次数，次数越少越偏向绿色（0则为灰色），次数越多越偏向红色（最大值50）；（B）比例尺表示全基因组大小（单位Mb）DIPs出现的次数，次数越少越偏向绿色（0则为灰色），次数越多越偏向红色（最大值50）；（C）梨火疫菌株种群系统发育分析，颜色背景代表菌株来源国家/地区。利用OrthoFinder 单拷贝同源基因分析13株中国梨火疫菌株与全球其他132梨火疫菌株以及10亚洲梨火疫菌株的进化关系；（D）矩形进化树由43梨火疫菌株基因组采用最大似然法（ML）构建（bootstrap test 基于1000个重复）；（E）测序15梨火疫菌株接种梨幼果症状分析。阴性对照用ddH2O灭菌水处理；（F）测序15梨火疫菌株接种后3、6、9d定量测定梨幼果表面坏死区域比率（计算侵染面积与总面积的比值）。所有数值均使用 ImageJ 通过颜色阈值量化得出。

中国梨火疫菌株在梨幼果中呈现致病性差异

已测序15株梨火疫菌株接种梨幼果后致病性存在显著差异（图2E；S9）。其中菌株KL17-17和KL18-6的致病性较弱，接种9天后果面坏死率分别为0.58%和4.75%，而菌株1506-4和31070-3的致病性较强，坏死率分别为67.31%和88.43%（图2F；S10）。然而，KL17-17和KL18-6菌株的致病性下降的原因和机制还有待进一步探究。本课题组后期将采用转录组、代谢组等多组学方法来揭示决定这些强致病株和弱致病株之间差异的关键因素。

方  法

菌株材料

本研究涉及的梨火疫菌株详细信息见表S1。各菌株的培养采用常规NB培养基（蛋白胨 5 g，蔗糖 10 g，酵母粉 1 g，牛肉浸膏 3 g，去离子水 1000 mL，调节 pH 6.8-7.2）。固体培养基 NA 另加 15 g/L 的琼脂粉，121 ℃灭菌 20 min。

全基因组测序和分析

梨火疫菌全基因组测序利用Genedenovo 生物技术公司（中国广州）提供的第三代PacBio和第二代Illumina 两种技术相结合的方式。依赖三代测序长读长的特点，保证更完整的基因组组装，并使用二代测序数据校正，保证组装结果更精确可靠。为保证测序质量，本研究采用细菌DNA提取试剂盒（中国天根）从菌株中提取高质量总DNA。以Qubit为标准，以Nanodrop为辅助检测（OD值260/280在1.8-2.0之间）检测每个样品的DNA质量。

用于PacBio测序的仪器主要是RSII和Sequel（美国太平洋生物科学）。利用G-tubes法构建基因组文库，将基因组DNA加工成8-10k片段。使用BluePippin（BP）核酸片段回收系统对文库进行片段选择，以去除短片段的文库分子。构建好文库后，使用Qubit进行质量检测，并利用Agilent 2100评估插入片段大小，随后用PacBio平台进行测序。基因组成分分析、功能注释、比较基因组等分析均基于基因组信息修正后的装配结果进行。然后利用dimand和COG（同源蛋白群聚类）进行数据库比较，得到相应的基因注释结果。根据注释结果对蛋白进行功能分类。

比较基因组学分析

为了比较中国梨火疫菌株与国外典型菌株之间的进化距离差异度，本试验选取CFBP 1430（isolated from France）作为参考基因组，利用共线性分析方法，使用MUMmer软件对目标基因组和参考基因组进行比对，确定了基因组间的大范围共线性关系。之后使用SyRI对区域间进行比对，确认局部位置排列关系，并从中查找易位（Translocation/Trans），倒置（Inversion/Inv）和易位+倒置（Trans+Inv）的区域。

ANI平均核苷酸同一性分析

在共线性比对完成的基础上，采用pyani法计算目标基因组与参考基因组的两两比对区域平均核苷酸同源性值（Average Nucleotide Identity，ANI）。使用JSpemies模型计算基于BLASTN的ANI值，以评估基因组序列之间的相似性。最后通过热图对相似值矩阵进行聚类和可视化。ANI值以95%作为分类阈值区分不同物种。

测序梨火疫菌株基因组结构变异检测与分析

检测了ZY20-3-2、KL17-17、KL20-28等15菌株的基因组结构变异，包括SNPs（单核苷酸多态性）和DIPs（短缺失、插入等多态性）。利用MUMmer软件检测、比较各菌株组装基因组中的SNPs和DIPs，并对其中SNPs和DIPs位置和突变结果进行统计学分析。以CFBP 1430基因组序列为参照，利用全基因组比对和TRAMS软件对SNPs、DIPs和其他多态性进行了探索。随后，使用对变异位置进行注释，并预测序列变化对相关基因功能的影响。

质粒鉴定与多样性分析

质粒测序与基因组测序同时进行。在获得各菌株中质粒的存在性和序列信息后，比较其与GenBank数据库中质粒的相似性。具体内容包括：组装不能与整个基因组连接的DNA序列序列；分析封闭的环状组装体的数量和覆盖范围；在GenBank中进行比较分析。如果在梨火疫中发现一个新质粒，则根据质粒鉴定结果进行命名。

基于梨火疫全基因组单拷贝基因构建系统发育树

对97株美国、9株亚洲、22株欧洲、2株非洲、2株大洋洲E. amylovora菌和10株E. pyrifoliae（亚洲梨火疫病）；与本试验测序的15菌株进行全基因组单拷贝基因的序列比对（表S13-S14）。除本试验15株菌株外，所有菌株的基因组均从NCBI（https://www. ncbi.nlm.nih.gov/datasets/genome）。之后利用Orthofinder2和ITOL（https:// itol.embl.de/）构建了具有1 000个bootstraps的ML系统发育树。进一步用上述相同方法对已测序15菌株与全球255菌株之间构建系统发育树并进行分析（图S8）。

E. amylovora接种梨幼果致病性分析

对未成熟的梨果实进行了致病性试验，具体方法如下：1、挑取单菌落接种到4mL NB液体培养基中，过夜摇菌；2、6000rpm，离心3 min，弃上清，用含无菌ddH2O悬浮沉淀的菌体，清洗2-3次（每次重悬后6000rpm，离心3 min）。3、梨幼果（Pyrus. cv. ‘Zaosu’）用75%酒精擦拭表面进行消毒处理，用无菌针头刺破表皮，注射 10 μL OD600 = 0.001 的菌悬液（OD600值在0.1，之后稀释100倍）。每个菌株接种3个果实。接种后的梨放置在 28°C、相对湿度为100% 的培养箱内，发病症状在第3、5、7、9天进行拍照记录。4、致病性通过坏死百分比来定量确定，即梨的表面积于坏死组织的表面积之比。