基于高通量测序的阿尔泰蝠蛾(鳞翅目:蝙蝠蛾科)幼虫转录组分析

孙 涛，张示渊，王 岩,陈创夫

(1. 石河子大学医学院；2. 石河子大学第三附属医院(石河子市人民医院)；3. 石河子大学动物科技学院；4. 石河子大学畜牧学博士后流动站，新疆维吾尔自治区石河子市 832000)

转录组测序(RNA-sequencing)技术的广泛运用，可以获得大量的转录本信息，从中发掘出有重要功能的基因，尤其适用于基因信息匮乏的物种，很大程度上促进了非模式生物的功能基因组分析(Nagalakshmietal., 2008)，也为研究鳞翅目昆虫分类、进化、发育、与寄主互作等问题提供了技术支持(杨帆等，2014)。与全基因组相比，转录组测序成本低、周期短，具有很高的性价比。对于目前还没测定全基因组的非模式昆虫来说，先测转录组是深入研究基因发掘、分子标记开发和功能鉴定的最佳选择(张棋麟和袁明龙，2013)。目前，已有甜菜夜蛾Spodopteraexigua(Pascualetal., 2012)、小菜蛾Plutellaxylostella(Youetal., 2013)、沙棘蠹木蛾Eogystiahippophaecolus(Cuietal., 2017)等鳞翅目非模式昆虫开展了高通量测序。

阿尔泰蝠蛾HepialusaltaicolaWang属于鳞翅目Lepidoptera蝙蝠蛾科Hepialidae蝠蛾属Hepialus，是新疆特有昆虫，仅分布于海拔1 000～3 000 m处的新疆阿尔泰山。其幼虫是新疆虫草菌Ophiocordycepsgracilis的寄主昆虫，受真菌感染后形成的虫菌复合体即新疆虫草，具有很高的药用价值和经济意义(索菲娅等，2008)。阿尔泰蝠蛾幼虫栖居于地下5～20 cm的土壤层中，以采食新疆芍药Paeoniasinjiangensis和新疆藜芦Veratrumlobelianum等植物嫩根部分为生，所栖息的土壤微环境的日平均温度在-5～15℃，每年还要经历长达半年的冻土期。另外，通过野外观察和室内饲养发现，幼虫在0～5℃还能维持取食和运动能力，具有很强的低温环境适应性。目前，对阿尔泰蝠蛾的研究主要集中在生物学特性方面，主要包括生殖规律的观察(赵恒等，1998a)、幼虫的饲养(赵恒等，1998b)、活蛹的无损鉴别(王岩等，2018)和食性的调查研究(王岩等，2020)等，而关于蝠蛾基因方面的研究相对较少。温度在决定昆虫的生存、分布和数量方面起着至关重要的作用，同时也影响昆虫生理生化过程、代谢活动以及生长发育状况(Sinclairetal., 2003)。阿尔泰蝠蛾幼虫独特的生境(低温高海拔地区)及低温下仍能维持取食和运动的能力，为研究昆虫低温适应性提供了良好的实验材料。

本研究采用Illumina HiSeqTM2500测序平台对阿尔泰蝠蛾幼虫进行转录组测序，并对获得的Unigenes进行基因功能注释、分类及代谢途径的分析，以期获得阿尔泰蝠蛾更多的转录本信息，从基因组水平上发掘蝠蛾应激代谢相关的功能基因，为研究蝠蛾低温适应相关分子和代谢调节机制提供依据，也为进一步相关基因的克隆与表达奠定基础。

1 材料与方法

1.1 供试昆虫

阿尔泰蝠蛾幼虫采自新疆阿勒泰地区布尔津县的阿尔泰山上(海拔1 260 m)，收集幼虫栖身土壤作为饲养土带回实验室。选择个体大小相近的健康活虫用于转录组研究，饲养于光照培养箱中(MGC-250P型，上海一恒)，饲养温度为12℃±1℃，相对湿度50%～60%，光周期16 L ∶8 D，饲养一周后用液氮速冻并保存于-80℃冰箱备用。

1.2 RNA提取、cDNA文库构建和Illumina测序

从冻存的蝠蛾幼虫中选取大小相近的3头，用TRIzol试剂提取RNA。用1%琼脂糖凝胶电泳分析总RNA降解程度及是否有污染，Nanodrop 2000分光光度计(IMPLEN, CA, USA)检测RNA的纯度(OD260/OD280)，并用Qubit RNA Assay Kit对RNA浓度进行定量，Agilent Bioanalyzer 2100(Agilent Technologies, CA, USA)检测RNA的完整性。

将检测合格的总RNA样品用带有Oligo(dT)的磁珠富集mRNA，加入fragmentation buffer使其随机打成短片段，再以片段化的mRNA为模板，用6碱基随机引物(random hexamers)合成cDNA第一链，然后加入缓冲液、dNTPs、RNaseH和DNA聚合酶I合成cDNA第二链，并使用AMPure XP beads对其进行纯化。洗脱纯化后的双链cDNA先进行末端修复、加A尾并连接测序接头，再用AMPure XP beads进行片段大小选择，最后通过PCR扩增得到cDNA文库。采用Illumina HiSeqTM2500(Illumina, San Diego, CA, USA)系统进行测序，测序读长为PE125。cDNA文库构建及测序工作由北京百迈克生物科技有限公司协助完成。

1.3 转录组测序数据组装和拼接

对测序得到的原始数据(raw reads)进行过滤，去除接头污染，去除未知碱基N大于10%的reads，去除低质量reads(质量值Q phred <=20的碱基数占整个reads的50%以上的reads)，得到过滤后的数据(clean reads)，再进行后续分析。在高质量的clean reads基础上计算Q20(Phred>20碱基占总碱基的百分比)、Q30(Phred>30碱基占总体碱基的百分比)、GC(碱基G和C的总和占总的碱基数的百分比)含量和重复序列水平。利用Trinity软件(Grabherretal, 2011)对clean reads进行拼接组装，使用Tgicl对组装结果进行聚类去冗余得到最终的单基因簇Unigenes用于后续分析，命名为“All-unigene”。

1.4 基因功能注释和分类

为获得全面的基因功能信息，使用Blast通过序列相似性将Unigenes与7个主要功能数据库(Nr，KOG/COG，Swiss-Prot，Pfam，KEGG和GO)进行比对，得到基因的功能注释信息。然后用Blast2GO和WEGO进行GO(Gene Ontology)条目和功能分类统计(Conesaetal., 2005; Yeetal., 2006)。利用KOBAS 2.0(Xieetal., 2011)获得KEGG(Kyoto Encyclopedia of Genes and Genomes)中Unigenes的直系同源结果。用HMMER(Eddy, 1998)将Unigene序列与Pfam数据库(Finnetal., 2014)进行比较，获得Unigenes的注释信息。

1.5 CDS预测

根据功能注释结果，将Unigenes按照Nr、Swissprot、KEGG、COG蛋白库的优先级顺序进行比对，挑选Unigene的最佳比对片段作为该Unigene的CDS(按照5′>3′的顺序)；没有比对上的序列或比对上未预测出结果的序列，则采用ESTScan软件预测，得到这部分基因编码的核酸序列和氨基酸序列。

1.6 与环境适应和抗逆性相关基因的推测

根据获得的阿尔泰蝠蛾基因功能注释信息，结合已有相关文献报道(Luetal., 2014; Cuietal., 2017; Torsonetal., 2017)，从中推测出与环境应激相关的代谢调节候选基因。采用FPKM (fragments per kilobase of transcript per million fragments mapped)值表示对应Unigenes的表达丰度。

2 结果与分析

2.1 阿尔泰蝠蛾转录组数据的组装

转录组样本数据量为21.31 G，经过分析，共组装获得144 329个转录本，在此基础上进行聚类和组装分析得到100 133个Unigenes，总长度86 319 112 bp，平均长度862 bp，N50(覆盖50%所有核苷酸的最大序列重叠群长度)长度为1 628 bp，其中长度在1 kb以上的有23 747条，占全部Unigenes的23.72%，长度为200～300 bp的Unigene所占比例最大，占总体的33.6%(表1)，阿尔泰蝠蛾转录组组装后得到Unigenes长度分布(图1)。序列的质量参数Q20的百分比在97%以上；Q30的百分比在92%以上(表2)，可对转录组质量进行评价。

表1 阿尔泰蝠蛾转录组组装结果

表2 阿尔泰蝠蛾cDNA 样品测序数据评估统计表Table 2 Summary for raw reads of cDNA samples of Hepialus altaicola

图1 阿尔泰蝠蛾转录组组装后得到Unigenes长度分布Fig.1 Length distribution of all assembled Unigenes of Hepialus altaicola

2.2 Unigene的功能注释、分类和代谢途径分析

2.2.1Unigene在各数据库中的比对

将获得的100 133个Unigenes分别在7个主要功能数据库中进行比对，取阈值e≤10-10，共有38 198条Unigenes获得注释，占38.15%；其中占比最高的是Nr数据库，注释32 381条Unigenes(32.34%)，与其他生物的已知蛋白具有不同程度的同源性。COG数据库注释基因共有11 084条，占11.03%；KOG数据库注释基因共有20 669条，占20.64%；Swiss-Prot数据库注释基因共有16 690条，占16.67%；Pfam数据库注释基因共有23 816条，占23.78%；GO数据库注释基因共有13 261条，占13.24%；KEGG数据库注释基因共有15 058条，占14.78%；此外，通过其中5种主要数据库(COG、NR、GO、Swiss-Prot和Pfam)维恩图可以更清晰了解各数据库交叉注释情况(图2)。

从所注释到NR数据库中的物种分布来看，有1 634个(占5.36%)基因与斜纹夜蛾Spodopteralitura相似基因序列最多；其次为脐橙螟蛾Amyeloistransitella，有1 503个(占4.93%)；有1 376个(占4.51%)基因与一种名为Pristionchuspacificus的寄生虫线虫同源；其它相似性序列数量大于2%的物种有小菜蛾(4.13%)，家蚕Bombyxmori(3.89%)，偏瞳蔽眼蝶Bicyclusanynana(3.63%)，棉铃虫Helicoverpaarmigera(3.59%)，烟芽夜蛾Heliothisvirescens(3.50%)，柑橘凤蝶Papilioxuthus(2.87%)，其它物种占63.59%(图3)。

图2 COG, GO, NR, Pfam及Swissprot功能注释维恩图Fig.2 A Venn diagram shows commonality and difference of annotation based on COG, GO, NR, Pfam and Swissprot

图3 阿尔泰蝠蛾Unigenes在Nr数据库中的物种分布图Fig.3 Species distribution of Unigenes of Hepialus altaicola in Nr database

2.2.2Unigene的COG分类注释

经过COG数据库功能预测和分类，共有11 084 个Unigenes被注释，根据其功能可分为26类(图4)。其中，一般功能预测类(General function prediction only，1 432个)基因最多；其次为转录、核糖体结构和生物合成(Translation, ribosomal structure and biogenesis，1 428个)；翻译后修饰，蛋白质折叠和分子伴侣(Posttranslational modification, protein turnover, chaperones，1 175个)、碳水化合物运输和代谢(Carbohydrate transport and metabolism，1 059个)，氨基酸转运和代谢(Amino acid transport and metabolism，826个)，其他类别的基因表达丰度各不相同。COG分类揭示了可能参与调节阿尔泰蝠蛾幼虫在冷适应过程中基因的特定反应和功能。

2.2.3Unigene的GO分类注释

共有13 261个Unigenes在GO数据库中3大类57个功能组中找到对应，分别是分子功能(Molecular function)分为15个亚类，细胞成分(Cellular component)的19个亚类，生物学过程(Biological process)23个亚类。分子功能类注释到15 023条，其中催化活性(Catalytic activity)和结合活性(Binging)的Unigenes数量较多，分别是6 871和5 786条，其余均在1 000条以下。细胞成分类注释到22 008条，其中细胞(Cell)和细胞部分(Cell part)的Unigenes数量较多，分别是4 529和4 519条；生物学过程类注释到23 963条，排在前三位的是代谢过程(Metabolic process)、细胞过程(Cellular process)和单一有机体过程(Single-organism process)，分别是6 257、6 217和3 581条，其余的均在2 000条以下(图5)。

2.2.4Unigene的KEGG分类注释

使用KEGG数据库进行Unigenes代谢途径分析。结果显示，共有15 058个Unigenes参与了环境信息处理(Environmental information processing)、遗传信息处理(Genetic information processing)、细胞过程(Cellular processes)、代谢(Metabolism)和有机体系统(Organismal systems)5类代谢分支，归属于305条通路。其中注释最多的5个代谢通路依次为核糖体(678)、碳代谢(431)、剪接体(414)、内质网蛋白加工(395)、嘌呤代谢(371)(表3)。

图4 阿尔泰蝠蛾Unigenes的COG分类图Fig.4 COG functional categories of Hepialus altaicola Unigenes注：图中括号内数字代表Unigenes的数量，百分数代表其所占比例。Note:The numbers in brackets in the figure represent the number of unigenes, and the percentage represents its propotion.

图5 Unigenes的GO分类图Fig.5 GO functional categories of Unigenes

表3 阿尔泰蝠蛾的KEGG通路总结

2.3 CDS预测

编码序列(CDS)预测结果显示，共54 002条序列可被编码，占全部基因的53.93%(54 002/100 133)。长度在300 nt以上的Unigene为40 237条，占74.51%。其中0～100 nt的最多，有23 591条；其次是100～200 nt(11 072条)及200～300 nt(5 574条)，其余长度区间的CDS都在5 000条以下(图6)。

图6 编码序列(CDS)长度分布图Fig.6 Coding sequence (CDS) length distribution

2.4 与环境应激相关的代谢调节基因的分析

对阿尔泰蝠蛾基因表达量进行分析，已有基因功能注释信息中的大多数Unigenes与生物过程中的新陈代谢有关。结合已有的相关文献报道，在转录组数据中筛选了311个可能与冷适应相关的蛋白和酶的代谢调节基因，例如昆虫表皮蛋白(Insect cuticle protein, ICP)，热休克蛋白(Heat shock proteins, HSPs)，海藻糖转运蛋白(Trehalose transporter 1, TRET1)，谷胱甘肽S-转移酶(Glutathione S-transferase, GST)和几丁质酶Chitinase(表4)。其中大多数与HSP相关的Unigenes被预测为编码HSP70(31条)和HSP90(18条)型，其余还包括HSP40、HSP60及分子量小于30 kDa的小分子热休克蛋白(small Heat Shock Proteins,sHSPs)。这些注释信息能够为阿尔泰蝠蛾幼虫环境适应性中发挥潜在作用的基因功能研究提供新的线索。

表4 阿尔泰蝠蛾抗逆性相关基因的推测

表5 阿尔泰蝠蛾转录组FPKM统计

3 结论与讨论

随着测序技术的发展，转录组分析逐渐成为非模式生物发掘功能基因的有效手段(Veraetal., 2008)。本研究应用Illumina高通量无参考基因组转录组测序，获得了阿尔泰蝠蛾的转录组信息，通过数据库分析和基因功能注释，从中发掘与应激代谢相关的抗逆性基因。通过序列的处理、Trinity拼接、转录本功能注释等步骤，得到数据量为21.31 G的原始转录组信息，共144 391个转录本，进一步分析得到100 133个Unigenes，总长度86 319 112 bp，平均长度和N50长度分别为862 bp和1 628 bp。同为蝙蝠蛾科的小金蝠蛾Hepialusxiaojinensis转录组测序得到50 164条Unigenes，N50为1 357 nt(Mengetal., 2015)。一般认为N50值越大，反映组装得到的长片段越多，组装效果越好。从获得基因的数量、平均长度和N50值等方面均显示本次阿尔泰蝠蛾转录组的拼接效果良好；碱基Q30值>92%，说明碱基识别率较高，Q30值在80%以上认为测序质量可靠(王兴春等，2015)。上述结果表明本研究测序和组装结果较好，其质量和长度满足转录组分析的基本要求。

图7 阿尔泰蝠蛾转录组中表达丰度较高的50个基因Fig.7 Top 50 most abundant Unigenes in the antennal transcriptome of Hepialus altaicola

基因数据库中参考序列的数量是决定转录组数据注释比例的主要因素，例如家蚕(Xiaetal., 2004)、大蜡螟Galleriamellonella(Vogeletal., 2011)等研究鳞翅目昆虫的模式生物，其全基因组测序为鳞翅目昆虫生物学和免疫学等基因功能的研究奠定了基础。本研究共计有38 198条(占38.15%)Unigenes被Nr、Swiss-Prot、Pfam、GO、COG/KOG和KEGG等数据库注释，但仍有61 935条(61.85%)序列无法成功比对，这可能与拼接的序列片段过短而缺乏注释信息有关(魏利斌等，2012)，也可能本身是非编码序列或者是新的某种基因(Houetal., 2011)，而与已报道其他昆虫缺少同源性。

GO是国际标准的基因功能分类体系，能够全面描述生物体基因和基因产物的属性，从宏观上认识某一物种的基因功能分布特征(Botsteinetal., 2000)。阿尔泰蝠蛾转录组Unigenes在GO数据库中共得到的13 261条功能注释，其中参与代谢、催化与结合活性的基因功能注释较多，这与其它昆虫的转录组研究结果相一致(Zhangetal., 2015; 孟翔等, 2016)。不同昆虫的GO二级功能注释有所差异，如意大利蝗Calliptamusitalicus二级功能注释有52个(向敏，2017)；大蜡螟有55个(杨爽，2019)；麦红吸浆虫Sitodiplosismosellana有50个(蒋月丽等，2020)；而阿尔泰蝠蛾GO二级功能注释有57个。阿尔泰蝠蛾幼虫转录组在KEGG数据库注释到305条代谢通路，主要涉及核糖体、碳代谢等，其中参与核糖体代谢(ko03010)的基因数量最多，达678条，该代谢通路由核糖体蛋白主导，其基因对核糖体蛋白的合成起重要作用，能够影响生物体的生命活动(Uechietal., 2001)。

在GO功能富集中，注释到与生物学过程中的代谢过程相关的Unigenes占比最高(占26.11%)。共分析得到311个与阿尔泰蝠蛾幼虫应激代谢调节相关的Unigenes，分属于11个不同的类群(表4)。这些基因被归类为以下功能：蛋白质折叠(GO：0006457)，对压力的反应(GO：0006950)，葡萄糖转运(GO：0015758)，跨膜转运蛋白活性(GO：0022857)，热激蛋白结合(GO：0031072)，ATP代谢过程(GO：0046034)，超氧化物代谢过程(GO：0006801)，水解酶活性(GO：0016787)，三羧酸循环(GO：0006099)等。响应冷适应调节的基因“热休克蛋白”(90个)，参与催化活性的“ATP合酶”(62个)的基因表达最多。热休克蛋白是重要的分子伴侣，它的表达和调控是物种响应环境胁迫的物质基础，在昆虫的应激响应中发挥重要作用(Wangetal., 2019)，对耐受低温环境得以生存至关重要(Rinehartetal., 2007)。在所注释的热休克蛋白中，HSP70和HSP90基因数量最多，这与分布于低温高原地区的竹蝗Ceracriskiangsu在15℃冷驯化下表达量类似(Lietal., 2019)。而热休克蛋白的合成，消耗了生物体大量的能量，因此，参与ATP合成的基因在低温胁迫后大量表达。在阿尔泰蝠蛾转录组中有22个Unigenes被注释为昆虫表皮蛋白，且表达量较高，这些表皮蛋白在昆虫表皮形成，它的产生是一种细胞保护策略，以应对环境挑战(Zhangetal., 2008)。海藻糖酶可以水解昆虫血淋巴中的海藻糖，产生2分子葡萄糖，能够为昆虫生长发育、代谢、大分子生物合成等提供所需能量，在昆虫蜕皮过程中的几丁质合成、耐寒性方面也发挥重要作用(Shietal., 2016)。有49个Unigenes被注释为海藻糖转运蛋白，帮助昆虫脂肪体中合成的海藻糖进入胞外，通过血淋巴运送到各组织，为生命活动提供能量(Kikawadaetal., 2007)。此外，有14个Unigenes被注释为水通道蛋白，参与体内渗透压的调节，其表达量与抗冻能力具有相关性(Philipetal., 2008)；有14个Unigenes被注释为超氧化物歧化酶(Superoxide dismutase, SOD)，这与耐寒昆虫准噶尔小胸鳖甲Microderapunctipennis在环境胁迫适应中SOD基因的表达相一致(Luetal., 2014)。阿尔泰蝠蛾幼虫栖息的低温高海拔环境，全年仅6-8月土壤层的日平均温度可达到10℃以上，冻土期还限制了幼虫的取食和运动。因此，这些Unigenes在野外采集的蝠蛾幼虫中表达，表明低温下的代谢补偿可能是其适应寒冷环境的重要策略。阿尔泰蝠蛾幼虫在环境变化过程中通过控制体内小分子运输，来降低应激条件下细胞内渗透压对刺激的反应能力，有助于阐明阿尔泰蝠蛾幼虫对低温的响应与应激代谢调节之间的相互关系。

温度是影响阿尔泰蝠蛾物种分布的重要因子，也是生态学研究的关键因素之一。生活于低温高海拔地区的阿尔泰蝠蛾幼虫，如何适应低温环境是一个值得探究的问题。本研究获得的转录组数据和分析结果，有助于快速发现阿尔泰蝠蛾广泛多样的候选基因，为进一步研究蝠蛾幼虫所栖居的独特生境下，冷适应相关机制和代谢平衡的调节提供分子水平依据，同时也丰富了鳞翅目蝙蝠蛾科昆虫的基因数据库。另外，不同昆虫参与冷适应相关的代谢调节基因的数目和种类也有所不同，被注释的这些基因的组织表达特异性和功能还需进一步分析验证。