抗氟苯虫酰胺小菜蛾差异表达基因及其通路

王成花，孙诗晴，徐巨龙，赵小龙，薛超彬

（山东农业大学植物保护学院/农药毒理与应用技术省级重点实验室，山东泰安 271018）

0 引言

【研究意义】小菜蛾（Plutella xylostella）属鳞翅目菜蛾科，主要危害甘蓝等十字花科植物，幼虫取食叶片，苗期集中危害心叶，严重影响蔬菜的产量及品质，成为世界性的蔬菜害虫。氟苯虫酰胺（flubendiamide）是由日本农药株式会社和德国拜耳公司联合开发的一种邻苯二甲酰胺类杀虫剂，其主要作用靶标为昆虫鱼尼丁受体（ryanodine receptor，RyR），该药剂已广泛应用于鳞翅目害虫的防治，效果优良，但近年来小菜蛾等鳞翅目害虫对该药剂的抗性越发突出。因此，明确小菜蛾对氟苯虫酰胺的抗性机理对合理控制小菜蛾及其抗药性发展具有重要意义。【前人研究进展】2008—2011年，泰国Tha Muang地区小菜蛾对氟苯虫酰胺的抗性由1.5倍发展到4 817倍[1]；2012年，泰国Sai Noi地区小菜蛾对氟苯虫酰胺的抗性倍数超过750倍；菲律宾 Sudlon地区小菜蛾田间种群对氟苯虫酰胺的抗性超过1 300倍[2]，2015年该Sudlon种群对氟苯虫酰胺的抗性超过10 000倍[3]。中国广东增城地区的小菜蛾田间种群对氟苯虫酰胺的抗性高达1 779倍[4]。除此之外，其他鳞翅目害虫也对氟苯虫酰胺产生了不同程度的抗性，如茶小卷叶蛾（Adoxophyes honmai）对氟苯虫酰胺的抗性达 105倍[5]；番茄潜叶蛾（Tuta absoluta）对氟苯虫酰胺的抗性超过1 000倍[6]。抗性机制研究发现，泰国和菲律宾的小菜蛾田间高抗种群中鱼尼丁受体基因（PxRyR）的4 946位均存在一个甘氨酸（G）到谷氨酸（E）的突变[7]，其位于受体基因羧基末端的跨膜区附近，且在所有昆虫鱼尼丁受体中高度保守，被认为是导致小菜蛾对氯虫苯甲酰胺产生抗性的主要原因。此外，解毒酶活性的提高，如细胞色素 P450单加氧酶活性的提高（P450）、羧酸酯酶（CarE）和谷胱甘肽S-转移酶（GSTs）活性的提高[8-9]；PxRyR中G4946E、E1338D、Q4594L和I4790M等4个位点的协同突变[10]和PxRyRmRNA转录表达变化等[4,11-12]，也是小菜蛾对氯虫苯甲酰胺产生抗性的重要机制。【本研究切入点】在小菜蛾对双酰胺类杀虫剂抗性机制的研究中，多数报道均以氯虫苯甲酰胺为药剂对象，但小菜蛾对氟苯虫酰胺与氯虫苯甲酰胺的抗性机制是否一致，目前还未有明确的认识。【拟解决的关键问题】以小菜蛾抗氟苯虫酰胺品系为试材，从转录组水平研究小菜蛾对氟苯虫酰胺的抗性机制，进一步揭示其抗性机理。

1 材料与方法

试验于 2016年在山东农业大学和深圳华大基因科技服务有限公司完成。

1.1 供试材料

小菜蛾敏感品系（S）于2006年采自山东农业大学南校区实验站园，在室内不接触任何药剂，用甘蓝苗长期饲养且稳定繁殖；抗氟苯虫酰胺品系（Rh）是敏感品系小菜蛾在氟苯虫酰胺药剂的持续选择压力下筛选至28代（Rh28），其抗性倍数为397.08倍；田间抗性种群（Rz）于2012年10月采集于广州增城市周边蔬菜田，在室内无农药暴露情况下繁殖到第 36代（Rz36），与 S品系相比，该种群对氟苯虫酰胺药剂的抗性倍数为98.12倍。

1.2 RNA提取及检测

选取标准一致的小菜蛾3龄幼虫，迅速用液氮冷冻，用RNAsimp Total RNA Kit试剂盒提取RNA，并利用 1.0%琼脂糖凝胶电泳检测，将得到的质量较高RNA样品送至深圳华大基因科技服务有限公司，利用Illumina HiSeqTM2000进行高通量测序。每个样品提供3次重复。并对质检合格的数据进行对比，包括抗氟苯虫酰胺品系与敏感品系对比（Rh28-vs-S）、抗氟苯虫酰胺田间抗性种群与敏感品系对比（Rz36-vs-S），及抗氟苯虫酰胺品系与抗氟苯虫酰胺田间抗性种群对比（Rh28-vs-Rz36），并进行相关生物信息分析。

1.3 基因定量方法

将3组样品间两两比较，用RSEM工具[13]进行基因表达定量。表达定量的结果以FPKM为单位，具体计算公式如下：

式中，FPKM（A）表示为基因A的表达丰度，C为唯一比对到基因A的fragments数，N为唯一比对到参考基因的总fragments数，L为基因A编码区的碱基数。FPKM法能够消除由于基因长度和测序质量差异的影响，使得到的基因量更加准确地比较不同样品间的表达差异。

1.4 明确差异基因

为了挖掘样本间的显著差异基因，基于泊松分布的分析方法，并参照 AUDIC等[14]的方法，假设观测到基因A对应的reads数为x，已知在一个大文库中，每个基因的表达量只占所有基因表达量的一小部分，在这种情况下，x的分布服从泊松分布：

之后，对差异检验的P值作多重假设检验校正，通过控制FDR（False Discovery Rate）来决定P值的域值[15]。在得到差异检验的FDR值同时，根据基因的表达量（FPKM值）计算该基因在不同样本间的差异表达倍数。FDR值越小，差异倍数越大，则表明表达差异越显著。本研究中，显著差异表达基因定义为FDR≤0.001且倍数差异在2倍以上的基因。

1.5 GO显著性富集分析

基因本体（Gene Ontology，GO）共有3个本体，分别描述基因的分子功能（molecular function）、细胞成分（cellular component）和参与的生物过程（biological process）。GO分析给出与参考基因组背景相比，在差异表达基因中显著富集的GO功能条目，从而筛选出差异表达基因与哪些生物学功能显著相关。本研究首先把所有差异表达基因向Gene Ontology数据库（http://www.geneontology.org/）的各个term映射，计算每个term的基因数目，然后应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著富集的GO条目。通过GO功能显著性富集分析能确定差异表达基因行使的主要生物学功能。

1.6 KEGG显著性富集分析

京都基因与基因组百科全书数据库（KEGG）是关于通路（Pathway）的主要公共数据库，Pathway显著性富集分析应用超几何检验，该分析的计算公式同GO功能显著性富集分析。在差异表达基因中只有满足Q值≤0.05才属于显著富集的 Pathway。通过Pathway显著性富集分析，能够确定差异表达基因参与的最主要的生化代谢途径以及信号转导途径。

2 结果

2.1 差异表达基因

通过将抗氟苯虫酰胺品系（Rh28）与敏感品系（S）的RNA-seq数据对比发现，Rh28中共有540个显著差异表达基因，其中上调表达的基因411个，下调表达的有129个（附表1）；同样的，将田间抗性种群（Rz36）与S品系的RNA-seq数据对比发现，Rz36中共有439个显著差异表达基因，其中上调表达的基因299个，下调表达的有140个；进一步将Rh28品系与Rz36种群数据对比发现，Rh28中共有75个显著差异表达基因，其中上调表达的48个，下调表达的有27个。

2.2 显著差异表达基因的基因本体（GO）分析

在得到每个差异基因的GO注释后，用WEGO软件[16]对 3个品系/种群间的对比（Rh28-vs-S、Rz36-vs-S和 Rh28-vs-Rz36）得到的差异基因 DEGs从生物过程/细胞组件和分子功能3大类进行分类统计，从宏观上认识差异基因的功能分布特征（附图1）。其中生物黏附（biological adhesion）、生物调节（biological regulation）、代谢进程（metabolic process）、免疫系统进程（immune system process）、对刺激物的反应（response to stimulus）、催化活性（catalytic activity）、酶调节活性（enzyme regulator activity）、结合因子（binding）、受体激活（receptor activity）、转运活力（transporter activity）、核酸绑定的转录因子激活（nucleic acid binding transcription factor activity）等功能都均与小菜蛾的抗性及免疫系统具有相关性。

2.3 显著差异表达基因的通路（pathway）分析

利用 KEGG数据库，对 Rh28-vs-S、Rz36-vs-S和Rh28-vs-Rz363组样品的差异表达基因可能参与或涉及的通路（pathway）进行分析，结果发现，Rh28与S品系之间共有14 527个被注释到188条pathway上，其中包含基因数量最多的一条是代谢通路（metabolic pathway），共有2 591个基因，占总数的17.84%，与抗性相关的 pathway还包括药物代谢——细胞色素P450（drug metabolism-cytochrome P450）、酪氨酸代谢（tyrosine metabolism）、谷胱甘肽代谢（glutathione metabolism）、促分裂素原活化蛋白激酶信号通路（MAPK signaling pathway）、外来物质代谢的细胞色素 P450（metabolism of xenobiotics by cytochrome P450）、心肌收缩（cardiac muscle contraction）、GABA能突触（GABAergic synapse）、咖啡因代谢（caffeine metabolism）、黑色素生成（melanogenesis）等。

Rh28与S品系之间540个显著差异的基因有289个被注释到通路当中，其中显著性富集（Q值＜0.05）的pathway只有1个，即过氧化物酶体增殖物激活受体信号通路（PPAR signaling pathway）。Rz36与S品系之间的差异基因中有237个被注释到了KEGG数据库中，得到172条pathway，其中显著富集的也只有1个，即 ubiquitin mediated proteolysis。而 Rh28与 Rz36品系之间的75个显著差异性基因有56个注释到了通路当中，得到71条pathway，没有显著富集的通路。KEGG富集排名前20的pathway散点图见附图2。

2.4 抗性显著上调基因

统计Rh28-vs-S，Rz36-vs-S和Rh28-vs-Rz363组间显著上调的基因数据发现，与S品系相比，Rh28品系显著上调的411个基因中，有218个在Rz36种群中同样显著上调表达。Rh28相对于Rz36显著上调的48个基因中，有 25个是Rh28品系特有的上调表达基因，而在Rz36种群中表现不显著（图1）。针对上述218及25个上调基因展开进一步研究。

2.5 抗性显著上调基因的GO分析

图1 上调表达的基因在小菜蛾不同种群中的分布Fig. 1 Distribution of up-regulated express genes in different populations of P. xylostella

Rh28-vs-S和Rz36-vs-S显著上调的218个交集基因，在GO 3大本体所占的比例及数目分布情况如图2-a所示；Rh28-vs-S和Rh28-vs-Rz36显著上调的25个交集基因如图2-b所示。结果表明，这部分交集差异基因在GO本体中的分布情况与前述的整体的分布情况基本一致，集中在代谢进程、对刺激物的反应、生物调节、催化活性、结合因子等方面。数据分析表明，RNA-Seq获得的小菜蛾上调的DEGs中，参与生物过程的基因占比较高，达50%以上。

对Rh28-vs-S与Rz36-vs-S交集的218个上调基因，注释到GO数据库BP中，将具有相同生物学功能的基因整合到一起（附表2）。部分差异基因在BP中所占的比例及数目见图3，其中主要涉及脂质代谢、蛋白质代谢、氨基酸转运和代谢、能量代谢、对环境压力的反应、对刺激物的反应、防御反应机制及转运等过程，附表3列举了Rh28-vs-S与Rh28-vs-Rz36交集上调的25个基因所参与的BP。

2.6 抗性显著上调基因的KEGG分析

为了较为准确地定位显著差异表达上调基因参与的生物学功能，针对其所参与的 pathway，通过查找统计KEGG数据库，发现Rh28-vs-S和Rz36-vs-S显著上调的218个交集基因中，能够注释到数据库的基因仅有 108个，主要集中在代谢通路上，有 110个（50.45%）基因在KEGG数据库中没有详细的代谢通路图。小菜蛾响应氟苯虫酰胺药物选择压力下差异显著的表达基因富集到的通路有：代谢通路、PPAR信号通路、MAPK信号通路、Wnt信号通路等。同一个基因可能同时注释到多条pathway之中，经统计，列举了部分可能与抗性相关的通路（图4）。

图2 上调交集基因在GO分析中的分布Fig. 2 Distribution of up-regulated intersection genes in GO analysis

图3 上调交集基因部分BP注释Fig. 3 Up-regulated intersection genes annotated to the corresponding BP

图4 上调交集基因中与小菜蛾抗氟苯虫酰胺相关的部分通路Fig. 4 Partial pathway associated with the resistance to flubendiamide in P. xylostella

2.7 抗性相关基因的聚类分析

利用高通量数据，对与杀虫剂抗性相关的基因进行了动态的监测，并从中筛选出可能与抗药性相关的通路基因，绘制了聚类热图（图5）。通过聚类热图的红色区域，可以看出Rh28品系与Rz36种群显著表达量上调的基因，很多集中在内质网上蛋白质的加工、酪氨酸代谢、咖啡因代谢、Wnt信号通路。

3 讨论

随着高通量测序技术日益成熟，为小菜蛾抗药性研究提供了新方法，通过进行小菜蛾基因的筛选与鉴定，进而分析基因的差异表达，对深入阐述小菜蛾的抗性机制具有重要意义。本文通过基因序列筛选，从小菜蛾S、Rh28和Rz36品系/种群的RNA-Seq序列信息中鉴定出多个与杀虫剂靶标及解毒代谢相关的基因，对这些差异表达基因在GO以及KEGG数据库进行功能预测与分类注释，从获得的2万多条差异表达基因序列中，筛选出满足差异显著性表达条件的只有几百条，这部分显著差异的序列，在GO数据库中差异表达基因显著富集的生物过程均在应激反应（response to stress）和防御反应（defense response）中，导致小菜蛾对氟苯虫酰胺抗性产生的基因主要集中在这部分生物过程中，其生物学功能可能与小菜蛾的抗性之间存在一定的相关性。

通过对差异基因的Pathway统计，信号通路主要集中于PPAR信号通路、蛋白质在内质网上的加工、甘油磷脂代谢、解毒酶代谢等。PPAR是一类配体激活的核转录因子超家族成员，包括PPAR-α、PPAR-β和PPAR-γ3种表型，其中以PPAR-γ的研究最为深入。现有的研究已证实 PPAR-γ主要参与调控脂肪和糖类的代谢、能量平衡等多种生物学功能的调节，小菜蛾的体内脂肪酸 LPL、胆固醇、甘油磷脂等各种脂类的代谢相关以及信号受体激活等一系列连锁反应[17]。蛋白质在内质网上的加工，主要参与协助一些蛋白质装配或者折叠，这类蛋白主要是Sec23蛋白和Hsps分子伴侣蛋白，它们与能量代谢、次级代谢途径密切相关。参与各种氨基酸、嘌呤，代谢产物（萜类、半萜类）代谢及生物合成、糖酵解；能量代谢过程中K/Na离子通道ATPase表达上调。研究表明，氯虫苯甲酰胺能显著影响甜菜夜蛾幼虫中枢神经元的兴奋性，使电压门控钠通道电流密度和通道功能发生改变，试验结果说明钠离子通道亦是氯虫苯甲酰胺的作用靶标之一[18]。本文通过聚类热图分析发现，上调表达的基因很多集中在内质网上蛋白质的加工、酪氨酸代谢和 Wnt信号通路中，与前人研究结果基本一致。

图5 小菜蛾抗氟苯虫酰胺相关基因表达分析Fig. 5 Expression profiling of genes associated with the resistance to flubendiamide in P. xylostella

谷胱甘肽-S-转移酶（GSTs）作为一种超基因家族酶类，可催化昆虫因各种化学药剂产生的有害物质的代谢。依据GSTs所处的细胞位置和功能，主要分为微粒体、线粒体和胞质型，而在昆虫中这种酶主要属于胞质型，其家族中的Delta家族（昆虫中特有的GSTs家族之一）[19]，在生物体内起到的作用主要包括催化活性、谷胱甘肽过氧化物酶活性、谷胱甘肽转移酶活性及裂解酶活性等，这些均与昆虫的代谢活性相关。研究发现，抗氯虫苯甲酰胺品系小菜蛾的GSTs 活性是敏感品系的3.34倍，因此GSTs可能在小菜蛾对氯虫苯甲酰胺的代谢抗性中发挥着主要作用[9]，且本文的研究也表明主要富集于谷胱甘肽代谢的 2个基因（XM_011569702.1、XM_011557174.1）表现为显著表达上调。

细胞色素 P450在昆虫生长发育及繁殖等方面起着重要作用，并且在次生物质和杀虫剂代谢方面也有重要作用，在COG数据库中，细胞色素P450被注释到了次生代谢产物的生物合成（second-metabolites biosynthesis）、运输和代谢（transport and catabolism）功能中，这也充分说明了细胞色素 P450是昆虫体内主要的代谢解毒酶之一[20]。其中编码 CYP2的基因（XM_011556598.1）表达量上调显著，而且已有研究报道CYP家族的细胞色素P450解毒酶在一些昆虫体内与抗药性密切相关[21]。

经过部分筛选分析，从小菜蛾中鉴定出多个热激蛋白（Hsp），发现该抗性小菜蛾体内的 Hsp的多个家族（如 Hsp40、Hsp70、Hsp90）表达量显著上调，并且参与了多个生物工程，笔者认为Hsp很可能与昆虫的抗性以及对环境的应激性有关，而且在GO数据库中，这些Hsp家族差异基因均注释进入了“对刺激物的反应”和“细胞蛋白质代谢过程（cellular protein metabolic process）”功能中，同时“对刺激物的反应”的功能注释更加表明了 Hsp可能与小菜蛾抗药性的产生有着密切的联系。Hsp基因家族可被压力诱导，在环境压力因素下，如重金属离子、极端温度、营养缺失、氧自由基及细菌或者病毒感染都有可能诱导热激蛋白的过量表达[22]。在粉虱类的抗逆基因研究中已经有过类似的报道，可能由于昆虫的这种蛋白可以应对环境中存在的极端胁迫压力[23-24]，导致一些入侵性的有害生物在分布地域方面不断扩大[25]。已有研究证实在高温、低温、辐射、干旱和农药等的环境压力下，都能够诱导昆虫体内不同Hsp的表达，这证实了Hsp可以提高昆虫对不良环境的耐受能力，从而保护昆虫免受和少受胁迫的伤害。如温度驯化的南美斑潜蝇（Lirionyza huidobrebsis），其体内Hsp70基因大量表达，在性状方面耐热能力也显著提高，表明Hsp表达和昆虫耐热性的获得有着密切关系，而且耐热性的诱导和消退动力学与Hsp的诱导和降解存在着平行关系[26]。抗有机磷类杀虫剂的摇蚊（Chiromomus yoshimatsui）品系Hsp的表达水平与敏感品系相比，是敏感品系的2—3倍，表明Hsp上调增强了昆虫的抗药能力[27]。本研究发现，小菜蛾抗氟苯虫酰胺品系的谷胱甘肽-S-转移酶、细胞色素 P450解毒酶及参与昆虫多种生理反应的热激蛋白超基因家族等基因的表达量均显著上调，与前人研究基本一致。因此，今后需要对这些基因及其所在的通路进行深入研究，进一步明确其抗性机制。

昆虫靶标基因的mRNA表达量变化也是昆虫抗药性产生的原因之一。本研究中小菜蛾田间抗性种群鱼尼丁受体相关8个基因表达量均下调，而且有6个显著下调。LIN等[28]通过比较转录组测序发现，对氯虫苯甲酰胺产生不同抗性水平小菜蛾种群的鱼尼丁受体mRNA的表达量均显著下调；WAN等[29]利用RNAi干扰马铃薯甲虫（Leptinotarsa decemlineata）鱼尼丁受体的表达后，再用氯虫苯甲酰胺处理，马铃薯甲虫的死亡率明显降低；YANG等[30]在白背飞虱（Sogatella furcifera）的研究中也发现了类似结果。这些结果均表明鱼尼丁受体基因表达量的下调可能会增强昆虫对双酰胺类药剂的抗性。笔者实验室正在开展鱼尼丁受体表达差异基因的验证工作，以明确该类基因在小菜蛾对双酰胺类杀虫剂抗性中所起的作用。

4 结论

通过转录组测序（RNA-Seq）获得了小菜蛾抗氟苯虫酰胺差异表达基因及抗性显著上调表达基因，其主要富集于代谢过程、应激反应及对刺激的反应等条目中，这些基因的相互协同调控是小菜蛾对氟苯虫酰胺产生抗性的重要机制。

[1]SUKONTHABHIROM S, DUMRONGSAK D, JUMROON S,SAROCH T, CHAWENG A, TANAKA T. Update on DBM diamide resistance from the Thailand: situation and causal factors//SRINIVASAN R, SHELTON A M, COLLINS H L.The SixthInternational Workshop on Management of the Diamondback Moth and Other Crucifer Insect Pests. AVRDC-The World Vegetable Center,2011: 202-206.

[2]TROCZKA B, ZIMMER C T, ELIAS J, SCHORN C, BASS C,DAVIES T G E, FIELD L M, WILLIAMSON M S, SLATER R,NAUEN R. Resistance to diamide insecticides in diamondback moth,Plutella xylostella(Lepidoptera: Plutellidae) is associated with a mutation in the membrane-spanning domain of the ryanodine receptor.Insect Biochemistry and Molecular Biology, 2012, 42(11):873-880.

[3]STEINBACH D, GUTBROD O, LÜMMEN P, MATTHIESEN S,SCHORN C, NAUEN R. Geographic spread, genetics and functional characteristics of ryanodine receptor based target-site resistance to diamide insecticides in diamondback moth,Plutella xylostella.Insect Biochemistry and Molecular Biology, 2015, 63: 14-22.

[4]YAN H H, XUE C B, LI G Y, ZHAO X L, CHE X Z, WANG L L.Flubendiamide resistance and Bi-PASA detection of ryanodine receptor G4946E mutation in the diamondback moth (Plutella xylostellaL.).Pesticide Biochemistry and Physiology, 2014, 115:73-77.

[5]UCHIYAMA T, OZAWA A. Rapid development of resistance to diamide insecticides in the smaller tea tortrix,Adoxophyes honmai(Lepidoptera: Tortricidae), in the tea fields of Shizuoka Prefecture,Japan.Applied Entomology and Zoology, 2014, 49(4): 529-534.

[6]RODITAKIS E, VASAKIS E, GRISPOU M, STAVRAKAKI M,NAUEN R, GRAVOUIL M, BASSI A. First report ofTuta absolutaresistance to diamide insecticides.Journal of Pest Science, 2015,88(1): 9-16.

[7]TROCZKA B J, WILLIAMS A J, WILLIAMSON M S, FIELD L M,LÜMMEN P, DAVIES T G E. Stable expression and functional characterization of the diamondback moth ryanodine receptor G4946E variant conferring resistance to diamide insecticides.Scientific Reports, 2015, 5: 14680.

[8]WANG X L, WU Y D. High levels of resistance to chlorantraniliprole evolved in field populations ofPlutella xylostella.Journal of Economic Entomology, 2012, 105(3): 1019-1023.

[9]HU Z D, FENG X, LIN Q S, CHEN H Y, Li Z Y, YIN F, LIANG P,GAO X W. Biochemical mechanism of chlorantraniliprole resistance in the diamondback moth,Plutella xylostellaLinnaeus.Journal of Integrative Agriculture, 2014, 13(11): 2452-2459.

[10]GUO L, LIANG P, ZHOU X G, GAO X W. Novel mutations and mutation combinations of ryanodine receptor in a chlorantraniliprole resistant population ofPlutella xylostella(L.).Scientific Reports, 2014,4: 6924.

[11]SUN L N, CUI L, RUI C H, YAN X J, YANG D B, YUAN H Z.Modulation of the expression of ryanodine receptor mRNA fromPlutella xylostellaas a result of diamide insecticide application.Gene,2012, 511: 265-273.

[12]LI X X, GUO L, ZHOU X G, GAO X W, LIANG P. miRNAs regulated overexpression of ryanodine receptor is involved in chlorantraniliprole resistance inPlutella xylostella(L.).Scientific Reports, 2015, 5: 14095.

[13]LI B, DEWEY C N. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome.BMC Bioinformatics,2011, 12(1): 323.

[14]AUDIC S, CLAVERIE J M. The significance of digital gene expression profiles.Genome Research, 1997, 7(10): 986-995.

[15]BENJAMINI Y, YEKUTIELI D. The control of the false discovery rate in multiple testing under dependency.The Annals of Statistics,2001, 29(4): 1165-1188.

[16]YE J, FANG L, ZHENG H K, ZHANG Y, CHEN J, ZHANG Z J,WANG J, LI S T, LI R Q, BOLUND L, WANG J. WEGO: a web tool for plotting GO annotations.Nucleic Acids Research, 2006, 34(2):W293-W297.

[17]马晶晶, 章涛. PPARγ功能与疾病关系研究进展. 中国药理学通报,2012, 28(5): 601-604.MA J J, ZHANG T. Function of PPARγ and its relationship with diseases: a present review.Chinese Pharmacological Bulletin, 2012,28(5): 601-604. (in Chinese)

[18]郭直岳, 王国红, 甄静, 贺秉军. 氯虫酰胺对甜菜夜蛾神经元钠离子通道电流的影响. 河南农业科学, 2013, 42(10): 79-83.GUO Z Y, WANG G H, ZHEN J, HE B J. Effects of chlorantraniliprole on sodium channel currents of central neurons fromLaphygma exigue.Journal of Henan Agricultural Sciences, 2013,42(10): 79-83. (in Chinese)

[19]HAYES J D, FLANAGAN J U, JOWSEY I R. Glutathione transferases.Annual Review of Pharmacology and Toxicology,2005,45: 51-88.

[20]FEYEREISEN R. Insect P450 enzymes.Annual Reviews of Entomology, 1999, 44: 507-533.

[21]ZHANG L, GAO X, LIANG P. Beta-cypermethrin resistance associated with high carboxylesterase activities in a strain of house fly,Musca domestica(Diptera: Muscidae).Pesticide Biochemistry and Physiology, 2007, 89: 65-72.

[22]GEHRMANN M, BRUNNER M, PFISTER K, REICHLE A,KREMMER E, MULTHOFF G. Differential up-regulation of cytosolic and membrane-bound heat shock protein 70 in tumor cells by anti-inflammatory drugs.Clinical Cancer Research, 2004, 10:3354-3364.

[23]KIM K K, KIM R, KIM S H. Crystal structure of a small heat-shock protein.Nature, 1998, 394: 595-599.

[24]FRANCK E, MADSEN O, VAN RHEEDE T, RICARD G,HUYNEN M A, DE JONG W W. Evolutionary diversity of vertebrate small heat shock proteins.Journal of Molecular Evolution, 2004, 59: 792-805.

[25]QIN W, NEAL S J, ROBERTSON R M, WESTWOOD J T,WALKER V K. Cold hardening and transcriptional change inDrosophila melanogaster.Insect Molecular Biology, 2005, 14(6):607-613.

[26]HUANG L H, KANG L. Cloning and interspecific altered expression of heat shock protein genes in two leafminer species in response to thermal stress.Insect Molecular Biology, 2007, 16(4): 491-500.

[27]YOSHIMI T, MINOWA K, KAROUNA-RENIER N K,WATANABE C, SUGAYA Y, MIURA T. Activation of a stressinduced gene by insecticides in the midge,Chironomus yoshimatsui.Journal of Biochemical and Molecular Toxicology, 2002, 16(1):10-17.

[28]LIN Q S, JIN F L, HU Z D, CHEN H Y, YIN F, LI Z Y, DONG X L,ZHANG D Y, REN S X, FENG X. Transcriptome analysis of chlorantraniliprole resistance development in the diamondback mothPlutella xylostella.PLoS ONE, 2013, 8(8): e72314.

[29]WAN P J, GUO W Y, YANG Y, LU F G, LU W P, LI G Q. RNAi suppression of the ryanodine receptor gene results in decreased susceptibility to chlorantraniliprole in Colorado potato beetleLeptinotarsa decemlineata.Journal of Insect Physiology, 2014, 63(1):48-55.

[30]YANG Y, WAN P J, HU X X, LI G Q. RNAi mediated knockdown of the ryanodine receptor gene decreases chlorantraniliprole susceptibility inSogatella furcifera.Pesticide Biochemistry and Physiology, 2014, 108: 58-65.