基于PacBio测序数据的蜜蜂球囊菌转录因子、融合基因及RNA编辑事件的鉴定

许雅静, 吴 鹰, 余岢骏, 孙明会, 刘佳美, 郭意龙, 徐细建,鲍佳益, 康育欣, 陈大福,2, 郭 睿,2,*, 付中民,2,*

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

蜜蜂球囊菌Ascosphaeraapis专性侵染蜜蜂幼虫而引发白垩病，该病可导致蜂群群势和生产力的急剧下降，给养蜂业造成较大损失(郭睿等, 2017a)。此前，由于参考基因组的长期缺失，蜜蜂球囊菌的组学和分子生物学研究举步维艰。Shang等(2016)利用二代测序技术测序和组装了蜜蜂球囊菌的参考基因组，并公布了完整的基因序列和功能注释，为其分子生物学和组学研究奠定了基础。笔者所在团队前期利用Illumina测序技术对蜜蜂球囊菌开展了较为系统的转录组研究(陈大福等, 2017; 郭睿等, 2017b; 张曌楠等, 2017)。

转录因子(transcription factor, TF)是指能直接或间接与真核生物基因的启动子区域中的顺式作用元件发生特异性相互作用，从而抑制或激活基因表达的一种DNA结合蛋白。作为转录水平调控基因表达的主要方式，TF可根据细胞类型、发育阶段和疾病状态来调节转录起始(Fengetal., 2020; Kimetal., 2020)。同时转录因子是由多个基因家族编码，大大增加了它们的数量并使转录调控机制复杂化，Lamber等(2018)根据DNA结合区的特征将TF分为71个家族。所谓融合基因(fusion gene)，是指将两个或多个基因的编码区首尾相连并置于同一套调控序列控制之下，从而构成的嵌合基因。基因融合可通过转录介导来实现，例如相邻基因的反式或顺式剪接、读取，其中顺式剪接发生在相同的前体mRNA分子内，而反式剪接可使两个单独的前体mRNA分子形成嵌合的非共线RNA，该RNA可以编码新蛋白质或充当非编码或调控RNA(Leietal., 2016)，也被称为转录衍生基因融合；此外还能通过各种结构重排来实现，例如基因复制或染色体易位、缺失、插入和倒位，也被称为DNA介导的基因融合(McCartneyetal., 2019)。目前，融合基因的相关研究多集中在人类重大疾病方面(Dupainetal., 2017; Yoonetal., 2019)，也有少量研究涉及细菌(Baietal., 2006)和病毒(Valencia-Herreraetal., 2019)。但真菌的融合基因研究未见报道。RNA编辑(RNA editing)指在RNA分子中特定核苷酸序列发生变化的过程，包括核苷酸的插入、缺失和替换等，这种转录后修饰不仅可以增加基因产物的多样性，还直接或间接地参与了基因表达调控，在诸多生命活动中发挥重要作用(Leongetal., 2019)。人们在锥虫的线粒体中首次鉴定到RNA编辑事件(Benneetal., 1986)，随后在人类肠道的载脂蛋白基因转录本、小麦Triticumaestivum的线粒体中也相继发现了RNA编辑现象(Powelletal., 1987; Covello and Gray, 1989)。随着高通量测序技术的日趋成熟，为RNA编辑事件的大规模挖掘提供了技术手段，研究表明RNA编辑现象在动物、植物和微生物中广泛存在(Bar-Yaacovetal., 2018; Porathetal., 2019; Tangetal., 2019)。

PacBio单分子实时(single molecule real-time, SMRT)测序技术是近年来兴起的新一代测序技术，与一代和二代测序技术相比，PacBio SMRT测序技术具有超长读长及可直接检测碱基上的化学修饰等显著优势，已成功应用于跳镰猛蚁Harpegnathossaltator(Shieldsetal., 2018)、卫氏并殖吸虫Paragonimuswestermani(Oeyetal., 2019)和小麦(Moolhuijzenetal., 2020)等动植物的相关研究。近期，笔者所在团队利用PacBio SMRT测序技术对蜜蜂球囊菌的纯化菌丝样品进行测序，构建和注释了蜜蜂球囊菌的首个高质量全长转录组，解析了基因的可变剪接与可变腺苷酸化，并鉴定和分析了长链非编码RNA(long non-coding RNA, lncRNA)和TF(Chenetal., 2020)。目前，蜜蜂球囊菌的TF相关信息十分有限，融合基因和RNA编辑的相关研究仍然缺失。本研究拟利用已获得的蜜蜂球囊菌菌丝和孢子的PacBio SMRT测序数据对TF、融合基因和RNA编辑事件进行系统鉴定和分析，以期丰富蜜蜂球囊菌的相关信息，并为深入探究它们的功能提供依据和基础。

1 材料与方法

1.1 蜜蜂球囊菌菌丝和孢子的PacBio SMRT测序数据来源

前期研究中，笔者所在团队已在实验室条件下获得蜜蜂球囊菌的纯培养，并制备了纯化的菌丝样品(AaM)和孢子样品(AaS)；此外，已利用PacBio SMRT测序技术对AaM和AaS进行测序，AaM共测得13 302 489条subreads，平均长度和居中长度(N50)分别为1 802 bp和3 077 bp；检测到464 043条环形一致性序列(circular consensus sequence, CCS)，平均长度为2 970 bp；经严格校正最终得到174 095条全长转录本，平均长度和N50分别为2 728 bp和3 543 bp(Chenetal., 2020)。AaS共测得9 911 345条subreads，平均长度和N50分别为1 742 bp和2 731 bp；检测到315 135条CCS，平均长度为2 733 bp；经校正最终得到103 845条全长转录本，平均长度和N50分别为2 502 bp和3 262 bp(未发表数据)。高质量的PacBio SMRT测序可为本研究中TF、融合基因和RNA编辑事件的鉴定与分析提供可靠的数据基础。

1.2 转录因子的鉴定及分析

利用BLASTx工具(http:∥www.ncbi.nlm.nih.gov/BLAST/)将AaM和AaS中鉴定到的全长转录本序列比对到Nr(http:∥www.ncbi.nlm.nih.gov), Swiss-Prot(http:∥us.expasy.org/sprot/)和KEGG(http:∥www.genome.jp/kegg)数据库，获得一致性最高的蛋白序列。再利用hmmscan软件将上述蛋白序列比对到Plant TFdb数据库(http:∥planttfdb.cbi.pku.edu.cn/)，从而获得TF的分类及注释信息。基因A为TOFU软件分析过程中对组成该融合基因中一个基因的编号，基因B为TOFU软件分析过程中对组成该融合基因中另一个基因的编号。

1.3 融合基因的鉴定及分析

采用TOFU分析套件中的fusion_finder.py程序进行融合基因的预测，参数设置为：dun-merge-5-shorter: turned off, min_locus_coverage: 0.05, min_locus_coverage_bp: 1, min_total_coverage: 0.99, min_dist_between_loci: 10 000。然后根据预测结果分析融合基因的序列和位置信息。

1.4 RNA编辑事件的预测及分析

首先使用SAMtools工具(Lietal., 2009)预测AaM和AaS中的RNA编辑事件。然后利用ANNOVAR软件(Wangetal., 2010)对RNA编辑事件进行注释，注释结果包括同义单核苷酸突变(synonymous single nucleotide mutation)、非同义单核苷酸突变(nonsynonymous single nucleotide mutation)和终止子获得(stop-gain)3种功能类型。最后利用基迪奥在线云平台(www.omicshare.com)的相关软件对RNA编辑位点基因进行功能和通路注释。

2 结果

2.1 蜜蜂球囊菌孢子中的转录因子

AaM中的转录因子的相关分析结果前期已另文发表(Chenetal., 2020)。在AaS中共鉴定到来源于17个TF家族的213个TF；其中来源于C2H2家族的TF数量最多，达到72个；其次为来源于bZIP和bHLH家族，TF数量均为25个；来源于TALE和Trilhelix家族的TF数量最少，均仅为1个(图1)。

图1 基于PacBio测序数据蜜蜂球囊菌孢子中鉴定到的转录因子

2.2 蜜蜂球囊菌菌丝和孢子中的融合基因

在AaM中共鉴定到921个融合基因，基因A来源于正链和负链的数量和占比分别为315个(34.2%)和606个(65.8%)，基因B来源于正链和负链的数量和占比分别为369个(40.1%)和552个(59.9%)。此外，基因A分布于39条scaffold，其中分布基因A数量最多的scaffold为AZGZ01000008.1(320个, 占34.7%)，其次为AZGZ01000001.1(108个, 占11.7%)和AZGZ01000002.1(84个, 占9.1%)；分布基因A数量最少的scaffold为AZGZ01000038.1(1个, 占0.1%)；基因B分布于57条scaffold，其中分布基因B数量最多的scaffold为AZGZ01000008.1(202个, 占21.9%)，其次为scaffold AZGZ01000014.1(37个, 占4.0%)和scaffold AZGZ01000011.1(36个, 占3.9%)，分布基因B数量最少的scaffold为AZGZ01000049.1(1个, 占0.1%)。在AaS中共鉴定到510个融合基因，基因A来源于正链和负链的数量和占比分别为160个(31.4%)和350个(68.6%)，基因B来源于正链和负链的数量和占比分别为214个(42.0%)和296个(58.0%)。此外，基因A分布于33条scaffold，其中分布基因A数量最多的scaffold为AZGZ01000008.1(200个, 占39.2%)，其次为scaffold AZGZ01000001.1(51个, 占10.0%)和scaffold AZGZ01000002.1(43个, 占8.4%)，分布基因A数量最少的scaffold为AZGZ01000034.1(1个, 占0.2%)；基因B分布位于54条scaffold，其中分布基因B数量最多的scaffold为AZGZ01000008.1(118个, 占23.1%)，其次为scaffold AZGZ01000029.1(27个, 占5.3%)和AZGZ01000013.1(15个, 占2.9%)，分布基因B数量最少的scaffold为AZGZ01000054.1(1个, 占0.2%)。进一步分析发现，AaM和AaS共有的融合基因为510个，AaM特有的融合基因为411个，而AaS无特有的融合基因。

2.3 蜜蜂球囊菌菌丝和孢子中的RNA编辑事件的鉴定与比较

在蜜蜂球囊菌中共鉴定到738次RNA编辑事件，涉及3种功能类型(图2: A)。在AaM中鉴定到的RNA编辑事件为547次，其中同义单核苷酸突变的数量最多，达到360次(占65.8%)，其次为非同义单核苷酸突变(171次, 占31.3%)和终止子获得(16次, 占2.9%)；在AaS中鉴定到的RNA编辑事件为191次，其中非同义单核苷酸突变的数量最多，为119次(占62.3%)，同义单核苷酸突变(63次, 占33%)和终止子获得(9次, 占4.7%)次之(图2: B)。

图2 基于PacBio测序数据蜜蜂球囊菌菌丝和孢子中RNA编辑事件的功能类型统计

在蜜蜂球囊菌中共鉴定到738次RNA编辑事件，涉及12种碱基替换类型(图3: A)。在AaM中鉴定到12种碱基替换类型，其中发生C->T的RNA编辑事件数量最多(158次，占28.8%)，其次为T->C(89次, 占16.2%)和C->A(58次, 占10.6%)。在AaS中鉴定到9种碱基替换类型，其中发生C->T和G->T的RNA编辑事件数量最多，均为42次(占21.9%)，其次为C->A(39次, 占20.4%)，A->T(27次, 占14.1%)和T->C(16次, 占8.3%)(图3: B)。

图3 基于PacBio测序数据蜜蜂球囊菌菌丝和孢子中RNA编辑事件碱基替换类型统计

GO数据库注释结果显示，AaM中RNA编辑位点基因可注释到生物学进程大类相关的8个功能条目，包括代谢进程(10)、细胞进程(8)和单细胞进程(6)等；分子功能大类相关的4个功能条目，包括结合(4)、催化活性(4)和结构分子活性(1)等；细胞组分大类相关的7个功能条目，包括细胞(5)、细胞部分(5)和细胞器(4)等(图5)。AaS中RNA编辑位点基因可注释到细胞进程(12)和代谢进程(11)等11个生物学进程大类相关功能条目，催化活性(13)和结合(8)等5个分子功能大类相关功能条目，细胞(10)和细胞部分(10)等8个细胞组分大类相关功能条目(图4)。括号内的数字代表注释在该条目的基因数量。

图4 基于PacBio测序数据蜜蜂球囊菌菌丝(AaM)和孢子(AaS)中RNA编辑位点基因的GO数据库注释

KEGG数据库注释结果显示，AaM中RNA编辑位点基因可注释到11条通路，包括内吞作用(2)、代谢通路(2)、磷脂酰肌醇信号系统(1)、核糖体(1)和次级代谢产物的生物合成(1)等(图5: A)；AaS中RNA编辑位点基因可注释到20条通路，包括代谢通路(3)、次级代谢产物的生物合成(2)、内吞作用(2)、丙酮酸代谢(1)和RNA转运(1)等(图5: B)。括号内的数字代表注释在该通路的基因数量。

图5 基于PacBio测序数据蜜蜂球囊菌菌丝(A)和孢子(B)中RNA编辑位点基因KEGG数据库注释

3 讨论

本研究中，在蜜蜂球囊菌孢子中鉴定到来源于17个TF家族的213个成员，总数少于菌丝中的TF；其中包含TF最多的家族也为C2H2(72)，与菌丝中的情况一致；但包含TF最少的家族为TALE(1)和Trihelix(1)(图1)，与菌丝中的情况存在差异；另外，菌丝中C2H2家族包含的TF数量多于孢子(Chenetal., 2020)。上述结果说明C2H2家族作为蜜蜂球囊菌最大的转录因子家族，其成员具有潜在的重要性。在动植物中，C2H2家族也是最重要的TF家族之一，Wu等(2019)研究阐明BmBlimp-1是家蚕Bombyxmori翅发育的重要调节因子。有研究表明该家族参与植物细胞生理状态变化期间细胞活动的调节(Maternaetal., 2006; Kametal., 2008)。棉花黄萎病菌VerticilliumdahliaeVdMsn2与酵母C2H2家族的转录因子Msn2具有同源性，有研究显示VdMsn2的缺失能够导致棉花黄萎病菌菌丝生长缓慢，隔膜和菌丝分枝增加(Tianetal., 2017)。因此，推测转录因子C2H2家族可能与蜜蜂球囊菌菌丝和孢子的生长发育和细胞活动等生物学过程密切相关。

基因的融合现象在人类疾病中广泛存在。此外，融合基因也是许多癌症的重要驱动因素，因而可作为抗癌治疗中潜在的诊断标记和治疗靶点(Yakushinaetal., 2018)。研究表明BRAF和KIAA1549基因的融合存在于80%的细胞性星形细胞瘤(Jeuken and Wesseling, 2010)；融合基因TRIM52-RACK1可促进口腔鳞状细胞癌(OSCC)细胞的增殖、迁移和侵袭，因此该融合基因有望成为治疗OSCC的有效靶点(Panetal., 2020)。相比于人类等少数哺乳动物，真菌的融合基因研究严重滞后。在蜜蜂病原中还没有相关研究报道。本研究综合蜜蜂球囊菌菌丝和孢子的PacBio SMRT测序数据共鉴定到921个融合基因，其中菌丝和孢子共有的融合基因为510个，菌丝和孢子特有的融合基因分别为411和0个。推测上述共有融合基因在蜜蜂球囊菌的不同形态中发挥重要作用，而菌丝的特有融合基因在菌丝的生长发育中扮演特殊角色，但仍需要进一步探究。进一步分析发现蜜蜂球囊菌菌丝和孢子中的融合基因均主要分布在负链，且主要由染色体AZGZ01000008.1产生，体现出菌丝和孢子中融合基因的共性。

目前仅有禾谷镰孢菌Fusariumgraminearum(Liuetal., 2016)、大孢指疫霉菌Sclerophthoramacrospora(Teichertetal., 2017)和粗糙链孢菌Neurosporacrassa(Liuetal., 2017)等少数几种真菌有RNA编辑的研究报道。在芽殖酵母Saccharomycescastellii、灵芝Ganodermalucidum和松生拟层孔菌Fomitopsispinicola中，大部分RNA编辑事件的碱基替换类型为T->C, G->A, C->T, A->G(Zhuetal., 2014; Wangetal., 2016; Wuetal., 2018)。本研究在蜜蜂球囊菌中共鉴定到738次RNA编辑事件，涉及12种碱基替换类型(图3: A)。其中最常见的碱基替换类型有4种(C->T, T->C, G->T, C->A)(图3: B)，与上述其他真菌的编辑类型并不完全一致，表明不同物种中RNA编辑事件的编辑类型具有物种特异性(马艳莉和俞嘉宁, 2009)，且RNA编辑酶复合物的识别具有特异性(Bahnetal., 2012)。此外，还发现蜜蜂球囊菌RNA编辑的功能类型中同义单核苷酸突变最为丰富(图2: A)，但在其他已报道真菌中非同义单核苷酸突变最为常见(Bianetal., 2019)，说明蜜蜂球囊菌在进化过程中需要保证基因编码蛋白的稳定性，从而利于生存和繁衍。进一步的比较分析发现，蜜蜂球囊菌菌丝中的不同类型RNA编辑事件数量普遍多于孢子中的(图2: B)。此前有研究表明RNA编辑具有阶段特异性(Liuetal., 2017)。在禾谷镰孢菌(Liuetal., 2016)和大孢指疫霉菌(Teichertetal., 2017)中，RNA编辑被证实参与了子实体的发育。在子囊菌中，RNA编辑可能是形成有性孢子所必需的(Liuetal., 2016)。推测一方面是由于孢子是一种休眠态，其新陈代谢较菌丝缓慢(Guoetal., 2018)，另一方面是由于在不同的真菌形态中RNA编辑发挥的功能不同(王淮和杨健康, 2020)。此外，AaM中C->T类型的RNA编辑事件数量最多，而AaS中G->T类型最为丰富(图3: B)。这可能是AaM中催化C->T突变的PPR蛋白含量较高，但背后的分子机理尚未明确(Chu and Wei, 2020)。本研究发现，AaM中RNA编辑位点基因可注释到代谢进程、细胞进程和催化活性等19个GO功能条目，以及内吞作用、核糖体和次级代谢产物的生物合成等11条KEGG通路；AaS中RNA编辑位点基因也可注释到细胞部分等24个功能条目以及RNA转运等20条KEGG通路(图4和5)。上述结果说明RNA编辑与蜜蜂球囊菌菌丝和孢子的生长发育、物质和能量代谢的潜在关联，值得进一步研究。