叉角厉蝽唾液腺转录组及差异分析*

2021-01-22廖贤斌李丽芳兰明先吴国星

云南农业大学学报(自然科学) 2021年1期

高平，廖贤斌，李丽芳，兰明先，赵航，陈斌，吴国星，高熹

(云南农业大学植物保护学院，云南昆明 650201)

可持续农业的发展，更加重视害虫的综合治理(integrated pest management)，利用天敌对作物害虫进行生物控制是其中重要的环节。常见的天敌有捕食性昆虫、寄生性昆虫和昆虫病原体。在作物生态系统中，捕食性昆虫和寄生蜂是大多数害虫的天然生物控制因子[1]。叉角厉蝽[Eocanthecona furcellata(Wolff)]属于半翅目(Hemiptera)蝽科(Pentatomidae) 益蝽亚科(Asopinae) 昆虫，其若虫和成虫的捕食量大，捕食的范围很广，尤其对鳞翅目的害虫具有较强的捕食能力，是一种重要的捕食性天敌[2-5]。叉角厉蝽易在实验室饲养形成种群，它的捕食潜能在东南亚、日本、印度以及中国台湾地区等都有报道[6-7]。

近年来，高通量测序技术已应用于昆虫唾液的研究与鉴定。通过高通量测序技术可以获得唾液腺的转录组，从基因层面研究唾液腺的特征与功能，并从转录水平预测唾液中的潜在分泌蛋白。昆虫的唾液主要由唾液腺分泌，通过分析唾液腺的表达基因来预测唾液蛋白成分是研究昆虫唾液的有效手段之一。目前，有报道采用转录组和蛋白组等手段分析烟粉虱、褐飞虱、灰飞虱和白背飞虱等唾液腺和唾液蛋白组分[8-11]。在这些昆虫唾液腺中发现很多与蛋白合成、食物消化和解毒等相关的基因。然而，在捕食性昆虫中，有学者对小花蝽(Orius laevigatus)唾液腺中差异基因表达谱和猎蝽科昆虫有毒唾液的新型肽和唾液中的毒素进行了研究[12-14]。

目前，国内对该蝽的研究主要集中在食性、生物学特性、人工饲料和捕食功能研究方面[15-17]，国外对该蝽的研究主要集中在生物防治、生活史和捕食功能等方面[18]。近年来，针对与叉角厉蝽同属于一个亚科的躅蝽(Arma chinensis)研究报道较多，王燕等[19]研究了躅蝽对草地贪夜蛾(Spodoptera frugiperda)的生防潜力，苗少明[20]对躅蝽进行了基因组测序，并对躅蝽和茶翅蝽(Halyomorpha halys)的唾液腺进行了转录组测序。益蝽亚科是蝽科中较为特化的类群，是由植食性种类分化形成的捕食性种类，虽然叉角厉蝽和躅蝽同属于益蝽亚科，但是近几年来国内尚未有关于叉角厉蝽唾液腺的研究报道。

捕食性蝽类将刺吸式口器插入猎物体内并分泌唾液，导致猎物快速麻痹和死亡[21]，之后再吸食猎物[22-23]。引起猎物麻痹和死亡的原因是唾液腺产生的分泌物释放到猎物内部产生的作用[24]。据报道组织蛋白酶B (Cathepsin B)和神经内分泌蛋白7b2 (Neuroendocrine protein 7b2)基因在小花蝽的捕食中起重要作用[14]。然而，叉角厉蝽唾液分泌物中引起猎物麻痹和死亡的唾液化合物仍然是未知的。因此，对唾液腺中基因表达的情况进行深入的研究，将有助于阐明唾液毒素分泌的分子机制。本研究的目的是对叉角厉蝽的唾液腺基因进行系统的研究，并确定与唾液腺毒素分泌有关的潜在基因，探究唾液腺毒素的分泌机制。期待发现对昆虫取食与猎物相互作用必不可少的唾液腺相关基因，为后续的基因调控研究和生物源农药开发打下基础。

本研究利用Illumina 测序技术对叉角厉蝽的三龄若虫和成虫唾液腺进行高通量转录组测序，将获得的序列信息进行拼接和组装，进行基因功能注释和功能分类，并对叉角厉蝽三龄若虫和成虫唾液腺的基因差异表达进行分析，为进一步挖掘基因、探究基因功能和分析代谢途径等奠定基础。

1 材料与方法

1.1 昆虫采集与饲养

叉角厉蝽采集于云南省元江县。带回室内采用笼罩式饲养[(27±1) ℃，RH 60%～80%，光周期14L∶10D)，每日饲喂黄粉虫(Tenebrio molitor)，建立实验室种群。

1.2 总RNA 提取与质检

将叉角厉蝽三龄若虫和成虫(雄性∶雌性=1∶1)置于−4 ℃冰箱中10 min，降低虫体活动能力。用剪刀把叉角厉蝽所有触角、足、口器和翅都剪除后，立即放入置于解剖镜下的培养皿，在培养皿中加入PBS 缓冲液，以没过虫体为宜。使用解剖镊和解剖针在显微镜下解剖获得叉角厉蝽唾液腺，将其移至盛有Trizol 的1.5 mL 离心管中，研磨并置于−80 ℃的冰箱中保存。三龄若虫和成虫各设3 次重复，每个重复解剖80 个唾液腺。采用Trizol Reagent 的方法提取得到上述样品的总RNA，RNA 提取后分别用琼脂糖凝胶电泳和NanoPhotometer spectrophotometer 检测RNA的纯度，检测合格之后寄送至北京诺禾致源生物信息科技有限公司进行转录组文库构建和测序。

1.3 cDNA 文库构建与测序

取1.5 μg 质检合格的RNA 样品进行文库构建，然后使用Agilent 2100 bioanalyzer 对文库的insert size 进行检测。库检合格后，把不同文库按照有效浓度及目标下机数据量的需求pooling 后进行Illumina 测序。

1.4 转录组拼接和功能基因注释

测序得到的原始数据要去除序列接头、ploy-N 和低质量序列，才能获得高质量的Clean Data，同时计算Q20、Q30 和GC 含量。获得高质量的测序数据再利用Trinity 软件进行组装[25]。将组装得到的unigenes 使用BLAST 软件在七大数据库中进行比对，获得功能基因的注释信息。

1.5 关键基因鉴定

为了鉴定叉角厉蝽唾液腺中引起猎物死亡的关键基因，在叉角厉蝽唾液腺NR 注释数据库里使用基因全称搜索BAEK 等[14]报道的组织蛋白酶B 和神经内分泌蛋白7b2 基因，并在转录组组装(transcriptome assembly)数据库里查找编码这些基因的序列。

2 结果与分析

2.1 转录组序列分析和组装

利用Illumina 测序平台对叉角厉蝽唾液腺转录组进行测序，每个样本都生成超过7.80 GB 的高质量数据，样本的GC 含量不低于43.71%，测得数据Q20 均大于97.00%，碱基Q30 百分比不少于92.00%，表明测序结果高度准确，可以用于后期的拼接组装(表1)。

获得的clean reads 进行组装(表2)，共得到51 845 个转录本(transcripts)，序列信息达到81 004 577 bp，平均长度为1 562 bp，N50 长度为2 615 bp。在转录本基础上进一步组装获得26 022条功能基因(unigenes)，长度为35 520 646 bp，N50 长度为2 298 bp；其中，长度超过1 kb 的unigenes 有10 634 条，占40.87%。

2.2 基因功能注释

unigenes 在七大数据库中比对得到功能信息(表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 数据库中注释的最少，为5 084 条，所占比例为19.53%。

表1 叉角厉蝽转录组测序产量统计Tab.1 Output statistics of E.furcellata transcriptome sequencing

表2 叉角厉蝽转录组组装质量Tab.2 Quality for de novo E.furcellata transcriptome assembly

E值等于0 说明该基因的注释可信。结果显示：所有匹配结果的E值均大于1e-5，E值大于1e-5 的占86.3% (图1a)。从图1b 可见：Nr 注释到的unigenes 基因中，相似性在95%以上的基因占15.5%。从匹配的物种来源分析(图1c)，叉角厉蝽唾液腺转录组测序拼接的unigenes 与茶翅蝽(Halyomorpha halys)的unigenes 相似性最高，达到76.0%。其次与温带臭虫(Cimex lectularius)、Cryptotermes secundus、褐飞虱(Nilaparvata lugens)和Lasius niger等的相似性分别为3.7%、2.7%、1.6%和1.0%；另外与其他物种同源的基因占15.0%。

表3 基因注释成功率统计表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 基因功能分类

叉角厉蝽唾液腺转录组有23 132 条(47.78%)unigenes 被GO 注释至生物学过程(biological process)的26 个功能亚类；注释至分子功能(molecular function)的unigenes 有10 660 条(22.02%)；14 621 条(30.20%) unigenes 归属于细胞组分(cellular component)的10 个功能亚类(图2)。在生物学过程中，注释至细胞过程(cellular process，GO：0009987)的unigenes 数量最多，达到4 859 条(21.00%)；注释到代谢过程(metabolic process，GO：0044237)的unigenes 数量较多，达到4 496条(19.44%)；注释到单生物过程(single-organism process，GO ：0044699) 的unigenes 达到3 791 条(16.39%)。其余均在3 000 条以下，其中注释至行为(behavior，GO：0007610)、节律性过程(rhythmic process，GO：0048511) 和解毒过程(detoxification，GO：0098754) 的unigenes 仅有5、5 和6 条。分子功能类型中，注释至蛋白结合活性(binding，GO：0005488) 和催化活性(catalytic activity，GO：0003824) 的unigenes 分别有5 096 条(47.80%) 和3 790 条(35.55%)；注释至转运活性(transporter activity，GO：0005215)的unigenes 数量达到677 条。其余功能亚类均在300 条以下，其中被注释至金属伴侣蛋白活性(metallochaperone activity，GO：0016530) 的unigenes 数量仅有4 条序列。细胞组分功能类型中，注释至细胞(cell，GO：0005623) 和细胞部分(cell part，GO：0044464) 的unigenes 都达到2 793 条(19.10%)；注释至细胞器(organelle，GO：0043226)、大分子复合物(macromolecular complex，GO：0032991)、细胞膜(membrane，GO：0016020)、细胞膜部分(membrane part，GO：004 4425)和细胞器部分(organelle part，GO：00444 22)的unigenes 数量也分别达到1 984 条(13.57%)、1 702 条(11.64%)、1 630 条(11.15%)、1 531 条(10.47%)和1 045 条(7.15%)；其余被注释的功能亚类均在1 000 条以下。

所有unigenes 经过KOG 数据库功能预测和分类，叉角厉蝽唾液腺转录组共有5 423 个unigenes 得到注释，并被分为26 个KOG 类(图3)。其中，共有4 个KOG 亚类包含450 条以上unigenes，依次是：一般功能预测(general function prediction only，R) 913 条unigenes，占比最大(16.84%)；信号转导机制(signal transduction mechanisms，T) 675 条(12.45%) unigenes；蛋白质翻译后修饰、折叠和分子伴侣(posttranslational mofdiication，protein turnover，chaperones) 570 条(10.51%) unigenes；未知功能(function unknown，S) 496 条(9.15%) unigenes。而细胞动力(cell motility，N)、核结构(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 数据库。在获得注释的unigenes 中，25.41%涉及有机系统(E，organismal systems)，其中主要参与内分泌系统(endocrine system，325 条)、消化系统(digestive system，209 条) 和免疫系统(immune system，205 条)，其余亚类均在200 条以下。24.63%涉及新陈代谢(D，metabolism)，主要参与糖类(carbohydrate metabolism，204 条)、脂质(lipid metabolism，182 条) 和蛋白质(amino acid metabolism，157 条)三大营养物质的代谢以及核苷酸代谢(nucleotide metabolism，117 条)等。除此以外，这些unigenes 还广泛参与到信号转导(signal transduction，548 条)通路中(图4)。

2.4 差异表达基因的筛选

通过FPKM 方法评估本次测序叉角厉蝽唾液腺样品中的基因表达情况，采用DESeq 进行若虫和成虫唾液腺样品组间的差异表达分析。在差异表达基因的筛选过程中，将Padj＜0.05 作为差异基因的筛选标准，比较2 组测序结果，筛选出大量差异表达的基因。共鉴定出3 437 个差异表达基因(图5)。

2.5 差异表达基因GO 富集分析

GO 富集柱状图能表现出差异表达基因在GO term 上富集的分布情况(图6)。将筛选到的3 437个差异表达基因(DEGs)分类到3 个ontology (本体)中，差异表达基因在分子功能和生物学过程中分别包含了15 和12 个功能小类。在生物过程这一类中，差异表达基因在蛋白水解(proteolysis)和DNA 复制(DNA replication)2 个功能小类中所占比例最高；在分子功能这一类中，差异表达基因在水解酶活性(hydrodase activity)和肽酶活性(peptidase activity) 2 个功能小类中占比最高。

2.6 差异表达基因的KEGG 注释

KEGG 代谢通路分析共得到1 509 个差异表达基因，这些差异表达基因参与261 个不同的代谢通路。表4 是含有最多差异基因的前20 个通路，在所有通路中代谢途径是含有差异基因序列最多的。KEGG 注释信息主要富集在嘌呤代谢(purine metabolism)、嘧啶代谢(pyrimidine metabolism)、半胱氨酸和蛋氨酸的代谢(cysteine and methionine metabolism)以及组氨酸代谢(histidine metabolism)等。在KEGG 代谢通路中对与分泌相关的差异表达基因数量进行统计(表5)，共得到58 个差异表达基因。

2.7 关键基因的鉴定

在叉角厉蝽唾液腺转录组数据库中鉴定了3 个编码组织蛋白酶B 的基因。组织蛋白酶B 是溶酶体半胱氨酸蛋白酶，由二硫键连接的重链和轻链的二聚体组成。同时，还鉴定了3 个编码神经内分泌蛋白7b2 的基因。神经内分泌蛋白7b2 在大多数神经元和内分泌细胞中起着至关重要的作用。

图4 叉角厉蝽转录组unigenes 的KEGG 通路分析Fig.4 Classification of E.furcellata transcriptome based on KEGG

图5 差异基因火山图Fig.5 Differential gene volcano map

3 讨论

转录组测序以N50 和Q30 来评估转录组数据的质量。一般认为N50 值越大就表明得到的长片段越多[26]，且Q30 值大于85% 说明测序结果可靠[27]。通过对叉角厉蝽成虫和三龄若虫唾液腺进行转录组测序，共获得唾液腺unigenes 序列26 022 条，平均长度1 365 bp，N50 长度为2 298 bp，Q30 值都大于92%。而茶翅蝽和躅蝽的唾液腺转录组测序数据为53 618 条和92 468条，平均长度为736 和729 bp，N50 长度仅为1 439 和1 311 bp[13]。从N50 长度和Q30 值等方面均显示本次叉角厉蝽唾液腺转录组的拼接结果良好。

将全部unigenes 序列与七大数据库进行功能注释，获得14 057 条unigenes (54.02%)成功注释，仍有11 965 条unigenes (45.98%)未被注释，这可能是unigenes 较短，所以未与公共数据库中的序列比对上[28]，也有可能是存在很多功能未知的新基因[29]。通过与NR 数据库进行比对，结果表明叉角厉蝽与茶翅蝽的同源序列最多，茶翅蝽虽是植食性害虫，主要危害蔬菜和水果[30]，但这可能与茶翅蝽和叉角厉蝽同属于蝽科，亲缘关系较近有关。然而本次同源序列比对并没有注释到与躅蝽的相似性，这可能是由于NR 数据库中暂时还没有收录到躅蝽的注释信息。

图6 GO 富集柱状图Fig.6 GO enrichment histogram

表4 差异基因KEGG 富集列表(前20 个)Tab.4 Differential gene KEGG enrichment list (top 20)

表5 与分泌相关差异表达基因的统计Tab.5 Statistics of differentially expressed genes in secretion

GO 功能富集结果显示：叉角厉蝽唾液腺参与生物进程的unigenes 最多，其中前3 个类别分别是细胞进程、代谢进程和单生物进程，由此可以看出唾液腺中的细胞新陈代谢十分活跃。在分子功能中，unigenes 数量最多的3 个类别是结合、催化活性和转运活性。这与白背飞虱(Sogatella furcifera)唾液腺[31]和烟粉虱(Bemisia tabaci)唾液腺[10]的研究结果较为一致。分泌唾液是昆虫唾液腺的主要功能，注释结果显示有很多基因参与结合和转运。在细胞组分中，前3 个是细胞、细胞部分和细胞器，由此可见细胞、细胞器和细胞分子是唾液腺的重要组成部分[32]。叉角厉蝽唾液腺的这些注释结果表明：在唾液的合成和分泌过程中需要活跃的代谢活动和高效的蛋白质合成与转运机制。

通过对叉角厉蝽成虫和三龄若虫唾液腺的差异表达分析，发现叉角厉蝽唾液腺的26 022 个unigenes 中有3 437 个unigenes 在成虫和三龄若虫中的表达水平不同。数据统计发现有1 871 个表达水平上调的unigenes，1 566 个表达水平下调的unigenes。GO 富集分析中，差异基因主要富集在水解酶活性、肽酶活性和蛋白水解，表明三龄若虫唾液腺还处在快速发育阶段。KEGG 代谢通路分析有1 509 个差异表达基因参与261 个不同的代谢或信号通路，同时发现13 个差异表达基因参与了唾液分泌，这些差异基因的功能有待进一步研究。

BAEK 等[14]报道组织蛋白酶B 可能参与了特定的蛋白水解，从而影响猎物的凝血或免疫力，最终导致杀虫作用，同时KUTSUKAKA 等[33]发现在Tuberaphis styraci中组织蛋白酶B 是毒素的主要成分。神经内分泌蛋白在激素原转化酶2(prohormone convertase-2) 的运输和激活中起关键作用，激素原转化酶2 负责处理神经肽和肽激素前体[34]。另外，神经内分泌蛋白7b2 可能调节分泌颗粒的形成和分泌[14]，并可能通过控制激素原转化酶2 的激活和转运参与唾液肽的成熟。组织蛋白酶B 和神经内分泌蛋白7b2 基因可能在唾液的功能中发挥重要的作用。