基于轉錄組分析椰心葉甲嚙小蜂復壯和保種種群差異表達基因

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

(1.中國熱帶農業科學院椰子研究所，海南省院士團隊創新中心，海南省熱帶油料作物生物學重點實驗室，海南文昌 571399；2.海南大學植物保護學院，海南海口 570228；3.海南大學熱帶作物學院，海南海口 570228)

重大危險性外來有害生物椰心葉甲Brontispa longissima原產于印度尼西亞和巴布亞新幾內亞，主要危害植物心葉組織，造成心葉大量死亡，嚴重影響棕櫚植物的景觀度和植株的健康生長[1]。2002年首次在海南省海口市椰子樹上發現椰心葉甲[2]，隨后短短幾年便傳播至海南全省[3]。近年來，椰心葉甲在海南省各地檳榔園頻繁發生，擴散蔓延快[4－5]，對椰心葉甲防治研究仍需加強。椰心葉甲嚙小峰Tetrastichus brontispae是椰心葉甲的蛹期寄生蜂[6]，國內外學者已對椰心葉甲嚙小蜂的成蟲習性、寄主適應性、生物生態學特性，以及室內大量繁殖和野外釋放等方面進行了相關研究[6－10]。飼養并釋放椰心葉甲嚙小蜂可作為防治椰心葉甲的有效手段。

自然界中的個體在種群進化的過程中，為了適應生存環境的變化會產生一系列的變異，當這種變異朝著各項生理指標下降的方向發展，就稱之為退化[11]。室內飼養的昆蟲退化往往表現在繁殖能力下降、種群整齊度下降、壽命縮短和抗病能力下降等方面。昆蟲退化現象是普遍存在的，目前關于黑腹果蠅Drosophila melanogaster[12]、煙蚜繭蜂Aphidius gifuensis[13]和管氏腫腿蜂Scleroderma guani[14]種群的退化和復壯已經進行了相關研究，都表明長時間的近親繁殖會導致昆蟲種群的退化，通過回交和外源種群雜交的方法可以有效復壯。為探究椰心葉甲嚙小蜂種群退化的分子機制，筆者基于椰心葉甲嚙小蜂2個種群的轉錄組數據，對差異基因進行分析。

1 材料與方法

1.1 供試蜂種 選用中國熱帶農業科學院椰子研究所天敵工廠種群(每年復壯1次)和中國熱帶農業科學院環境與植物保護研究所天敵工廠種群(自2004年引進后隔離繁殖，未經復壯)，各設置3次重復，每重復20頭，均為1～2日齡的椰心葉甲嚙小蜂(不分雌雄)。

1.2 轉錄組數據 用于對比分析的轉錄組數據已提交至NCBI SRA數據庫中，Genbank登錄號:PRJNA678031。原始數據質量分析詳見劉華偉等[15]。測序的質量較好，可用于后續差異基因的分析。

1.3 差異基因表達分析及功能富集分析 皮爾遜相關系數(pearson′s correlation coefficient，r2)作為樣品間相關性的評估指標[16]，r2越接近1，說明兩個樣品之間的相關性越強。通過聚類分析判斷差異表達基因在不同實驗條件下的表達模式，在樣品間表達量相關性高的基因被歸為一類，使用R語言Pheatmap軟件包對所有比較組的差異基因的并集和樣品進行雙向聚類分析，根據同一基因在不同樣品中的表達水平和同一樣品中不同基因的表達模式進行聚類，采用Euclidean方法計算距離，采用層次聚類最長距離法(Complete Linkage)進行聚類。通過DESeq(http://www.bioconductor.org/packages/release/bioc/html/DESeq.html)軟件對復壯種群(WC)和保種種群(DZ)進行差異表達分析，篩選差異表達基因條件為:表達差異倍數log2FoldChange＞1，顯著性 P-value＜0.05。 使用topGO進行 GO富集分析，顯著富集的標準為P-value＜0.05，確定差異基因行使的主要生物學功能。將差異基因與KEGG數據庫比對，顯著富集的標準為P-value＜0.05，獲得差異基因顯著性富集的通路注釋信息。

1.4 RT-qPCR驗證 選擇與椰心葉甲嚙小蜂各種行為相關的 7個基因(編號為 DN935，DN4191，DN4534，DN26397，DN14402，DN535，DN1980)進行熒光定量PCR，以驗證轉錄組數據的可靠性。

2 結果與分析

2.1 椰心葉甲嚙小蜂差異基因表達量分析 通過對復壯種群3個重復(WCA、WCB和WCC)與保種種群3個重復(DZA、DZB和DZC)間的差異基因表達量進行分析，各樣本間的相關系數r2均大于0.95，表示種群內各樣本間及復壯與保種種群之間的相關性較強(表1)。通過聚類分析發現，在復壯種群和保種種群的6個樣本中，高表達基因和低表達基因的分布在不同種群之間出現分離，在同一種群的不同生物學重復中聚集在一起。表明所用樣本在組內的基因表達模式相似，生物學重復性較好，且樣本分組較合理(圖1)。

表1 樣品基因表達量相關性Tab.1 Correlation of sequenced samples based on gene expression level

圖1 差異表達基因表達模式聚類圖Fig.1 Expression pattern clustering of DEGs

相較于復壯種群，保種種群(DZ/WC)共有553個差異表達基因，其中上調基因206個、下調基因347個。兩組樣本間基因差異表達的分布情況見圖2。其中，有3個基因的注釋結果為蜂王漿主蛋白(major royal jelly proteins，MRJPs)，3 個基因的注釋結果為細胞色素P450(cytochrome P450，CYP450)，它們在保種種群中下調表達；此外，還發現幾種嗅覺相關蛋白。

圖2 基因差異表達火山圖分析Fig.2 Volcano plot of DEGs

2.2 差異表達基因的GO富集分析和KEGG通路分析 差異表達基因共有14個基因顯著富集到分子功能(molecular function，MF)中，并且全部下調，富集到單氧酶活動(monooxygenase activity)的差異基因數量最多，占供試基因總數的28.57%。有25個基因顯著富集到細胞組分(cellular component，CC)中，其中10個基因上調，15個基因下調，富集到細胞膜(cell membrane)中的差異基因數量最多，占供試基因總數的68%。有25個基因顯著富集到生物過程(biological process，BP)中，其中13個基因上調，12個基因下調，此分類中各項顯著富集到的差異基因都較少，均為1～2個。在顯著富集的差異基因中，包括幾種氣味受體和細胞色素 P450(圖3)。

圖3 差異表達基因GO分類Fig.3 GO classification of DEGs

差異表達基因共有4條顯著富集的代謝通路，分別為Hh信號通路(Hedgehog signaling pathway)、Hh信號通路－fly(Hedgehog signaling pathwayfly)、胰腺分泌(pancreatic secretion)和神經活性配體受體相互作用(neuroactive ligand-receptor interaction)。因4條通路涉及的基因有重合，所以共計10個差異基因(表2)。

表2 差異表達基因的顯著富集代謝通路Tab.2 KEGG analysis of DEGs

2.3 RT-qPCR驗證轉錄組數據 在椰心葉甲嚙小蜂保種和復壯2個種群中，7個基因的表達變化趨勢與轉錄組的測序結果是一致的，可以認為轉錄組的數據是可靠準確的(圖4)。

圖4 各基因轉錄組及RT-qPCR統計圖Fig.4 Transcriptome and RT-qPCR of each gene

3 結論與討論

有研究證明胰蛋白酶(TRINITY＿DN334＿c0＿g1，TRINITY＿DN7050＿c0＿g1)與自身免疫防御有關[17]，還發現離子型谷氨酸受體(TRINITY＿DN29200＿c0＿g1)屬于嗅覺相關蛋白中的離子型受體(IRs)，參與化學感受基因感覺[18]。在差異表達基因的KEGG通路富集分析中發現有2條通路中的基因全部下調，表示保種種群的椰心葉甲嚙小蜂的自身免疫防御和嗅覺系統可能出現退化。

蜂王漿主蛋白是蜂王漿中的水溶性蛋白，與蜂王漿的生物學功能和醫療保健作用可能有密切的關系，比如蜂王漿主蛋白參與雌蜂的級型分化，為幼蟲發育提供可利用的氮元素并可能在蜜蜂行為調控中起作用。另外還發現蜂王漿主蛋白在其他動物中具有重要的醫療保健作用，主要表現為高效的抗菌活性、較強的免疫活性和促進特殊細胞增殖等[19]。在椰心葉甲嚙小蜂轉錄組數據中，共檢測出25個蜂王漿主蛋白基因，其中3個基因差異表達且全部下調，可以說明保種種群椰心葉甲嚙小蜂的抗菌活性、免疫活性等出現退化。

細胞色素 P450(cytochrome P450，CYP450)是一個超家族酶系，幾乎存在于昆蟲的所有組織中，參與外源化合物的代謝、蛻皮激素和保幼激素的合成與降解等過程，還可以調節昆蟲的生長和發育[20]。在椰心葉甲嚙小蜂轉錄組中發現124個P450基因，其中有3個基因差異表達，全部下調。在對飛蝗CYP450的研究中發現，CYP303A1可能與蛻皮激素代謝調控和影響新表皮增厚有關[21]，進而搜索椰心葉甲嚙小蜂轉錄組數據時發現一個幼蟲表皮蛋白也同樣下調，表明保種種群的椰心葉甲嚙小蜂的幼蟲發育過程可能受阻。但在果蠅的研究中表示，CYP303A1與果蠅剛毛結構和機械感覺有關[22]，在茶小綠葉蟬 Empoasca flavescens[23]中還發現CYP303A1參與藥物代謝解毒過程，出現這種情況的原因可能是物種的親緣關系較遠，同一個基因在不同物種間可能存在差異。

嗅覺系統是幫助昆蟲識別各種化學氣味分子，引導昆蟲進行覓食、交配、選擇產卵地點和躲避天敵等生命活動的重要系統[24－25]。研究表明，有多種嗅覺相關蛋白參與了嗅覺識別過程，主要有氣味結合蛋白(odorant binding proteins，OBPs)、氣味降解酶(odor degrading enzymes，ODEs)、嗅覺受體(olfactory receptors，OLRs)、氣味受體(odorant receptors，ORs)、化學感受蛋白(chemosensory proteins，CSPs)、感覺神經元膜蛋白(sensory neuron membrane proteins，SNMPs)及離子型受體(ionotropic receptors，IRs)等[26]。在椰心葉甲嚙小蜂轉錄組數據中，共發現153個 ORs基因，394個OBPs基因，18個SNMPs基因和16個IRs基因，并未發現ODEs、OLRs、CSPs相關基因。在這些基因中，共有7個基因在兩個種群中差異表達，包括4個ORs基因上調，2個OBPs基因和1個IRs基因下調。一般認為OBPs和ORs都是在嗅覺識別進程中起作用，氣味結合蛋白結合和轉運氣味分子，氣味受體完成信號的激活和轉導，但在本研究中ORs(上調)和OBPs(下調)的表達模式相反。在對果蠅的研究中發現，OR82a基因會抑制果蠅的求偶行為[27]，但其他ORs的功能仍不明確，需要繼續研究。

通過對椰心葉甲嚙小蜂2個種群的轉錄組數據進行差異基因分析，發現保種種群在一些方面已經出現了明顯的退化，這也表明在室內大量繁殖過程中，合理的繁殖技術和適時的復壯是非常有必要的，后續應針對退化規律、復壯時間等展開研究，以期為更新現有的繁蜂技術提供數據支撐。