基于PacBio测序数据的蜜蜂球囊菌SNP与InDel位点发掘及分析

蔡宗兵, 张文德, 隆 琦, 余岢骏, 孙明会, 吴 鹰, 许雅静,刘佳美, 郭意龙, 徐细建, 陈大福,2,*, 郭 睿,2,*

(1. 福建农林大学动物科学学院(蜂学学院), 福州 350002; 2. 福建农林大学蜂疗研究所, 福州 350002;3. 江西省养蜂研究所, 南昌 330000)

蜜蜂是一种具有代表性的真社会性昆虫和不可替代的授粉昆虫，但因群居性而易遭受病原微生物的侵袭(梁勤和陈大福, 2009)。蜜蜂球囊菌Ascosphaeraapis是一种专性侵染蜜蜂幼虫的真菌病原，导致的白垩病可引起蜂群中成年工蜂数量和群势的骤降，据统计，白垩病可致使蜂蜜产量下降5%～37%，每年给养蜂业造成较大损失(Zaghlouletal., 2005; Aronstein and Murray, 2010)。Shang等(2016)利用二代测序技术组装和注释了蜜蜂球囊菌ARSEF 7405菌株的参考基因组，为深入开展组学和分子生物学相关研究奠定了基础。

单核苷酸多态性(single nucleotide polymorphism, SNP)和插入缺失(insertion-deletion, InDel)突变作为新型分子标记被广泛应用于动植物和微生物的基因图位克隆、基因型分型、遗传图谱绘制及遗传多样性等研究(Rafalski, 2002; Haraksingh and Snyder, 2013)。在球孢白僵菌Beauveriabassiana的研究中，Coates等(2002)鉴定了核小亚基中的InDel突变位点并将其应用于昆虫宿主中分离菌株的鉴定。在HapMap计划中，研究者通过对人类HomosapiensSNP位点进行基因分型绘制单倍型图，为精确定位复杂疾病的致病基因提供了重要信息(Gonzaga-Jaureguietal., 2012)。此外，通过SNP和InDel位点与疾病或个体敏感性、耐药性的相关性研究可指导遗传育种和药物设计，例如在结核分枝杆菌Mycobacteriumtuberculosis的研究中，Li等(2018)通过研究与利福平和异烟肼耐药性相关的SNP和InDel位点信息检测到大量耐药性相关基因。

相比于基因组测序，转录组测序的成本较低、周期较短，为高通量鉴定SNP和InDel提供了方便快捷的工具(张倩倩等, 2019)。目前，以Illumina为代表的二代测序技术已应用于人类(姜玥, 2013)、尼罗罗非鱼Oreochromisniloticus(Yáezetal., 2020)和小麦Triticumaestivum(陈广凤等, 2015)等物种的SNP和InDel研究。但蜜蜂球囊菌的分子标记相关研究较为滞后。笔者所在团队前期基于高质量的Illumina测序数据组装和注释了蜜蜂球囊菌的参考转录组，并对其SSR位点进行了大规模鉴定和分析(张曌楠等, 2017)。随着高通量测序技术的不断革新与发展，以PacBio单分子实时(single molecule real-time, SMRT)测序技术和纳米孔(nanopore)长读段测序技术为代表的第三代测序技术日趋成熟，已成功应用于中华章鱼Octopussinensis(Lietal., 2020)、拟南芥Arabidopsisthaliana(Cuiet

al., 2020)和东方蜜蜂微孢子虫Nosemaceranae(陈华枝等, 2020, 2021)等动植物和微生物的全长转录组研究。相比于二代测序，三代测序具有超长读长和直接读取碱基修饰等显著优势。目前，蜜蜂病原的三代组学研究匮乏。前期研究中，笔者所在团队利用PacBio SMRT测序技术对纯化的蜜蜂球囊菌菌丝和孢子分别进行测序，基于高质量的长读段测序数据构建和注释了蜜蜂球囊菌的首个全长转录组(杜宇等, 2021)。目前，对于包括蜜蜂球囊菌在内的蜜蜂病原，尚没有利用三代测序数据开发分子标记的研究报道。

本研究拟利用已获得的蜜蜂球囊菌菌丝PacBio SMRT测序数据发掘蜜蜂球囊菌的SNP和InDel位点，并分析其突变类型、基因组功能元件分布和密码子突变类型，进而通过功能和通路注释探讨SNP和InDel位点基因的潜在功能，以期丰富蜜蜂球囊菌的SNP和InDel位点信息，并为新型分子标记的开发和应用提供基础。

1 材料与方法

1.1 蜜蜂球囊菌的PacBio SMRT测序数据来源

利用PacBio SMRT测序技术对蜜蜂球囊菌的纯化菌丝样品进行测序，获得了高质量的全长转录组数据(Chenetal., 2020)，共测得13 302 489条subreads(约23.97 Gb)，平均读长和居中长度(N50)分别为1 802和3 077 bp；经多轮校正共得到174 095条高质量全长转录本，总碱基数为474 928 820，平均读长和N50分别为2 728和3 543 bp(Chenetal., 2020)。经严格质控的数据用于后续的SNP和InDel位点鉴定与分析。

1.2 SNP位点和InDel位点的鉴定与分析

基于蜜蜂球囊菌菌丝的PacBio SMRT测序数据，参照Li等(2009)的方法，采用SAMtools软件进行蜜蜂球囊菌菌丝中全长转录本识别。参照Wang等(2010)的方法，利用ANNOVAR软件将识别的全长转录本比对到蜜蜂球囊菌参考基因组(assembly AAP 1.0)以检测SNP位点和InDel位点，检测内容包括：(1)SNP位点的突变的类型和数量，发生转换和颠换的SNP位点的数量与比率；(2)SNP位点和InDel位点在参考基因组7种功能元件(外显子区、内含子区、基因上游区、基因下游区、基因上下游重叠区、基因间隔区、剪接区)上的分布数量和占比；(3)SNP位点和InDel位点中不同类型密码子突变的数量和占比。

1.3 SNP位点和InDel位点基因的数据库注释

使用基迪奥生物云平台(https:∥www.omicshare.com/tools/home/index/index)的相关生物信息学工具对蜜蜂球囊菌菌丝中SNP位点和InDel位点基因进行GO数据库(https:∥www.omicshare.com/tools/Home/Soft/gogsea)和KEGG数据库(https:∥www.omicshare.com/tools/Home/Soft/pathwaygsea)比对，从而获得相应的功能条目(term)和通路(pathway)的注释。

2 结果

2.1 蜜蜂球囊菌SNP位点的鉴定与分析

共鉴定到蜜蜂球囊菌的6 743个SNP位点，包括6 091个纯合位点和652个杂合位点。SNP位点的详细信息如表1所示。

表1 基于PacBio SMRT测序数据的蜜蜂球囊菌SNP位点的详细信息(仅展示10个)Table 1 Detailed information of SNP sites in Ascosphaera apis (only 10 presented) based on PacBio SMRT sequencing data

上述蜜蜂球囊菌的SNP位点的碱基突变类型共分为12种，包括C/T, G/A, A/G, T/C, G/T, A/T, C/A, T/A, T/G, A/C, C/G和G/C(图1: A)；其中最丰富的突变类型为C/T型(1 296个)，最少的突变类型为G/C(173个)(图1: A)。进一步分析发现，发生转换和颠换的SNP位点数量分别为4 887和1 856个；虽然发生转换的SNP数量较颠换更多，但是后者的SNP类型更为丰富，达到前者的2倍(图1: A)。基因组功能元件分布的统计和分析结果显示，分布在外显子区的SNP位点最多，共计3 860个(占57.24%)(图1: B)；其次是分布在基因间隔区(1 117个，占16.57%)、基因下游区(781个，占11.58%)、基因上游区(614个，占9.11%)、基因上下游重叠区(207个，占3.07%)和内含子区(160个，占2.37%)(图1: B)；分布SNP位点数量最少的是剪接区，仅有4个(占0.06%)(图1: B)。进一步对SNP位点涉及的密码子突变类型进行统计和分析，结果显示发生同义单核苷酸突变的SNP位点数量最多，达到2 892个(占74.92%)，其次是非同义单核苷酸突变(950个，占24.61%)和终止子增加(17个，占0.44%)；发生终止子减少的SNP数量最少，仅为1个(占0.03%)(图1: C)。

图1 基于PacBio SMRT测序数据的蜜蜂球囊菌SNP的突变类型(A)及基因组位置分布(B)和功能类型(C)Fig. 1 Mutation type (A), location distribution in genome (B) and functional type (C) of SNPsin Ascosphaera apis based on PacBio SMRT sequencing dataSNV: 单核苷酸突变Single nucleotide mutation.

2.2 蜜蜂球囊菌SNP位点所在基因的数据库注释

GO数据库注释结果显示，蜜蜂球囊菌SNP位点基因可注释到生物学进程、细胞组分和分子功能大类相关的34个功能条目，涉及细胞进程和代谢进程等13个生物学进程相关条目，细胞和细胞部分等8个细胞组分相关条目，催化活性和结合等4个分子功能相关条目(图2: A)。KEGG数据库注释结果显示，SNP位点基因可注释到细胞进程、环境信息处理、遗传信息处理和代谢相关的76条通路，其中注释基因数量最多的通路是代谢通路、次生代谢物的生物合成、线粒体自噬-酵母、泛素介导的蛋白水解及内吞作用等(图2: B)。

图2 基于PacBio SMRT测序数据的蜜蜂球囊菌SNP位点基因的GO(A)和KEGG(B)数据库注释Fig. 2 GO (A) and KEGG (B) database annotation of genes with SNP site in Ascosphaera apisbased on PacBio SMRT sequencing data

2.3 蜜蜂球囊菌的InDel位点的鉴定与分析

共鉴定到蜜蜂球囊菌的597个InDel位点，包括349个纯合位点和248个杂合位点。InDel位点的详细信息如表2所示。

此外，InDel位点分布最多的为基因下游区(182个，占30.49%)，其次是基因上游区(146个，占24.45%)、外显子区(109个，占18.26%)、基因间隔区(63个，占10.55%)、基因上下游重叠区(60个，占10.05%)和内含子区(37个，占6.20%)(图3: A)。进一步分析发现，最丰富密码子突变类型为非移码插入(37个，占33.94%)，其次是非移码缺失(33个，占30.28%)、移码缺失(22个，占20.18%)和移码插入(16个，占14.68%)，最少的终止子增加仅有1个，占0.92%(图3: B)。

表2 基于PacBio SMRT测序数据的蜜蜂球囊菌InDel位点的详细信息(仅展示10个)Table 2 Detailed information of InDel sites in Ascosphaera apis (only 10 presented)based on PacBio SMRT sequencing data

图3 基于PacBio SMRT测序数据的蜜蜂球囊菌中InDel位点所在基因组位置分布(A)和功能类型(B)Fig. 3 Location distribution in genome (A) and functional type (B) of genes with InDel sites in Ascosphaera apisbased on PacBio SMRT sequencing data

2.4 蜜蜂球囊菌的InDel位点基因的数据库注释

GO数据库注释结果显示，上述InDel位点基因可注释到细胞进程和代谢进程等17个生物学进程相关条目，细胞和细胞部分等12个细胞组分相关条目，催化活性和结合等10个分子功能相关条目(图4: A)。KEGG数据库注释结果显示，上述InDel位点基因还能注释到细胞进程、环境信息处理、遗传信息处理和代谢相关的87条通路，包括碳代谢、氨基酸的生物合成、细胞周期-酵母、真核生物中的核糖体生物发生和RNA转运等(图4: B)。

图4 基于PacBio SMRT测序数据的蜜蜂球囊菌InDel位点所在基因的GO(A)和KEGG(B)数据库注释Fig. 4 GO (A) and KEGG (B) database annotation of genes with InDel site in Ascosphaera apisbased on PacBio SMRT sequencing data

3 讨论

本研究利用前期已获得的蜜蜂球囊菌菌丝的PacBio SMRT测序数据，共鉴定到6 743个SNP位点，包含12种碱基突变类型，其中最丰富的类型为C/T型；发生转换和颠换的SNP分别为4 887和1 856个；SNP位点在蜜蜂球囊菌参考基因组的7种功能元件上均有分布，分布数量最多的元件为外显子；SNP位点涉及8种密码子突变类型，以同义单核苷酸突变最为丰富(图1)。此外还鉴定到597个InDel位点，在基因下游区的分布数量最多；InDel位点涉及5种密码子突变类型，其中最丰富的类型为非移码插入突变(图3)。

本研究检测到的SNP颠换类型共有8种，为转换类型的2倍，推测这是由于基因组中不同碱基的突变频率不同所致。另外，在SNP的12种碱基突变类型中，发生频率最高的为C/T突变，这可能是由于CpG二核苷酸上的胞嘧啶残基极易发生甲基化突变，进而自发脱去氨基而形成胸腺嘧啶。上述结果与蚕豆Viciafaba、水稻和椰心叶甲啮小蜂Tetrastichusbrontispae等物种的研究结论(Temnykhetal., 2001; Ocaaetal., 2015; 张洁慧等, 2021)相似。不同物种的SNP发生频率存在差异，例如在小麦条锈菌Pucciniastriiformis基因组中为1个/0.67 kb(Kiranetal., 2017)，在向日葵锈菌Pucciniahelianthi基因组中为1个/2.8 kb(王妍等, 2018)。本研究中，蜜蜂球囊菌的SNP发生频率约为1个/70 kb，远低于上述物种。Saranathan等(2017)在研究鲍氏不动杆菌Acinetobacterbaumannii时发现更高的SNP发生频率有利于菌株的快速进化，并能导致菌株产生更强的毒力和耐药性。这说明SNP发生频率的高低可能会影响病原菌对宿主的侵染力及环境适应性。需要强调的是，测序深度和检测方法等因素也可能对SNP发生频率的检出产生一定影响。

外显子与氨基酸翻译过程直接相关，其碱基的改变极易引起编码氨基酸的改变，最终影响蛋白质的功能。本研究中，蜜蜂球囊菌的SNP位点主要分布在基因组的外显子区(3 860个，占57.2%)，而在其他区域的分布相对较少；进一步分析发现对于蜜蜂球囊菌的SNP位点，同义单核苷酸突变的数量最多(2 892个，占74.9%)(图1: C)，鉴于同义单核苷酸突变并不会造成氨基酸序列的变化，上述结果表明蜜蜂球囊菌在长期的进化过程中需要保持基因组编码蛋白质的稳定，这对蜜蜂球囊菌的存活与进化至关重要。InDel可能发生移码突变，致使mRNA在翻译时遇到错误的终止密码子。Zhang ZK等(2020)将17株来自不同宿主的球孢白僵菌Beauveriabassiana与参考基因组比较，共鉴定出10 098个非同义突变基因，功能注释结果表明其中的大部分都涉及毒力蛋白的生物学功能。Tambong等(2014)以采自加拿大安大略省南部查尔斯顿湖附近林地的黏质沙雷氏菌Serratiamarcescens为研究材料，通过与黏质沙雷氏菌E-15参考菌株的金属蛋白酶基因比较后鉴定到位于该基因上的72个SNP位点和3个InDel位点，进而鉴定得到8个非同义核苷酸突变，对蛋白质变异影响的分析表明，由蛋白质结构中的非同义核苷酸突变产生的新天冬氨酸残基可能对其生物学功能产生最显著的影响。鉴于此，SNP和InDel位点对蜜蜂球囊菌的生长发育具有潜在影响。

Chen等(2021)大规模开发了玉米Zeamays中的SNP标记，GO数据库注释结果显示SNP基因可注释到大部分初级和次级代谢通路。Zhang Y等(2020)在全基因组范围鉴定了凤头鸭的SNP和InDel位点，使用GO和KEGG数据库注释后发现SNP和InDel位点基因可富集到包括骨化、软骨发育、大分子生物合成过程等条目和Hedgehog信号通路、甘油脂代谢、磷脂酰肌醇信号系统等通路。Calarco等(2018)曾检测了犬新孢子Neosporacaninum中的SNP和InDel位点，鉴定并验证了核基因组编码区的SNP热点区域，进一步的GO数据库注释结果显示SNP和InDel基因可注释到与蛋白质结合、水解酶活性、转录和翻译等相关的GO条目中。本研究中，SNP/InDel位点基因可注释到生殖进程和发育进程等生长发育相关的功能条目；此外还可注释到精氨酸生物合成和赖氨酸降解等15条氨基酸代谢相关通路，磷酸戊糖通路和半乳糖代谢等8条碳水化合物代谢相关通路，甘油磷脂代谢和甘油脂代谢等7条脂类代谢相关通路，以及三羧酸循环和氧化磷酸化等2条能量代谢相关通路(图2和图4)。以上结果表明SNP/InDel位点潜在与蜜蜂球囊菌的生长、发育、生殖及物质和能量代谢密切相关。

丝裂原活化蛋白激酶(mitogen-activated protein kinase, MAPK)级联反应在真菌的生长发育、致病性、繁殖等方面发挥关键的调控作用(Igbariaetal., 2008; 张楠等, 2017)。Di Pietro等(2003)靶向敲除了尖孢镰刀菌FusariumoxysporumMAPK基因簇中的Fmk1基因后，发现产生的突变体可以在人工培养基上生长而不能侵入番茄Lycopersiconesculentum根部，导致其丧失了对番茄的致病能力。张楠等(2017)敲除了胶孢炭疽菌ColletotrichumgloeosporioidesMAPK信号通路上的CgSho1基因，通过与野生型相比发现其营养生长缓慢，菌丝稀疏且产孢量下降，致病力也明显减弱。Jin等(2014)通过基因敲除技术对蝗绿僵菌Metarhiziumacridum中MAPK信号通路相关基因MaMk1进行靶向敲除后，发现其不能穿透蝗虫的角质层，完全丧失了致病性。本研究中，共有4个SNP位点基因和6个InDel位点基因可注释到MAPK信号通路，说明相关SNP/InDel位点与蜜蜂球囊菌的MAPK信号通路具有潜在关联，未来可参考上述前人已报道的方法尝试对注释到MAPK信号通路的SNP和InDel位点基因进行敲除以揭示其功能，进而为白垩病的治疗提供分子靶点。