南亚实蝇Bactrocera tau (Walker),属双翅目Diptera实蝇科Tephritidae,果实蝇属Bactrocera Macquart[1,2,3],为害丝瓜、南瓜、苦瓜等13种瓜类,柑桔、芒果等共计70多种热带和亚热带瓜果。中国主要分布于广东、广西、海南、福建、云南、贵州、四川、湖北、浙江、江西等地[4]。因其抗逆性强,个体食量大,发育历期短,产卵期长,在很多地区其危害甚至超过橘小实蝇B. dorsalis (Hendel)[5]。
目前,对南亚实蝇的防治措施主要借鉴橘小实蝇B. dorsalis和瓜实蝇B. cucuribitae的4步防控技术[6,7],其中性诱灭雄技术[指施用甲基丁香酚（methyl eugenol,简称ME）和诱蝇酮（cuelure,简称CUE）等对相关种类实蝇雄性具有强吸引作用的副信息素进行诱杀而降低雄虫数量]从实蝇防治初期就一直发挥着重要作用[8]。这2种物质对很多种类实蝇雄性都具有很强的引诱活性,但世界上没有一类昆虫像实蝇科昆虫那样对副信息素类物质的反应存在如此大差异[9,10]。就南亚实蝇而言,对ME和CUE都有趋性,但ME/CUE对其引诱活性远不如ME对橘小实蝇、CUE对瓜实蝇那样强烈,这就导致了性诱灭雄技术在南亚实蝇上效果非常不理想,近年来,随着高通量测序技术的发展,转录组分析逐渐成为昆虫中寻求新的防治措施的有效手段[11]。本研究以HiSeq2500高通量测序技术,首次对南亚实蝇进行了转录组测序研究,并以该转录本为依据对其嗅觉基因进行了挖掘,对OBPs进行了进化分析,为南亚实蝇嗅觉机理的探讨和防治方法的寻求提供了重要数据。
1 材料与方法
1.1 材料
1.1.2 测序样本 正常饲养、辐照处理、CUE/ME诱导的南亚实蝇在1日龄（对ME有触角电位反应,但没趋性[6]）、5日龄（对ME/CUE开始有趋性反应,但性未成熟）、10日龄（性成熟、交配前后）、17日龄（产卵后）、19日龄（部分实蝇开始出现2次交配）雌雄各取5只,平行取3管,液氮速冻后送至北京百迈客生物科技有限公司进行测序和分析,测序系统为HiSeq2500PE125测序仪。
1.2 方法
1.2.1 RNA的提取 样品充分研磨后用Trizol Reagent提取RNA,用Nanodrop分光光度计（德国Implen公司）、Qubit 2.0 荧光定量仪（美国Life Technologies公司）、Agilent 2100系统（美国Agilent Technologies公司）方法检测RNA样品的纯度、浓度和完整性。
1.2.2 文库的构建、测序及数据的拼接、组装 RNA样品检测合格后,按照Preparing Samples for Sequencing of mRNA试剂盒说明书（Illumina）构建测序文库。当库检合格后,用Illumina HiSeq 2500 的PE125 bp双向测序技术对文库进行测序。将测序得到碱基质量打分能达到或超过Q30高质量的Raw data进行数据过滤,去除其中接头序列及低质量Reads,获得高质量的Clean data。然后利用短reads组装软件Trinity将具有重叠区域的Clean data进行序列组装,获得Unigene库。
1.2.3 Unigene功能注释 使用BLAST软件将Unigene序列与NR（Non-redundant protein database,非冗余的蛋白质序列数据库）、Swiss-Prot（A manually annotated and reviewed protein sequence database,人工注释和审查蛋白质序列数据库）、GO（Gene ontology,基因本体）、COG（Cluster of orthologous groups,蛋白相邻类的聚簇）、KOG（Eukaryotic orthologous groups,真核生物同源组）、Pfam（Pfam database,同源蛋白家族数据库）、KEGG（Kyoto Encyclopedia of Genes and Genomes,京都基因与基因组百科全书）数据库比对,通过基因的相似性进行功能注释。
2 结果与分析
2.1 RNA的提取与检测
总RNA提取完整,OD260/280为2.09,OD260/230为0.91,28S/18S为1.4,浓度分别为2873.4 ng/μL,RNA纯度和质量符合后续实验要求。
2.2 转录组测序及序列拼接组装
对测序数据进行产出统计,共得到10.97 Gb Clean Data,Q30比例（测序错误率<0.1%）为90.05%（>90%）,GC含量为43.58%。说明此样品测序的数据符合要求,可以保证后续序列拼接组装质量。
采用Trinity软件把获得的高质量测序数据进行序列组装,取每条基因中最长的转录本作为Unigene,得到36 109条Unigenes（表1）。其中N50和平均长度分别为1897和955.14 bp,说明此转录本文库的测序和组装结果较好,能够进行后续生物信息学分析。
表1 组装结果统计表Tab. 1 The statistical table of assemble result |
| 长度范围 Length range/bp | 转绿本 Transcript | 独立基因 Unigene | ||
|---|---|---|---|---|
| 数量 Number | 百分比 Percentage/% | 数量 Number | 百分比 Percentage/% | |
| 200~300 | 14 306 | 23.28 | 11 402 | 31.58 |
| 300~500 | 11 761 | 19.14 | 8 271 | 22.91 |
| 500~1000 | 11 941 | 19.43 | 6 562 | 18.17 |
| 1000~2000 | 11 036 | 17.96 | 5 201 | 14.40 |
| >2000 | 12 412 | 20.20 | 4 673 | 12.94 |
| 总数 Total number | 61 456 | 36 109 | ||
| 总长度 Total length | 80 188 825 | 34 489 176 | ||
| N50长度 N50 length/bp | 2 523 | 1 897 | ||
| 平均长度 Mean length/bp | 1 304.82 | 955.14 | ||
Note: N50, length of the smallest contig which when added to a set of larger contigs yields at least 50% of the genome. | |
注：N50为覆盖50%所有核苷酸的最大Unigene长度。 |
2.3 Unigenes的功能注释
使用BLAST软件将36 109条Unigenes与NR、Swiss-Prot、GO、COG、KOG、Pfam、KEGG数据库比对,获得Unigene的注释信息。最终获得21 127条有注释信息的Unigenes,占总Unigenes的58.51%,其中在NR数据库中注释成功的Unigenes数目最多,20 948条,占总数的58.01%;其后依次是Swiss-Prot有14 103条（39.06%）,GO有14 029条（38.85%）,KOG有13 479条（37.33%）,Pfam有12 670条（35.09%）,KEGG有6442条（17.84%）和COG有5984条（16.57%）（表2）。Unigene注释到NR数据库中的物种分布图显示（图1）,地中海实蝇C. capitata最多为13 526条,其后依次为黑腹果蝇D. melanogaster（3112条）、拟果蝇D. simulans（302条）、果蝇D. sechellia（196条）、中华支睾吸虫Clonorchis sinensis（192条）、家蚕Bombyx mori（176条）、黑果蝇D. virilis（173条）、非洲果蝇D. yakuba（159条）、罗阿丝虫Loa loa（154条）、南美热带果蝇D. willistoni（153条）;其他物种2718条。
表2 独立基因注释统计表Tab. 2 The statistical table of unigene annotation |
| 注释数据库 Annotated databases | 独立基因Unigenes | 注释百分比 Percentage/% |
|---|---|---|
| COG | 5 984 | 16.57 |
| GO | 14 029 | 38.85 |
| KEGG | 6 442 | 17.84 |
| KOG | 13 479 | 37.33 |
| Pfam | 12 670 | 35.09 |
| Swiss-Prot | 14 103 | 39.06 |
| NR | 20 948 | 58.01 |
| 至少一个数据库注释成功 Annotated in at least one database | 21 127 | 58.51 |
| 总独立基因数Total unigenes | 36 109 |
2.4 Unigenes的功能分类
对南亚实蝇转录组的Unigenes进行GO分析发现,在GO注释的14 029条Unigenes中,注释到生物学过程（Biological process）的基因最多（52 151条）,其次是分子功能（Molecular functiion）（25 047条）,细胞组分（Cellular component）的最少（18 959条）。在57个亚类中,以细胞过程、单生物过程、代谢过程在生物学过程中的比例较高,细胞和细胞部分在细胞组分中所占比例较高,蛋白结合和催化活性在分子功能分类中占有较高比例（图2）。对南亚实蝇的Unigenes基于GO数据库进行功能分类,可以从宏观上认识南亚实蝇基因功能分布特征,为下一步目标基因的挖掘提供了便利。
将组装得到的Unigenes与KOG数据库比对,共有13 479个Unigenes基因被注释,占37.33%,分为25个类别,最多的是一般功能预测基因（1681个）,占12.47%;其次,是翻译、核糖体结构与合成基因（711个）,占5.27%;第三,是复制、重组和修复基因（629个）,占4.67%;第四,是翻译后的修饰、蛋白质转化、伴侣基因（584个）,占4.33%;最少的是核酸结构基因（2个）,占0.01%,细胞外结构基因（0个）（图3）,说明了南亚实蝇Unigenes功能分类的总体情况。
KEGG数据库对南亚实蝇的Unigenes可能参与或涉及的代谢途径进行分析表明,6442个Unigenes参与了细胞过程（cellular processing）、环境信息处理（environmental information processing）、遗传信息处理（genetic information processing）、代谢（metabolism）和生物系统（organismal systems）5大类共152个代谢通路,其中注释基因数最多的10个代谢通路依次为：核糖体（276个）、氧化磷酸化（268个）、内质网蛋白加工（246个）、RNA运输（233个）、嘌呤代谢（181个）、剪接体（174个）、细胞内吞作用（146个）、泛素介导的蛋白质降解（143个）、溶酶体（141个）和吞噬体（133个）（表3）。
表3 注释在京都基因与基因组百科全书上独立基因最多的10个代谢通路Tab. 3 Ten unigene pathway hierarchy with the most annonated unigenes of Bactrocera tau on KEGG |
| ID | 代谢通路 Pathway hierarchy | 基因数 Number of genes |
|---|---|---|
| ko03010 | 核糖体 Ribosome | 276 |
| ko00190 | 氧化磷酸化 Oxidative phosphorylation | 268 |
| ko04141 | 内质网蛋白加工 Protein processing in endoplasmic reticulum | 246 |
| ko03013 | RNA运输 RNA transport | 233 |
| ko00230 | 嘌呤代谢 Purine metabolism | 181 |
| ko03040 | 剪接体 Spliceosome | 174 |
| ko04144 | 细胞内吞作用 Endocytosis | 146 |
| ko04120 | 泛素介导的蛋白质降解 Ubiquitin mediated proteolysis | 143 |
| ko04142 | 溶酶体 Lysosome | 141 |
| ko04145 | 吞噬体 Phagosome | 133 |
2.5 嗅觉相关基因分析
通过对7大公共数据库同源比对和注释分析,528个嗅觉相关基因在南亚实蝇转录组中获得鉴定。其中包括34个OBPs、3个化学感受蛋白基因[Chemosensory protein (CSP) genes]、20个嗅觉受体基因[Olfactory receptor (Or) genes]、43个离子型受体基因[Ionotropic receptor (IR) genes]、204个昆虫触角G蛋白基因（GTP binding protein genes）、5个昆虫感觉神经元膜蛋白[Sensory neuron membrane proteins (SNMPs)]基因和219个气味降解酶基因[Odor degrading enzymes (ODE) genes]（表4）。对OBPs的分析发现,28条具有完整的开放阅读框,进一步氨基酸序列分析表明,2条属于Minus-OBPs,1条为Dimer-OBPs,1条为“Plus-C”OBP,24条为典型OBPs;对30条结构较完整的南亚实蝇OBPs序列与黑腹果蝇D. melanogaster（58条）、瓜实蝇B. cucuribitae（34条）、地中海实蝇C. capitata（28条）、橘小实蝇B. dorsalis（47条）OBPs序列进行进化分析表明（图4）,B. tau OBP多数至少与一种实蝇OBP同源物聚类在一支,且拥有80%以上的序列一致性,仅有极个别B. tau OBP与同一支上的实蝇OBP序列一致性较低,如C29372.graph_c0与BdorOBP（AKM45842.1）序列一致性仅为44%且在进化树中形成一个独立的小分支,表明不同种类的实蝇OBPs出现了各自的分化。
图4 MEGA5.0构建的南亚实蝇与黑腹果蝇、地中海实蝇、橘小实蝇、瓜实蝇OBPs的系统发生关系○：南亚实蝇候选气味结合蛋白;▲：DemlOBPs,黑腹果蝇气味结合蛋白;♦:CcapOBPS,地中海实蝇气味结合蛋白;：BdorOBPs,橘小实蝇气味结合蛋白;：BcucOBPs,瓜实蝇气味结合蛋白。 Fig. 4 Neighbor-joining tree of odorant binding protein (OBP) sequenees from B. Tau, D.melanogaster, C. capitata, B. dorsalis and B. cucuribitae, created using MEGA5.0 ○: Candidate odorant binding protein of B. tau; ▲: Odorant binding protein of D.melanogaster; ♦: Odorant binding protein of C. capitata; : Odorant binding protein of B. dorsalis; : Odorant binding protein of B. cucuribitae. |
表4 南亚实蝇转录组中与嗅觉相关基因统计Tab. 4 Statistics of unique transcripts associated with olfactory of Bactrocera tau |
| 嗅觉相关基因 Olfaction-related genes | 被注释基因数 Number of annotated unigenes | NR 注释数 NR annotation number |
|---|---|---|
| 气味结合蛋白基因 Odorant binding protein (OBP) genes | 33 | 32 |
| 化学感受蛋白基因 Chemosensory protein (CSP) genes | 3 | 3 |
| 嗅觉受体基因 Olfactory receptor (Or) genes | 20 | 17 |
| 离子型受体基因 Ionotropic receptor (IR) genes | 43 | 41 |
| 昆虫触角G蛋白基因 GTP binding protein genes | 204 | 121 |
| 昆虫感觉神经元膜蛋白 Sensory neuron membrane protein (SNMPs) genes | 5 | 3 |
| 气味降解酶基因 Odor degrading enzymes (ODE) genes | 219 | 217 |
3 讨论
副信息素类引诱剂从发现其对实蝇相关种类雄性具有强引诱活性开始,就一直在实蝇综合防治和监控上发挥着重要作用[8],在很多地区甚至是唯一的防治措施。但在施用过程中发现：相关实蝇种类的雄虫取食少量的引诱剂不但不会死亡还可大幅度提高其交配竞争能力以及对雌性的吸引力[6, 16];实蝇经过一段时间的ME处理后不但有助于性成熟[17],且其当代和后代雄性对ME的敏感性下降,有些学者据此担心副信息素类诱剂不合理的使用会导致田间出现其抗性种群[8];实蝇辐照成为不育雄蝇后,对引诱剂的反应活性降低且开始反应的日龄推迟,饲喂一定量的引诱剂后交配能力得到一定程度的恢复[18],据此副信息素类引诱剂被广泛应用到实蝇SIT项目中[19]。为探讨副信息素类引诱剂使用对南亚实蝇嗅觉相关基因的影响,进而指导其在实蝇监测、防治中的合理应用,本研究综合考虑转录组测序的特点[20]和南亚实蝇当时的发育历期[21,22],对正常饲养、辐照处理、副信息素类引诱剂诱导下的虫源在不同日龄进行了混合取样和测序。
该样品通过组装共获得36 109条Unigenes,其中N50长度为1897 bp,获得了测序和组装结果较好,能满足后续对功能基因挖掘要求的转录本文库。进一步将得到的Unigenes通过Blast与7大公共数据库进行同源性搜索获得21 127条Unigenes（58.51%）的注释,其中在NR注释成功的Unigenes数目最多,20 948条（58.01%）,且在NR数据库中,注释到地中海实蝇C. capitata上的Unigenes最多,为12 526条,占在NR数据库中总注释Unigenes的58.80%;注释到黑腹果蝇D. melanogaster上的Unigenes为3112条,占在NR数据库中总注释Unigenes的14.86%,这不但说明了地中海实蝇与南亚实蝇的亲缘关系较近,也进一步说明了此转录本文库组装和注释的正确,还有,这2种昆虫都有基因组数据,这为下一步南亚实蝇功能基因的挖掘和生物信息学分析提供了便利。同时,利用GO、KOG数据库对南亚实蝇Unigenes进行了基因功能分类和功能预测,用KEGG数据库对南亚实蝇6442条Unigenes进行了代谢通路初步分析,这为下一步挖掘功能基因及开展靶标基因的克隆及功能验证提供了基础数据。
嗅觉在昆虫寄主定位与选择、躲避天敌、寻找配偶、产卵行为等生命活动中发挥着感知寄主挥发物、昆虫利它素、昆虫信息素等气味物质的重要作用,昆虫在长期进化过程中形成了一套高度专一、极其灵敏、复杂的嗅觉系统,系统中涉及OBP、CSP、OR、SNMP、IR、G-Protein和ODE等多种蛋白家族[23]。其中OBP的功能是将外界疏水性的气味分子运载到嗅觉神经树突膜上的感觉受体（OR）,是昆虫实现嗅觉信号转导的第一步,也是目前研究较多、机理了解最多的一步。在该转录本中,我们初步筛选出了528个嗅觉相关基因,其中包括34个OBPs、3个CSPs、20个Ors、43个IRs、204个G蛋白基因、5个SNMPs和219个ODEs。204个G蛋白基因中包括13个小G蛋白基因;ODEs包括细胞色素P450基因150个,谷胱甘肽S-转移酶基因（Glutathione S-transferase）22个,醛脱氢酶基因（Aldehyde dehydrogenase,ALDHs）25个、醇脱氢酶基因（Alcohol dehydrogenase,ADHs）22个。34个OBPs中有28条具有完整的开放阅读框,氨基酸序列大多符合Cys-X20-66-Cys-X3-Cys-X21-43-Cys- X8-14-Cys-X8-Cys、Cys-X15-39-Cys-X3-Cys-X21-44Cys-X7-12-Cys-X8-Cys或X22-68-Cys-X25-68-Cys-X3- Cys-X31-46-Cys-X8-29-Cys-X8-9-Cys-X5-71典型OBPs的结构特点[24,25],与黑腹果蝇D. melanogaster、瓜实蝇B. cucuribitae、地中海实蝇C. capitata、橘小实蝇B. dorsalis在NCBI上已登录的所有OBPs中至少有一条同源物聚类在同一支,且拥有80%以上的序列一致性,这为下一步借鉴这些实蝇在基因组方面的研究成果进一步挖掘研究南亚实蝇OBPs及其他功能基因提供了依据,也为下一步以该蛋白为基础,探讨南亚实蝇对副信息素的引诱机理并筛选新的引诱成分提供了基础数据。