二化螟幼虫热胁迫响应的转录组分析

朱 宇，刘 洋

(温州科技职业学院，浙江 温州 325006)

温度是影响昆虫生长、发育、繁殖、存活、分布范围和多样性的重要非生物因子[1-2]。昆虫是变温动物，对温度的变化敏感[3]。在高温条件下，昆虫发育速度变得更快，温度升高会使昆虫向高纬度地区扩张迁移[1]。高温会导致昆虫种群动态的改变、物种丰富度的下降，以及昆虫存活率的降低[3-5]。昆虫也可以通过调整自身的生理活动来应对高温环境，包括热激蛋白的表达、丙醛浓度的增高以及过氧化氢酶、谷胱甘肽S-转移酶和超氧化物歧化酶等抗氧化物质的提高[6-8]。随着全球气候变暖，越来越多的研究者开始关注昆虫对高温或热胁迫的适应机制[9]。

二化螟是水稻和茭白等作物上的重要害虫，每年给农业生产造成巨大的经济损失[2]。二化螟对温度的适应性强，其越冬幼虫能在-20 ℃下存活[10]。罗举等[11]发现，每天在40 ℃下4 h的二化螟幼虫，也能发育为成虫。二化螟发育的最适温度在27 ℃左右，发育温度的最低阈值在10～15 ℃，最高阈值在35～40 ℃，在40 ℃温度下孵出的幼虫不能存活[11-12]。在高温胁迫下，二化螟幼虫体内会产生一系列的生理反应来消除氧化应激反应对身体的伤害，包括热激蛋白的产生、抗氧化酶的增多等[13-18]。目前，有关二化螟高温胁迫响应的研究主要集中在单个基因的功能，尚未通过转录组学的方法对二化螟全基因表达的影响进行分析。

随着第二代测序技术的成本降低，转录组高通量测序方法已经被应用到多种生物学试验中[19]。通过高通量测序，我们可以同时测定几万个基因的表达水平，并在全转录组范围内找出不同试验处理下的差异表达基因[20]。前人已经有许多关于昆虫热胁迫响应的比较转录组学研究[21-23]。在这些研究中，亚洲柑橘木虱(Diaphorinacitri)在热胁迫下，它的分子伴侣、蛋白生物合成和抗氧化相关基因被诱导表达[21]。高温也可以诱导桑绢丝野螟(GlyphodespyloalisWalker)热激蛋白、解毒酶和氧化还原酶的表达，但糖酵解、糖异生和脂肪酸的生物合成等代谢活动被抑制[23]。东方粘虫(MythimnaseparataWalker)在30 ℃热胁迫下，其热激蛋白和泛素介导的蛋白酶体的表达会显著上调[22]。本研究通过比较转录组学的方法，分析二化螟幼虫应对热胁迫的差异表达基因及相关代谢通路，从而更深入地了解二化螟幼虫对高温胁迫的适应机理。

1 材料与方法

1.1 供试昆虫与样品准备

二化螟幼虫从浙江省温州市种子种苗基地水稻田采集，并在室内用稻苗饲养法进行饲养，温度设为(27±1)℃，光照黑暗时间比16 h/8 h[24]。试验时选取大小基本一致的二化螟四龄幼虫。将热胁迫组(heat stress，HS)和对照组(Control)二化螟幼虫分别置于38 ℃和26 ℃中2 h，随后，从对照组和热胁迫组中分别取2头二化螟放入离心管中，每组处理做3次生物学重复，加入TRIzol以备提取RNA。

1.2 cDNA文库构建与测序

RNA的提取、RNA质量检测、cDNA文库的构建及上机测序由武汉未来组生物科技有限公司(武汉)完成，测序平台为Illumina NovaSeq 6000。

1.3 转录组拼装与注释

由Illumina平台测序得到的原始图像数据经过碱基识别转化为原始数据(raw data)，使用Trimmomatic软件对原始数据进行过滤得到clean data[25]。然后使用FastQC软件对clean data进行质控[26]。使用HISAT2和StringTie软件，将获得的clean data与二化螟参考基因组(下载地址ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/004/000/445/GCA_004000445.1_ASM400044v1)进行比对拼装得到转录本，并计算出每个基因在不同处理中的FPKM值[27-28]。

1.4 差异基因表达分析及功能富集分析

通过DESeq2软件对热胁迫组和对照组进行差异表达分析[29]。对差异检验的P值作多重假设检验校正，通过控制FDR(false discovery rate)来决定P值的阈值[30]。在本次分析中，差异表达基因的阈值为FDR≤0.05。利用NCBI的BLASTN软件、GO(gene ontology)和KEGG(Kyoto encyclopedia of genes and genomes)数据库对筛选出的差异表达基因进行功能注释[31-33]。利用显著上调的差异表达基因和KEGG富集的结果，应用超几何检验，找出显著富集的KEGG 代谢途径[32]。

2 结果与分析

2.1 二化螟转录组测序与序列组装

我们应用RNA-seq测序技术对两组不同温度下二化螟幼虫的基因表达进行了测定。转录组测序的原始数据已经提交到NCBI，BioProject登录号为PRJNA578414。通过测序，我们共得到34.2 GB 原始数据，425 986 598个clean reads(表1)。6个样品(HS_1、HS_2、HS_3、Control_1、Control_2、Control_3) clean data的Q30值均大于93%(表1)。使用HISAT2软件将6个样品产生的测序序列与二化螟参考基因组进行比对，能够比对到参考基因组上的序列占总体的百分比均超过89%，其中，在参考基因组中有唯一比对位置的测序序列占总体的百分比均高于82%，有多个比对位置的序列占总体的百分比均不超过9%(表1)。

2.2 二化螟在热胁迫下的差异表达基因

通过比较两组不同温度处理后的二化螟转录组，得到387个差异表达基因，其中221个基因为上调表达，166个基因为下调表达(图1)。通过BLASTN对这些差异表达基因进行注释，其中，有11个基因的注释结果为热激蛋白(heat shock protein)，它们在热胁迫下上调表达。这些热激蛋白中，有1个属于热激蛋白90家族，有4个属于热激蛋白70家族，2个属于热激蛋白40家族，还有4个是小分子热激蛋白(图2)。此外，在上调表达的差异表达基因中，有3个基因为细胞色素P450(CYP，cytochrome P450)(图3)。除了这两类蛋白外，在上调表达的差异表达基因中，未找到其他与抗氧化相关的酶。

2.3 差异表达基因的GO注释和KEGG通路分析

利用GO数据库对所获得的差异表达基因进行功能注释，并对差异表达基因进行功能分类统计(图4)。在上调的差异表达基因中，有136个基因比对到生物过程(biological process)数据库，其中，代谢过程(metabolic process)这一条目中的差异基因数量最多，其次是单组织过程(single-organism process)和细胞过程(cellular process)。有90个上调差异表达基因比对到细胞组分(cellular component)数据库中，其中，细胞膜(membrane)中的差异基因数量最多，其次是细胞(cell)和细胞部分(cell part)。有105个上调差异表达基因比对到分子功能(molecular function)数据库中，其中，结合(binding)这一条目中的差异基因数量最多，其次是催化活性(catalytic activity)。

利用KEGG数据库对差异表达基因进行注释，上调的221个差异表达基因分布在69条代谢通路中(图5)。应用超几何检验，找到上调表达基因显著富集的代谢通路4条，分别为内质网内蛋白质加工、寿命调节通路-多物种、植物病原体相互作用和抗原处理及呈现(表2)。在内质网内蛋白质加工代谢通路的9个基因中， 8个为热激蛋白，1个为热激蛋白结合蛋白。

表1 二化螟转录组组装结果统计

Table 1 Summary statistics of assembly results ofChilosuppressalistranscriptome

Q30比率表示质量值大于等于30的碱基所占百分比。

Q30 rate， The percentage of bases with a quality score of 30 or higher.

红点表示上调表达的差异表达基因；绿点表示下调表达的差异表达基因；蓝点表达差异不显著的基因；Fold change，差异倍数；padj，经过FDR校正的P值。Red dot indicated up-regulated DEGs (differentially expressed genes); Green dot indicated down-regulated DEGs; Blue dot indicated non-DEGs; Fold change， The ratio of heat stress group versus control group; padj， FDR-adjusted P value.图1 差异基因火山图Fig.1 Volcano map of DEGs

HS，热胁迫组；Control，对照组；GeneID，二化螟参考基因组注释文件中的基因ID；HSP，热激蛋白；FPKM，每1百万个匹配上的读长中匹配到外显子的每1 000个碱基上的片段个数。HS， Heat stress group; Control， control group; GeneID， Gene ID from reference annotation file of Chilo suppressalis; HSP， Heat shock protein; FPKM， Fragments per kilobase of exon model per million reads mapped.图2 二化螟幼虫在热胁迫下差异表达的热激蛋白Fig.2 Differentially expressed genes encoding for heat shock proteins in Chilo suppressalis larvae under high temperature stress

HS，热胁迫组；Control，对照组；GeneID，二化螟参考基因组注释文件中的基因ID；CYP，细胞色素P450；FPKM，每1百万个匹配上的读长中匹配到外显子的每一千个碱基上的片段个数。HS， Heat stress group; Control， control group; GeneID， Gene ID from reference annotation file of Chilo suppressalis; CYP， Cytochrome P450; FPKM， Fragments per kilobase of exon model per million reads mapped.图3 二化螟幼虫在热胁迫下差异表达的细胞色素P450Fig.3 Differentially expressed genes encoding for CYPs in Chilo suppressalis larvae under high temperature stress

Up，上调表达的差异表达基因；Down，下调表达的差异表达基因。Up， Up-regulated DEGs; Down， Down-regulated DEGs.图4 差异表达基因GO分类Fig.4 GO classification of DEGs

3 讨论

本研究通过高通量测序技术对分别置于38 ℃和26 ℃中2 h的二化螟幼虫进行转录组测序和拼装，并进行了差异表达分析。测序结果显示，6个试验样品的Q30值均大于93%，说明转录组测序的质量较高。将转录组数据与二化螟参考基因组进行比对拼接，比对到参考基因组上的序列占总体的百分比均超过89%，在参考基因组中有唯一比对位置的测序序列占总体的百分比均高于82%，有多个比对位置的序列占总体的百分比均不超过9%，说明转录组拼接的质量较高。通过比较两组不同温度处理后的二化螟转录组，我们得到了387个差异表达基因(221个上调表达、166个下调表达)，部分差异表达基因与二化螟幼虫的解毒、抗氧化应激、免疫以及代谢相关。通过对上调表达的差异表达基因进行代谢通路富集分析，发现内质网内蛋白质加工、寿命调节通路-多物种、植物病原体相互作用和抗原处理及呈现这4条通路在二化螟热胁迫时显著富集，说明这些代谢通路在二化螟热胁迫响应的过程中发挥着重要的作用。

Rich factor表示富集上的基因占注释到的基因的百分比；纵坐标表示富集上的条目；Gene number表示富集上的基因个数；颜色表示的Q value值，越低越显著。Rich factor indicated the percentage of the gene on the enrichment to the annotated gene; the ordinate indicated the term of the enrichment; Gene number indicated the number of genes on the enrichment， and the Q value represented by the color， the lower the more significant.图5 上调差异表达基因前20条富集的KEGG代谢通路Fig.5 The top 20 enriched KEGG pathways of the up-regulated differentially expressed genes

表2 二化螟上调表达差异基因的显著富集代谢通路

Table 2 Summary of significantly enriched pathway terms associated with up-regulated differentially expressed genes

通路Pathway基因数Genenumber参考基因数Reference genenumber校正后的P值Corrected P valueMap登录号Map ID内质网内蛋白质加工Protein processing in endoplasmic reticulum91300.000416map04141寿命调节通路-多物种Longevity regulating pathway-multiple species5580.008634map04213植物病原体相互作用Plant-pathogen interaction3140.008634map04626抗原处理及呈现Antigen processing and presentation3220.025823map04612

在昆虫体内，热激蛋白通常和分子伴侣或辅助蛋白一起协作，通过促进变形蛋白修复或参与蛋白的合成、折叠和组装等生命活动来保持昆虫细胞内环境的稳定[34]。热激蛋白可以在昆虫的正常生长发育过程中表达，也会在高温、寒冷、种群密度过高、缺氧以及农药等胁迫环境下诱导表达[35-39]。根据相对分子质量大小，热激蛋白被划分为许多家族，包括小分子热激蛋白、热激蛋白40、热激蛋白60、热激蛋白70和热激蛋白90[34]。热激蛋白能在高温胁迫下诱导表达，已经被许多文献报道。在热胁迫下，葱蝇(Deliaantiqua)非滞育蛹的热激蛋白90表达量增高[40]。肉蝇(Sarcophagacrassipalpis)非滞育蛹中的热激蛋白70在高温胁迫下上调表达[41]。在麦红吸浆虫(Sitodiplosismosellana)中也有类似结果，麦红吸浆虫滞育幼虫的热激蛋白70和热激蛋白90分别在35 ℃和30 ℃表达量最高[42]。本研究中，1个二化螟热激蛋白90和4个热激蛋白70在高温胁迫下上调表达。小分子热激蛋白是一种十分普遍的蛋白，它能够阻止蛋白聚集，帮助细胞适应外界不良环境[43]。本研究中，4个二化螟小分子热激蛋白在高温条件下大量表达。上述结果说明，二化螟热激蛋白70、热激蛋白90和小分子热激蛋白能够提高二化螟幼虫对高温环境的适应。

细胞色素P450是一类包含血红蛋白的家族蛋白[44]。昆虫的细胞色素P450基因家族有4个进化枝，分别为CYP2家族、CYP3家族、CYP4家族和线粒体P450家族[45]。昆虫细胞色素P450基因能够参与昆虫的多种生理活动，包括昆虫的激素代谢、植物化感物质的解毒、昆虫的抗药性以及昆虫的抗氧化应激反应等[46-48]。细胞色素P450参与昆虫的抗氧化应激作用，已经被许多文献报道。在埃及伊蚊(Aedesaegypti)体内发现了4个细胞色素P450基因，参与抗氧化应激作用[48]。在中华蜜蜂(Apisceranacerana)体内，也发现1个细胞色素P450基因AccCYP336A1，具有抗氧化应激作用[49]。昆虫在高温环境下，会产生活性氧和其他有毒物质，对昆虫细胞造成损伤[50]。在本研究中，有3个细胞色素P450基因在38 ℃处理2 h后显著上调表达，说明细胞色素P450可能在二化螟热胁迫响应和抗氧化应激过程中发挥作用。

本研究利用转录组测序技术对分别置于38 ℃和26 ℃中2 h的二化螟幼虫进行转录组测序和分析，从转录组水平揭示短时高温胁迫对二化螟幼虫的影响。结果表明，二化螟幼虫在高温胁迫过程中，热激蛋白和细胞色素P450这两类与热胁迫响应相关的基因表达量发生了显著变化。这些结果可为二化螟热胁迫响应相关的试验提供参考。