基于蜜蜂球囊菌納米孔測序數據的基因非翻譯區延長、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長讀段測序數據對現有的蜜蜂球囊菌參考基因組的序列和功能注釋進行了完善，為相關組學及分子生物學研究的深入開展提供了重要的參考信息，也為其他物種的基因組完善提供了方法借鑒。