徐志文 任雪敏 王 俊 張枝全 劉乃勇 閻 偉 朱家穎

(1. 西南林業(yè)大學(xué)云南省森林災(zāi)害預(yù)警與控制重點實驗室，云南 昆明 650224；2. 中國熱帶農(nóng)業(yè)科學(xué)院椰子研究所，海南 文昌 571339)

椰子織蛾 (Opisinaarenosella) 又名椰子木蛾和椰蛀蛾等，為鱗翅目 (Lepidoptera)、織蛾科 (Oecophoridae) 昆蟲。它的原產(chǎn)地為南亞的印度和斯里蘭卡[1]，后傳播到南亞的其他地區(qū)和東南亞地區(qū)，包括巴基斯坦、泰國、馬來西亞、緬甸等多個國家和地區(qū)[2]，并造成嚴重危害。隨著頻繁的國際貿(mào)易發(fā)展，該蟲首次于2013年8月在我國海南省萬寧市被發(fā)現(xiàn)，并在該地區(qū)造成一定危害[3]。我國于2014年將其列為檢疫性入侵害蟲，現(xiàn)主要分布于我國的海南、廣西、廣東和福建[4]，其他地區(qū)還未見報道。由于椰子織蛾的環(huán)境適應(yīng)性強，加之寄主種類繁多(主要為棕櫚科植物)[5]，因此該害蟲在我國的潛在分布范圍廣[6]。一旦椰子織蛾在我國大面積傳播擴散，將會造成難以預(yù)估的經(jīng)濟損失。

目前，國內(nèi)外針對椰子織蛾的研究主要限于形態(tài)特征、生物學(xué)特性、危害和防治等方面[2,4,7]，甚少有涉及分子層面的研究報道。截至目前，僅Rajesh等對來自5個地區(qū)的18份椰子織蛾幼蟲線粒體細胞色素氧化酶亞基Ⅰ基因進行測序，分析揭示了它們的遺傳多樣性[8]。在NCBI中，現(xiàn)僅檢索到椰子織蛾的基因序列不到30條，也未見轉(zhuǎn)錄組學(xué)或其他組學(xué)手段應(yīng)用于該害蟲的研究中。隨著快速、經(jīng)濟、高效的二代高通量測序技術(shù)的發(fā)展和成熟，且轉(zhuǎn)錄組具有不需要基因組為背景、準確的大數(shù)據(jù)量、低冗余性和檢測低豐富度基因等優(yōu)點[9]。轉(zhuǎn)錄組學(xué)手段被廣泛的應(yīng)用于昆蟲分子生物學(xué)研究領(lǐng)域，其在挖掘基因、探究基因功能和基因差異表達、分析代謝途徑、開發(fā)單核苷酸多態(tài)性和微衛(wèi)星分子標記等方面具有其特有的優(yōu)勢[10-12]。據(jù)此，為了得到比較全面的椰子織蛾轉(zhuǎn)錄組數(shù)據(jù)，本研究基于Illumina HiSeqTM2000測序平臺對該害蟲幼蟲和成蟲的混合樣品開展轉(zhuǎn)錄組測序及分析，測序獲得的轉(zhuǎn)錄組數(shù)據(jù)填補了椰子織蛾轉(zhuǎn)錄信息的空白，將為后期挖掘分子標記、研究椰子織蛾的遺傳多樣性及其功能基因提供參考。

1 材料與方法

1.1 昆蟲材料

本實驗的供試蟲源椰子織蛾幼蟲、雌雄成蟲，均由中國海南熱帶農(nóng)業(yè)科學(xué)院椰子研究所提供。室內(nèi)飼養(yǎng)條件為：溫度 (28 ± 1)℃，飽和濕度，光照12 h/d。選擇中齡幼蟲、雌雄成蟲各2頭，混合研磨后轉(zhuǎn)移至盛有Trizol試劑的離心管中，保存在-80 ℃超低溫冰箱用于轉(zhuǎn)錄組測序。

1.2 樣品總RNA提取與質(zhì)檢

參照Trizol Reagent說明書提取樣品總RNA，分別利用1%瓊脂糖凝膠電泳、Agilent 2100 (Agilent Technologies, CA, USA)、Nanodrop分光光度計 (IMPLEN, CA, USA) 和Qubit RNA Assay Kit (Life Technologies, CA, USA) 檢測分析總RNA質(zhì)量、完整性和純度 (OD260與OD280比值)，并精確定量。

1.3 構(gòu)建cDNA文庫和上機測序

用帶有Oligo (dT) 的磁珠富集mRNA。在高溫條件下，加入fragmentation buffer將mRNA打斷成短片段，以此為目的模板，用六堿基隨機引物 (random hexamers) 和逆轉(zhuǎn)錄酶 (RNase H-) 合成一鏈cDNA，之后加入緩沖液、dNTPs、DNA polymerase I和RNase H合成二鏈cDNA，再用AMPure XP beads純化雙鏈cDNA。純化的雙鏈cDNA先進行末端修復(fù)、加A尾并連接測序接頭，再用AMPure XP beads選擇大小在150～200 bp的片段。最后進行PCR擴增，后用AMPure XP beads純化PCR產(chǎn)物，得到最終的文庫。文庫質(zhì)檢合格和定量后，將cDNA文庫置于Illumina HiSeqTM2000平臺上機測序。

1.4 轉(zhuǎn)錄組拼接和功能基因注釋

測序的原始數(shù)據(jù) (raw read) 經(jīng)質(zhì)量評估合格后，去除序列接頭、poly-N的比例大于10%和低質(zhì)量的read，得到高質(zhì)量序列 (clean read)，并計算Q20、Q30、GC含量和序列重復(fù)水平。利用Trinity軟件將高質(zhì)量的clean read進行拼接[13]，將拼接得到的unigene運用BLAST軟件在七大數(shù)據(jù)庫：Nr (NCBI non-redundant protein sequences)；Nt (NCBI non-redundant nucleotide sequences)；Pfam (Protein family)；KOG/COG (Clusters of Orthologous Groups)；Swiss-Prot (A manually annotated and reviewed protein sequence database)； KEGG (Kyoto Encyclopedia of Genes and Genomes)；GO (Gene Ontology) 中進行比對 (E-value < 1.0 × 10-5)，獲得功能基因注釋信息。

1.5 微衛(wèi)星鑒定

基于拼接得到的椰子織蛾unigene序列，利用MISA軟件 (http://pgrc.ipk-gatersleben.de/misa/misa.html) 進行微衛(wèi)星鑒定。搜索參數(shù)設(shè)置為：單堿基重復(fù)次數(shù) ≥10，二堿基重復(fù)次數(shù) ≥6，三堿基重復(fù)次數(shù) ≥5，四堿基重復(fù)次數(shù) ≥5，五堿基重復(fù)次數(shù) ≥5，六堿基重復(fù)次數(shù) ≥5，相鄰2個微衛(wèi)星重復(fù)單元間的最大間隔長度為100 bp。

2 結(jié)果與分析

2.1 測序質(zhì)控分析和轉(zhuǎn)錄組組裝

測序獲得54 571 146條原始序列，經(jīng)過濾去除低質(zhì)量序列后得到50 606 604條clean reads，共計5.06 Gb。樣本的GC含量為49.85%，Q20為96.86%，堿基Q30為89.44%。經(jīng)Trinity軟件組裝后 (表1)，共得到51 578條轉(zhuǎn)錄本 (transcript) 和37 416條unigenes。轉(zhuǎn)錄本平均長度為1 035 bp，N50為1 939 bp，在200～300 bp長度之間的轉(zhuǎn)錄本數(shù)量最大，占總比約26.99%，而長度大于2 000 bp的占比最小，占13.84%。Unigene平均長度為803 bp，N50為1 415 bp，在200～300 bp長度之間的unigene數(shù)量最多，占總比約32.91%，長度大于2 000 bp的數(shù)量最少，占8.91%。

表1 椰子織蛾轉(zhuǎn)錄組組裝結(jié)果統(tǒng)計Table 1 Statistics of transcriptome assembly for O.arenosella

2.2 Unigene功能注釋結(jié)果

Unigene在七大數(shù)據(jù)庫中Blast比對得到功能注釋信息 (表2)。組裝獲得的37 416條unigenes被成功注釋19 154條，注釋率為51.19%。在每一個數(shù)據(jù)庫中都被注釋的有1 977條，占5.28%。其中，在Nr中被注釋的unigene最多 (17 428條，46.57%)，之后則依次是GO (14 131條，37.76%)、PFAM (12 217條，32.65%)、SwissProt (11 561條，30.89%)、KOG (8 606條，23%) 和KEGG (6 171條，16.49%)；在Nt數(shù)據(jù)庫中被注釋的最少，為4 144條，占11.07%。從在Nr數(shù)據(jù)庫中注釋的物種相似分布可知 (圖1)，與黑脈金斑蝶 (Danausplexippus) 的相似基因最多 (68.45%)，之后依次是家蠶 (Bombyxmori) (6.63%)、柑橘鳳蝶 (Papilioxuthus) (3.66%)、赤擬谷盜 (Triboliumcastaneum) (2.13%)和草地貪夜蛾 (Spodopterafrugiperda) (1.41%)，而與其他物種的達到17.72%。

表2 椰子織蛾轉(zhuǎn)錄組功能注釋結(jié)果統(tǒng)計Table 2 Statistics of transcriptome functional annotation for O.arenosella

圖1椰子織蛾unigene與Nr數(shù)據(jù)庫其他物種相似比對
Fig.1 Similarity comparison of unigene fromO.arenosellawith other species in Nr database

2.3 Unigene功能分類

在GO數(shù)據(jù)庫中，被注釋到的14 131條unigenes共得到了83 418次功能注釋，其被分為3個類別 (圖2)：生物學(xué)過程、細胞組分和分子功能，3個功能組又被劃分為54個第2層級分類。注釋到生物學(xué)過程的最多，為38 686條，占比46.38%；細胞組分次之，為26 506條，占比31.77%；分子功能最少，為18 226條，占比21.85%。其生物學(xué)過程中細胞過程 (Cellular process，8 378條) 和代謝過程 (Metabolic process，7 550條) 的數(shù)量最多，節(jié)律進程 (Rhythmic process，6條) 的最少；細胞組分中細胞部分 (Cell part，5 301條) 和細胞 (Cell，5 301條) 的最多，核狀小體 (Nucleoid，4條) 的最少；分子功能中結(jié)合活性 (Binding，8 314條) 和催化活性 (Catalytic activity，6 248條) 的最多，調(diào)節(jié)受體活性 (Receptor regulator activity，3條) 的最少。

圖2椰子織蛾unigene的GO分類
Fig.2 GO categories ofO.arenosellaunigene

轉(zhuǎn)錄組拼接得到的unigene在COG數(shù)據(jù)庫中比對分析，被成功注釋的8 606條unigenes共得到9 637個功能注釋，注釋到26個分類 (圖3)。其中R類 (具有一般功能) 注釋到的數(shù)量最多，有1 847條 (19.17%)；其次是T類 (信號轉(zhuǎn)導(dǎo)機制) 1 234條 (12.8%)，O類 (翻譯后修飾、蛋白轉(zhuǎn)運及蛋白伴侶) 828條 (8.59%)，S類 (未知功能) 552條 (5.73%)，J類 (參與翻譯、核糖體結(jié)構(gòu)及生物合成) 502條 (5.21%)；其余的都小于500條，最少的是X類 (未命名蛋白) 2條 (0.02%)。

進行KEGG注釋到的6 171條unigenes參與5大類通路分支 (圖4)：細胞過程 (Cellular processes)，環(huán)境信息處理 (Environmental information processing)，遺傳信息處理 (Genetic information processing)，代謝 (Metabolism) 和有機系統(tǒng) (Organismal systems)，其又歸屬于信號轉(zhuǎn)導(dǎo) (Signal transduction)、翻譯 (Translation) 和內(nèi)分泌系統(tǒng) (Endocrine system) 等32類亞通路。將KEGG數(shù)據(jù)庫作為參考，最后將unigene定位于260條通路 (數(shù)據(jù)未展示)，有16條與脂類代謝 (Lipid metabolism) 相關(guān)的通路，14條與異物生物降解與代謝 (Xenobiotics biodegradation and metabolism) 相關(guān)的通路。其中有大于200條unigenes被注釋到的代謝途徑只有核糖體 (ribosome，220條)，在100～200條unigenes之間的有18條生化代謝途徑，其余代謝途徑上注釋到的unigene數(shù)量都小于100條。

圖3椰子織蛾unigene的COG分類
Fig.3 COG categories ofO.arenosellaunigene

A. 細胞過程Cellular processes；B. 環(huán)境信息處理Environmental information processing；C. 遺傳信息處理Genetic information processing；D. 代謝Metabolism；E. 有機系統(tǒng)Organismal systems。

圖4椰子織蛾unigene的KEGG分類
Fig.4 KEGG categories ofO.arenosellaunigene

2.4 微衛(wèi)星位點分析

在37 416條unigenes中，2 849條unigenes中存在微衛(wèi)星位點，含1個以上微衛(wèi)星位點的序列數(shù)為317條。在鑒定到的3 248個微衛(wèi)星位點中 (表3)，有2 282個單堿基重復(fù)單元，重復(fù)次數(shù)最多的是10次；312個二堿基重復(fù)單元，重復(fù)次數(shù)最多的是6次；592個三堿基重復(fù)單元，重復(fù)次數(shù)最多的是5次；56個四堿基重復(fù)單元，重復(fù)次數(shù)最多的是5次；4個五堿基重復(fù)單元，重復(fù)次數(shù)最多的是5次；2個六堿基重復(fù)單元，重復(fù)次數(shù)僅為6次。

表3 椰子織蛾轉(zhuǎn)錄組的微衛(wèi)星重復(fù)單元分析統(tǒng)計Table 3 Summary of repeat units of microsatellite in O.arenosella transcriptome

椰子織蛾轉(zhuǎn)錄組中不同微衛(wèi)星重復(fù)類型分布見表4。在單堿基重復(fù)單元中，1 176個的A重復(fù)類型最多，占比為36.21%；在二堿基重復(fù)單元中，67個的AT重復(fù)類型最多，占比為2.06%；在三堿基重復(fù)單元中，50個的GGC重復(fù)類型最多，占比為1.54%；在四堿基重復(fù)單元中，各為5個的AATA和ATGG重復(fù)類型最多，占比為0.16%；在五堿基重復(fù)單元中，其重復(fù)類型僅為4類，分別是CCTCG、 CTGCA、 TACAA和TGGGT；在六堿基重復(fù)單元中，僅含有一個GATGAC和一個TGCCCA重復(fù)類型。

表4 椰子織蛾轉(zhuǎn)錄組中22個不同微衛(wèi)星重復(fù)類型統(tǒng)計Table 4 Statistics of repeat types of 22 microsatellites in O.arenosella transcriptome

3 結(jié)論與討論

轉(zhuǎn)錄組測序去除低質(zhì)量的數(shù)據(jù)后，通常以N50和Q30來評估轉(zhuǎn)錄組數(shù)據(jù)質(zhì)量。據(jù)目前轉(zhuǎn)錄組研究報道，認為N50值越大得到的長片段就越多，拼接效果越好，且當(dāng)Q30值大于80%時，說明其測序質(zhì)量可靠[14]，但也有學(xué)者認為Q30值大于85%測序質(zhì)量才可靠[15]。雖然本次轉(zhuǎn)錄組測序的N50的值比目前報道的大墊尖翅蝗 (Epacromiuscoerulipes)[16]、桑絹絲野螟 (Glyphodespyloalis)[17]等部分昆蟲轉(zhuǎn)錄組測序得到的N50小，但卻大于小菜蛾 (Plutellaxylostella)[18]、日本蚱 (Tetrixjaponica)[19]、亞洲玉米螟 (Ostriniafurnacalis)[20]等昆蟲轉(zhuǎn)錄組測序的N50。并且，本研究椰子織蛾轉(zhuǎn)錄組測序得到的Q30值大于85%，表明本次測序質(zhì)量可信度較高[18,21]。因此，本研究獲得的轉(zhuǎn)錄組數(shù)據(jù)質(zhì)量高，可用于功能注釋及新基因挖掘等后續(xù)研究。

去除低質(zhì)量序列拼接得到的unigene在Nr、Nt、Pfam、Swiss-prot、GO、KEGG和KOG數(shù)據(jù)庫中進行比對分析，預(yù)測其可能的功能分類和分子代謝通路。其中，未被成功注釋的有18 262條 (48.81%)，表明得到的unigene中存在很多功能未知的新基因[22]。另外，未被注釋unigene的存在與拼接的片段長度短、缺少基因組信息和目前仍有大量的基因功能信息缺乏有關(guān)[23-24]。在Nr數(shù)據(jù)庫中的同源比對搜索結(jié)果顯示，椰子織蛾unigene注釋到黑脈金斑蝶和家蠶上的序列最多，而注釋到赤擬谷盜和麗蠅蛹集金小蜂 (Nasoniavitripennis) 等昆蟲上的卻較少。這是由于黑脈金斑蝶和家蠶的基因組已被測定，且與椰子織蛾都屬于鱗翅目昆蟲所致。

通過GO這個國際標準化的基因功能分類體系，明確了椰子織蛾的GO功能分類，其第1層級基因注釋分類比例與鱗翅目的柞蠶 (Antheraeapernyi)[25]、菜粉蝶 (Pierisrapae)[24]等的相似。在KEGG分析中，發(fā)現(xiàn)了一些椰子織蛾unigene被注釋到脂類代謝和異物生物降解與代謝兩類代謝通路上。這兩類代謝通路中包含與昆蟲抗藥性密切相關(guān)的代謝途徑，如注釋到50條unigenes的細胞色素P450外源物質(zhì)代謝、注釋到42條unigenes的細胞色素P450藥物代謝和注釋到48條unigenes的其他酶的藥物代謝途徑，這為后續(xù)研究椰子織蛾的抗藥性機制做了鋪墊。但是，直接通過生物信息學(xué)方法得到的基因GO和KEGG注釋結(jié)果會存在假陽性，這已經(jīng)在黑腹果蠅 (Drosophilamelanogaster)、布氏錐蟲 (Trypanosomabrucei) 等中通過實驗驗證具體的基因功能得到了證實[26]。因此，該研究得到的GO和KEGG注釋結(jié)果也會存在假陽性。GO和KEGG注釋結(jié)果雖然能為挖掘研究椰子織蛾的基因功能提供了方便，但還需依據(jù)得到的功能分類或通路結(jié)果進行進一步的試驗驗證才能真正揭示相關(guān)基因的生理功能。從轉(zhuǎn)錄組中鑒定得到的3 248個微衛(wèi)星，其單堿基、二堿基和三堿基重復(fù)較為豐富，其他重復(fù)數(shù)量較少，這與黃粉甲 (Tenebriomolitor)、大墊尖翅蝗等[12,16]昆蟲的微衛(wèi)星研究結(jié)果類似。本研究獲得的椰子織蛾轉(zhuǎn)錄組數(shù)據(jù)為挖掘功能基因從分子層面揭示該害蟲的生物學(xué)和生理學(xué)相關(guān)機制奠定基礎(chǔ)，鑒定獲得的微衛(wèi)星有助于今后開發(fā)多態(tài)性微衛(wèi)星引物來研究該害蟲的種群遺傳結(jié)構(gòu)及其分子機制。