叉角厲蝽唾液腺轉(zhuǎn)錄組及差異分析*

高 平，廖賢斌，李麗芳，蘭明先，趙 航，陳 斌，吳國星，高 熹

(云南農(nóng)業(yè)大學(xué) 植物保護(hù)學(xué)院，云南 昆明 650201)

可持續(xù)農(nóng)業(yè)的發(fā)展，更加重視害蟲的綜合治理(integrated pest management)，利用天敵對作物害蟲進(jìn)行生物控制是其中重要的環(huán)節(jié)。常見的天敵有捕食性昆蟲、寄生性昆蟲和昆蟲病原體。在作物生態(tài)系統(tǒng)中，捕食性昆蟲和寄生蜂是大多數(shù)害蟲的天然生物控制因子[1]。叉角厲蝽[Eocanthecona furcellata(Wolff)]屬于半翅目(Hemiptera)蝽科(Pentatomidae) 益蝽亞科(Asopinae) 昆蟲，其若蟲和成蟲的捕食量大，捕食的范圍很廣，尤其對鱗翅目的害蟲具有較強(qiáng)的捕食能力，是一種重要的捕食性天敵[2-5]。叉角厲蝽易在實(shí)驗(yàn)室飼養(yǎng)形成種群，它的捕食潛能在東南亞、日本、印度以及中國臺灣地區(qū)等都有報(bào)道[6-7]。

近年來，高通量測序技術(shù)已應(yīng)用于昆蟲唾液的研究與鑒定。通過高通量測序技術(shù)可以獲得唾液腺的轉(zhuǎn)錄組，從基因?qū)用嫜芯客僖合俚奶卣髋c功能，并從轉(zhuǎn)錄水平預(yù)測唾液中的潛在分泌蛋白。昆蟲的唾液主要由唾液腺分泌，通過分析唾液腺的表達(dá)基因來預(yù)測唾液蛋白成分是研究昆蟲唾液的有效手段之一。目前，有報(bào)道采用轉(zhuǎn)錄組和蛋白組等手段分析煙粉虱、褐飛虱、灰飛虱和白背飛虱等唾液腺和唾液蛋白組分[8-11]。在這些昆蟲唾液腺中發(fā)現(xiàn)很多與蛋白合成、食物消化和解毒等相關(guān)的基因。然而，在捕食性昆蟲中，有學(xué)者對小花蝽(Orius laevigatus)唾液腺中差異基因表達(dá)譜和獵蝽科昆蟲有毒唾液的新型肽和唾液中的毒素進(jìn)行了研究[12-14]。

目前，國內(nèi)對該蝽的研究主要集中在食性、生物學(xué)特性、人工飼料和捕食功能研究方面[15-17]，國外對該蝽的研究主要集中在生物防治、生活史和捕食功能等方面[18]。近年來，針對與叉角厲蝽同屬于一個(gè)亞科的躅蝽(Arma chinensis)研究報(bào)道較多，王燕等[19]研究了躅蝽對草地貪夜蛾(Spodoptera frugiperda)的生防潛力，苗少明[20]對躅蝽進(jìn)行了基因組測序，并對躅蝽和茶翅蝽(Halyomorpha halys)的唾液腺進(jìn)行了轉(zhuǎn)錄組測序。益蝽亞科是蝽科中較為特化的類群，是由植食性種類分化形成的捕食性種類，雖然叉角厲蝽和躅蝽同屬于益蝽亞科，但是近幾年來國內(nèi)尚未有關(guān)于叉角厲蝽唾液腺的研究報(bào)道。

捕食性蝽類將刺吸式口器插入獵物體內(nèi)并分泌唾液，導(dǎo)致獵物快速麻痹和死亡[21]，之后再吸食獵物[22-23]。引起獵物麻痹和死亡的原因是唾液腺產(chǎn)生的分泌物釋放到獵物內(nèi)部產(chǎn)生的作用[24]。據(jù)報(bào)道組織蛋白酶B (Cathepsin B)和神經(jīng)內(nèi)分泌蛋白7b2 (Neuroendocrine protein 7b2)基因在小花蝽的捕食中起重要作用[14]。然而，叉角厲蝽唾液分泌物中引起獵物麻痹和死亡的唾液化合物仍然是未知的。因此，對唾液腺中基因表達(dá)的情況進(jìn)行深入的研究，將有助于闡明唾液毒素分泌的分子機(jī)制。本研究的目的是對叉角厲蝽的唾液腺基因進(jìn)行系統(tǒng)的研究，并確定與唾液腺毒素分泌有關(guān)的潛在基因，探究唾液腺毒素的分泌機(jī)制。期待發(fā)現(xiàn)對昆蟲取食與獵物相互作用必不可少的唾液腺相關(guān)基因，為后續(xù)的基因調(diào)控研究和生物源農(nóng)藥開發(fā)打下基礎(chǔ)。

本研究利用Illumina 測序技術(shù)對叉角厲蝽的三齡若蟲和成蟲唾液腺進(jìn)行高通量轉(zhuǎn)錄組測序，將獲得的序列信息進(jìn)行拼接和組裝，進(jìn)行基因功能注釋和功能分類，并對叉角厲蝽三齡若蟲和成蟲唾液腺的基因差異表達(dá)進(jìn)行分析，為進(jìn)一步挖掘基因、探究基因功能和分析代謝途徑等奠定基礎(chǔ)。

1 材料與方法

1.1 昆蟲采集與飼養(yǎng)

叉角厲蝽采集于云南省元江縣。帶回室內(nèi)采用籠罩式飼養(yǎng)[(27±1) ℃，RH 60%～80%，光周期14L∶10D)，每日飼喂黃粉蟲(Tenebrio molitor)，建立實(shí)驗(yàn)室種群。

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

將叉角厲蝽三齡若蟲和成蟲(雄性∶雌性=1∶1)置于?4 ℃冰箱中10 min，降低蟲體活動能力。用剪刀把叉角厲蝽所有觸角、足、口器和翅都剪除后，立即放入置于解剖鏡下的培養(yǎng)皿，在培養(yǎng)皿中加入PBS 緩沖液，以沒過蟲體為宜。使用解剖鑷和解剖針在顯微鏡下解剖獲得叉角厲蝽唾液腺，將其移至盛有Trizol 的1.5 mL 離心管中，研磨并置于?80 ℃的冰箱中保存。三齡若蟲和成蟲各設(shè)3 次重復(fù)，每個(gè)重復(fù)解剖80 個(gè)唾液腺。采用Trizol Reagent 的方法提取得到上述樣品的總RNA，RNA 提取后分別用瓊脂糖凝膠電泳和NanoPhotometer spectrophotometer 檢測RNA的純度，檢測合格之后寄送至北京諾禾致源生物信息科技有限公司進(jìn)行轉(zhuǎn)錄組文庫構(gòu)建和測序。

1.3 cDNA 文庫構(gòu)建與測序

取1.5 μg 質(zhì)檢合格的RNA 樣品進(jìn)行文庫構(gòu)建，然后使用Agilent 2100 bioanalyzer 對文庫的insert size 進(jìn)行檢測。庫檢合格后，把不同文庫按照有效濃度及目標(biāo)下機(jī)數(shù)據(jù)量的需求pooling 后進(jìn)行Illumina 測序。

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

測序得到的原始數(shù)據(jù)要去除序列接頭、ploy-N 和低質(zhì)量序列，才能獲得高質(zhì)量的Clean Data，同時(shí)計(jì)算Q20、Q30 和GC 含量。獲得高質(zhì)量的測序數(shù)據(jù)再利用Trinity 軟件進(jìn)行組裝[25]。將組裝得到的unigenes 使用BLAST 軟件在七大數(shù)據(jù)庫中進(jìn)行比對，獲得功能基因的注釋信息。

1.5 關(guān)鍵基因鑒定

為了鑒定叉角厲蝽唾液腺中引起獵物死亡的關(guān)鍵基因，在叉角厲蝽唾液腺NR 注釋數(shù)據(jù)庫里使用基因全稱搜索BAEK 等[14]報(bào)道的組織蛋白酶B 和神經(jīng)內(nèi)分泌蛋白7b2 基因，并在轉(zhuǎn)錄組組裝(transcriptome assembly)數(shù)據(jù)庫里查找編碼這些基因的序列。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄組序列分析和組裝

利用Illumina 測序平臺對叉角厲蝽唾液腺轉(zhuǎn)錄組進(jìn)行測序，每個(gè)樣本都生成超過7.80 GB 的高質(zhì)量數(shù)據(jù)，樣本的GC 含量不低于43.71%，測得數(shù)據(jù)Q20 均大于97.00%，堿基Q30 百分比不少于92.00%，表明測序結(jié)果高度準(zhǔn)確，可以用于后期的拼接組裝(表1)。

獲得的clean reads 進(jìn)行組裝(表2)，共得到51 845 個(gè)轉(zhuǎn)錄本(transcripts)，序列信息達(dá)到81 004 577 bp，平均長度為1 562 bp，N50 長度為2 615 bp。在轉(zhuǎn)錄本基礎(chǔ)上進(jìn)一步組裝獲得26 022條功能基因(unigenes)，長度為35 520 646 bp，N50 長度為2 298 bp；其中，長度超過1 kb 的unigenes 有10 634 條，占40.87%。

2.2 基因功能注釋

unigenes 在七大數(shù)據(jù)庫中比對得到功能信息(表3)。在組裝獲得的26 022 條 unigenes 中成功注釋14 057 條，注釋率為54.01%。其中，NR 注釋到的unigenes 最多(12 514，48.09%)，之后依次為Pfam (9 025，34.68%)、GO (9 025，34.68%)、NT (8 643，33.21%)、SwissProt (8 573，32.94%)和KOG (5 423，20.84%)；在KO 數(shù)據(jù)庫中注釋的最少，為5 084 條，所占比例為19.53%。

表1 叉角厲蝽轉(zhuǎn)錄組測序產(chǎn)量統(tǒng)計(jì)Tab.1 Output statistics of E.furcellata transcriptome sequencing

表2 叉角厲蝽轉(zhuǎn)錄組組裝質(zhì)量Tab.2 Quality for de novo E.furcellata transcriptome assembly

E值等于0 說明該基因的注釋可信。結(jié)果顯示：所有匹配結(jié)果的E值均大于1e-5，E值大于1e-5 的占86.3% (圖1a)。從圖1b 可見：Nr 注釋到的unigenes 基因中，相似性在95%以上的基因占15.5%。從匹配的物種來源分析(圖1c)，叉角厲蝽唾液腺轉(zhuǎn)錄組測序拼接的unigenes 與茶翅蝽(Halyomorpha halys)的unigenes 相似性最高，達(dá)到76.0%。其次與溫帶臭蟲(Cimex lectularius)、Cryptotermes secundus、褐飛虱(Nilaparvata lugens)和Lasius niger等的相似性分別為3.7%、2.7%、1.6%和1.0%；另外與其他物種同源的基因占15.0%。

表3 基因注釋成功率統(tǒng)計(jì)表Tab.3 Success rate statistics of unigene annotation

圖1 NR 注釋的E 值分布和物種分布Fig.1 E-value distribution and species classification of sequences matched to the NR database

2.3 基因功能分類

叉角厲蝽唾液腺轉(zhuǎn)錄組有23 132 條(47.78%)unigenes 被GO 注釋至生物學(xué)過程(biological process)的26 個(gè)功能亞類；注釋至分子功能(molecular function)的unigenes 有10 660 條(22.02%)；14 621 條(30.20%) unigenes 歸屬于細(xì)胞組分(cellular component)的10 個(gè)功能亞類(圖2)。在生物學(xué)過程中，注釋至細(xì)胞過程(cellular process，GO：0009987)的unigenes 數(shù)量最多，達(dá)到4 859 條(21.00%)；注釋到代謝過程(metabolic process，GO：0044237)的unigenes 數(shù)量較多，達(dá)到4 496條(19.44%)；注釋到單生物過程(single-organism process，GO ：0044699) 的unigenes 達(dá)到3 791 條(16.39%)。其余均在3 000 條以下，其中注釋至行為(behavior，GO：0007610)、節(jié)律性過程(rhythmic process，GO：0048511) 和解毒過程(detoxification，GO：0098754) 的unigenes 僅有5、5 和6 條。分子功能類型中，注釋至蛋白結(jié)合活性(binding，GO：0005488) 和催化活性(catalytic activity，GO：0003824) 的unigenes 分別有5 096 條(47.80%) 和3 790 條(35.55%)；注釋至轉(zhuǎn)運(yùn)活性(transporter activity，GO：0005215)的unigenes 數(shù)量達(dá)到677 條。其余功能亞類均在300 條以下，其中被注釋至金屬伴侶蛋白活性(metallochaperone activity，GO：0016530) 的unigenes 數(shù)量僅有4 條序列。細(xì)胞組分功能類型中，注釋至細(xì)胞(cell，GO：0005623) 和細(xì)胞部分(cell part，GO：0044464) 的unigenes 都達(dá)到2 793 條(19.10%)；注釋至細(xì)胞器(organelle，GO：0043226)、大分子復(fù)合物(macromolecular complex，GO：0032991)、細(xì)胞膜(membrane，GO：0016020)、細(xì)胞膜部分(membrane part，GO：004 4425)和細(xì)胞器部分(organelle part，GO：00444 22)的unigenes 數(shù)量也分別達(dá)到1 984 條(13.57%)、1 702 條(11.64%)、1 630 條(11.15%)、1 531 條(10.47%)和1 045 條(7.15%)；其余被注釋的功能亞類均在1 000 條以下。

所有unigenes 經(jīng)過KOG 數(shù)據(jù)庫功能預(yù)測和分類，叉角厲蝽唾液腺轉(zhuǎn)錄組共有5 423 個(gè)unigenes 得到注釋，并被分為26 個(gè)KOG 類(圖3)。其中，共有4 個(gè)KOG 亞類包含450 條以上unigenes，依次是：一般功能預(yù)測(general function prediction only，R) 913 條unigenes，占比最大(16.84%)；信號轉(zhuǎn)導(dǎo)機(jī)制(signal transduction mechanisms，T) 675 條(12.45%) unigenes；蛋白質(zhì)翻譯后修飾、折疊和分子伴侶(posttranslational mofdiication，protein turnover，chaperones) 570 條(10.51%) unigenes；未知功能(function unknown，S) 496 條(9.15%) unigenes。而細(xì)胞動力(cell motility，N)、核結(jié)構(gòu)(nuclear structure，Y) 和未命名蛋白(unnamed protein，X)的KOG 亞類均包含不足30 條unigenes。

圖2 叉角厲蝽unigenes 的GO 分類Fig.2 Gene ontology (GO) classification of assembled E.furcellata unigenes

圖3 叉角厲蝽unigenes 的KOG 分類Fig.3 Karyotic orthologous groups (KOG) classification of E.furcellata unigenes

共有5 084 條unigenes 比對到KEGG 數(shù)據(jù)庫。在獲得注釋的unigenes 中，25.41%涉及有機(jī)系統(tǒng)(E，organismal systems)，其中主要參與內(nèi)分泌系統(tǒng)(endocrine system，325 條)、消化系統(tǒng)(digestive system，209 條) 和免疫系統(tǒng)(immune system，205 條)，其余亞類均在200 條以下。24.63%涉及新陳代謝(D，metabolism)，主要參與糖類(carbohydrate metabolism，204 條)、脂質(zhì)(lipid metabolism，182 條) 和蛋白質(zhì)(amino acid metabolism，157 條)三大營養(yǎng)物質(zhì)的代謝以及核苷酸代謝(nucleotide metabolism，117 條)等。除此以外，這些unigenes 還廣泛參與到信號轉(zhuǎn)導(dǎo)(signal transduction，548 條)通路中(圖4)。

2.4 差異表達(dá)基因的篩選

通過FPKM 方法評估本次測序叉角厲蝽唾液腺樣品中的基因表達(dá)情況，采用DESeq 進(jìn)行若蟲和成蟲唾液腺樣品組間的差異表達(dá)分析。在差異表達(dá)基因的篩選過程中，將Padj＜0.05 作為差異基因的篩選標(biāo)準(zhǔn)，比較2 組測序結(jié)果，篩選出大量差異表達(dá)的基因。共鑒定出3 437 個(gè)差異表達(dá)基因(圖5)。

2.5 差異表達(dá)基因GO 富集分析

GO 富集柱狀圖能表現(xiàn)出差異表達(dá)基因在GO term 上富集的分布情況(圖6)。將篩選到的3 437個(gè)差異表達(dá)基因(DEGs)分類到3 個(gè)ontology (本體)中，差異表達(dá)基因在分子功能和生物學(xué)過程中分別包含了15 和12 個(gè)功能小類。在生物過程這一類中，差異表達(dá)基因在蛋白水解(proteolysis)和DNA 復(fù)制(DNA replication)2 個(gè)功能小類中所占比例最高；在分子功能這一類中，差異表達(dá)基因在水解酶活性(hydrodase activity)和肽酶活性(peptidase activity) 2 個(gè)功能小類中占比最高。

2.6 差異表達(dá)基因的KEGG 注釋

KEGG 代謝通路分析共得到1 509 個(gè)差異表達(dá)基因，這些差異表達(dá)基因參與261 個(gè)不同的代謝通路。表4 是含有最多差異基因的前20 個(gè)通路，在所有通路中代謝途徑是含有差異基因序列最多的。KEGG 注釋信息主要富集在嘌呤代謝(purine metabolism)、嘧啶代謝(pyrimidine metabolism)、半胱氨酸和蛋氨酸的代謝(cysteine and methionine metabolism)以及組氨酸代謝(histidine metabolism)等。在KEGG 代謝通路中對與分泌相關(guān)的差異表達(dá)基因數(shù)量進(jìn)行統(tǒng)計(jì)(表5)，共得到58 個(gè)差異表達(dá)基因。

2.7 關(guān)鍵基因的鑒定

在叉角厲蝽唾液腺轉(zhuǎn)錄組數(shù)據(jù)庫中鑒定了3 個(gè)編碼組織蛋白酶B 的基因。組織蛋白酶B 是溶酶體半胱氨酸蛋白酶，由二硫鍵連接的重鏈和輕鏈的二聚體組成。同時(shí)，還鑒定了3 個(gè)編碼神經(jīng)內(nèi)分泌蛋白7b2 的基因。神經(jīng)內(nèi)分泌蛋白7b2 在大多數(shù)神經(jīng)元和內(nèi)分泌細(xì)胞中起著至關(guān)重要的作用。

圖4 叉角厲蝽轉(zhuǎn)錄組unigenes 的KEGG 通路分析Fig.4 Classification of E.furcellata transcriptome based on KEGG

圖5 差異基因火山圖Fig.5 Differential gene volcano map

3 討論

轉(zhuǎn)錄組測序以N50 和Q30 來評估轉(zhuǎn)錄組數(shù)據(jù)的質(zhì)量。一般認(rèn)為N50 值越大就表明得到的長片段越多[26]，且Q30 值大于85% 說明測序結(jié)果可靠[27]。通過對叉角厲蝽成蟲和三齡若蟲唾液腺進(jìn)行轉(zhuǎn)錄組測序，共獲得唾液腺unigenes 序列26 022 條，平均長度1 365 bp，N50 長度為2 298 bp，Q30 值都大于92%。而茶翅蝽和躅蝽的唾液腺轉(zhuǎn)錄組測序數(shù)據(jù)為53 618 條和92 468條，平均長度為736 和729 bp，N50 長度僅為1 439 和1 311 bp[13]。從N50 長度和Q30 值等方面均顯示本次叉角厲蝽唾液腺轉(zhuǎn)錄組的拼接結(jié)果良好。

將全部unigenes 序列與七大數(shù)據(jù)庫進(jìn)行功能注釋，獲得14 057 條unigenes (54.02%)成功注釋，仍有11 965 條unigenes (45.98%)未被注釋，這可能是unigenes 較短，所以未與公共數(shù)據(jù)庫中的序列比對上[28]，也有可能是存在很多功能未知的新基因[29]。通過與NR 數(shù)據(jù)庫進(jìn)行比對，結(jié)果表明叉角厲蝽與茶翅蝽的同源序列最多，茶翅蝽雖是植食性害蟲，主要危害蔬菜和水果[30]，但這可能與茶翅蝽和叉角厲蝽同屬于蝽科，親緣關(guān)系較近有關(guān)。然而本次同源序列比對并沒有注釋到與躅蝽的相似性，這可能是由于NR 數(shù)據(jù)庫中暫時(shí)還沒有收錄到躅蝽的注釋信息。

圖6 GO 富集柱狀圖Fig.6 GO enrichment histogram

表4 差異基因KEGG 富集列表(前20 個(gè))Tab.4 Differential gene KEGG enrichment list (top 20)

表5 與分泌相關(guān)差異表達(dá)基因的統(tǒng)計(jì)Tab.5 Statistics of differentially expressed genes in secretion

GO 功能富集結(jié)果顯示：叉角厲蝽唾液腺參與生物進(jìn)程的unigenes 最多，其中前3 個(gè)類別分別是細(xì)胞進(jìn)程、代謝進(jìn)程和單生物進(jìn)程，由此可以看出唾液腺中的細(xì)胞新陳代謝十分活躍。在分子功能中，unigenes 數(shù)量最多的3 個(gè)類別是結(jié)合、催化活性和轉(zhuǎn)運(yùn)活性。這與白背飛虱(Sogatella furcifera)唾液腺[31]和煙粉虱(Bemisia tabaci)唾液腺[10]的研究結(jié)果較為一致。分泌唾液是昆蟲唾液腺的主要功能，注釋結(jié)果顯示有很多基因參與結(jié)合和轉(zhuǎn)運(yùn)。在細(xì)胞組分中，前3 個(gè)是細(xì)胞、細(xì)胞部分和細(xì)胞器，由此可見細(xì)胞、細(xì)胞器和細(xì)胞分子是唾液腺的重要組成部分[32]。叉角厲蝽唾液腺的這些注釋結(jié)果表明：在唾液的合成和分泌過程中需要活躍的代謝活動和高效的蛋白質(zhì)合成與轉(zhuǎn)運(yùn)機(jī)制。

通過對叉角厲蝽成蟲和三齡若蟲唾液腺的差異表達(dá)分析，發(fā)現(xiàn)叉角厲蝽唾液腺的26 022 個(gè)unigenes 中有3 437 個(gè)unigenes 在成蟲和三齡若蟲中的表達(dá)水平不同。數(shù)據(jù)統(tǒng)計(jì)發(fā)現(xiàn)有1 871 個(gè)表達(dá)水平上調(diào)的unigenes，1 566 個(gè)表達(dá)水平下調(diào)的unigenes。GO 富集分析中，差異基因主要富集在水解酶活性、肽酶活性和蛋白水解，表明三齡若蟲唾液腺還處在快速發(fā)育階段。KEGG 代謝通路分析有1 509 個(gè)差異表達(dá)基因參與261 個(gè)不同的代謝或信號通路，同時(shí)發(fā)現(xiàn)13 個(gè)差異表達(dá)基因參與了唾液分泌，這些差異基因的功能有待進(jìn)一步研究。

BAEK 等[14]報(bào)道組織蛋白酶B 可能參與了特定的蛋白水解，從而影響獵物的凝血或免疫力，最終導(dǎo)致殺蟲作用，同時(shí)KUTSUKAKA 等[33]發(fā)現(xiàn)在Tuberaphis styraci中組織蛋白酶B 是毒素的主要成分。神經(jīng)內(nèi)分泌蛋白在激素原轉(zhuǎn)化酶2(prohormone convertase-2) 的運(yùn)輸和激活中起關(guān)鍵作用，激素原轉(zhuǎn)化酶2 負(fù)責(zé)處理神經(jīng)肽和肽激素前體[34]。另外，神經(jīng)內(nèi)分泌蛋白7b2 可能調(diào)節(jié)分泌顆粒的形成和分泌[14]，并可能通過控制激素原轉(zhuǎn)化酶2 的激活和轉(zhuǎn)運(yùn)參與唾液肽的成熟。組織蛋白酶B 和神經(jīng)內(nèi)分泌蛋白7b2 基因可能在唾液的功能中發(fā)揮重要的作用。

4 結(jié)論

唾液腺是昆蟲攝食的重要器官[35]，研究叉角厲蝽唾液腺轉(zhuǎn)錄組意義重大。這些轉(zhuǎn)錄組數(shù)據(jù)將為叉角厲蝽生命活動中功能基因的預(yù)測、分子標(biāo)記的開發(fā)、昆蟲與其他生物相互作用和信號通路的研究提供有用的分子資源，也可為其他半翅目昆蟲提供數(shù)據(jù)參考。