基于高通量测序的玉米象转录组分析

李媛媛，赵金红

(皖南医学院医学寄生虫学教研室，安徽芜湖 241002)

玉米象Sitophiluszeamais属鞘翅目Coleoptera象甲科Curculionidae，别名米牛、铁嘴，在全世界范围内均有分布，是储粮的世界性头号害虫。玉米象可孳生于多种谷物及加工品、干果、豆类、药材、油料等，尤其是阴暗潮湿的仓库为最为严重(郑旭等，2014；姜洪等，2018)，不仅可引起储粮中谷物数量减少，而且还可改变谷物及其产品的质量，导致种子活力的退化(Haqetal.，2000；王鹏杰等，2020)。目前对于玉米象的研究主要集中在生物、生态学特性以及防治研究方面，对其使用化学合成熏蒸剂进行防治，但长期使用化学制剂可导致玉米象产生严重抗药性(吕建华等，2010；侯昌亮等，2016)。目前，由于缺乏在基因组和转录组方面研究，有效的分子标记尚未被发现，导致玉米象的抗药性和种群遗传结构的研究尚未开展。因此，本文通过获得大量转录组信息，以期为功能基因的挖掘提供信息资源，为玉米象种群遗传结构、抗药性监测以及防治奠定理论基础。

转录组广义上是指生物体细胞或组织在特定状态下所转录出来的所有RNA的总和(张春兰等，2012)。转录组学是从整体转录水平系统研究基因转录图谱，并揭示复杂生物学通路和性状调控网络分子机制的学科(崔凯等，2019)。随着测序技术的发展，核酸的检测和定量更加便捷和准确，越来越多的研究学者将高通量测序技术应用到转录组研究中。对于危害农作物的昆虫转录组研究工作也有相关报道，如大垫尖翅蝗Epacromiuscoerulipes(金永玲等，2015)、东方粘虫Mythimnaseparata(胡艳华等，2015；李微等，2017)、沟眶象Eucryptorrhynchuschinensis(武政梅，2016)、疣蝗Trilophidiaannulata(邱忠营等，2016)、印度谷螟Plodiainterpunctella(李慧等，2019)、麦红吸浆虫Sitodiplosismosellana(蒋月丽等，2020)、红棕象甲Rhynchophorusferrugineus(Yangetal.，2020)。本研究采用高通量测序技术对玉米象进行转录组测序、拼接与组装，对所获得的unigenes进行基因功能注释、代谢途径及SSR检测等分析，为玉米象进一步的分子标记开发和基因功能研究提供分子基础。

1 材料和方法

1.1 试虫

供试的玉米象成虫采集于安徽省芜湖市粮库，按照文献(邓树华等，2011)培养条件进行培养。恒温培养箱：温度27℃±2℃，相对湿度75%±5%，光周期L ∶D=16 h ∶8 h，小麦饲养7 d后，筛去成虫，待新一代成虫羽化大量出现一周左右，选取2头羽化后7～14 d的健壮活泼玉米象成虫作为研究对象，用无菌水清洗虫体，液氮速冻后保存于-80℃冰箱备用。

1.2 方法

1.2.1总RNA提取

采用Total RNA Extractor试剂盒提取液氮冷冻保存的玉米象总RNA，将提取的RNA进行电泳检测，Nanodrp检测RNA的纯度(OD260/OD280)，利用Qubit对RNA浓度进行精确定量，RNA的完整性用Agilent 2100进行检测。将满足测序要求的RNA样本进行转录组测序。

1.2.2文库构建和测序

样品检测合格后，先采用Oligo(dT)方法分离并纯化出mRNA，随后将mRNA打断为短片段，用随机引物合成一链cDNA(互补DNA)和二链cDNA，对双链cDNA进行纯化、末端修饰、片段大小的选择，最后进行PCR富集得到最终的cDNA文库。使用Agilent 2100对文库进行准确定量，以保证文库的质量。库检合格后进行Illumina Hiseq测序。

1.2.3数据组装与功能注释

采用Illumina Hiseq测序得到原始图像数据，经碱基识别转化为序列数据，分别去除带接头的序列、无法确定的碱基数量(N%)>10%的序列和低质量(质量值≤Q20的碱基数占50%以上)的序列。利用Trinity软件对样品数据进行组装拼接后，分别与7个数据库(Nr、GO、Swissprot、KOG、Pfam、KEGG、TrEMBL)中的蛋白质序列进行BLAST比对，从而对其功能进行注释。

1.2.4简单重复序列(SSR)及单核苷酸多态性(SNP)分析

利用MISA软件对玉米象中筛选得到的1 kb以上的Unigene进行SSR位点分析。通过Bcftools软件对单核苷酸多态性(SNP)进行鉴定。

2 结果与分析

2.1 玉米象转录组测序质控与组装

对获得的原始数据进行统计，共产生44 277 838条原始reads，总碱基数为6.64 G。再消除带接头和低质量的reads后，获得pair-end Reads总数为43 331 830，总碱基数为6.34 G，GC含量为38.66%，Q20碱基比例平均超过99.27%。经过组装，共获得64 358条unigenes，总长度39 481 752 bp，最短201 bp，最长29 046 bp，平均长度613 bp，N50长度为871 bp。(图1)

图1 组装序列长度分布图Fig.1 Size distribution of transcripts and unigenes

2.2 功能注释

分别在7个数据库中对组装得到的64 358条unigenes进行注释(表1)。不同数据库中注释成功的unigenes基因数及其所占比例有所差别。有915条unigenes同时在7个数据库中注释成功，占总数的1.42%；在这7个数据库中至少有一个数据库注释成功的有26 580条，占总数的41.3%。在Nr数据库当中，有25 271条unigenes能够寻觅到近似序列，约占unigenes总数的39.27%；在GO数据库当中，有20 747条unigenes获取了注释，约为总数的32.24%；在Pfam和KEGG数据库汇中获得注释的unigenes数量都10 000条以下，分别为9 471(占总体数的14.78%)和3 370(占总体数的5.24%)。

2.2.1Unigenes的GO分类注释

按照GO功能分类方式，将GO注释的20 747条unigenes分为生物过程(biological process)、细胞组分(cellular component)和分子功能(molecular function)3大类(图2)。

所有大类可细分为68个二级分类。对unigenes在二级分类中的分布情况进行统计分析，“生物学过程”包含26个不同的亚类，它也是3大类中最多的亚类，其中细胞过程(cellular process)和生物调节(biological regulation)所占比例较高，分别为13 426条和7 814条，分别占GO注释信息总数的67.71%和37.66%；在“细胞组分”类别中有22个亚类，其中细胞(cell)和细胞部分(cell part)类型的unigenes较高，分别为13 710(66.08%)和13 664(65.86%)；在“分子功能”类别中有20个亚类，其中连接(binding)和催化活性(catalytic activity)的unigenes较高，分别为12 704(61.23%)和9 356(45.1%)。

2.2.2Unigenes的KOG分类注释

在KOG数据库中比对玉米象转录组unigenes，结果显示共12 060条unigenes获得注释信息，根据其功能大致可分为25类(图3)。其中，信号传导机制(Signal transduction mechanisms)、一般功能预测基因(General function prediction only)、翻译后修饰、蛋白质折叠和分子伴侣(Posttranslational modification, protein turnover， chaperones)获得注释最多，分别为1 962条、1 648条、1 135条，而细胞活性(Cell motility)获得注释最少，仅20条。

图2 玉米象GO功能分类图Fig.2 Statistical diagram of GO function of Sitophilus zeamais注：1,抗氧化功能；2,连接功能；3,催化功能；4,趋化活性；5,化学抗拮；6,电子转移活性；7,金属伴侣活性；8,分子功能调节；9,分子传导活性；10,形态发生活性；11,营养库活性；12,蛋白标志物；13,信号转导活性；14,结构分子活性；15,转录因子活性，蛋白质结合；16,翻译调节功能；17,载体活性；18,通道调节活性；19,酶调节功能；20,受体调节功能；21,生物习性；22,生物黏附；23,生物相；24,生物调节；25,细胞聚集；26,细胞杀伤；27,组织或生物源的细胞组分；28,细胞过程；29,解毒；30,发育过程；31,生长；32,免疫反应过程；33,定位；34,移动；35,新陈代谢过程；36,多有机体过程；37,多细胞生物过程；38,生物负调控过程；39,生物正调控过程；40,生物过程调控；41,生殖；42,生殖过程；43,刺激反应；44,节律进程；45,信号；46,建立定位；47,细胞；48,细胞结合；49,细胞成分；50,细胞外基质；51,细胞外基质成分；52,细胞外区域；53,细胞外区域成分；54,蛋白质复合物；55,细胞膜；56,细胞膜成分；57,膜包围的内腔；58,内核；59,细胞器；60,细胞器成分；61,其他有机体；62,其他有机体成分；63,超分子纤维；64,共质体；65,突触；66,突触成分；67,病毒；68,病毒部分。

图3 玉米象KOG功能统计图Fig.3 Statistical diagram of KOG function of Sitophilus zeamais

2.2.3Unigenes的KEGG分类注释

使用KEGG数据库对玉米象的unigenes可能参与或涉及的代谢途径进行了统计分析，在组装得到的unigenes中，有3 370条得到注释。这些通路信息分可为5大类，分别为细胞进程(Cellular Processes)、环境信息处理(Environmental Information Processing)、遗传信息处理(Genetic Information Processing)、新陈代谢(Metabolism)和有机系统(Organismal Systems)，这5大类可以进一步分为33个亚类274个功能通路(图4)。其中信号转导(Signal transduction)获得注释信息最多，有724条，参与较多的几条信号转导通路分别是PI3K-Akt信号通路(102条)、MAPK信号通路(89条)、mTOR信号通路(89条)、cMAP信号通路(79条)，这些代谢通路与环境信息大类中信号转导有关；翻译(Translation)获得的注释信息有435条，位居第2位；膜运输(Membrane transport)和其它此生代谢产物的生物合成(Biosynthesis of other secondary metabolites)的注释最少，仅有24条和26条。

图4 玉米象KEGG通路分类统计图Fig.4 Statistical diagram of KEGG Pathway of Sitophilus zeamais

2.3 玉米象SNP特征分析

对玉米象序列进行检测，共发现146 081个单核苷酸多态性(SNP，Single Nucleotide Polymorphsims)位点。对这些SNP位点进行类型统计，结果显示转换突变类型的SNP有102 850个(占70.41%)，颠换突变类型的SNP有43 231个(占29.59%)。在转换突变类型中，由胞嘧啶转换为胸腺嘧啶的突变最多(25 808个)，其次为鸟嘌呤转换为腺嘌呤的突变(25 786个)；在颠换突变类型中，由腺嘌呤颠换为胸腺嘧啶的突变最多(7 152个)，而鸟嘌呤颠换为胞嘧啶的突变最少(4 084个)(图5)。

图5 玉米象SNP突变类型统计图Fig.5 Statistical diagram of SNP mutation of Sitophilus zeamais

2.4 玉米象微卫星位点(SSR)分析

使用软件MISA对玉米象转录组中的64 358条unigenes进行搜索，共发现5 002个SSR位点，出现频率为7.772%(检测出的SSR数量与总unigenes数目的比值)，分布于4 606条unigenes中，发生频率为7.157%(含有SSR的unigenes的数目与总unigenes数目的比值)。玉米象SSR重复类型丰富，一至六核苷酸重复均有出现，但是数量分布差异较大。单核苷酸重复所占比例最大，为84.09%；其次为三核苷酸和二核苷酸重复，分别所占比例为9.42%和5.12%；四核苷酸重复、五核苷酸重复和六核苷酸重复分别占1.10%、0.26%和0.01%(表2)。其中，发生从分布情况来看，玉米象转录组序列中平均每11.805 kb(序列总长度与SSR总数目的比值)出现一个SSR，表明该物种转录组SSR数量较为丰富。

表2 玉米象转录组中SSR不同重复类型数量和分布

3 结论与讨论

转录组学能从整体水平上研究细胞中基因转录的情况及其转录调控规律(Shendure，2008；Grabherr，2011)，通过高通量测序技术和生物信息分析手段对玉米象的RNA进行测序和分析，获得的转录组数据可为玉米象分子生物学研究提供基因组数据来源。

本研究在无参考基因组的情况下，采用Illumina HiSeqTM测序技术对玉米象进行转录组测序，分析玉米象的转录组特征，共获得6.64 G原始数据，组装后得到64 358条unigenes，平均长度分别613 bp，N50长度为871 bp。结合生物信息学分析方法，对玉米象unigenes与7个公共数据库(Nr、GO、Swissprot、KOG、Pfam、KEGG、TrEMBL)进行比对分析，分别有25 271、20 747、16 477、12 060、9 471、3 370、25 500条unigenes获得注释，在7个数据库中同时得到注释的unigenes有26 580条，占全部unigenes总数的41.3%，但仍有37 778条(58.7%)未得到准确的定位，出现这一现象的原因，可能与目前测序技术的局限性组装得到的unigenes片段长度太短、数据库基因注释信息不足以及该物种存在新基因等因素有关。

其中，与Nr数据库进行比对时，注释的25 271条unigenes与中欧山松大小蠹Dendroctonusponderosae匹配率最高，为51.40%，这个结果和之前学者对沟眶象E.chinensis(武政梅，2016)的转录组进行研究时，得出的结论相同，其对获得的65 186条unigenes进行功能注释，有21 742条注释到Nr数据库，与中欧山松大小蠹的匹配率最高达到49.41%。

在GO注释的unigenes中，与细胞过程和生物调节、细胞组件及连接催化活性相关的基因最多，这些可能与玉米象细胞的增殖和分化，自身的生长发育及生命活动过程有着密切联系。通过比对KOG数据库，获得12 060个unigenes注释信息，对玉米象的unigenes的功能分布状况有了初步的了解。

KEGG通路注释分类中，共有3 370条玉米象unigenes参与到了274个代谢通路中，其中信号转导获得注释最多，参与PI3K-Ak(102条)、MAPK(89条)、mTOR(89条)信号通路的unigenes数量最多，这些信号通路在细胞生长、分化和凋亡方面都起到重要的作用，目前上述信号通路作为肿瘤领域的研究热点，被证明在多种人类肿瘤中异常激活，通路中的某些组分的突变可引起细胞转化、调控肿瘤细胞的存活、增殖及迁移，同时抑制自噬(舒婷等，2016；柳望舒等，2016)。在现今大多数昆虫抗药性研究中，功能基因组技术已被广泛用于候选基因调控抗药性机制的研究(陈皓等，2020)；随着昆虫基因组学、蛋白质组学、遗传学与分子生物学的应用，昆虫抗药性机制的研究已取得突破性进展(吴有刚等，2019)。通过对玉米象转录组的注释分析，可以从代谢调控中找到相关功能基因，并了解基因在代谢通路中的调控位置，对后续玉米象的抗药性、防治研究方面提供了基础数据。

SSR分子标记具有操作简便、多态性好、共显性以及数量丰富等特点(何海等，2015)。目前高通量测序可以从上万条基因序列中获得大量的微卫星。本研究对拼接得到的64 358条unigenes进行分析，发现5 002个SSR位点，统计发现单核苷酸和三核苷酸重复所占比例最大，分别为84.09%和9.42%。在之前的研究中，如温带臭虫Cimeslectularius(李敏等，2019)的SSR、东方粘虫(李微等，2017)的SSR、星天牛Anoplophorachinensis(韩小红等，2019)的SSR、大猿叶甲Colaphellusbowringi(沙君雪等，2018)以单核苷酸为主要类型，其次为三核苷酸重复，与本研究相符；疣蝗(邱忠营等，2016)的SSR主要以二核苷酸为主要重复类型；与玉米象同属于象甲科的沟眶象(武政梅等，2016)的SSR以三核苷酸为主。不同种属的昆虫SSR的主要重复类别不同，可以通过对SSR的密度、分布特点对不同的种进行区分。

目前对于象甲科昆虫分子生物学的研究大多集中在线粒体全基因组序列的主要结构特征特点分析上，转录组的相关报道并不多见。近期，有学者对和玉米象同属于隐颏象亚科Dryophthorinae的红棕象甲的转录组进行了研究，获得了蛹、幼虫和成虫3个时期大量的转录组数据，对昆虫发育相关基因转录组数据的分析将有助于害虫防治。本研究获得的数据可为玉米象后续的多态性监测、种群遗传结构提供基础资料，对玉米象的生物防治和其抗药性的研究提供参考依据。