美洲大蠊microRNA的初步研究

宝钲， 杨茗羽， 张修月， 岳碧松， ， ， 范振鑫*

(1.生物资源与生态环境教育部重点实验室，四川大学生命科学学院，成都610065；2.四川药用动物工程技术研究中心，四川西昌615000； 3. 药用美洲大蠊四川省重点实验室，成都610031)

microRNA(miRNA)是一类长度18～25 nt的非编码小RNA，广泛存在于动物、植物、微生物中，在转录后水平调控基因表达和翻译(Bartel，2004；Bartel & Chen，2004)。miRNA在生物中的调控方式一般有3种：一种是与靶mRNA的3’UTR完全互补配对，导致mRNA直接被降解(Llaveetal.，2002)；一种是与靶mRNA的3’UTR不完全互补配对，不影响mRNA的稳定性，只是抑制其翻译(Slacketal.，2000)；还有一种是同时具有上述2种方式(Lee & Ambros，2001)。目前大量研究已揭示miRNA在昆虫多种生理活动中均有重要功能，如发育(Cristinoetal.，2011；Surridgeetal.，2011；Chenetal.，2012)、激素信号传导(Rebijithetal.，2016)、繁殖(Pauloetal.，2017)等。在蜚蠊目Blattodea中，目前仅有德国小蠊BlattellagermanicamiRNA的相关报道(Cristinoetal.，2011)。

美洲大蠊Periplanetaamericana属节肢动物门Arthropoda昆虫纲Insecta蜚蠊目，是蜚蠊科Blattidae中体型最大的一种，也是世界上最常见的室内蟑螂之一。美洲大蠊对各类环境适应性强，同时也传播病菌和寄生虫，被认为是世界性卫生害虫(Chenetal.，2015；Salama，2015；Zhangetal.，2016；Lietal.，2017)。但美洲大蠊也是传统的药用昆虫，其提取物具有提高免疫、抗肿瘤、抗炎镇痛、保肝、组织修复等效果(肖小芹等，2007；Wangetal.，2011；Jiangetal.，2012；王鹏飞等，2015；张丹等，2015；吴道勋等，2016)。目前有关美洲大蠊的研究多集中在基因组、繁殖、消化特性和转录组等(Wicheretal.，2006；Nishinoetal.，2010；Chenetal.，2013，2015；Tamakietal.，2014；Kimetal.，2016；Zhangetal.，2016；晋家正等，2018；Lietal.，2018)，美洲大蠊的miRNA及其在基因表达中的调控机制迄今尚未见报道。尤其是雌雄美洲大蠊在形态和寿命上存在明显的差异，雌虫比雄虫体型大，雄虫比雌虫寿命短(陈梦林，2002)，这种差异是否与雌雄之间miRNA调控的差异有关？本研究对实验室饲养的雌雄美洲大蠊进行小RNA的高通量测序，去除低质量和受污染的序列后，将所得序列与数据库(NCBI、Rfam)进行比对注释，旨在初步了解美洲大蠊小RNA的组成，包括tRNA、rRNA、snRNA、snoRNA等；将非miRNA的小RNA过滤后，检测已知的miRNA、预测新miRNA并寻找雌雄之间的差异表达miRNA。

1 材料与方法

1.1 样品准备与小RNA的文库构建和测序

美洲大蠊成虫和饲料均由四川好医生攀西药业有限责任公司提供。取饲养于本实验室(温度25 ℃，湿度60%)的雌雄成虫各3只。为避免消化道中细菌对测序结果的干扰，将去除消化道的美洲大蠊搅碎后加入磷酸盐缓冲液离心，加入异丙醇提取总RNA。检测合格的总RNA由北京诺禾致源科技股份有限公司进行建库和高通量测序，测序平台为Illumina HiSeq 2000，测序长度为50 bp。

1.2 数据处理

测序完成后，去除低质量和被污染的序列，剩下的序列(干净序列)用于后续分析。

为了防止非miRNA(包括tRNA、rRNA、scRNA、snRNA、snoRNA等)序列影响分析，将得到的干净序列用Blastn与数据库(NCBI、Rfam)中非miRNA进行比较，再利用Repeatmasker找出小RNA中的重复序列，然后注释这些非miRNA并在后续分析中去除。

剩下的未注释序列利用miRDeep2(Friedländeretal.，2011)检测已知的miRNA和预测潜在的新miRNA。因为miRBase 21数据库(ftp：//mirbase.org/pub/mirbase/21)中目前还没有美洲大蠊的miRNA数据，但miRNA在物种间具有保守性，所以通过miRBase 21数据库中其他昆虫已知的miRNA对美洲大蠊进行检测。分析所需要的基因组数据来自本实验室测序和组装注释的美洲大蠊基因组(晋家正等，2018)。

1.3 预测靶基因及KEGG通路分析

利用miRanda(Enrightetal.，2003)和RNAhybrid(Krüger & Rehmsmeier，2006)预测美洲大蠊miRNA的靶基因。miRanda的阈值设定为：分数≥155，自由能ΔG≤-20 kcal/mol。RNAhybrid的阈值设定为：ΔG≤-20 kcal/mol，其余参数为默认值(Enrightetal.，2003；Krüger & Rehmsmeier，2006)。为使预测结果更加可信，将2个软件的预测结果取交集进行后续的分析。

通过KEGG通路分析对miRNA的靶基因进行注释。利用KAAS (Moriyaetal.，2007)和KEGG mapper software (Kanehisaetal.，2012)对靶基因进行KEGG注释。

1.4 已知miRNA的差异表达

将雄性成虫和雌性成虫的miRNA分别进行统计分析以鉴定雌雄之间的差异表达。先将miRNA的表达量归一化，然后再计算P值。

归一化公式：归一化表达量(NE)=实际miRNA数量/clean reads总数×1 000 000。

P值公式：

式中，N1和x属于相同的组，N1代表clean reads的总数，x代表归一化的表达水平。N2和y属于相同的组，N2代表clean reads的总数，y代表归一化的表达水平(Zhangetal.，2013)。

2 结果

2.1 测序数据、长度分布和其他小RNA的注释

通过测序，在雄性成虫和雌性成虫中分别得到17 259 315条和17 249 947条小RNA序列，剔除质量低和被污染的序列后，分别得到12 155 616条和9 847 263条。其中，共有序列占总序列的84.08%，雄性成虫特有的序列占总序列的8.81%，雌性成虫特有的序列占总序列的7.11%(图1：a)。雌雄成虫中所测得的所有小RNA序列均主要分布于18～23 nt，在22 nt和29 nt处有2个峰值(图1：b)。在这些序列中，将相同的序列统计为同一种得到唯一序列，最终在雄性成虫和雌性成虫中分别得到2 247 751条和2 051 813条小RNA序列。

通过对小RNA注释，发现许多非miRNA存在，包括rRNA、tRNA、snRNA、snoRNA、scRNA、短小的重复序列、内含子、外显子等(表1)。为了防止这些序列干扰后续分析，去除注释为非miRNA的序列。最终在雄性成虫和雌性成虫中分别获得11 507 387条和9 207 534条可能为miRNA的序列，这些序列除部分序列没有注释外，绝大多数被注释为miRNA。

2.2 miRNA及靶基因预测和KEGG通路分析

在雄性成虫和雌性成虫中分别发现了57种和53种已知的miRNA，去除二者共有的miRNA，共发现59种已知的miRNA。表达量前10的miRNA见表2，其中，miR-1、miR-10-5p、miR-8-5p、miR-31a等广泛参与了各种组织和细胞的发育、生长和代谢等重要生理过程。

图1 美洲大蠊雌雄成虫相同序列的小RNA数量和长度分布Fig. 1 Venn diagram of the same sequence of small RNAs and their length distributions in female and male Periplaneta americana

组别雄性成虫雌性成虫唯一序列数量/条(占总小RNA的比例)全序列数量/条(占总小RNA的比例)唯一序列数量/条(占总小RNA的比例)全序列数量/条(占总小RNA的比例)总小RNA2 247 751(100%)12 155 616(100%)2 051 813(100%)9 847 263(100%)rRNA47 172(2.10%)384 358(3.16%)57 580(2.81%)392 256(3.98%)tRNA29 363(1.31%)227 840(1.87%)28 868(1.41%)218 086(2.21%)snoRNA3 145(0.14%)9 626(0.08%)2 489(0.12%)6 692(0.07%)snRNA3 688(0.16%)11 157(0.09%)2 796(0.14%)9 860(0.10%)piRNA247(0.01%)3 866(0.03%)183(0.01%)2 106(0.02%)scRNA0(0%)0(0%)1(0%)1(0%)外显子23(0%)80(0%)15(0%)54(0%)内含子1 440(0.06%)3 611(0.03%)1 958(0.10%)5 208(0.05%)重复序列3 187(0.14%)7 691(0.06%)2 401(0.12%)5 466(0.06%)未注释序列2 159 486(96.07%)11 507 387(94.67%)1 955 522(95.31%)9 207 534(93.50%)

表2 美洲大蠊雌雄成虫表达量前10的miRNATable 2 Top 10 highly expressed miRNAs in female and male Periplaneta americana

同时，在雄性成虫和雌性成虫中预测出152种和94种潜在的新miRNA，去除共有的miRNA，共发现152种新miRNA。这些新miRNA在miRBase 21数据库中尚未报道，可能是美洲大蠊所特有的miRNA。

进一步对预测的潜在的新miRNA的靶基因进行预测，将miRanda和RNAhybrid预测结果整合后，共计在雄性成虫和雌性成虫中预测了5 389个和3 823个靶基因，去掉共有的，共5 822个靶基因。KEGG通路分析(图2)显示了这些靶基因的功能，其中，信号转导通路注释到的靶基因最多，且雄性成虫比雌性成虫多，说明美洲大蠊miRNA调控的靶基因主要是参与信号转导通路的基因。

2.3 已知的miRNA在雌雄美洲大蠊中的差异表达

雌雄美洲大蠊成虫差异表达的miRNA的散点图显示仅有1条显著差异表达的miRNA：miR-750，其在雌性成虫中的表达显著高于雄性成虫(图3)。

3 讨论

miRNA越来越受到研究者的关注，这是由于miRNA被发现在越来越多生物过程中发挥着重要作用，如发育(Reinhartetal.，2000)、细胞分化(Kawasaki & Taira，2003)、细胞增殖(Brenneckeetal.，2003)等。很多昆虫的miRNA已被报道(Weietal.，2009；Jagadeeswaranetal.，2010；Skalskyetal.，2010；Zhangetal.，2012；Huangetal.，2014；Heetal.，2016)。本研究是首次对美洲大蠊的miRNA进行报道，在雄性和雌性成虫中分别发现了2 247 751条和2 051 813条miRNA序列。在miRNA序列长度分布上，显示出2个峰值，在21～23 nt的峰代表的是miRNA，在28～30 nt的峰代表的可能是Piwi蛋白互作RNA(piRNA)，这是一种与miRNA完全不同的RNA沉默因子(Kawaokaetal.，2011)。这一现象在其他昆虫，包括飞蛾Manducasexta(Zhangetal.，2012)、桔小实蝇Bactroceradorsalis(Huangetal.，2014)等都存在相似的结果。与德国小蠊的miRNA(Cristinoetal.，2011)相比，美洲大蠊表达量最高的前10个miRNA中，有4个与德国小蠊相同：miR-1、miR-8-5p、miR-10-5p和miR-31a。在德国小蠊中表达、在美洲大蠊中不表达的miRNA有miR-184、miR-276a、miR-279b、miR-375和miR-iab-4as-5p。然而，更多的miRNA在美洲大蠊中表达而未在德国小蠊中发现，如miR-133、miR-277-3p、miR-278、miR-2765-5p、miR-281-5p、miR-316和miR-9b-5p。其中，miR-1在心脏和骨骼肌祖细胞的发育上发挥重要作用(Zhao & Srivastava，2007)；miR-10-5p的靶基因Hox基因决定体轴的正确前后样式(Lund，2010)；miR-8-5p参与众多生物学过程，包括发育(Bellesetal.，2012)、神经肌肉协调(Karresetal.，2007；Rubioetal.，2013)、卵及卵巢形成(Lucasetal.，2015)等；miR-31a可能在蜕皮过程中调控上皮代谢(Yuetal.，2008)；miR-276通过调控1个转录共活化基因调控孵化同步性(Heetal.，2016)；miR-133的其中1个靶基因是结缔组织生长因子(CTGF)，其在纤维化过程中发挥重要作用(Duistersetal.，2009)；miR-281的靶基因是登革热病毒血清2型(DENV-2)，增强在白纹伊蚊Aedesalbopictus中DENV-2的病毒复制(Zhouetal.，2014)。本研究结果显示，美洲大蠊的miRNA广泛参与了昆虫组织和细胞的生长、发育、代谢等多种重要的调控机制。此外，美洲大蠊和德国小蠊都存在一些特有的miRNA调控的基因，显示2种蜚蠊目昆虫存在不同的调控机制，这或许是造成二者形态和生理差异的原因之一。

图2 KEGG通路中美洲大蠊雌雄成虫miRNA的靶基因数量
Fig. 2 Numbers of target genes of miRNA in female and malePeriplanetaamericanaannotated by KEGG

图3 美洲大蠊雌雄成虫已知miRNA的差异表达散点图Fig. 3 The scatter figure of differentially expressed known miRNAs between female and male Periplaneta americana

差异表达分析显示，美洲大蠊雌雄成虫miRNA在组成和表达量上几乎没有差异，只有miR-750在雌性成虫中高表达。miR-750在被白斑综合症病毒WSSV感染的斑节对虾Penaeusmonodon中表达量显著降低，其预测的靶基因是ATG14基因，但随后的实验证实ATG14基因不是miR-750的靶基因(Kaewkascholkuletal.，2016)。Rebijith等(2016)通过RNAi实验在棕榈蓟马Thripspalmi中证明JH候选受体——过剩气门蛋白是miR-750的靶基因，表明miR-750可能通过调节卵黄蛋白原基因参与免疫、应激反应和激素信号传导。卵黄蛋白原主要存在于非哺乳类性成熟的卵生雌性动物的血淋巴、脂肪体和卵中，所以miR-750在雌性成虫中高表达。总之，对雌雄美洲大蠊进行了小RNA的测序，检测出了57种已知的miRNA，并预测出152种潜在的新miRNA，但这些数据是用软件预测出来的，这些潜在的新miRNA是否具有功能及其相应的靶基因，还需要进一步的验证。这些miRNA的发现为将来美洲大蠊中miRNA的研究提供了基础，其表达量和功能的研究将推进对于美洲大蠊调控网络的理解，并为美洲大蠊的虫害防治和在中药中的应用提供新的视角。