花椒窄吉丁触角转录组及嗅觉相关基因分析

巩雪芳,谢寿安,杨 平,郭 莉,陈 迪,车显荣,王延来,吕淑杰

(西北农林科技大学林学院,陕西杨陵 712100)

花椒窄吉丁Agriluszanthoxylumi属鞘翅目(Coleoptera)吉丁科(Buprestidae)窄吉丁属Agrilus，是我国花椒林中最主要的毁灭性蛀干害虫(吴海,2006;权丽君,2016)。1987年在陕西省韩城市椒园首次发现花椒窄吉丁是致使花椒大量死亡的主要原因(薛耀宗,1989)。目前，主要分布于陕西、甘肃、山东、山西等地(党国刚等,2017)。幼虫主要蛀食枝干形成层，引起树势衰弱，产量品质下降(李辉,2013)，主要危害特征是受害部位虫疤松软湿润，有无色或淡褐色液体渗出，俗称“流油”(焦浩,2011)。成虫取食叶片补充营养，造成叶片的缺刻、空洞，树皮流胶、软化、腐烂、干枯、龟裂最后脱落至整株死亡(杜平,2012;王燕等,2012)。花椒窄吉丁的研究主要集中于虫情调查、种群动态、发生规律(李孟楼等,1990)及其化学防治(康克功等,2007)等方面。近年来，在花椒窄吉丁与寄主植物互作的分子机制方面，对寄主挥发物、后肠与粪便提取物对其引诱和驱避作用进行了初步研究(袁丽芳等,2016;刘绥鹏等,2017)；对花椒窄吉丁雌、雄成虫触角差异表达基因的进行了筛选(巩雪芳等,2019)；对气味结合蛋白基因AzanOBP3的克隆、原核表达及组织表达谱分析也进行了初探(巩雪芳等,2020)，但是对其作用的分子靶标尚不清楚。

触角是昆虫感受外界信号的主要器官，在长期进化过程中形成了能够感知各种化学信息物的高度灵敏且特异的嗅觉感受系统(王桂荣等,2002;修伟明等,2005;游灵等,2012;Britoetal.,2016)，在定位寄主、寻求配偶、觅食或产卵等种群生存及繁衍过程中均发挥重要作用(Zhangetal.,2015)。嗅觉识别过程需要多种蛋白的参与，主要有气味结合蛋白(odorant binding proteins,OBPs)、化学感受蛋白(chemosensory proteins,CSPs)、气味受体(odorant receptors,ORs)、离子型受体(ionotropic receptors,IRs)、感受神经元膜蛋白(sensory neuron membrane proteins,SNMPs)、气味降解酶(odorant-degrading enzymes,ODEs)等(詹文会,2018)。昆虫感受器淋巴液内发现两类作为环境信息分子结合载体的水溶性蛋白，即气味结合蛋白和化学感受蛋白(Angelietal.,1999;Yangetal.,2011)。

花椒窄吉丁是如何感知寄主气味物质并准确定位到唯一的寄主植物，在寻找寄主过程中嗅觉是如何发挥作用的？深入研究花椒窄吉丁嗅觉系统有助于阐明其定位寄主及寻找配偶时识别信息物质及其化学通讯与感受机制，而获得大量嗅觉相关基因是进行分子机制研究的前提。为此，本研究借助高通量测序对花椒窄吉丁成虫触角转录组进行生物信息学分析，获得了参与嗅觉感受的一些重要化学感受基因OBP,CSP,OR,IR和SNMP家族基因，为今后花椒窄吉丁嗅觉感受蛋白的作用机制提供理论支撑，对开展害虫行为调控剂的筛选及绿色防控技术研究具有重要的意义。

1 材料与方法

1.1 供试昆虫

2018年5-7月于陕西省蓝田县普化镇花椒园采集花椒窄吉丁成虫后带回实验室于25℃下饲养。选取1 000头活动能力较强的成虫，用灭菌三蒸水冲洗虫体表面，70%的酒精表面消毒30 s，灭菌三蒸水冲洗虫体2次后解剖镜下于冰上迅速切下触角雌雄分开收集，液氮冷冻保存于-80℃备用。

1.2 雌雄成虫触角总RNA的提取与质量检测

采用Trizol法分别提取花椒窄吉丁雌、雄成虫触角总RNA。利用1.0%琼脂糖凝胶电泳和凯奥K5500分光光度计检测RNA纯度(OD260/OD280比值)，安捷伦2100 RNA Nano 6000 Assay Kit检测完整性和浓度，质量合格的RNA样品保存于-80℃备用或者立即进行反转录。

1.3 文库构建与测序

1.2节样品RNA质量及纯度经检测合格后，用带有Oligo(dT)的磁珠富集mRNA，加入片段缓冲液到mRNA中使其片段化成为短片段，以片段后的mRNA为模板，用六碱基随机引物合成cDNA第1链，并加入缓冲液、dNTPs、RNaseH和DNA PolymeraseⅠ继续合成cDNA第2链，采用QIA Quick PCR试剂盒纯化并加EB缓冲液洗脱。对洗脱纯化后得到的双链cDNA再进行末端修复、加碱基A、加测序接头处理，回收目的大小片段进行PCR扩增。在获得clean reads后需采用Trinity软件对clean reads拼接和组装，获取分析所需的参考序列。转录本具有可变剪接的特点，组装出全长转录本，对测序得到的原始数据(raw data)过滤，获得高质量的数据(clean data)，在所有的转录本中，取最长的contigs作为unigenes。

1.4 转录组数据分析

高通量测序中，Q20值是指测序时碱基识别过程中，对所识别的碱基给出的错误概率，错误识别概率是1%，正确率99%；质量值Q30，错误识别概率0.1%，正确率99.9%；N50即覆盖50%核苷酸的最大unigene长度或最大序列重叠群长度。unigenes序列与NR,NT,Swiss-Prot,GO,KEGG,BLASTX,eggNOG,Pfam,TmHMM,SignalP,KO,Map,BLASTP和RNAMMER公共数据库比对，对获得的所有unigenes进行功能注释及结果统计。

1.5 嗅觉相关基因的鉴定与序列分析

BLAST同源搜索从触角转录组中鉴定OBPs,CSPs,Ors,IRs和SNMPs嗅觉相关基因。从NCBI数据库中下载苹果小吉丁A.mali和白蜡窄吉丁A.planipennis等其他鞘翅目昆虫的嗅觉相关基因核酸序列，采用ClustalX进行多序列比对，利用MEGA软件邻接法(neighbour-joining method,NJ method)构建系统发育树，对结果进行1 000次bootstrap验证。

1.6 嗅觉相关基因的表达分析

衡量基因表达量时，单纯以map到的read数计算基因的表达量，统计上是不合理的。本研究利用每百万条reads中来自于某基因每千碱基长度的reads数(reads per kilobase per million mapped reads,RPKM)(Mortazavietal.,2008)，对嗅觉相关基因表达量进行分析。

2 结果

2.1 测序结果及从头组装

花椒窄吉丁雌雄成虫触角转录组数据(SUB6796283)经过滤后分别获得43 124 393和45 854 648条clean read，组装得到36 209条基因和90 982条转录本，转录本和基因的N50分别为2 523和2 103 bp (图1)。样品clean data均达到96%以上，Q20均为100%，Q30均不少于92%，组装完整性较高(表1)。

图1 花椒窄吉丁雌雄成虫触角转录组unigenes长度分布Fig.1 Length distribution of unigenes in the antennal transcriptome of female and male adults of Agrilus zanthoxylumi

表1 花椒窄吉丁雌雄成虫触角转录组组装结果统计Table 1 Summary of assembly for the antennal transcriptome of female and male adults of Agrilus zanthoxylumi

2.2 基因功能注释

使用BLAST软件将获得的花椒窄吉丁雌雄成虫触角转录组unigenes 序列与NR,Swiss-Prot，GO,BLASTX和BLASTP等公共数据库比对(表2)，并对注释结果做比较韦恩图(图2)，共有36 209条unigenes得到注释，其中包括29 927条(54.21%)雌成虫触角的和26 706条(50.78%)雄成虫触角的。NR数据库注释到的unigene最多(15 070条,41.62%)，同时注释到NT,NR,BLASTX和BLASTP四大数据库的有2 682条。同源蛋白的物种来源分析结果表明，与赤拟谷盗的相似基因最多(19%)(图3)。

表2 花椒窄吉丁雌雄成虫触角转录组基因功能注释结果统计Table 2 Statistics of gene function annotation in the antennal transcriptome of female and male adults of Agrilus zanthoxylumi

图2 花椒窄吉丁雌雄成虫触角转录组基因功能注释结果韦恩图Fig.2 Venn diagram of gene annotation results of the antennal transcriptome of female and male adults of Agrilus zanthoxylumi

图3 NR数据库中注释的花椒窄吉丁雌雄成虫触角unigenes的物种分布Fig.3 Species distribution of unigenes in the antennal transcriptome of female and male adults of Agrilus zanthoxylumi annotated in the NR database

2.3 GO功能分类

GO功能分类统计(图4)表明，比对到11 614条unigenes，分为细胞组分(cellular component)、分子功能(molecular function)与生物学进程(biological process)三大类共57个分支。其中，生物学进程大类中参与细胞进程(cellular process)基因最多(72.55%)，其次是代谢进程(metabolic process)(58.15%)、生物学调控(biological regulation)(42.91%)，其余均在40%以下；细胞组分大类中细胞部分(cell part)基因最多(78.54%)，其次是细胞器(organelle)(49.19%)、细胞器部分(organelle part)(40.69%)、膜(membrane)(30.48%)与膜部分(membrane part)(27.00%)、大分子复合物(macromolecular complex)(26.47%)，其余均低于10%；分子功能大类中结合(binding)(70.57%)与催化活性(catalytic activity)(45.51%)的基因最多，其余均在10%以下。

图4 花椒窄吉丁雌雄成虫触角unigenes的GO分类Fig.4 GO classification of unigenes of female and male antennae of Agrilus zanthoxylumi

2.4 KEGG代谢通路分析

通过整合数据库，分析unigenes可能参与或涉及的代谢途径的注释结果(图5)发现有7 427条unigenes参与了5类代谢通路：细胞进程(cellular process)，环境信息处理(environmental information processing)，遗传信息处理(genetic information processing)，新陈代谢(metabolism)和有机系统(organismal system)。其中unigenes注释最多的代谢通路为信号转导(signal transduction)(815条)，翻译(translation)(586条)，内分泌系统(endocrine system)(466条)，运输和分解代谢(transport and catabolism)(440条)，折叠、分拣和降解(folding,sorting and degradation)(417条),碳水化合物代谢(carbohydrate metabolism)(397条)。

图5 花椒窄吉丁雌雄成虫触角unigenes的KEGG代谢通路Fig.5 KEGG metabolic pathways of the antennal unigenes of female and male adults of Agrilus zanthoxylumiA:细胞进程Cellular process;B:遗传信息处理Genetic information processing;C:环境信息处理Environmental information processing;D:新陈代谢Metabolism;E:有机体系统Organismal system.

2.5 嗅觉相关基因的筛选及进化分析

初步鉴定到7个OBP基因(AzanOBP1-7)，其中只有5个OBP基因具有全长开放阅读框，AzanOBP1,AzanOBP2,AzanOBP3,AzanOBP4和AzanOBP5具有6个保守的半胱氨酸位点，符合气味结合蛋白的特征(Cys-X15-39-Cys-X3-Cys-X21-44-Cys-X7-12-Cys-X8-Cys，X代表任意一种氨基酸)，为典型的OBP。5个具有全长开放阅读框的CSP基因(AzanCSP1-5)，3个OR基因(AzanOR1-3)，5个IR基因(AzanIR1-5)和5个SNMP基因(AzanSNMP1-5)。基于AzanOBP和AzanCSP基因以及其他类群昆虫OBP和CSP基因核苷酸序列构建系统发育树，发现花椒窄吉丁的7个OBP基因分支聚类较为分散(图6)，但都至少与一个鞘翅目昆虫的OBP同源基因聚类在一支，且序列一致性很高。AzanOBP1和AzanOBP3的序列一致性较高，且分别与白蜡窄吉丁AgrilusplanipennisAplaGOBP83a (GenBank登录号:XM_025974145.1)和苹果小吉丁AgrilusmaliAmalOBP2 (GenBank登录号:MF615491.1)聚到一支；AzanOBP2和AzanOBP4均与苹果小吉丁的AmalOBP3 (GenBank登录号:MF615492.1)亲缘关系近，三者聚到同一支。

图6 应用邻接法基于核苷酸序列构建的花椒窄吉丁雌雄成虫触角转录组中7个AzanOBP基因与其他鞘翅目昆虫OBP基因的系统发育树(1 000次重复)Fig.6 Phylogenetic tree of seven AzanOBP genes in the antennal transcriptome of female and male adults of Agrilus zanthoxylumi and OBP genes from other coleopteran insect species based on nucleotide sequences by using neighbor-joining method (1 000 replicates)基因来源物种Origin species of genes:AzanOBP:花椒窄吉丁Agrilus zanthoxylumi;AmalOBP:苹果小吉丁Agrilus mali;AplaOBP:白蜡窄吉丁Agrilus planipennis;DhelOBP:花绒寄甲Dastarcus helophoroides;HparOBP:暗黑鳃金龟Holotrichia parallela;PaenOBP:榆绿毛萤叶甲Pyrrhalta aenescens;PmacOBP:榆黄毛萤叶甲Pyrrhalta maculicollis;AcorOBP:铜绿丽金龟Anomala corpulenta;AchiOBP:星天牛Anoplophora chinensis;AglaOBP:光肩星天牛Anoplophora glabripennis;MaltOBP:松墨天牛Monochamus alternatus;TcasOBP:赤拟谷盗Tribolium castaneum;RdomOBP:谷蠹Rhyzopertha dominica;DarmOBP:华山松大小蠹Dendroctonus armandi.

花椒窄吉丁5个CSP基因聚类较为集中(图7)；AzanCSP1,AzanCSP2和AzanCSP3聚集到同一分支，5个CSP基因均与苹果小吉丁的AmalCSP亲缘关系最近。对于鞘翅目其他不同种的昆虫的CSP基因，也出现一定的聚类关系分别在进化树中形成多个独立的小分支；同一种昆虫的不同CSP基因聚类关系较分散也形成了多个独立的小分支。

图7 应用邻接法基于核苷酸序列构建的花椒窄吉丁雌雄成虫触角转录组中5个AzanCSP基因与其他鞘翅目昆虫CSP基因的系统发育树(1 000次重复)Fig.7 Phylogenetic tree of five AzanCSP genes in the antennal transcriptome of female and male adults of Agrilus zanthoxylumi and CSP genes from other coleopteran insect species based on nucleotide sequences by using neighbor-joining method (1 000 replicates)CSP基因来源物种Origin species of CSP genes:AzanCSP:花椒窄吉丁Agrilus zanthoxylumi;AmalCSP:苹果小吉丁Agrilus mali;AplaCSP:白蜡窄吉丁Agrilus planipennis;TmolCSP:黄粉虫Tenebrio molitor；AchiOBP:星天牛Anoplophora chinensis;TcasCSP:赤拟谷盗Tribolium castaneum;RdomCSP:谷蠹Rhyzopertha dominica;DarmCSP:华山松大小蠹 Dendroctonus armandi；DponCSP:中欧山松大小蠹Dendroctonus ponderosae;PmacCSP:榆黄毛萤叶甲Pyrrhalta maculicollis;PaenCSP:榆蓝叶甲Pyrrhalta aenescens；AglaCSP:光肩星天牛Anoplophora glabripennis;CbowCSP:大猿叶甲Colaphellus bowringi;MaltCSP:松墨天牛Monochamus alternatus;GdauCSP:沙葱萤叶甲Galeruca daurica;HparCSP:暗黑鳃金龟Holotrichia parallela;AcorCSP:铜绿丽金龟Anomala corpulenta;AgraCSP:棉铃象甲Anthonomus grandis;DhelCSP:花绒寄甲Dastarcus helophoroides.

2.6 嗅觉相关基因的表达量

利用RPKM值评估不同嗅觉相关基因在雌雄成虫触角中的表达量，结果表明(图8)，AzanOBP1和AzanOBP2在雌成虫触角中不表达，在雄成虫触角中微量表达；AzanOBP3在雄成虫触角中高丰度表达(RPKM>1 000)，在雌成虫触角中中量表达。AzanOBP4在雄成虫触角中中量表达(RPKM>400)，AzanOBP基因在雄成虫触角中的表达量排序为:AzanOBP3 (RPKM=1 160.57)>AzanOBP4(RPKM=244.93)>AzanOBP6 (RPKM=145.09)>AzanOBP7(RPKM=44.78)>AzanOBP5(RPKM=42.36)>AzanOBP1(RPKM=2.30)和AzanOBP2(RPKM=1.79)，在雌成虫触角中的表达量排序为:AzanOBP3(RPKM=205.25)>AzanOBP6(RPKM=42.94)>AzanOBP4(RPKM=39.69)>AzanOBP7(RPKM=30.66)>AzanOBP5(RPKM=14.04)>AzanOBP1(RPKM=0)和AzanOBP2(RPKM=0)(图8:A)。除AzanCSP1(RPKM>200)在雄成虫触角中中量表达外，其他AzanCSP基因在雌雄成虫触角中表达量均较低(RPKM<200);AzanCSP2在雌成虫触角中不表达(RPKM=0)，在雄成虫触角中微量表达(RPKM=3.064);AzanCSP基因在雄成虫触角中的表达量排序为:AzanCSP1(RPKM=302.16>AzanCSP5(RPKM=64.58)>AzanCSP4(RPKM=32.07)>AzanCSP3(RPKM=8.64)>AzanCSP2(RPKM=3.06)，AzanCSP基因在雌成虫触角中的表达量排序为：AzanCSP1(RPKM=207.85)>AzanCSP5(RPKM=79.65)>AzanCSP4(RPKM=27.21)>AzanCSP3(RPKM=4.24)>AzanCSP2(RPKM=0)(图8:B)。气味受体基因AzanOR基因在雌成虫触角中的表达量排序为:AzanOR1(RPKM=3.89)>AzanOR2(RPKM=0.28)>AzanOR3(RPKM=0.18)，AzanOR基因在雄成虫触角中的表达量排序为：AzanOR1(RPKM=18.13)>AzanOR3(RPKM=2.81)>AzanOR2(RPKM=1.33)(图8:C)。离子型受体基因AzanIR基因在雌成虫触角中的表达量排序为:AzanIR4(RPKM=2.86)>AzanIR2(RPKM=1.89)>AzanIR3(RPKM=1.73)>AzanIR1(RPKM=0.95)>AzanIR5(RPKM=0.14)，在雄成虫触角中的表达量排序为:AzanIR2(RPKM=6.08)>AzanIR4(RPKM=4.21)>AzanIR1(RPKM=3.70)>AzanIR3(RPKM=2.08)>AzanIR5(RPKM=0.48)(图8:D)。神经元感受膜蛋白基因AzanSNMP5在雌成虫触角中不表达(RPKM=0)，在雄成虫触角中微量表达(RPKM=0.24)，AzanSNMP基因在雌成虫触角中的表达量排序为:AzanSNMP2(RPKM=19.90)>AzanSNMP1(RPKM=0.42)>AzanSNMP3(RPKM=0.05)>AzanSNMP4(RPKM=0.02)>AzanSNMP5(RPKM=0)，在雄成虫触角中的表达量排序为:AzanSNMP2(RPKM=105.17)>AzanSNMP4(RPKM=0.46)>AzanSNMP1(RPKM=0.32)>AzanSNMP5(RPKM=0.24)>AzanSNMP3(RPKM=0.06)(图8:E)。

图8 花椒窄吉丁雌雄成虫触角转录组中嗅觉相关基因的RPKM值Fig.8 RPKM values of olfaction-related genes in the antenna of female and male adults of Agrilus zanthoxylumiA:AzanOBP;B:AzanCSP;C:AzanOR;D:AzanIR;E:AzanSNMP.AnF:雌成虫触角Antenna of female adult;AnM:雄成虫触角Antenna of male adult.

3 讨论与结论

嗅觉相关基因的鉴定可为基于反向化学生态学的研究提供潜在分子靶标(Zhuetal.,2017)，利用昆虫行为调节剂和分子干扰技术调控昆虫个体间、昆虫与寄主间的化学通讯来防控害虫是国内外研究的热点(杨海博等,2018)。RNA-Seq技术已成为转录组研究的重要手段，成功运用于非模式昆虫功能基因挖掘(郑海霞等,2018)，为深入研究昆虫嗅觉相关基因功能奠定了坚实的基础。本研究对花椒窄吉丁雌雄成虫触角转录组进行测序，共获得90 982条转录本。GO功能分类表明，分子功能分类中参与结合的基因占比最多(70.57%)，猜测花椒窄吉丁成虫触角中存在大量结合化学信息物质的蛋白。仍有59%的unigenes未被注释到公共数据库，这些序列可能是一些非编码序列或者新发现的基因序列，也可能是因为某些unigene长度过短并未能与公共数据库中的已有基因序列比对到，此外，也与公共数据库中现有的鞘翅目昆虫吉丁科的参考序列的数量不够丰富有关。本研究发现少数仅在花椒窄吉丁触角中存在的气味结合蛋白基因，这些特异表达的基因可能与花椒窄吉丁物种生境中所特有的化学物质有关，与它的强隐蔽性和对寄主专一性也紧密关联。KEGG数据库对花椒窄吉丁成虫触角的unigenes可能参与或涉及的代谢途径的注释结果分析表明，花椒窄吉丁信号感受与传递的主要场所是触角。前期我们已经对雌雄成虫触角中这些靶标基因的表达差异进行了初步探讨(巩雪芳等,2019)。本研究中AzanOBP3在触角中高丰度表达(图8:A)，与前期的研究结果(巩雪芳等,2020)不一致，主要是因为本研究对象是花椒窄吉丁触角(不包含头部)，而前期目的基因(AzanOBP3)在触角中的表达研究中的研究对象是花椒窄吉丁头部(包含触角)，且根据系统发育分析发现AzanOBP3与鞘翅目昆虫的性信息素结合蛋白聚类到同一分支，初步猜测并推断其很可能是与花椒窄吉丁性信息素相关的气味结合蛋白(巩雪芳等,2020)，所以在雌雄成虫触角中的表达量差异是显著的，在雌雄成虫头部的表达量差异是不显著的。

筛选到7个候选OBP基因，相比于其他鞘翅目昆虫的触角转录组研究，本研究所获得的嗅觉基因数量较偏少，除了受NCBI基因数据库(https:∥www.ncbi.nih.gov/)中参考序列的数量的影响外，还可能是因为化学感受基因中的气味结合蛋白和化学感受蛋白两个家族进化速率较高，这些基因进化很可能与识别特定植物化合物有关，本研究中的研究对象花椒窄吉丁为单食性昆虫，目前发现寄主植物只有花椒一种，这可能造成其OBP和CSP基因不如其他昆虫的数量多。发现花椒窄吉丁OBP并非都是典型的OBP(含有6个保守的半胱氨酸)(巩雪芳等,2019)，这7个OBP基因中，只有5个OBP基因具有全长开放阅读框；AzanOBP1,AzanOBP2,AzanOBP3,AzanOBP4和AzanOBP5的编码蛋白具有6个保守的半胱氨酸位点，符合气味结合蛋白的特征(Cys-X15-39-Cys-X3-Cys-X21-44-Cys-X7-12-Cys-X8-Cys，X代表任意一种氨基酸)，为典型的OBP；AzanOBP6的编码蛋白序列中总共有10个半胱氨酸位点，是非典型的气味结合蛋白(atypical OBPs)；AzanOBP7的编码蛋白具有4个保守的半胱氨酸位点，初步鉴定为DimerOBP类。

利用RPKM值对基因表达量的评估结果表明(图8)，AzanOBP1和AzanOBP2在雌成虫触角中不表达(RPKM=0)，在雄成虫触角中微量表达；AzanOBP3在雄成虫触角中高丰度表达(RPKM>1 000)，在雌成虫触角中中量表达(RPKM=205.25)，推测其在触角中高丰度表达的OBPs功能与雄成虫识别同类异性释放的信息素或寄主植物释放的挥发物有关。

本研究首次构建了花椒蛀干害虫花椒窄吉丁的雌雄成虫触角转录组数据库，并成功筛选得到了部分嗅觉相关基因，为将来利用嗅觉靶标无公害防治农林害虫提供了重要分子信息；基于KEGG数据库把从已经完整测序的花椒窄吉丁触角转录组中得到的基因目录与更高级别的细胞、物种和生态系统水平的系统功能关联起来。研究结果为鞘翅目其他近缘物种昆虫的触角转录组的研究提供了丰富的转录本基础数据信息，对鞘翅目吉丁科非模式昆虫的研究有重要的推动意义。