杜宇，祝智威，王杰，王秀娜，蔣海賓，范元嬋，范小雪，陳華枝，隆琦，蔡宗兵，熊翠玲,2，鄭燕珍，付中民,2，陳大福,2，郭睿,2

利用第三代納米孔長讀段測序技術構建和注釋蜜蜂球囊菌的全長轉錄組

杜宇1，祝智威1，王杰1，王秀娜3,4，蔣海賓1，范元嬋1，范小雪1，陳華枝1，隆琦1，蔡宗兵1，熊翠玲1,2，鄭燕珍1，付中民1,2，陳大福1,2，郭睿1,2

1福建農林大學動物科學學院（蜂學學院），福州 350002；2福建農林大學蜂療研究所，福州 350002；3福建農林大學生命科學學院，福州 350002；4福建省病原真菌與真菌毒素重點實驗室（福建農林大學），福州 350002

【】利用第三代納米孔（nanopore）長讀段測序技術對蜜蜂球囊菌（，簡稱球囊菌）的純化菌絲（Aam）和孢子（Aas）進行測序，構建和注釋球囊菌的高質量全長轉錄組。通過Oxford Nanopore PromethION平臺對Aam和Aas進行測序。利用Guppy軟件對原始讀段（raw reads）進行堿基識別（base calling），通過過濾短片段和低質量原始讀段得到有效讀段（clean reads）。通過識別兩端引物鑒定全長轉錄本序列。通過比對Nr、Swissprot、KOG、eggNOG、Pfam、GO和KEGG數據庫獲得全長轉錄本的注釋信息。分別利用CPC、CNCI、CPAT、Pfam 4種方法對長鏈非編碼RNA（long non-coding RNA，lncRNA）進行預測，取四者的交集作為高可信度的lncRNA。Aam和Aas的納米孔測序分別測得6 321 704和6 259 727條原始讀段，經質控得到5 669 436和6 233 159條有效讀段，其中包含的全長有效讀段分別為4 497 102（79.32%）和4 963 101（79.62%）條。共鑒定到9 859和16 795條非冗余全長轉錄本，N50分別為1 482和1 658 bp，平均長度分別為1 187和1 303 bp，最大長度分別為6 472和6 815 bp。Venn分析結果顯示有6 512條非冗余全長轉錄本為菌絲和孢子所共有，分別有3 347和10 283個非冗余全長轉錄本為二者特有。此外，在球囊菌菌絲和孢子中共鑒定到20 142條全長轉錄本，其中分別有20 809、11 151、17 723、12 164、11 340和9 833條全長轉錄本可注釋到Nr、KOG、eggNOG、Pfam、GO和KEGG數據庫。注釋全長轉錄本數量最多的物種是球囊菌、和莢膜組織胞漿菌（）。GO數據庫注釋結果顯示，上述全長轉錄本可注釋到45個功能條目，涉及細胞組件、細胞和細胞器等細胞組分相關條目；催化活性、結合和轉運器活性等分子功能相關條目；以及細胞進程、代謝進程和單一組織進程等生物學進程相關條目。KEGG數據庫注釋結果顯示，上述全長轉錄本還可注釋到抗生素的生物合成、核糖體、氨基酸的生物合成、碳代謝和剪接體等49條通路。此外，鑒定到648條高可信度的lncRNA，包含480條基因間區(qū)lncRNA、119條反義鏈lncRNA和49條正義鏈lncRNA。構建和注釋了球囊菌的首個高質量全長轉錄組，為探究球囊菌轉錄組的復雜性，完善參考基因組的序列和功能注釋信息以及深入開展球囊菌可變剪接體的功能研究提供了關鍵依據。

第三代高通量測序技術；納米孔測序；全長轉錄本；參考轉錄組；蜜蜂；蜜蜂球囊菌

0 引言

【研究意義】蜜蜂球囊菌（，簡稱球囊菌）是專性侵染蜜蜂幼蟲的致死性真菌病原，引發(fā)的白堊病是長期危害養(yǎng)蜂生產的頑疾，不僅可導致蜜蜂幼蟲的大量死亡，還能導致成年蜜蜂數量的銳減以及蜂群群勢和蜂產品產量的驟降[1-2]。目前，球囊菌的基因組注釋信息尚不完善，高質量參考轉錄組匱乏，嚴重限制了球囊菌的組學和分子生物學研究。利用納米孔（nanopore）長讀段測序技術構建和注釋球囊菌的全長轉錄組，有利于完善球囊菌的基因組注釋，揭示其轉錄組的復雜性，并為其組學和分子生物學研究提供重要的參考信息?！厩叭搜芯窟M展】Qin等[3]曾利用Sanger測序法對球囊菌0.5-1 A和A10菌株的菌絲進行測序，基于156 Mb的測序數據拼接出8 092條contig（總長度約為21.57 Mb），進一步組裝成1 627條scafford（總長度約為21.28 Mb）；但作者當時僅公布了基因序列信息，沒有同時公布基因功能注釋信息，導致無法利用版本基因組開展球囊菌的組學研究。直到2016年，Shang等[4]測序并公布了球囊菌ARSEF 7405菌株的完整基因組信息（assembly AAP 1.0），為其組學和分子生物學研究打下了基礎。筆者所在團隊前期對球囊菌開展了較為系統(tǒng)的轉錄組研究[5-13]。例如，基于Illimina短讀段測序得到的146 135 308條短讀段組裝出球囊菌的42 609個unigene，其中有29 316個unigene在Nr、Swiss-prot、KOG和KEGG數據庫中具有功能和通路注釋信息[5]；并在此基礎上通過比較轉錄組分析初步揭示了球囊菌對意大利蜜蜂（）幼蟲和中華蜜蜂（）幼蟲的侵染機制[6-7]。第二代測序技術雖具有通量高和準確性較高的優(yōu)勢，但因得到的讀段較短（＜300 bp），需要對短讀段進行拼接得到轉錄本，無法直接獲得轉錄本的全長信息。近年來，以PacBio單分子實時（single molecule real time，SMRT）測序和Oxford Nanopore測序為代表第三代高通量測序技術因具有超長讀長、較短測序周期及直接讀取核酸修飾等優(yōu)勢[14]，已成功應用于動植物的全長轉錄組研究[15-20]。但微生物的全長轉錄組研究相對滯后，有限的研究多集中在病毒[21-22]?！颈狙芯壳腥朦c】目前，利用納米孔測序技術對球囊菌的轉錄組研究十分滯后，高質量的球囊菌全長轉錄組缺失?！緮M解決的關鍵問題】利用納米孔長讀段測序技術對球囊菌的純化菌絲（Aam）和純化孢子（Aas）分別進行測序，將高質量的三代測序數據混合后用于構建球囊菌的全長轉錄組，并通過比對主流數據庫進行功能注釋，同時對球囊菌的長鏈非編碼RNA（long non-coding RNA，lncRNA）進行鑒定和分析。

1 材料與方法

試驗于2019年8月至2020年1月在福建農林大學動物科學學院（蜂學學院）蜜蜂保護實驗室完成。

1.1 供試生物材料

球囊菌菌株[1,5]由福建農林大學動物科學學院（蜂學學院）蜜蜂保護實驗室分離、純化和保存。參照筆者所在實驗室前期已建立的技術流程[23-24]進行球囊菌活化及菌絲和孢子純化。純化得到的純凈菌絲樣品和孢子樣品經液氮速凍后迅速轉移到-80℃超低溫冰箱保存?zhèn)溆谩?/p>

1.2 RNA提取、cDNA文庫構建及納米孔測序

（1）參照說明書步驟，利用TRizol試劑盒（Thermo Fisher公司，美國）分別提取Aam和Aas的總RNA；（2）引物退火，利用Maxima H Minus Reverse Transcriptase試劑盒（Thermo Fisher公司，美國）進行反轉錄，得到的cDNA添加switch oligo，再合成互補鏈；（3）對DNA進行損傷修復和末端修復，再利用磁珠對cDNA進行純化；（4）委托北京百邁克生物科技有限公司對上述構建好的cDNA文庫進行全長轉錄組測序，測序平臺為PromethION（Oxford Nanopore Technologies公司，英國）。

1.3 數據質控及全長轉錄本鑒定

參照Chen等[25-26]的方法進行數據質控：（1）因PromethION測序下機的原始讀段格式為二代FAST5格式，包含所有原始測序信號，故利用MinKNOW2.2軟件包中的Guppy軟件[18]對原始讀段進行堿基識別（base calling），將數據轉換為FASTQ格式；（2）進一步過濾短片段和低質量的原始讀段，得到高質量的有效讀段；（3）根據納米孔cDNA測序原理[27-28]，對有效讀段的兩端進行引物識別，兩端均識別到引物則判定其為全長轉錄本序列。

1.4 全長轉錄本的數據庫注釋

利用Blast工具將上述所有全長轉錄本比對Nr[29]、Swissprot[30]、KOG[31]、eggNOG[32]、Pfam[33]、GO（Gene Ontology）[34]和KEGG（Kyoto Encyclopedia of Genes and Genomes）[35]數據庫，獲得相應的功能和通路注釋信息。

1.5 LncRNA的鑒定和分析

通常認為lncRNA不具備蛋白編碼能力，因而可通過對轉錄本進行編碼潛能篩選判定該轉錄本是否為lncRNA[36]。分別利用CPC[37]、CNCI[38]、CPAT[39]、Pfam[33]蛋白結構域分析4種方法對上述轉錄本進行l(wèi)ncRNA的鑒定，取四者的交集作為高可信度的結果。

2 結果

2.1 球囊菌菌絲和孢子的納米孔測序數據質控

球囊菌菌絲和孢子的納米孔測序分別得到6 321 704和6 259 727條原始讀段，N50分別達到1 094和1 157 bp，平均長度分別為992和1 047 bp，最大長度分別為9 421和13 060 bp（表1）。來源于Aam和Aas的原始讀段的長度分布介于1—10 kb以上，其中分布reads數最多的長度均為1 kb（圖1-A、1-B）；原始讀段的Q值分布介于Q6—Q15，分布reads數最多的質量值分別為Q9和Q11（圖1-C、1-D）。

表1 納米孔長讀段測序產生的原始讀段信息概要

A：球囊菌菌絲測序產生的原始讀段的長度分布Length distribution of raw reads produced from sequencing of Aam；B：球囊菌孢子測序產生的原始讀段的長度分布Length distribution of raw reads produced from sequencing of Aas；C：球囊菌菌絲測序產生的原始讀段的質量值分布Quality distribution of raw reads produced from sequencing of Aam；D：球囊菌孢子測序產生的原始讀段的質量值分布Quality distribution of raw reads produced from sequencing of Aas

2.2 球囊菌全長轉錄本的鑒定和分析

對Aam和Aas樣品測序產生的原始讀段進行質控，分別得到5 669 436和6 233 159條有效讀段，其中包含的全長有效讀段分別為4 497 102（79.32%）和4 963 101（79.62%）條（表2）；全長有效讀段的長度介于1—9 kb，其中分布在1 kb的全長有效讀段數最多（圖2-A、2-B）。

表2 全長有效讀段的信息概覽

A：球囊菌菌絲測序產生的全長有效讀段Full-length clean reads yielded from sequencing of Aam；B：球囊菌孢子測序產生的全長有效讀段Full-length clean reads yielded from sequencing of Aas；C：球囊菌菌絲測序產生的全長轉錄本Full-length transcripts yielded from sequencing of Aam；D：球囊菌孢子測序產生的全長轉錄本 Full-length transcripts yielded from sequencing of Aas

進一步過濾冗余全長有效讀段，分別得到9 859和16 795條非冗余全長轉錄本，N50分別達到1 482和1 658 bp，平均長度分別達到1 187和1 303 bp，最大長度分別為6 472和6 815 bp（表3）；上述非冗余全長轉錄本的長度介于1—7 kb，其中分布在1 kb的全長轉錄本數最多（圖2-C、2-D）。進一步對Aam和Aas的非冗余全長轉錄本進行Venn分析，結果顯示有6 512個非冗余全長轉錄本為菌絲和孢子所共有，分別有3 347和10 283個非冗余全長轉錄本為二者特有（圖3）。

圖3 球囊菌菌絲和孢子全長轉錄本的Venn分析

2.3 球囊菌全長轉錄本的數據庫注釋

在球囊菌菌絲和孢子中共鑒定出20 142條全長轉錄本，數據庫注釋結果顯示，分別有20 809、11 151、17 723、12 164、11 340和9 833全長轉錄本可注釋到Nr、KOG、eggNOG、Pfam、GO和KEGG數據庫。注釋全長轉錄本數量最多的物種是球囊菌、和莢膜組織胞漿菌（）（圖4-A）。球囊菌的全長轉錄本可注釋到KOG數據庫的24個功能分類，注釋數量最多的是一般功能預測（1 658，13.29%），其次是翻譯后修飾、蛋白質翻轉和分子伴侶（1 180，9.46%），以及翻譯、核糖體結構及生物合成（1 147，9.19%）（圖4-B）。上述全長轉錄本還能注釋到eggNOG數據庫的25個功能類別，注釋數量最多的前3位分別是未知功能（7 721，42.33%），翻譯、核糖體結構和生物合成（1 285，7.04%），以及翻譯后修飾、蛋白質翻轉和分子伴侶（1 265，6.94%）（圖4-C）。

GO數據庫注釋結果顯示，上述全長轉錄本可注釋到45個功能條目，其中細胞組分大類中注釋數量最多的是細胞組件（6 423，56.64%）、細胞（6 349，55.99%）和細胞器（4 769，42.04%），分子功能大類中注釋數量最多的是催化活性（5 178，45.66%）、結合（4 315，38.05%）和轉運器活性（714，6.30%），生物學進程大類中注釋數量最多的是細胞進程（6 138，54.13%）、代謝進程（5 853，51.35%）和單一組織進程（3 769，33.24%）（圖5）。

圖4 球囊菌全長轉錄本的Nr（A）、KOG（B）和eggNOG（C）數據庫注釋

KEGG數據庫注釋結果顯示，上述全長轉錄本還能注釋到的49條通路，其中注釋數最多的通路是抗生素的生物合成（760，13.15%）、核糖體（589，10.19%）、氨基酸的生物合成（395，6.84%）、碳代謝（372，6.44%）和剪接體（335，5.80%）（圖6）。

表3 過濾掉冗余全長有效讀段的全長轉錄本概覽

1：胞外區(qū)Extracellular region；2：細胞Cell；3：擬核Nucleoid；4：細胞膜Membrane；5：病毒Virion；6：細胞膜內腔Membrane-enclosed lumen；7：大分子復合物Macromolecular complex；8：細胞器Organelle；9：胞外區(qū)Extracellular region part；10：細胞器組件Organelle part；11：病毒體組件Virion part；12：細胞膜組件Membrane part；13：細胞組件Cell part；14：超分子復合物Supramolecular complex；15：轉錄因子活性，蛋白結合Transcription factor activity, protein binding；16：核酸結合轉錄因子活性Nucleic acid binding transcription factor activity；17：催化活性Catalytic activity；18：信號傳感器活性Signal transducer activity；19：結構分子活性Structural molecule activity；20：轉運器活性Transporter activity；21：結合Binding；22：電子載體活性Electron carrier activity；23：抗氧化活性Antioxidant activity；24：金屬伴侶活性 Metallochaperone activity；25：蛋白標簽Protein tag；26：翻譯常規(guī)活性 Translation regular activity；27：分子轉換器活性Molecular transducer activity；28：分子功能調節(jié)器Molecular function regulator；29：生殖Reproduction；30：免疫系統(tǒng)進程Immune system process；31：代謝進程Metabolic process；32：細胞進程Cellular process；33：生殖進程Reproductive process；34：生物黏附Biological adhesion；35：信號Signaling；36：多細胞組織進程 Multicellular organismal process；37：發(fā)育進程Developmental process；38：生長Growth；39：單一組織進程Single-organism process；40：應激反應Response to stimulus；41：定位Localization；42：多組織進程Multi-organism process；43：生物調控Biological regulation；44：細胞成分組織或生物合成Cellular component organization or biogenesis；45：解毒作用Detoxification

2.4 球囊菌lncRNA的鑒定及分析

利用CNCI、CPC、Pfam和CPAT 4種方法分別鑒定出750、1 906、648和1 682條lncRNA，四者的交集為648個（圖7-A）；其中基因間區(qū)lncRNA（long intergenic RNA，lincRNA）、反義鏈lncRNA（anti-sense lncRNA）和正義鏈lncRNA（sense lncRNA）的數量分別為480、119和49個（圖7-B）。

圖6 球囊菌全長轉錄本的KEGG數據庫注釋

圖7 球囊菌lncRNA的數量（A）和種類（B）

3 討論

對于一個物種，高質量的參考轉錄組可為其組學和分子生物學研究提供可靠的參考信息。此前，筆者所在課題組利用二代測序技術對球囊菌開展了一系列轉錄組研究[5-13]。但對于包括球囊菌在內的所有蜜蜂病原，迄今尚沒有基于納米孔測序技術的組學研究報道。本研究利用納米孔長讀段測序技術對球囊菌的純化菌絲和純化孢子樣品進行測序，分別測得6 321 704和6 259 727條原始讀段，質控后得到5 669 436和6 233 159條有效讀段；分別鑒定到9 859和16 795條非冗余的全長轉錄本；有20 809、11 151、17 723、12 164、11 340和9 833條全長轉錄本可注釋到Nr、KOG、eggNOG、Pfam、GO和KEGG數據庫。此前，筆者所在團隊利用PacBio SMRT測序技術對球囊菌的純化菌絲進行測序，鑒定出168 740條可比對上參考基因組的全長轉錄本，并通過比對Nr、KOG、GO和KEGG數據庫對這些全長轉錄本進行了功能注釋[40]。PacBio SMRT測序技術和納米孔長讀段測序技術為當前三代測序技術的兩大代表，前者的單堿基準確性更高，在轉錄本結構分析方面更具優(yōu)勢；而后者不僅可以進行轉錄本結構分析，同時還能進行轉錄本定量分析，此外因測序設備體型較小、便于攜帶而具有更廣泛的應用性[14,41]。本研究構建的球囊菌全長轉錄組版本與此前構建的全長轉錄組版本相互補充，同時提供兩套高質量的全長轉錄本集，可為球囊菌的基因組注釋完善、轉錄組分析和基因全長序列克隆提供寶貴資源。

第一代測序技術即Sanger測序技術的讀長最長可達約1 000 bp，具有準確性高的優(yōu)點，但高成本和低通量對其大規(guī)模應用產生了極大限制。以Illumina HiSeq技術為代表的第二代測序技術雖具有通量高和準確性較高的優(yōu)勢，但因得到的讀段較短（＜300 bp），轉錄本需要由短讀段拼接而來，難以獲得轉錄本的全長信息。筆者所在團隊前期基于Illimina HiSeq測序得到的短讀段組裝出42 609條球囊菌unigene，N50和平均長度分別為1 550和966 bp[5]。本研究中，球囊菌菌絲和孢子的全長轉錄本的N50分別達到1 482和1 658 bp，與前期研究結果相比有所提升；平均長度分別為1 187和1 303 bp，顯著優(yōu)于二代測序研究結果。Workman等[15]曾利用納米孔測序技術對人類B淋巴細胞GM12878細胞系進行測序和分析，鑒定到的全長轉錄本的N50和平均長度分別為1 334和771 bp，與本研究的結果相似。以上結果表明納米孔長讀段測序技術在鑒定全長轉錄本方面具有獨特優(yōu)勢。

菌絲和孢子是真菌生長發(fā)育階段的兩種不同形態(tài)，當外界環(huán)境適宜時真菌形成菌絲進行無性或有性繁殖，當外界環(huán)境不適宜時真菌形成休眠態(tài)孢子確保生存[42]。球囊菌孢子被蜜蜂幼蟲經口攝入后，在中腸低水平萌發(fā)，并伴有菌絲的少量生長，至預蛹期中腸和后腸隔膜消失、相互連通，孢子隨食物殘渣涌入后腸并在此劇烈生長，進而致死宿主[43]。本研究發(fā)現，有6 512條非冗余全長轉錄本為菌絲和孢子所共有，分別有3 347和10 283個非冗余全長轉錄本在菌絲和孢子中特異性表達。鑒于本研究的測序材料來源于球囊菌的純培養(yǎng)，不同于處于侵染和增殖狀態(tài)的球囊菌，推測二者共有的全長轉錄本及其編碼蛋白在球囊菌生長發(fā)育的不同階段都發(fā)揮必要功能；特異性表達的全長轉錄本及其編碼蛋白在球囊菌生長發(fā)育的不同階段發(fā)揮不同作用；這些共有和特有全長轉錄本與病原的毒力和致病性存在潛在聯系，未來可結合處于侵染和增殖狀態(tài)的球囊菌全長轉錄組數據進一步挖掘和驗證。進一步對球囊菌的全長轉錄本進行數據庫注釋，分別有20 809、11 151、17 723、12 164、11 340和9 833條全長轉錄本可注釋到Nr、KOG、eggNOG、Pfam、GO和KEGG數據庫。Nr數據庫是NCBI中的非冗余蛋白質數據庫，包含Swissprot、PIR（Protein Information Resource）、PRF（Protein Research Foundation）、PDB（Protein Data Bank）蛋白質數據庫及從GenBank和RefSeq的CDS數據翻譯過來的蛋白質數據信息。本研究中，共有多達19 275（92.63%）條全長轉錄本注釋到球囊菌（圖4-A），與實際情況相符。另外，共有5 135（24.68%）條全長轉錄本在Nr數據庫中注釋為假定蛋白（hypothetical protein），說明球囊菌的基因功能注釋信息還很不完善。究其原因，一是球囊菌目前的參考基因組版本公布時間較晚（2016年），而且由于該版本是基于二代測序短讀段組裝而成，基因組注釋質量仍有較大的提升空間；二是由于相關轉基因操作技術體系尚未完全建立，球囊菌的基因功能研究相對滯后且進展緩慢。鑒于此，若要進行球囊菌的基因功能研究，首先需要對其CDS序列進行克隆。本研究鑒定出的球囊菌全長轉錄本序列可為基因全長序列克隆提供關鍵的數據基礎。Tauber等[44]利用靶向-葡聚糖合成蛋白編碼基因以及Ras家族編碼基因的雙鏈RNA（dsRNA）處理球囊菌，發(fā)現外源遺傳物質可能在病原萌發(fā)初期被吸收，進而抑制相關轉錄本，導致孢子萌發(fā)率降低。該研究為球囊菌的基因功能研究提供了方法借鑒。

目前，球囊菌的lncRNA研究滯后，相關信息匱乏。前期研究中，筆者所在團隊基于球囊菌菌絲和孢子的二代測序得到的短讀段數據，鑒定出379條lncRNA，包括123條反義鏈lncRNA、242條lincRNA、13條正義鏈lncRNA和1條內含子lncRNA[12]。本研究鑒定到480條lincRNA、119條反義鏈lncRNA和49條正義鏈lncRNA，但沒有鑒定到內含子lncRNA。本研究鑒定到的lncRNA可進一步豐富球囊菌的lncRNA信息，為今后開展lncRNA的功能研究提供數據支撐。

4 結論

構建和注釋了球囊菌的首個高質量全長轉錄組，為探究球囊菌轉錄組的復雜性、完善參考基因組的序列和功能注釋信息以及深入開展球囊菌可變剪接體的功能研究提供了關鍵依據。

[1] CHEN D F, GUO R, XU X J, XIONG C L, LIANG Q, ZHENG Y Z, LUO Q, ZHANG Z N, HUANG Z J, KUMAR D, XI W J, ZOU X, LIU M. Uncovering the immune responses oflarval gut toinfection utilizing transcriptome sequencing.Gene, 2017, 621: 40-50.

[2] GUO R, CHEN D F, DIAO Q Y, XIONG C L, ZHENG Y Z, HOU C S. Transcriptomic investigation of immune responses of thelarval gut infected by.Journal of Invertebrate Pathology, 2019, 166: 107210.

[3] QIN X, EVANS J D, ARONSTEIN K A, MURRAY K D, WEINSTOCK G M. Genome sequences of the honey bee pathogensand. Insect Molecular Biology, 2006, 15(5): 715-718.

[4] SHANG Y F, XIAO G H, ZHENG P, CEN K, ZHAN S, WANG C S. Divergent and convergent evolution of fungal pathogenicity. Genome Biology and Evolution, 2016, 8(5): 1374-1387.

[5] 張曌楠, 熊翠玲, 徐細建, 黃枳腱, 鄭燕珍, 駱群, 劉敏, 李汶東, 童新宇, 張琦, 梁勤, 郭睿, 陳大福. 蜜蜂球囊菌的參考轉錄組組裝及SSR分子標記開發(fā). 昆蟲學報, 2017, 60(1): 34-44.

ZHANG Z N, XIONG C L, XU X J, HUANG Z J, ZHENG Y Z, LUO Q, LIU M, LI W D, TONG X Y, ZHANG Q, LIANG Q, GUO R, CHEN D F.assembly of a reference transcriptome and development of SSR markers forActa Entomologica Sinica, 2017, 60(1): 34-44. (in Chinese)

[6] 陳大福, 郭睿, 熊翠玲, 梁勤, 鄭燕珍, 徐細建, 黃枳腱, 張曌楠, 張璐, 李汶東, 童新宇, 席偉軍. 脅迫意大利蜜蜂幼蟲腸道的球囊菌的轉錄組分析. 昆蟲學報, 2017, 60(4): 401-411.

CHEN D F, GUO R, XIONG C L, LIANG Q, ZHENG Y Z, XU X J, HUANG Z J, ZHANG Z N, ZHANG L, LI W D, TONG X Y, XI W J.

Transcriptomic analysis ofstressing larval gut of(Hyemenoptera: Apidae).Acta Entomologica Sinica, 2017, 60(4): 401-411. (in Chinese)

[7] 郭睿, 陳大福, 黃枳腱, 梁勤, 熊翠玲, 徐細建, 鄭燕珍, 張曌楠, 解彥玲, 童新宇, 侯志賢, 江亮亮, 刀晨. 球囊菌脅迫中華蜜蜂幼蟲腸道過程中病原的轉錄組學研究. 微生物學報, 2017, 57(12): 1865-1878.

GUO R, CHEN D F, HUANG Z J, LIANG Q, XIONG C L, XU X J, ZHENG Y Z, ZHANG Z N, XIE Y L, TONG X Y, HOU Z X, JIANG L L, DAO C. Transcriptome analysis ofstressing larval gut of. Acta Microbiologica Sinica, 2017, 57(12): 1865-1878. (in Chinese)

[8] 郭睿, 李汶東, 陳大福, 熊翠玲, 鄭燕珍, 付中民, 徐細建, 黃枳腱, 駱群. 意大利蜜蜂幼蟲腸道內球囊菌及其純培養(yǎng)的高表達基因差異分析. 微生物學通報, 2018, 45(2): 368-375.

GUO R, LI W D, CHEN D F, XIONG C L, ZHENG Y Z, FU Z M, XU X J, HUANG Z J, LUO Q. Highly-expressed gene differences betweenstressing the larval gut ofand the pure culture of. Microbiology China, 2018, 45(2): 368-375. (in Chinese)

[9] 陳大福, 王鴻權, 李汶東, 熊翠玲, 鄭燕珍, 付中民, 徐細建, 黃枳腱, 郭睿. 脅迫中華蜜蜂幼蟲腸道的球囊菌及其體外培養(yǎng)的高表達基因分析. 福建農林大學學報(自然科學版), 2017, 46(5): 562-568.

CHEN D F, WANG H Q, LI W D, XIONG C L, ZHENG Y Z, FU Z M, XU X J, HUANG Z J, GUO R. Analysis of highly expressed genes ofinfecting the gut oflarvae and itsculture. Journal of Fujian Agriculture and Forestry University (Natural Science Edition),2017, 46(5): 562-568. (in Chinese)

[10] 郭睿, 陳華枝, 童新宇, 熊翠玲, 鄭燕珍, 付中民, 解彥玲, 王海朋, 趙紅霞, 陳大福. 蜜蜂球囊菌基因結構優(yōu)化及新基因鑒定. 中國農業(yè)大學學報, 2019, 24(1): 61-68.

GUO R, CHEN H Z, TONG X Y, XIONG C L, ZHENG Y Z, FU Z M, XIE Y L, WANG H P, ZHAO H X, CHEN D F. Structural optimization of annotated genes and identification of novel genes in.Journal of China Agricultural University,2019, 24(1): 61-68. (in Chinese)

[11] 郭睿, 王海朋, 陳華枝, 熊翠玲, 鄭燕珍, 付中民, 趙紅霞, 陳大福. 蜜蜂球囊菌的microRNA鑒定及其調控網絡分析. 微生物學報, 2018, 58(6): 1077-1089.

GUO R, WANG H P, CHEN H Z, XIONG C L, ZHENG Y Z, FU Z M, ZHAO H X, CHEN D F. Identification ofmicroRNAs and investigation of their regulation networks. Acta Microbiologica Sinica, 2018, 58(6): 1077-1089. (in Chinese)

[12] GUO R, CHEN D F, XIONG C L, HOU C S, ZHENG Y Z, FU Z M, DIAO Q Y, ZHANG L, WANG H Q, HOU Z X, LI W D, KUMAR D, LIANG Q. Identification of long non-coding RNAs in the chalkbrood disease pathogen. Journal of Invertebrate Pathology, 2018, 156: 1-5.

[13] GUO R, CHEN D F, CHEN H Z, FU Z M, XIONG C L, HOU C S, ZHENG Y Z, GUO Y L, WANG H P, DU Y, DIAO Q Y. Systematic investigation of circular RNAs in, a fungal pathogen of honeybee larvae. Gene, 2018, 678: 17-22.

[14] LU H Y, GIORDANO F, NING Z M. Oxford nanopore MinION sequencing and genome assembly. Genomics Proteomics and Bioinformatics,2016, 14(5): 265-279.

[15] WORKMAN R E, TANG A D, TANG P S, JAIN M, TYSON J R, RAZAGHI R, ZUZARTE P C, GILPATRICK T, PAYNE A, QUICK J,. Nanopore native RNA sequencing of a human poly (A) transcriptome. Nature Methods, 2019, 16(12): 1297-1305.

[16] LEA W A, PARNELL S C, WALLACE D P, CALVET J P, ZELENCHUK L V, ALVAREZ N S, WARD C J. Human-specific abnormal alternative splicing of wild-typeinduces premature termination of polycystin-1. Journal of the American Society of Nephrology, 2018, 29(10): 2482-2492.

[17] CHEN S Y, DENG F L, JIA X B, LI C, LAI S J. A transcriptome atlas of rabbit revealed by PacBio single-molecule long-read sequencing. Scientific Reports,2017, 7: 7648.

[18] BAYEGA A, OIKONOMOPOULOS S, ZORBAS E, WANG Y C, GREGORIOU M E, TSOUMANI K T, MATHIOPOULOS K D, RAGOUSSIS J. Transcriptome landscape of the developing olive fruit fly embryo delineated by Oxford Nanopore long-read RNA-Seq. bioRxiv, 2018. doi: https://doi.org/10.1101/478172.

[19] CHAO Q, GAO Z F, ZHANG D, ZHAO B G, DONG F Q, FU C X, LIU L J, WANG B C. The developmental dynamics of thestem transcriptome. Plant Biotechnology Journal, 2019, 17(1): 206-219.

[20] ZHU C H, LI X F, ZHENG J Y. Transcriptome profiling using Illumina- and SMRT-based RNA-seq of hot pepper for in-depth understanding of genes involved in CMV infection. Gene, 2018, 666: 123-133.

[21] TOMBáCZ D, BALáZS Z, CSABAI Z, MOLDOVáN N, SZ?CS A, SHARON D, SNYDER M, BOLDOGK?I Z. Characterization of the dynamic transcriptome of a herpesvirus with long-read single molecule real-time sequencing. Scientific Reports, 2017, 7: 43751.

[22] TOMBáCZ D, BALáZS Z, CSABAI Z, SNYDER M, BOLDOGKOI Z. Long-read sequencing revealed an extensive transcript complexity in herpesviruses. Frontiers in Genetics, 2018, 9: 259.

[23] 陳華枝, 祝智威, 蔣海賓, 王杰, 范元嬋, 范小雪, 萬潔琦, 盧家軒, 熊翠玲, 鄭燕珍, 付中民, 陳大福, 郭睿. 蜜蜂球囊菌菌絲和孢子中微小RNA及其靶mRNA的比較分析. 中國農業(yè)科學, 2020, 53(17): 3606-3619.

CHEN H Z, ZHU Z W, JIANG H B, WANG J, FAN Y C, FAN X X, WAN J Q, LU J X, XIONG C L, ZHENG Y Z, FU Z M, CHEN D F, GUO R.Comparative analysis of microRNAs and corresponding target mRNAs inmycelium and spore. Scientia Agricultura Sinica, 2020, 53(17): 3606-3619. (in Chinese)

[24] 陳華枝, 王杰, 祝智威, 蔣海賓, 范元嬋, 范小雪, 萬潔琦, 盧家軒, 鄭燕珍, 付中民, 徐國鈞, 陳大福, 郭睿. 蜜蜂球囊菌菌絲和孢子中長鏈非編碼RNA的比較及其潛在功能分析. 中國農業(yè)科學, 2021, 54(2): 435-448.

CHEN H Z, WANG J, ZHU Z W, JIANG H B, FAN Y C, FAN X X, WAN J Q, LU J X, ZHENG Y Z, FU Z M, XU G J, CHEN D F, GUO R. Comparison and potential functional analysis of long non-coding RNAs betweenmycelium and spore.Scientia Agricultura Sinica, 2021, 54(2): 435-448. (in Chinese)

[25] CHEN H Z, FAN X X, DU Y, FAN Y C, WANG J, JIANG H B, XIONG C L, ZHENG Y Z, CHEN D F, GUO R. Nanopore-based long-read transcriptome data of-infected and un-infected western honeybee workers’ midguts. bioRxiv, 2020. doi: https://doi.org/10.1101/2020.03.21.001958.

[26] DU Y, FAN Y C, CHEN H Z, WANG J, XIONG C L, ZHENG Y Z, CHEN D F, GUO R. A full-length transcriptome dataset of normal and-challenged midgut tissues of eastern honeybee workers. bioRxiv, 2020. doi: https://doi.org/10.1101/2020.03.18. 997981.

[27] JENJAROENPUN P, WONGSURAWAT T, PEREIRA R, PATUMCHAROENPOL P, USSERY D W, NIELSEN J, NOOKAEW I. Complete genomic and transcriptional landscape analysis using third-generation sequencing: a case study ofCEN.PK113-7D. Nucleic Acids Research, 2018, 46(7): e38.

[28] BOLDOGKOI Z, MOLDOVAN N, BALAZS Z, SNYDER M, TOMBACZ D. Long-read sequencing-A powerful tool in viraltranscriptome research. Trends in Microbiology, 2019, 27(7): 578-592.

[29] 鄧泱泱, 荔建琦, 吳松鋒, 朱云平, 陳耀文, 賀福初. nr數據庫分析及其本地化. 計算機工程, 2006, 32(5): 71-73, 76.

DENG Y Y, LI J Q, WU S F, ZHU Y P, CHEN Y W, HE F C. Integrated nr database in protein annotation system and its localization. Computer Engineering, 2006, 32(5): 71-73, 76. (in Chinese)

[30] The Uniprot Consortium. UniProt: the universal protein knowledgebase. Nucleic Acids Research, 2017, 45(D1): D158-D169.

[31] KOONIN E V, FEDOROVA N D, JACKSON J D, JACOBS A R, KRYLOV D M,MAKAROVA K S, MAZUMDER R, MEKHEDOV S L, NIKOLSKAYA A N, RAO B S,. A comprehensive evolutionary classification of proteins encoded in complete eukaryotic genomes. Genome Biology, 2004, 5(2): R7.

[32] POWELL S, FORSLUND K, SZKLARCZYK D, TRACHANA K, ROTH A, HUERTA-CEPAS J, GABALDóN T, RATTEI T, CREEVEY C, KUHN M, JENSEN L J, VON MERING C, BORK P. eggNOG v4.0: nested orthology inference across 3686 organisms. Nucleic Acids Research, 2014, 42(Database issue): D231-D239.

[33] FINN R D, BATEMAN A, CLEMENTS J, COGGILL P, EBERHARDT R Y, EDDY S R, HEGER A, HETHERINGTON K, HOLM L, MISTRY J, SONNHAMMER E L L, TATE J, PUNTA M. Pfam: the protein families database. Nucleic Acids Research, 2014, 42(Database issue): D222-D230.

[34] ASHBURNER M, BALL C A, BLAKE J A, BOTSTEIN D, BUTLER H, CHERRY J M, DAVIS A P, DOLINSKI K, DWIGHT S S, EPPIG J T,. Gene ontology: tool for the unification of biology. Nature Genetics, 2000, 25(1): 25-29.

[35] KANEHISA M, GOTO S, KAWASHIMA S, OKUNO Y, HATTORI M. The KEGG resource for deciphering the genome. Nucleic Acids Research, 2004, 32(Database issue): D277-D280.

[36] 熊翠玲, 耿四海, 王心蕊, 劉思亞, 陳大福, 鄭燕珍, 付中民, 杜宇, 王海朋, 陳華枝, 周丁丁, 郭睿. 意大利蜜蜂工蜂中腸的長鏈非編碼RNA的預測、分析及鑒定. 應用昆蟲學報, 2018, 55(6): 1034-1044.

XIONG C L, GENG S H, WANG X R, LIU S Y, CHEN D F, ZHENG Y Z, FU Z M, DU Y, WANG H P, CHEN H Z, ZHOU D D, GUO R.Prediction, analysis and identification of long non-coding RNA in the midguts ofworkers. Chinese Journal of Applied Entomology, 2018, 55(6): 1034-1044. (in Chinese)

[37] KONG L, ZHANG Y, YE Z Q, LIU X Q, ZHAO S Q, WEI L, GAO G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Research, 2007, 35(Web Server issue): W345-W349.

[38] SUN L, LUO H T, BU D C, ZHAO G G, YU K T, ZHANG C H, LIU Y N, CHEN R S, ZHAO Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Research, 2013, 41(17): e166.

[39] WANG L, PARK H J, DASARI S, WANG S, KOCHER J P, LI W. CPAT: coding-potential assessment tool using an alignment-free logistic regression model. Nucleic Acids Research, 2013, 41(6): e74.

[40] CHEN D F, DU Y, FAN X X, ZHU Z W, JIANG H B, WANG J, FAN Y C, CHEN H Z, ZHOU D D, XIONG C L, ZHENG Y Z, XU X J, LUO Q, GUO R. Reconstruction and functional annotation offull-length transcriptome via PacBio single-molecule long-read sequencing. bioRxiv, 2019. doi: https://doi.org/10.1101/770040.

[41] MAGI A, SEMERARO R, MINGRINO A, GIUSTI B, D’AURIZIO R. Nanopore sequencing data analysis: state of the art, applications and challenges. Briefings in Bioinformatics, 2018, 19(6): 1256-1272.

[42] Aronstein K A, Murray K D. Chalkbrood disease in honey bees. Journal of Invertebrate Pathology, 2010, 103(Suppl.1): S20-S29.

[43] 李江紅, 鄭志陽, 陳大福, 梁勤. 影響蜜蜂球囊菌侵染蜜蜂幼蟲的因素及侵染過程觀察. 昆蟲學報, 2012, 55(7): 790-797.

LI J H, ZHENG Z Y, CHEN D F, LIANG Q. Factors influencinginfection on honeybee larvae and observation on the infection process. Acta Entomologica Sinica, 2012, 55(7): 790-797. (in Chinese)

[44] TAUBER J P, EINSPANIER R, EVANS J D, MCMAHON D P. Co-incubation of dsRNA reduces proportion of viable spores of, a honey bee fungal pathogen. Journal of Apicultural Research, 2020, 59(5): 791-799.

Construction and Annotation ofFull-Length Transcriptome Utilizing Nanopore Third-Generation Long-Read Sequencing Technology

DU Yu1, ZHU ZhiWei1, WANG Jie1, WANG XiuNa3,4, JIANG HaiBin1, FAN YuanChan1, FAN XiaoXue1, CHEN HuaZhi1, LONG Qi1, CAI ZongBing1, XIONG CuiLing1,2, ZHENG YanZhen1, FU ZhongMin1,2, CHEN DaFu1,2, GUO Rui1,2

1College of Animal Sciences (College of Bee Science), Fujian Agriculture and Forestry University, Fuzhou 350002;2Apitherapy Research Institution, Fujian Agriculture and Forestry University, Fuzhou 350002;3College of Life Sciences, Fujian Agriculture and Forestry University, Fuzhou 350002;4Key Laboratory of Pathogenic Fungi and Mycotoxins of Fujian Province (Fujian Agriculture and Forestry University), Fuzhou 350002

【】Purified mycelia sample (Aam) and spore sample (Aas) were sequenced using third-generation nanopore long-read sequencing technology, followed by construction and annotation of high-quality full-length transcriptome.【】Aam and Aas were respectively sequenced using Oxford Nanopore PromethION platform. Guppy software was used to conduct base calling of raw reads. Clean reads were obtained after filtering out short fragments and low-quality raw reads. Full-length transcripts were identified by recognizing primers at both ends of clean reads. Full-length transcripts were aligned to Nr, Swissprot, KOG, eggNOG, Pfam, GO and KEGG databases to gain corresponding annotations. Four approaches such as CPC, CNCI, CPAT, and Pfam were used to predict lncRNAs, and the intersection was deemed to be high-reliability lncRNAs.【】In total, 6 321 704 and 6 259 727 raw reads were yielded from nanopore sequencing of Aam and Aas, and after quality control, 5 669 436 and 6 233 159 clean reads were obtained, including 4 497 102 (79.32%) and 4 963 101 (79.62%) full-length clean reads. Additionally, 9 859 and 16 795 non-redundant full-length transcripts were identified, with a N50 of 1 482and 1 658bp, an average length of 1 187 and 1 303 bp, and a maximum length of 6 472 and 6 815 bp, respectively. Venn analysis showed that 6 512 non-redundant full-length transcripts were shared by Aam and Aas, while 3 347 and 10 283 ones were specific for Aam and Aas, respectively.Besides, a total of 20 142 full-length transcripts were identified in Aam and Aas, among them 20 809, 11 151, 17 723, 12 164, 11 340 and 9 833 full-length transcripts could be annotated to Nr, KOG, eggNOG, Pfam, GO and KEGG databases, respectively. Most of full-length transcripts were annotated to,and. Moreover, GO database annotation demonstrated that the above-mentioned full-length transcripts could be annotated to 45 functional terms, involving in cell component-associated terms such as cell part, cell and organelle; molecular function-associated terms such as catalytic activity, binding and transporter activity; and biological process-associated terms such as cellular processes, metabolic processes and single-organism processes. KEGG database annotation indicated that these full-length transcripts could be annotated to 49 pathways, including biosynthesis of antibiotics, ribosome, biosynthesis of amino acid, carbon metabolism, spliceosome and so on. In addition, 648 lncRNAs were identified, including 480 long intergenic RNAs (lincRNAs), 119 anti-sense lncRNAs and 49 sense lncRNAs. 【】The first high-quality full-length transcriptome was constructed and annotated in this work, which offers a key basis for exploration of the complexity oftranscriptome, improvement of sequence and functional annotation of reference genome and further study on isoforms’ function of.

third-generation high-throughput sequencing technology; nanopore sequencing; full-length transcript; reference transcriptome; honeybee;

10.3864/j.issn.0578-1752.2021.04.017

2020-05-04；

2020-05-22

國家現代農業(yè)產業(yè)技術體系建設專項（CARS-44-KXJ7）、福建省自然科學基金（2018J05042）、福建農林大學杰出青年科研人才計劃（xjq201814）、福建省病原真菌與真菌毒素重點實驗室開放課題（郭睿）、江西省蜜蜂生物學與飼養(yǎng)重點實驗室開放基金（JXKLHBB-2020-04）、福建農林大學優(yōu)秀碩士學位論文資助基金（杜宇）

杜宇，E-mail：m18505700830@163.com。祝智威，E-mail：zzw15235470398@163.com。杜宇和祝智威為同等貢獻作者。通信作者陳大福，E-mail：dfchen826@fafu.edu.cn。通信作者郭睿，E-mail：ruiguo@fafu.edu.cn

（責任編輯 岳梅）