椰心葉甲嚙小蜂轉錄組分析及基因功能注釋

劉華偉，李朝緒，李 芬，呂朝軍，吳少英*，覃偉權*

（1.海南大學熱帶作物學院，?？?570228；2.中國熱帶農(nóng)業(yè)科學院椰子研究所，?？?571399；3.海南大學植物保護學院，?？?570228；4.海南省院士團隊創(chuàng)新中心，海口 571339；5.海南省熱帶油料作物生物學重點實驗室，海口 571399）

椰心葉甲Brontispalongissima(Gestro)屬鞘翅目Coleoptera鐵甲科Hispidae害蟲，其成蟲和幼蟲均以棕櫚科植物未展開的心葉為食，并且?guī)缀跛械淖貦翱浦参锒伎墒芷錇楹1,2]。自2002年6月首次在海南省發(fā)現(xiàn)椰心葉甲為害，至2006年僅四年的時間，已蔓延至海南全省[3]。椰心葉甲嚙小蜂TetrastichusbrontispaeFerrière是椰心葉甲的優(yōu)勢天敵寄生蜂。目前國內(nèi)外已經(jīng)對該蜂的生物生態(tài)學、室內(nèi)繁殖和野外釋放等方面進行了相關研究[4-8]，但由于缺乏遺傳信息，椰心葉甲嚙小蜂的遺傳多樣性、基因功能等的研究還較為滯后。二代測序技術的轉錄組測序技術（RNA-seq），具有低成本、快速和高準確性等優(yōu)勢，能夠在缺乏基因信息的條件下獲得物種的代謝和生長規(guī)律，并揭示其基因與生物學特性內(nèi)在關聯(lián)，同時還可獲得物種大多數(shù)的轉錄產(chǎn)物[9]。本研究通過RNA-seq技術對兩個種群椰心葉甲嚙小蜂進行轉錄組測序，并進行了功能基因注釋和分析，旨在為椰心葉甲嚙小蜂的功能基因挖掘等分子生物學研究中提供理論基礎。

1 材料與方法

1.1 材料

為后續(xù)基于轉錄組數(shù)據(jù)分析連續(xù)自交是否會引起椰心葉甲嚙小蜂的種群退化，退化程度以及退化對其影響。本試驗中所采用的兩個種群椰心葉甲嚙小蜂分別選用中國熱帶農(nóng)業(yè)科學院椰子研究所天敵工廠提供的經(jīng)過復壯（每年一次）的種群和中國熱帶農(nóng)業(yè)科學院環(huán)境與植物保護研究所天敵工廠提供的自引進后隔離繁殖未經(jīng)復壯蜂種。兩個種群各設置3次試驗重復，每個重復由20頭椰心葉甲嚙小蜂（雌雄蜂混合的1～2日齡成蜂）混合而成。

1.2 總RNA提取，文庫構建和轉錄組測序

本試驗將收集的兩個種群椰心葉甲嚙小蜂經(jīng)液氮快速冷凍后，用干冰將樣品送至上海派森諾生物科技股份有限公司進行總RNA提取、文庫構建和轉錄組測序。

1.3 序列組裝和轉錄本拼接

RNA-seq測序完成后，統(tǒng)計總reads的數(shù)量和長度、Q30值、模糊堿基（N）所占比例以及Q20和Q30所占比例等。經(jīng)數(shù)據(jù)過濾后獲得 clean reads，繼而統(tǒng)計 clean reads的數(shù)量、總長度和占比。通過 Trinity Software（http:///trinityrnaseq.Github.io/）軟件對clean reads進行拼裝獲得轉錄本（Transcript），選取每個基因最長的轉錄本作為代表序列（Unigene）[9,10]，并進一步統(tǒng)計Transcript和Unigene的總長度、序列總數(shù)、最大長度、平均長度、N50、N50%、N90、N90%以及GC含量等指標。

1.4 Unigene功能注釋

根據(jù)基因的相似性，通過 BLAST比對工具（參數(shù)設置：E≤1e-5）將椰心葉甲嚙小蜂轉錄組獲得的Unigene與 NCBI non-redundant protein sequences（NR）數(shù)據(jù)庫、Gene Ontology（GO）數(shù)據(jù)庫、Kyoto Encyclopedia of Genes and Genome（KEGG）數(shù)據(jù)庫和 evolutionary genealogy of genes: Non-supervised Orthologous Groups（eggNOG）數(shù)據(jù)庫比對，進行功能注釋。

2 結果與分析

2.1 椰心葉甲嚙小蜂轉錄組測序數(shù)據(jù)及轉錄本拼接質(zhì)量情況

轉錄組數(shù)據(jù)集已保存在NCBI SRA數(shù)據(jù)庫中，編號：PRJNA678031。對每個樣品的原始數(shù)據(jù)（Raw data）進行統(tǒng)計（表1），兩個種群的6個樣品分別得到48514092、43095558、44266534、49405442、47396830、47326810條reads，且每個樣本的堿基總量均在6.46 Gb以上；堿基百分比Q20均大于97.47%、堿基百分比Q30均大于93.44%，含量相近；模糊堿基占比在0.001485%～0.001512%的較低水平，通過數(shù)據(jù)過濾后統(tǒng)計發(fā)現(xiàn)（表2），Clean Reads占reads總數(shù)的93.94%～94.37%。對Transcript和Unigene序列進行統(tǒng)計（表3），共得到78930條Transcript和29535條Unigene，Unigene的長度分布統(tǒng)計圖見圖1，總長度為51330466 bp，平均長度為1737.95 bp，N50值為3547 bp。

圖1 椰心葉甲嚙小蜂轉錄組Unigene組裝長度分布統(tǒng)計Fig.1 Statistics of Unigene length distribution of T.brontispae transcriptome

表1 椰心葉甲嚙小蜂轉錄組測序數(shù)據(jù)及質(zhì)量情況統(tǒng)計Table 1 Transcriptome sequencing data and quality statistics of T.brontispae

表2 椰心葉甲嚙小蜂轉錄組高質(zhì)量Reads與堿基統(tǒng)計Table 2 High quality Reads and base statistics of T.brontispae transcriptome

表3 椰心葉甲嚙小蜂轉錄組Transcript和Unigene序列統(tǒng)計表Table 3 Statistics of transcript and Unigene sequences of T.brontispae transcriptome

以上數(shù)據(jù)表明樣本的測序質(zhì)量，文庫構建質(zhì)量和拼接的完整性都較好，可以用于后續(xù)的分析。

2.2 Unigene功能注釋情況

通過Unigene和四個公共數(shù)據(jù)庫的比對并進行功能注釋（表4）。共有13796條Unigene被注釋，其在4個數(shù)據(jù)庫中的注釋如下：13401條基因（45.37%）在NR數(shù)據(jù)庫中獲得了注釋，3834條基因（12.98%）在GO數(shù)據(jù)庫中獲得了注釋，12707條基因（43.02%）在eggNOG數(shù)據(jù)庫中獲得了注釋，5999條基因（20.31%）在KEGG數(shù)據(jù)庫中獲得了注釋。

表4 椰心葉甲嚙小蜂轉錄組注釋結果匯總Table 4 Summary of annotated results of T.brontispae transcriptome

2.3 椰心葉甲嚙小蜂轉錄組Unigene的NR功能分類

通過 Blastx將椰心葉甲嚙小蜂轉錄組中的Unigene與NR數(shù)據(jù)庫進行比對，注釋到13401條Unigene。統(tǒng)計注釋結果并繪制物種分布圖（圖2），結果顯示椰心葉甲嚙小蜂與蠅繭蜂Diachasmaalloeum、佛羅里達弓背蟻Camponotusfloridanus、黑褐毛蟻Lasiusniger、赤眼蜂Trichogrammapretiosum、多胚跳小蜂Copidosomafloridanum、榕小蜂Ceratosolensolmsimarchali、麗蠅蛹集金小蜂Nasoniavitripennis都有同源序列分布，其中與麗蠅蛹集金小蜂相似序列最多，占50.34%，與榕小蜂、佛羅里達跳小蜂和赤眼蜂分別有13.3%、7.42%和5.95%的相似序列，與黑蟻、佛羅里達弓背蟻和蠅繭蜂中的同源序列較少（0.98%～1.57%）。

圖2 椰心葉甲嚙小蜂轉錄組NR注釋結果統(tǒng)計圖Fig.2 NR annotated species distribution of T.brontispae transcriptome

2.4 椰心葉甲嚙小蜂轉錄組Unigene的GO分類注釋

對椰心葉甲嚙小蜂轉錄組Unigene進行GO功能分類，共有3834條（12.98%）Unigene獲得注釋。將注釋得到的Unigene劃分為三大類（生物過程、細胞成分和分子功能）67個分支（圖3），統(tǒng)計注釋到每一類的基因數(shù)量，發(fā)現(xiàn)在生物過程類中細胞過程（1527條Unigene）的占比最大（39.83%）；在細胞成分大類中膜（1324條Unigene）的占比最大（34.53%）；在分子功能大類中結合（1733條Unigene）的占比最大（45.20%）。

圖3 椰心葉甲嚙小蜂轉錄組Unigene的GO功能分類Fig.3 GO analysis of T.brontispae transcriptome Unigene

2.5 椰心葉甲嚙小蜂轉錄組Unigene的eggNOG功能注釋

將椰心葉甲嚙小蜂轉錄組的Unigene與eggNOG數(shù)據(jù)庫進行比對，注釋到12707條Unigene，根據(jù)功能可將其劃分為26個功能區(qū)域（圖4），統(tǒng)計注釋到各類功能的基因數(shù)量，發(fā)現(xiàn)無特征基因占比最多（63.06%），其次是新陳代謝類（33.49%），細胞過程和信號傳遞，信息存儲與處理類分別占 20.82%和19.67%。

圖4 椰心葉甲嚙小蜂轉錄組Unigene的eggNOG功能分類Fig.4 eggNOG analysis of T.brontispae transcriptome Unigene

2.6 椰心葉甲嚙小蜂轉錄組Unigene的KEGG代謝途徑分類分析

將椰心葉甲嚙小蜂轉錄組的Unigene與KEGG數(shù)據(jù)庫進行對比，共有5999條Unigene獲得注釋，涉及的代謝通路可歸為五個大類（新陳代謝、遺傳信息處理、環(huán)境信息處理、細胞進程和有機系統(tǒng)）35個子類（圖 5）。統(tǒng)計注釋到各類通路的基因數(shù)量，發(fā)現(xiàn)在新陳代謝一類中獲得注釋最多的是碳水化合物代謝（299條Unigene）；遺傳信息處理一類中翻譯獲得注釋最多（361條Unigene）；環(huán)境信息處理一類中獲得注釋最多的是信號轉導（735條Unigene），細胞進程和有機系統(tǒng)兩類中獲得注釋最多的分別為運輸和分解代謝（382條Unigene）和內(nèi)分泌系統(tǒng)（388條Unigene）。

圖5 椰心葉甲嚙小蜂轉錄組KEGG注釋統(tǒng)計圖Fig.5 KEGG annotation statistical chart of T.brontispae transcriptome

3 討論

第二代高通量測序技術（Next-Generation Sequencing，NGS）因測序時間短、成本低、高準確性和所獲得數(shù)據(jù)量大等優(yōu)點，被廣泛應用于非模式生物分子生物學研究中[9,10]。椰心葉甲嚙小蜂是椰心葉甲的優(yōu)勢寄生蜂，但是由于其遺傳信息數(shù)據(jù)的缺乏，對其分子生物學的研究仍然較少。本研究旨在通過該技術分別對經(jīng)過復壯和未經(jīng)復壯兩個種群的椰心葉甲嚙小蜂進行轉錄組測序和拼裝，揭示椰心葉甲嚙小蜂整體基因表達特征，并為后續(xù)基于轉錄組數(shù)據(jù)分析復壯對椰心葉甲嚙小蜂的影響，種群遺傳結構等研究提供數(shù)據(jù)支撐。

一般來說，Q30在80%以上，N50值不小于800 bp（N50值越大表示長片段越多）就可以認為測序質(zhì)量可靠，序列組裝的完整性較好[11-13]。本研究通過無參考基因組分析椰心葉甲嚙小蜂轉錄組的特異性，序列拼接后共獲得29535條Unigene，Q30和N50值分別為93.44%和3547 bp，可以認為本研究中所得到的測序結果質(zhì)量和序列拼裝的完整性都較好，可以滿足后續(xù)分析的基本要求。

對組裝后所獲得的29535條Unigene在NR、GO、KEGG和eggNOG數(shù)據(jù)庫進行基因功能注釋，共有13796條Unigenes被注釋，仍有15739條Unigenes未被注釋。這一結果在許多生物的轉錄組測序中都有存在，如植物的彩色馬蹄蓮[10,14]，動物中的云錦杜鵑[15]和昆蟲中的麗蠅蛹集金小蜂[16]都有出現(xiàn)，可能由于椰心葉甲嚙小蜂缺乏基因組方面研究的基礎資料，使得部分Unigene在數(shù)據(jù)庫中無法得以注釋[11]，也可能是由于部分Unigene片段太短或椰心葉甲嚙小蜂中存在新的功能基因而導致的[17]。

與Nr公共數(shù)據(jù)庫比對發(fā)現(xiàn)，有13401條Unigene獲得注釋（45.37%），其中與麗蠅蛹集金小蜂相似的序列最多（50.34%）。在椰心葉甲嚙小蜂的轉錄組數(shù)據(jù)中有19.36%的Unigene屬于其他序列，這在很多物種的轉錄組測序中也有出現(xiàn)[18-20]，可能屬于椰心葉甲嚙小蜂自身特有的與大多數(shù)物種不同的序列，也可能這些基因?qū)儆诜蔷幋aRNA或不是功能基因。此外，還有54.63%的Unigene未獲得相關注釋信息，此現(xiàn)象在其他昆蟲中也存在，如二化螟盤絨繭蜂，菜粉蝶等[21,22]，其原因可能是Unigene片段過短，基因信息的暫時缺乏，或者椰心葉甲嚙小蜂中存在新的功能基因[23]。在G0數(shù)據(jù)庫中有3834條Unigene獲得注釋，僅占總Unigene的12.98%，相對注釋率較低，這進一步說明了椰心葉甲嚙小蜂基因信息的缺乏[24]，也可能與G0數(shù)據(jù)庫信息不夠完善有關[25]。eggNOG數(shù)據(jù)庫中注釋到的Unigene也較多（43.02%），除無特征序列外，注釋到與新陳代謝有關的Unigene較多（24.44%），說明椰心葉甲嚙小蜂自身具有較強的代謝能力。在 KEGG數(shù)據(jù)庫中注釋到 5999條 Unigene，同樣與新陳代謝有關的代謝通路最多，有 1842條 Unigene（30.7%），與eggNOG分析結果相一致。這些注釋結果為在分子層面研究椰心葉甲嚙小蜂提供了充足的依據(jù)。

綜上，本研究利用轉錄組測序技術，對椰心葉甲嚙小蜂復壯和未復壯兩個種群進行了轉錄組測序，并利用四大數(shù)據(jù)庫對獲取的Unigene進行了功能注釋、分類和代謝途徑預測等，為下一步深入研究椰心葉甲嚙小蜂兩個種群的差異表達基因和其他分子生物學研究提供了參考。