殺蟲劑處理的梨小食心蟲轉錄組測序分析及差異表達解毒酶基因篩選

劉曉慶 ，龐欽瑋 ，冀佳悅 ，楊 春 ，高玲玲 ，郭艷瓊

（1.山西農業大學 植物保護學院，山西 太谷 030801；2.山西農業大學 果樹研究所，山西 太谷 030815；3.澳大利亞聯邦科學與工業研究組織農業與食品部，西澳 文布利 6913）

梨小食心蟲（Grapholita molesta）屬鱗翅目卷蛾科，簡稱“梨小”，廣泛分布于亞洲、歐洲和美洲，我國從南至北均有發生[1-2]，主要危害蘋果、梨、核桃和紅棗等仁果類和核果類果樹。目前，梨小食心蟲的防治手段仍以化學防治為主。梨小食心蟲幼蟲具有鉆蛀性，難以完全接觸藥劑，因此，果園噴施農藥通常是用于滅殺成蟲，阿維菌素、高效氯氟氰菊酯是防治梨小食心蟲的有效藥劑，吡蟲啉作為果園防治害蟲常用殺蟲劑，對梨小食心蟲也有一定的防治效果。

杜恩強等[3]試驗表明，阿維菌素能殺滅梨小成蟲、抑制成蟲產卵；庾琴等[4]研究表明，阿維菌素、高效氯氟氰菊酯等殺蟲劑對梨小初孵幼蟲有一定的防治效果；陳靜等[5]研究表明，高效氯氟氰菊酯、阿維菌素對梨小具有一定的防治效果；田間藥效試驗發現，4.5%高效氯氰菊酯EC 對梨小15 d 田間防效達98%[6]。然而，由于梨小的鉆蛀特性，使部分齡期難以完全接觸藥劑，且化學農藥的長期不合理使用，使梨小對不少農藥耐受性增加。害蟲對殺蟲劑產生耐藥性的過程，與解毒酶的活性密切相關。解毒酶代謝活性增加在害蟲產生耐藥性過程中起重要作用[7-9]。

細胞色素P450（Cytochrome P450，CYP）、羧酸酯酶（Carboxylesterases，CarEs）和谷胱甘肽硫轉移酶（GlutathioneS-transferases，GSTs）是昆蟲三大解毒酶系。細胞色素P450 可廣泛地代謝外源和內源化合物[10]，在解毒和激活方面具有雙重作用。近年來，對雙翅目、鱗翅目和鞘翅目等多種昆蟲進行了全基因組測序，與昆蟲解毒代謝相關的P450基因大部分集中在CYP4、CYP6 和CYP9 家族中，通過聚類分析發現，這些細胞色素基因分屬于CYP2 簇、CYP3 簇、CYP4 簇和線粒體簇中[11]。CarEs 主要是通過將包含酯鍵的內源性和外源性化合物分解為醇和酸[12]，達到解毒目的，其共有14 個進化支（A-N），分為3 個功能組，即與代謝解毒相關的酶類（A-C）、與激素和信息素降解相關的酶類（D-G）和與神經發育相關的酶類（H-N）。GSTs 通過催化谷胱甘肽與有毒的疏水親電物質發生軛合反應，使毒性降低，可溶性增強，進而排出體外，達到解毒的目的[13]，GSTs 在昆蟲中可分為七大類：sigma、zeta、omega、theta、microsomal、delta和epsilon[14]。這3 種解毒酶廣泛分布于動植物和微生物中[15]，可將進入昆蟲體內的有毒物質轉變為弱毒或無毒的形式[16]。

隨著轉錄組測序研究的深入，高通量雙端RNA測序被廣泛應用，它可以測得不同組織或不同條件下的轉錄組數據[17]，例如，人、鼠、昆蟲、微生物等的轉錄組[18-20]，此外，RNA 測序還能進行差異表達基因的分析[21]。CRAVA 等[22]對櫻桃果蠅Drosophila suzukii的轉錄組分析，篩選出了OR 基因；TIAN等[23]對桃蛀果蛾Carposina sasakii轉錄組進行分析，篩選出了化學感受蛋白基因；CHEN 等[24]對韭菜遲眼蕈蚊Bradysia odoriphaga轉錄組分析，篩選出4 個與吡蟲啉代謝有關的P450 基因。

本研究對3 種殺蟲劑不同濃度處理的梨小轉錄組文庫進行高通量測序，篩選可能與殺蟲劑代謝相關的解毒酶基因，旨在為合理科學使用殺蟲劑及深入開展研究梨小解毒酶基因的功能提供理論依據。

1 材料和方法

1.1 試驗材料

梨小食心蟲成蟲初始種群于2012 年采自山西省太谷縣西山底村桃林（37°21′42″N，112°34′19″E），將其帶回實驗室放置于溫度（26±1）℃、相對濕度70%±10%、光周期15L∶9D 的人工氣候培養箱中飼養。

RNAiso Plus，購買自TaKaRa；cDNA 第1 鏈合成試劑盒PrimeScript RT reagent kit with gDNA Eraser，購買自北京寶日醫生物技術有限公司；殺蟲劑為阿維菌素（avermectin）、吡蟲啉（imidacloprid）、高效氯氟氰菊酯（lambda-cyhalothrin），均購買自SIGMA-ALDRICH 公司。PRX-450C 人工氣候箱（寧波海曙賽福儀器廠）；HS-840 超凈工作臺（蘇州安泰空氣技術有限公司）。

1.2 試驗方法

1.2.1 3 種殺蟲劑室內毒力測定 用藥膜法[25]測定吡蟲啉、阿維菌素和高效氯氟氰菊酯對梨小食心蟲的室內毒力。以丙酮為溶劑將原藥配制成母液，并用丙酮稀釋成一系列質量分數的供試藥液。分別吸取1 mL 藥液加入到指形管（直徑16 mm，長度80 mm）中；立即滾動小瓶，使農藥均勻地分布在內表面，在丙酮揮發后用吸水棉蓋住。將7 只健康的梨小成蟲接入帶有殺蟲劑薄膜的指形管中，每組設置3 個生物學重復，原飼養條件下飼養24 h 觀察死亡情況并記錄死亡率，輕觸蟲體不動為死亡。

1.2.2 樣品制備 根據室內毒力測定結果，得到吡蟲啉、阿維菌素和高效氯氟氰菊酯亞致死濃度（LC10、LC30）、致死中濃度（LC50），分別吸取3 種殺蟲劑3 種濃度1 mL 藥液于指形管中，立即滾動小瓶，使農藥均勻地分布在內表面，并將7 只健康成蟲放入帶有殺蟲劑薄膜的小瓶24 h，然后收集存活的成蟲作為測序的一個樣本，每組設置3 個生物學重復。所有的樣品立即在液氮中冷凍，然后儲存在-80 ℃用于提取RNA。

1.2.3 RNA 提取 按照說明用Trizol RNA 提取試劑和Rnase 提取1.2.2 樣品中的總RNA，然后用Rnasy 離心柱純化。利用超微量分光光度計檢測RNA 的質量，用瓊脂糖凝膠電泳檢測RNA 的完整性。

1.2.4 轉錄組測序及數據組裝 采用Illumine HiSeq 2000 測序儀進行測序，對獲得的原始測序數據進行過濾，去除接頭和含未知堿基N 比例超過5%的reads 以及低質量測序序列，得到待分析數據（Clean reads）。使用Trinity 軟件對獲得的高質量測序數據進行組裝，然后使用Tgicl 軟件進行聚類去冗余，獲得Unigene。

1.2.5 功能注釋 為獲得全面的基因功能信息，將非重復序列基因通過BLAST 分別與NR、NT、SWISS-PROT、Pfam、GO、KOG、KEGG 數據庫進行比對，獲得相應的注釋信息。

1.2.6 差異表達基因分析 使用DESeq 軟件對差異表達基因進行分析，篩選差異表達基因的標準為誤報率（falsediscovery rate，FDR）小于0.05 且差異倍數大于2。通過MEGA 7.0 軟件使用P 距離模型1 000 次引導值的鄰接法進行系統發育樹分析[26]。

第三、加快推進社會體育改革，構建社區體育發展機制。經過十幾年體育改革發展，我國體育社會化進度在加快，體育體制開始分化，社會普及程度較高的足球、籃球、排球等項目實現職業化發展。計劃經濟體制下的“單位體育”逐漸向社區體育轉變，體育社團、體育協會、體育俱樂部等體育社會組織的發展，推動社區體育發展，社區體育成為全民健身運動發展基本載體。

2 結果與分析

2.1 3 種殺蟲劑對梨小的毒力測定

采用藥膜法進行吡蟲啉、阿維菌素和高效氯氟氰菊酯對梨小的室內毒性測定，計算出3 種殺蟲劑的毒力回歸方程以及LC10、LC30、LC50濃度（表1）。

表1 3 種殺蟲劑的室內毒力測定Tab.1 Indoor toxicity determination of three insecticidesμg/mL

2.2 梨小食心蟲轉錄組測序與組裝結果

梨小轉錄組測序數據中，Q20 堿基和Q30 堿基的比值分別＞96%和＞87%，說明RNA-Seq 數據質量較好，測序數據準確可靠（表2）。

表2 測序數據統計Tab.2 Sequencing data statistics

經過組裝后產生了101 873 條Unigenes，總長度、平均長度以及覆蓋50%所有核苷酸的最大序列重疊群長度（N50）分別為202 980 197、1 992、3 037 bp，GC 含量為40.64%（表3）。Unigenes 其不同長度的具體分布如圖1 所示，可看出＞3 kb 的Unigenes 最多，有22 913 條，說明組裝結果較好。

圖1 梨小食心蟲轉錄組拼接后的Unigenes 長度分布Fig.1 Length distribution of the assembled Unigenes in the G. molesta transcriptome

表3 組裝結果統計Tab.3 Summary of assembly results

2.3 梨小食心蟲轉錄組基因功能注釋

將梨小食心蟲轉錄組數據組裝后的Unigenes序列與NR、NT、KEGG、KOG、SWISS-PROT、Pfam和GO數據庫進行BLAST比對，對獲得的101 873條Unigenes進行功能注釋，其中通過NR、NT、SWISSPROT、Pfam、KEGG、KOG、GO 數據庫比對獲得注釋的Unigenes 分別有57 209、39 186、43 749、46 168、46 678、42 247、27 486 條，共有66 658 條Unigenes 成功得到注釋，占總Unigenes 的65.43%（表4）。根據NR 數據庫比對，序列同源性最高的物種是大蠟螟Galleria mellonella（9.35%），其次是粉紋夜蛾Trichoplusia ni（9.32%）、棉鈴蟲Helicoverpa armigera（9.28%）、海波斯莫科馬屬蛾Hyposmocoma kahamanoa（9.16%）和臍橙螟蛾Amyelois transitella（8.13%）（圖2）。

圖2 梨小食心蟲轉錄組Unigenes 在NR 數據庫中的物種分布Fig.2 Species distribution of Unigenes in the G. molesta transcriptome in NR database

表4 梨小食心蟲轉錄組Unigene 的功能注釋信息Tab.4 Functional annotation information of G. molesta transcriptome

使用基因本體（GO）方法，將注釋基因分為3 類（共24 個功能組）：生物過程、細胞組成和分子功能。在生物過程中，細胞過程包含了最多的Unigenes；在細胞組分中，膜部分占主導地位；在分子功能分類中，注釋基因數最多的是結合和催化活性（圖3）。與KEGG 數據庫比對獲得注釋的46 678 條Unigene主要涉及6 大類功能，包括細胞過程、環境信息處理、遺傳信息處理、人類疾病、新陳代謝和生物系統。最具代表性的途徑是信號轉導以及運輸和分解代謝（圖4）。

圖3 梨小食心蟲轉錄組Unigenes GO 功能分類Fig.3 GO functional classification of Unigenes in the G. molesta transcriptome

圖4 梨小食心蟲轉錄組Unigenes KEGG 功能分類Fig.4 KEGG functional classification of Unigenes in the G. molesta transcriptome

2.4 梨小食心蟲差異表達基因及功能注釋

采用生物信息學方法分別篩選3 種農藥處理的梨小食心蟲差異表達基因，阿維菌素處理梨小后篩選出825 個差異表達基因，其中639 個上調基因和186 個下調基因（圖5-A），并對其功能進行注釋，注釋到結合（binding）的Unigenes 最多（圖5-B）。

圖5 差異表達量倍數變化及功能注釋Fig.5 Differential expression fold change and functional annotation

吡蟲啉處理后篩選出641 個差異表達基因，包括333 個上調基因和308 個下調基因（圖5-C），功能注釋到結合（binding）的Unigenes 最多（圖5-D）。

高效氯氟氰菊酯處理后篩選出1 391 個差異表達基因，包括538 個上調基因和853 個下調基因（圖5-E），并對其功能進行注釋，注釋到結合與催化活性（binding and catalytic activity）的Unigenes 最多（圖5-F）。

2.5 梨小食心蟲差異表達解毒酶基因篩選

經阿維菌素3 種濃度處理梨小后的轉錄組數據中，差異表達基因主要為P450s 和GSTs，其中有5 個P450，4 個上調，1 個下調；2 個GST 基因，1 個上調，1 個下調；其中P450 基因Unigene16556_All在LC10和LC50濃度下均顯著上調；GST 基因CL996.Contig4_All 在LC10和LC50濃度下均顯著下調（表5）。

表5 阿維菌素處理下差異表達的解毒基因Tab.5 Detoxification genes differentially expressed under avermectin treatment

表6 吡蟲啉處理下差異表達的解毒基因Tab.6 Detoxification genes differentially expressed under imidacloprid treatment

高效氯氟氰菊酯3 種濃度處理梨小后的轉錄組數據庫中，3 種解毒酶基因均有差異表達，其中有6 個P450 基因，4 個上調，2 個下調；2 個GST 基因均上調；3 個CarE 基因均上調（表7）。

將3 種殺蟲劑處理梨小后的差異表達基因進行匯總，共篩選到解毒代謝相關基因20 個，包括11 個P450 基因、5 個GST 基因、4 個CarE 基因，20 個差異表達基因的表達量熱圖聚類分析如圖6所示。分析發現，2個P450基因（Unigene12076_All、Unigene18672_All）在阿維菌素、吡蟲啉處理下均顯著上調；P450 基因CL4554.Contig2_All 在吡蟲啉、高效氯氟氰菊酯處理下顯著下調；P450 基因Unigene15542_All、GST 基因CL2631.Contig4_All在吡蟲啉、高效氯氟氰菊酯處理下顯著上調；P450基因Unigene19792_All 在吡蟲啉處理下顯著下調，在高效氯氟氰菊酯處理下顯著上調；P450 基因Unigene16556_All 在3 種農藥處理下均顯著上調。

圖6 差異表達基因表達量聚類熱圖Fig.6 Clustering heat map of differentially expressed gene expression

將11 條P450 基因序列與已知鱗翅目昆蟲的CYP 蛋白序列比對分析并構建進化樹，其中有7 個基因屬于CYP6 家族，2 個基因屬于CYP4 家族，2 個基因屬于CYP12 家族（圖7）。將5 條GST 基因序列與已知鱗翅目昆蟲的GST 蛋白序列比對分析并構建進化樹，均屬于sigma 家族（圖8），將4 條CarE基因序列與已知鱗翅目昆蟲的CarE 蛋白序列比對分析并構建進化樹，其中3 個聚在A 簇，另一個聚在C 簇（圖9）。

圖7 梨小食心蟲細胞色素P450 系統發育分析Fig.7 Phylogenetic analysis of cytochrome P450 in G. molesta

圖8 梨小食心蟲GST 系統發育分析Fig.8 Phylogenetic analysis of GST in G. molesta

圖9 梨小食心蟲4 個CarEs 系統發育分析Fig.9 Phylogenetic analysis of four CarEs in G. molesta

3 結論與討論

梨小食心蟲是世界范圍內重要的果樹害蟲之一，幼蟲蛀食梨、蘋果和桃等的果實以及嫩梢，給林果產業帶來了威脅。近年來，許多學者對梨小食心蟲進行了深入研究，但針對殺蟲劑亞致死劑量下梨小體內差異表達解毒酶基因鮮有報道。本研究通過Illumina Hiseq 2000 測序儀對阿維菌素、高效氯氟氰菊酯、吡蟲啉3 種殺蟲劑處理后梨小食心蟲成蟲轉錄組進行測序和差異表達分析，得到解毒代謝相關基因，為下一步的分子研究奠定基礎。本次研究共獲得101 873 條Unigene，N50 為3 037 bp，Q20堿基和Q30 堿基的比值分別＞96%和＞87%。研究表明，Q30 在80%以上就認為測序質量可靠[27]。因此，本研究的測序數據組裝質量和長度滿足轉錄組分析的基本要求。

本研究共計有66 658 條Unigenes 被NR、NT、SWISS-PROT 瑞士蛋白質數據庫（包括Pfam）、GO、KOG、KEGG 數據庫注釋，但仍有35 215 條序列沒有成功比對，這可能與拼接的序列片段過短、缺乏注釋信息有關[28]，也可能本身是非編碼序列或者是新的某種基因[29]。在NR 數據庫中，梨小注釋到大蠟螟的序列信息最多，可能與2 種昆蟲均屬鱗翅目，親緣關系相對較近有關。梨小轉錄組Unigenes 在 GO 數據庫得到27 486 條功能注釋，其中參與催化活性與結合活性的基因功能注釋較多，這與其他昆蟲轉錄組結果一致[30]。在KEGG 數據庫中得到46 678 條功能注釋，其中注釋到信號轉導通路以及運輸和分解代謝通路的基因較多，反映出梨小成蟲有較強的代謝活力，為深入挖掘、鑒定解毒代謝基因提供基礎數據。

山西農業大學植物保護學院郭艷瓊教授課題組前期對未做任何處理的梨小成蟲進行了轉錄組測序，得到77 條CYP 基因、46 條CarEs 基因、28 條GSTs 基因[31]，為進一步研究梨小食心蟲對殺蟲劑的代謝作用提供了寶貴信息?；诖?，進行3 種殺蟲劑處理的梨小轉錄組測序分析及差異表達解毒酶基因篩選，得到11 條CYP 基因、4 條 CarEs 基因、5 條GSTs 基因。研究表明，外源性和內源性化合物可以在一些昆蟲物種中誘導一些解毒酶基因的表達[32]，如黑腹果蠅[33]、小菜蛾[34]、埃及伊蚊[35]等。與未經氯菊酯處理的家蠅相比，暴露于氯菊酯的家蠅多個P450 基因上調[36]。噻蟲嗪誘導Q 型煙粉虱GST14的過表達[37]。家蠅CarEs 在氯菊酯處理后可進一步升高[38]。有研究認為，P450s 在昆蟲體內的誘導參與了昆蟲對環境的適應和對殺蟲劑抗性的形成[39]。CYP6A1基因的過度表達使褐飛虱其對吡蟲啉的抗性水平提高[40]。吡蟲啉誘導后，麥二叉蚜中有6 個屬于CYP3 家族的P450 基因表達量升高[41]。阿維菌素作用后，小菜蛾CYP340W1 基因表達量升高，RNA 技術證明此基因參與阿維菌素的代謝[42]。高效氯氟氰菊酯誘導后，斜紋夜蛾CYP301A1、CYP301B1、CYP367A18和CYP304F1可能參與了對高效氯氟氰菊酯的解毒代謝過程[43]。GST 的過表達也增強了殺蟲劑的代謝，柑橘全爪螨Panonychus citri中的一個GST 基因（PcGSTm5）在阿維菌素處理下高表達，參與阿維菌素的解毒[44]。殺蟲劑對CarE 基因的誘導使岡比亞按蚊[45]、埃及伊蚊[46]、朱砂葉螨[47]等昆蟲代謝解毒能力增強。本研究經阿維菌素處理梨小后，有4 個上調的P450 基因、1 個上調的GST 基因；經吡蟲啉處理梨小后，有5 個上調的P450 基因、1 個上調的GST 基因；經高效氯氟氰菊酯處理梨小后，有4 個上調的P450 基因、2 個上調的GST 基因、3 個上調的CarEs基因。由此推測，篩選得到的解毒酶基因可能使昆蟲解毒能力增強，從而參與對3 種殺蟲劑的代謝。

本研究通過對阿維菌素、吡蟲啉、高效氯氟氰菊酯3 種殺蟲劑不同亞致死劑量處理的梨小食心蟲轉錄組的測序分析以及差異表達基因的篩選，獲得與3 種殺蟲劑解毒代謝相關基因的信息，為相關基因的克隆、表達、功能研究及進一步研究梨小食心蟲解毒代謝機制奠定基礎。