基于蜜蜂球囊菌纳米孔测序数据的基因非翻译区延长、SSR位点发掘及未注释基因和转录本鉴定

杜 宇, 付中民, 祝智威, 王 杰, 冯睿蓉, 王秀娜, 蒋海宾,范元婵, 范小雪, 熊翠玲, 郑燕珍, 徐国钧, 陈大福, 郭 睿,*

(1. 福建农林大学动物科学学院(蜂学学院), 福州 350002; 2. 福建农林大学生命科学学院, 福州 350002;3. 福建农林大学, 福建省病原真菌与真菌毒素重点实验室, 福州 350002)

蜜蜂是自然界最重要的授粉昆虫，在农业生产和生态维持方面发挥不可替代的作用(Montoya-Pfeifferetal., 2020)。此外，蜜蜂生产的蜂王浆、蜂蜜、蜂胶和蜂蜡等蜂产品具有重要的经济和药用价值(Ahmadetal., 2020)。但作为群居性昆虫，蜜蜂易遭受细菌、真菌和病毒等病原微生物的侵袭而罹患疾病。其中，蜜蜂白垩病是一种长期困扰养蜂生产的顽疾，由蜜蜂球囊菌Ascopshaeraapis侵染蜜蜂幼虫而引发(Jensenetal., 2013)。到目前为止，养蜂生产中对于白垩病仍缺乏有效的防治手段(陈大福等, 2017)。

Qin等(2006)通过对蜜蜂球囊菌0.5-1 A和A10菌株进行Sanger测序，组装了蜜蜂球囊菌的基因组草图，但作者当时仅公布了基因序列信息，并没有同时公布基因功能注释信息，导致该版本的基因组长期无法被有效利用，阻碍了蜜蜂球囊菌的进一步研究。Shang等(2016)运用二代测序技术对蜜蜂球囊菌ARSEF 7405菌株进行测序，重新组装和注释了scafford水平的蜜蜂球囊菌参考基因组(AAP 1.0)，同时公布了完整的基因序列和基因功能注释信息，为该真菌病原的组学和分子生物学研究奠定了基础。由于测序技术的限制，除人类(Audanoetal., 2019)、小鼠Musmusculus(Mouse Genome Sequencing Consortium, 2009)和黑腹果蝇Drosophilamelanogaster(Solaresetal., 2018)等极少数模式生物的基因组组装到染色体水平外，多数物种的基因组仅组装到contig或scafford水平，仍有较大的提升空间。近年来，以牛津纳米孔(Oxford Nanopore)长读段测序技术和PacBio单分子实时(single-molecule real-time, SMRT)测序技术为代表的三代测序技术逐渐兴起并快速发展。三代测序技术因具有超长读长的显著优势而能够轻松跨越重复序列，目前已成为基因组研究的利器(Luetal., 2016; Nakanoetal., 2017)。人们已利用纯三代测序或三代测序结合二代测序将人类(Pendletonetal., 2015)、跳镰猛蚁Harpegnathossaltator(Shieldsetal., 2018)和苹果Malusdomestica(Daccordetal., 2017)等物种的基因组组装到染色体水平。但目前基于三代测序技术的基因组测序成本较高，对一些基因组较大的物种进行基因组测序成本仍然高昂；对于一些经费有限的实验室，利用三代测序技术进行基因组测序还存在较大困难。与基于三代测序技术的基因组测序相比，通过三代测序技术进行转录组测序的周期较短且成本较低(Magrinietal., 2018)，因此利用三代全长转录组数据对现有的参考基因组注释进行完善是可行性较高的替代策略。近期，利用PacBio SMRT测序得到的全长转录组数据对锡兰勾虫Ancylostomaceylanicum(Magrinietal., 2018)和小麦Triticumaestivum(Dongetal., 2015)基因组注释进行完善的研究已见诸报道。然而，利用基于Nanopore测序得到的长读段数据对基因组注释进行完善的研究报道匮乏。

为开展蜜蜂球囊菌的全长转录组研究，笔者前期已利用Nanopore长读段测序技术对蜜蜂球囊菌的纯化菌丝(AaM)和纯化孢子(AaS)分别进行测序，基于高质量的测序数据构建和注释了蜜蜂球囊菌的首个全长转录组(未发表数据)；并对蜜蜂球囊菌基因的可变剪切和可变腺苷酸化进行了系统鉴定和分析(未发表数据)。本研究利用已获得的高质量Nanopore长读段测序对现有的蜜蜂球囊菌参考基因组中已注释基因进行结构优化，对未注释的简单重复序列(simple sequence repeat, SSR)位点进行鉴定，进而对未注释的新基因和新转录本进行鉴定和功能注释，并预测完整开放阅读框(open reading frame, ORF)。研究结果可为蜜蜂球囊菌参考基因组的序列和功能注释提供重要补充，也能为其他物种的基因组完善提供思路和方法借鉴。

1 材料与方法

1.1 长读段测序数据来源

前期已通过Oxford Nanopore技术对来源于纯培养的蜜蜂球囊菌AaM和AaS分别进行全长转录组测序，获得了高质量的长读段测序数据，分别测得6 321 704和6 259 727条原始读段(raw reads)，居中长度(N50)分别为1 094和1 157 bp，平均读长分别为992和1 047 bp，最大读长分别为9 421和13 060 bp；分别鉴定出9 859和16 795条非冗余全长转录本，N50分别达1 482和1 658 bp，平均长度分别为1 187和1 303 bp，最大长度分别为6 472和6 815 bp (未发表数据)。纳米孔测序原始数据已上传NCBI SRA数据库，获得BioProject号: PRJNA645872。

1.2 基因结构优化

由于软件和数据本身的局限性，导致多数基因组的基因结构信息不够精确，需要进一步优化。为最大限度对蜜蜂球囊菌的参考基因组注释进行完善，本研究将AaM和AaS的长读段测序数据混合后采用gffcompare软件(http:∥ccb.jhu.edu/software/stringtie/gffcompare.shtml)将鉴定到蜜蜂球囊菌的全长转录本与蜜蜂球囊菌参考基因组(AAP 1.0)注释的转录本进行比较，然后对基因组注释的基因结构信息进行优化。若在注释基因边界之外的区域有比对上的读段(mapped reads)支持，则将注释基因的非翻译区(untranslated region, UTR)向上游或下游延伸以修正注释基因的边界。

1.3 完整ORF的生物信息学预测

利用TransDecoder软件(http:∥transdecoder.sourceforge.net/)基于ORF长度、对数似然函数值、氨基酸序列及Pfam数据库蛋白质结构域序列的比对等信息，从蜜蜂球囊菌AaM和AaS的长读段测序混合数据鉴定到的新转录本序列中识别可靠的潜在编码区序列(coding sequence, CDS)及其对应氨基酸序列，同时预测包含起始密码子和终止密码子的完整ORF。

1.4 SSR位点的鉴定及分析

MISA软件(http: ∥pgrc.ipk-gatersleben.de/misa/)可以通过分析转录本序列鉴定出8种类型的SSR，包括单核苷酸重复(p1)、双核苷酸重复(p2)、三核苷酸重复(p3)、四核苷酸重复(p4)、五核苷酸重复(p5)、六核苷酸重复(p6)、混合SSR(c和c*)(即两个SSR之间的距离小于100 bp)，其中c类型的SSR重复序列之间包含若干个碱基，而c*类型的SSR重复序列之间没有或只有一个其他碱基(Thieletal., 2003)。从去冗余的蜜蜂球囊菌全长转录本中筛选长度在500 bp以上的全长转录本，利用MISA软件预测SSR位点，采用默认参数。

1.5 新基因和新转录本的鉴定及功能注释

通过将蜜蜂球囊菌的全长转录本与参考基因组注释的基因和转录本进行比较，鉴定现有参考基因组上未注释的新基因和新转录本。利用Blast工具将上述新基因和新转录本分别比对Nr, Swiss-Prot, Pfam, KOG, eggNOG, GO和KEGG数据库以获得相应的功能注释。

2 结果

2.1 蜜蜂球囊菌参考基因组已注释基因的5′UTR和3′UTR延长

共对蜜蜂球囊菌的9 481个基因的结构进行优化，其中5′UTR和3′UTR延长的基因分别有4 744和4 737个。部分蜜蜂球囊菌基因的结构优化信息如表1所示。

表1 蜜蜂球囊菌参考基因组已注释的10个基因的结构优化信息概要

2.2 蜜蜂球囊菌基因组中完整ORF预测

共预测出10 492个完整ORF，它们编码的氨基酸序列长度分布介于0～400 aa，其中分布在0～100 aa的ORF数量最多，为4 088个(占38.96%)；其次为分布在100～200, 200～300和300～400 aa的ORF，数量分别为3 872个(占36.90%), 1 525个(占14.53%)和595个(占5.67%)(图1)。

2.3 蜜蜂球囊菌参考基因组未注释SSR位点

本研究在24 294 167 bp的序列中共鉴定到5 286个SSR位点，含有SSR位点超过1个的基因数为1 004个，混合SSR位点有434个。此外，p1, p2, p3, p4, p5和p6的数量分别为1 870, 826, 2 398, 138, 43和11个(表2)。进一步分析发现，p3类型的SSR密度最大，达到83.72个/Mb，其次为p1, p2, c, p4, p5, c*和p6，分别达到65.20, 27.91, 15.77, 4.86, 1.48, 0.45和0.33个/Mb(图2)。

表2 蜜蜂球囊菌参考基因组中SSR位点的MISA软件分析结果

2.4 蜜蜂球囊菌参考基因组中未注释的新基因的鉴定及功能注释

图2 蜜蜂球囊菌参考基因组中不同类型SSR的密度统计

共鉴定到1 558个新基因，其中分别有1 556, 731, 330, 592, 1 177, 709和589个新基因可分别被注释到Nr, Swiss-Prot, Pfam, KOG, eggNOG, GO和KEGG数据库。Nr数据库中新基因注释数量最多的物种是蜜蜂球囊菌，其次为Polytolypahystricis和伊蒙微小菌Emmonsiaparva(图3: A)。新基因可注释到KOG数据库的25个功能类别，注释数量最多的是仅一般功能预测(general function prediction only)，其次是翻译后修饰、蛋白质转换和分子伴侣(posttranslational modification, protein turnover, chaperones)，氨基酸转运和代谢(amino acid transport and metabolism)，信号转导机制(signal transduction mechanisms)以及翻译、核糖体结构和生物合成(translation, ribosomal structure and biogenesis)等(图3: B)。此外，新基因可被注释到eggNOG数据库的25个功能类别，数量最多的为未知功能(function unknown)，其次为碳水化合物转运及代谢(carbohydrate transport and metabolism)，翻译后修饰、蛋白质转换和分子伴侣，细胞内移动、分泌和囊泡运输(intracellular trafficking, secretion, and vesicular transport)，转录(transcription)以及翻译、核糖体结构和生物合成等(图3: C)。

图3 蜜蜂球囊菌参考基因组中新基因的Nr(A)、KOG(B)和eggNOG(C)数据库注释

蜜蜂球囊菌的新基因还能被注释到GO数据库的37个功能条目，包括细胞组件(cell part)(347个)，细胞(cell)(340个)，细胞器(organelle)(262个)等细胞组分相关GO term；催化活性(catalytic activity)(328个)，结合(binding)(254个)等分子功能相关GO term；细胞进程(cellular process)(359个)，代谢进程(metabolism process)(340个)，单一组织进程(single-organism process)(245个)等生物学过程相关GO term(图4)。

此外，上述新基因还可被注释到KEGG数据库的101条通路，包括抗生素的生物合成(biosynthesis of antibiotics)(52个)，碳代谢(carbon metabolism)(29个)，氨基酸的生物合成(biosynthesis of amino acids)(27个)，剪接体(spliceosome)(23个)，糖酵解/糖异生(glycolysis/gluconeogenesis)(20个)，细胞周期-酵母(cell cycle-yeast)(20个)，核糖体(ribosome)(18个)，RNA转运(RNA transport)(18个)，泛素介导的蛋白水解(ubiquitin mediated proteolysis)(15个)以及嘌呤代谢(purine metabolism)(14个)等(图5)，条目或通路后的括号内数字代表注释的新基因占比。

图5 蜜蜂球囊菌参考基因组中新基因的KEGG数据库注释

2.5 蜜蜂球囊菌参考基因组中未注释的新转录本的鉴定及功能注释

共鉴定出14 403条新转录本，其中分别有14 376, 8 524, 7 276, 7 405, 12 035, 7 891和6 855条新转录本可被分别注释到Nr, Swiss-Prot, Pfam, KOG, eggNOG, GO和KEGG数据库。Nr数据库中新转录本注释数量最多的物种是蜜蜂球囊菌，其次为Polytolypahystricis和Helicocarpusgriseus(图6: A)。新转录本可被注释到KOG数据库的25个功能类别，包括仅一般功能预测，翻译、核糖体结构和生物合成，翻译后修饰、蛋白质转换和分子伴侣，信号转导机制，氨基酸转运和代谢，细胞内移动、分泌和囊泡运输，能量生产和转换(energy production and conversion)，RNA加工与修饰(RNA processing and modification)，未知功能以及碳水化合物转运及代谢等(图6: B)。此外，新转录本还可被注释到eggNOG数据库的25个功能类别，包括未知功能，翻译、核糖体结构和生物合成，翻译后修饰、蛋白质转换和分子伴侣，细胞内移动、分泌和囊泡运输，碳水化合物转运及代谢，氨基酸转运和代谢，转录，能量生产和转换，脂转运及代谢(lipid transport and metabolism)以及信号转导机制等(图6: C)。图6括号内数字代表注释到该条目或通路的新转录本数量和占比。

图6 蜜蜂球囊菌参考基因组中新转录本的Nr(A)、KOG(B)和eggNOG(C)数据库注释

上述新转录本还能被注释到GO数据库的44个功能条目，主要涉及细胞(4 494条)，细胞组件(4 448条)，细胞器(3 356条)，细胞膜(2 332条)，大分子复合物(macromolecular complex)(1 951条)等细胞组分相关GO term；催化活性(3 539条)，结合(2 976条)等分子功能相关GO term；细胞进程(4 281条)，代谢进程(4 055条)，单一组织进程(2 584条)等生物学过程相关GO term(图7)。

此外，这些新转录本还可被注释到KEGG数据库的119条通路，注释数量最多的是抗生素的生物合成(550条)，其次是核糖体(495条)，氨基酸的生物合成(284条)，碳代谢(275条)及剪接体(253条)等(图8)。

图8 蜜蜂球囊菌参考基因组中新转录本的KEGG数据库注释

3 讨论

目前，蜜蜂球囊菌的基因组尚未组装到染色体水平，其序列和功能注释信息仍需进一步优化完善。此前，笔者所在课题组利用Illumina测序得到的短读段数据对蜜蜂球囊菌的参考基因组注释进行完善，分别对51和50个已注释基因的5′UTR和3′UTR进行延长，鉴定出373个新基因并对部分新基因进行了功能注释(郭睿等, 2019)。Nanopore长读段测序技术作为当前主流的三代测序技术已成功应用于人类(Leaetal., 2018)、大豆Glycinemax(Flemingetal., 2018)和杆状病毒(Moldovánetal., 2018)等物种的全长转录组研究。然而对于绝大多数物种还没有基于Nanopore长读段测序数据完善基因组的研究报道。本研究利用前期已获得的Nanopore长读段测序数据对蜜蜂球囊菌的参考基因组注释进行完善，分别延长了4 744和4 737个已注释基因的5′UTR和3′UTR，数量远多于此前基于二代测序数据延长的注释基因数量，说明Nanopore长读段测序技术在优化基因结构方面具有显著优势。鉴于UTR与真核生物的基因表达调控存在密切关系(Barrettetal., 2012)，本研究中蜜蜂球囊菌基因的5′UTR和3′UTR的延长对于基因表达调控的深入研究具有重要意义。此外，本研究还预测出10 492个完整ORF，可为蜜蜂球囊菌基因全长序列的克隆及功能研究提供宝贵的参考信息。

第二代分子标记SSR是以1～6个核苷酸为重复单元组成的简单串联重复序列，具有实验操作易、重复性好和多态性高等优点(Jarne and Lagoda, 1996)。与传统方法相比，利用二代转录组数据开发SSR具有高通量的特点，使SSR的大规模开发成为现实(郭欢等, 2018; 黎东海和赵萍, 2019)。笔者所在课题组前期也基于RNA-seq数据大规模开发了中华蜜蜂Apisceranacerana(熊翠玲等, 2017)和意大利蜜蜂Apismelliferaligustica(郭睿等, 2018)的SSR。目前，已开发和利用的蜜蜂球囊菌SSR较为有限。笔者所在课题组前期利用蜜蜂球囊菌的Illumina测序数据大规模挖掘出7 968个SSR，最主要的SSR类型是三核苷酸重复(53.15%)，其次是二核苷酸重复(32.32%)和四核苷酸重复(8.46%)(李汶东等, 2017)。本研究共鉴定到5 286个SSR位点，其中最主要的类型同样为三核苷酸重复(45.37%)，其次为单核苷酸重复(35.38%)和二核苷酸重复(15.63%)，表明基于三代长读段数据和二代短读段数据开发出的SSR类型相似，但也存在一些差异。但基于三代长读段数据开发出的SSR总数明显少于基于二代短读段数据开发出的SSR总数，究其原因，可能是前期基于二代测序数据组装得到的unigene总数多达42 610个(李汶东等, 2017)，远多于蜜蜂球囊菌参考基因组包含的基因总数(6 442)，这是由于二代测序得到的片段较短(不超过300 bp)，需要利用生物信息学软件对短片段进行拼接。下一步将通过毛细管电泳和荧光标记对两种测序技术开发出的SSR进行有效性和多态性检测，进而明确何种测序技术在大规模开发SSR方面更胜一筹。

前期研究中，笔者所在课题组基于蜜蜂球囊菌的RNA-seq数据鉴定到373个新基因(郭睿等, 2019)。本研究中，共鉴定到现有参考基因组未注释的1 558个新基因，占注释基因总数的24.19%，说明基于Nanopore长读段测序数据较二代短读段测序数据在鉴定新基因方面具有显著优势。共有1 314个新基因注释到蜜蜂球囊菌，与实际情况相符；分别有11和10个新基因注释到P.hystricis和伊蒙微小菌(图3: A)，表明上述新基因在蜜蜂球囊菌与这两个物种之间具有一定的保守性。共有1 177个新基因可注释到eggNOG数据库，但注释到Swiss-Prot, Pfam, KOG, GO和KEGG数据库的新基因数量偏少，分别为731, 330, 592, 709和589个，说明这些数据库收录的蜜蜂球囊菌及近缘物种的蛋白功能注释信息较少。蜜蜂球囊菌的成熟转基因操作技术体系迄今尚未建立，导致蜜蜂球囊菌的基因功能研究严重滞后。近期，Tauber等(2019)通过体外转录合成β-葡聚糖合成蛋白编码基因以及Ras家族编码基因双链RNA(dsRNA)并处理蜜蜂球囊菌，结果显示上述dsRNA可能在蜜蜂球囊菌孢子萌发初期被吸收，相关转录本受到抑制，孢子萌发率也相应降低。该研究为蜜蜂球囊菌的基因功能研究提供了思路借鉴。现有的蜜蜂球囊菌参考基因组注释的转录本数量为6 442条，本研究鉴定到14 403条新转录本，高于注释转录本的数量，说明由于二代测序产生的短读段的限制，蜜蜂球囊菌和其他物种的大量转录本有待挖掘，Nanopore长读段测序技术在新转录本的鉴定方面大有作为。这些鉴定出的未注释的全长转录本可为基因全长序列克隆及功能研究提供可靠的数据基础。新转录本注释数量最多的物种同样是蜜蜂球囊菌，与现实情况相符，分别有70和58条新转录本注释到P.hystricis和H.griseus(图6: A)，与新基因的注释情况略有差异。此外，分别有14 376, 8 524, 7 276, 7 405, 12 035, 7 891和6 855条新转录本可被分别注释到Nr, Swiss-Prot, Pfam, KOG, eggNOG, GO和KEGG数据库，这些信息可进一步完善蜜蜂球囊菌的参考基因组注释。

综上所述，本研究利用高质量的Nanopore长读段测序数据对现有的蜜蜂球囊菌参考基因组的序列和功能注释进行了完善，为相关组学及分子生物学研究的深入开展提供了重要的参考信息，也为其他物种的基因组完善提供了方法借鉴。