寶鉦， 楊茗羽， 張修月， 岳碧松， ， ， 范振鑫*

(1.生物資源與生態(tài)環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，四川大學(xué)生命科學(xué)學(xué)院，成都610065；2.四川藥用動(dòng)物工程技術(shù)研究中心，四川西昌615000； 3. 藥用美洲大蠊四川省重點(diǎn)實(shí)驗(yàn)室，成都610031)

microRNA(miRNA)是一類長度18～25 nt的非編碼小RNA，廣泛存在于動(dòng)物、植物、微生物中，在轉(zhuǎn)錄后水平調(diào)控基因表達(dá)和翻譯(Bartel，2004；Bartel & Chen，2004)。miRNA在生物中的調(diào)控方式一般有3種：一種是與靶mRNA的3’UTR完全互補(bǔ)配對，導(dǎo)致mRNA直接被降解(Llaveetal.，2002)；一種是與靶mRNA的3’UTR不完全互補(bǔ)配對，不影響mRNA的穩(wěn)定性，只是抑制其翻譯(Slacketal.，2000)；還有一種是同時(shí)具有上述2種方式(Lee & Ambros，2001)。目前大量研究已揭示miRNA在昆蟲多種生理活動(dòng)中均有重要功能，如發(fā)育(Cristinoetal.，2011；Surridgeetal.，2011；Chenetal.，2012)、激素信號(hào)傳導(dǎo)(Rebijithetal.，2016)、繁殖(Pauloetal.，2017)等。在蜚蠊目Blattodea中，目前僅有德國小蠊BlattellagermanicamiRNA的相關(guān)報(bào)道(Cristinoetal.，2011)。

美洲大蠊Periplanetaamericana屬節(jié)肢動(dòng)物門Arthropoda昆蟲綱Insecta蜚蠊目，是蜚蠊科Blattidae中體型最大的一種，也是世界上最常見的室內(nèi)蟑螂之一。美洲大蠊對各類環(huán)境適應(yīng)性強(qiáng)，同時(shí)也傳播病菌和寄生蟲，被認(rèn)為是世界性衛(wèi)生害蟲(Chenetal.，2015；Salama，2015；Zhangetal.，2016；Lietal.，2017)。但美洲大蠊也是傳統(tǒng)的藥用昆蟲，其提取物具有提高免疫、抗腫瘤、抗炎鎮(zhèn)痛、保肝、組織修復(fù)等效果(肖小芹等，2007；Wangetal.，2011；Jiangetal.，2012；王鵬飛等，2015；張丹等，2015；吳道勛等，2016)。目前有關(guān)美洲大蠊的研究多集中在基因組、繁殖、消化特性和轉(zhuǎn)錄組等(Wicheretal.，2006；Nishinoetal.，2010；Chenetal.，2013，2015；Tamakietal.，2014；Kimetal.，2016；Zhangetal.，2016；晉家正等，2018；Lietal.，2018)，美洲大蠊的miRNA及其在基因表達(dá)中的調(diào)控機(jī)制迄今尚未見報(bào)道。尤其是雌雄美洲大蠊在形態(tài)和壽命上存在明顯的差異，雌蟲比雄蟲體型大，雄蟲比雌蟲壽命短(陳夢林，2002)，這種差異是否與雌雄之間miRNA調(diào)控的差異有關(guān)？本研究對實(shí)驗(yàn)室飼養(yǎng)的雌雄美洲大蠊進(jìn)行小RNA的高通量測序，去除低質(zhì)量和受污染的序列后，將所得序列與數(shù)據(jù)庫(NCBI、Rfam)進(jìn)行比對注釋，旨在初步了解美洲大蠊小RNA的組成，包括tRNA、rRNA、snRNA、snoRNA等；將非miRNA的小RNA過濾后，檢測已知的miRNA、預(yù)測新miRNA并尋找雌雄之間的差異表達(dá)miRNA。

1 材料與方法

1.1 樣品準(zhǔn)備與小RNA的文庫構(gòu)建和測序

美洲大蠊成蟲和飼料均由四川好醫(yī)生攀西藥業(yè)有限責(zé)任公司提供。取飼養(yǎng)于本實(shí)驗(yàn)室(溫度25 ℃，濕度60%)的雌雄成蟲各3只。為避免消化道中細(xì)菌對測序結(jié)果的干擾，將去除消化道的美洲大蠊攪碎后加入磷酸鹽緩沖液離心，加入異丙醇提取總RNA。檢測合格的總RNA由北京諾禾致源科技股份有限公司進(jìn)行建庫和高通量測序，測序平臺(tái)為Illumina HiSeq 2000，測序長度為50 bp。

1.2 數(shù)據(jù)處理

測序完成后，去除低質(zhì)量和被污染的序列，剩下的序列(干凈序列)用于后續(xù)分析。

為了防止非miRNA(包括tRNA、rRNA、scRNA、snRNA、snoRNA等)序列影響分析，將得到的干凈序列用Blastn與數(shù)據(jù)庫(NCBI、Rfam)中非miRNA進(jìn)行比較，再利用Repeatmasker找出小RNA中的重復(fù)序列，然后注釋這些非miRNA并在后續(xù)分析中去除。

剩下的未注釋序列利用miRDeep2(Friedl?nderetal.，2011)檢測已知的miRNA和預(yù)測潛在的新miRNA。因?yàn)閙iRBase 21數(shù)據(jù)庫(ftp：//mirbase.org/pub/mirbase/21)中目前還沒有美洲大蠊的miRNA數(shù)據(jù)，但miRNA在物種間具有保守性，所以通過miRBase 21數(shù)據(jù)庫中其他昆蟲已知的miRNA對美洲大蠊進(jìn)行檢測。分析所需要的基因組數(shù)據(jù)來自本實(shí)驗(yàn)室測序和組裝注釋的美洲大蠊基因組(晉家正等，2018)。

1.3 預(yù)測靶基因及KEGG通路分析

利用miRanda(Enrightetal.，2003)和RNAhybrid(Krüger & Rehmsmeier，2006)預(yù)測美洲大蠊miRNA的靶基因。miRanda的閾值設(shè)定為：分?jǐn)?shù)≥155，自由能ΔG≤-20 kcal/mol。RNAhybrid的閾值設(shè)定為：ΔG≤-20 kcal/mol，其余參數(shù)為默認(rèn)值(Enrightetal.，2003；Krüger & Rehmsmeier，2006)。為使預(yù)測結(jié)果更加可信，將2個(gè)軟件的預(yù)測結(jié)果取交集進(jìn)行后續(xù)的分析。

通過KEGG通路分析對miRNA的靶基因進(jìn)行注釋。利用KAAS (Moriyaetal.，2007)和KEGG mapper software (Kanehisaetal.，2012)對靶基因進(jìn)行KEGG注釋。

1.4 已知miRNA的差異表達(dá)

將雄性成蟲和雌性成蟲的miRNA分別進(jìn)行統(tǒng)計(jì)分析以鑒定雌雄之間的差異表達(dá)。先將miRNA的表達(dá)量歸一化，然后再計(jì)算P值。

歸一化公式：歸一化表達(dá)量(NE)=實(shí)際miRNA數(shù)量/clean reads總數(shù)×1 000 000。

P值公式：

式中，N1和x屬于相同的組，N1代表clean reads的總數(shù)，x代表歸一化的表達(dá)水平。N2和y屬于相同的組，N2代表clean reads的總數(shù)，y代表歸一化的表達(dá)水平(Zhangetal.，2013)。

2 結(jié)果

2.1 測序數(shù)據(jù)、長度分布和其他小RNA的注釋

通過測序，在雄性成蟲和雌性成蟲中分別得到17 259 315條和17 249 947條小RNA序列，剔除質(zhì)量低和被污染的序列后，分別得到12 155 616條和9 847 263條。其中，共有序列占總序列的84.08%，雄性成蟲特有的序列占總序列的8.81%，雌性成蟲特有的序列占總序列的7.11%(圖1：a)。雌雄成蟲中所測得的所有小RNA序列均主要分布于18～23 nt，在22 nt和29 nt處有2個(gè)峰值(圖1：b)。在這些序列中，將相同的序列統(tǒng)計(jì)為同一種得到唯一序列，最終在雄性成蟲和雌性成蟲中分別得到2 247 751條和2 051 813條小RNA序列。

通過對小RNA注釋，發(fā)現(xiàn)許多非miRNA存在，包括rRNA、tRNA、snRNA、snoRNA、scRNA、短小的重復(fù)序列、內(nèi)含子、外顯子等(表1)。為了防止這些序列干擾后續(xù)分析，去除注釋為非miRNA的序列。最終在雄性成蟲和雌性成蟲中分別獲得11 507 387條和9 207 534條可能為miRNA的序列，這些序列除部分序列沒有注釋外，絕大多數(shù)被注釋為miRNA。

2.2 miRNA及靶基因預(yù)測和KEGG通路分析

在雄性成蟲和雌性成蟲中分別發(fā)現(xiàn)了57種和53種已知的miRNA，去除二者共有的miRNA，共發(fā)現(xiàn)59種已知的miRNA。表達(dá)量前10的miRNA見表2，其中，miR-1、miR-10-5p、miR-8-5p、miR-31a等廣泛參與了各種組織和細(xì)胞的發(fā)育、生長和代謝等重要生理過程。

圖1 美洲大蠊雌雄成蟲相同序列的小RNA數(shù)量和長度分布Fig. 1 Venn diagram of the same sequence of small RNAs and their length distributions in female and male Periplaneta americana

組別雄性成蟲雌性成蟲唯一序列數(shù)量/條(占總小RNA的比例)全序列數(shù)量/條(占總小RNA的比例)唯一序列數(shù)量/條(占總小RNA的比例)全序列數(shù)量/條(占總小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%)內(nèi)含子1 440(0.06%)3 611(0.03%)1 958(0.10%)5 208(0.05%)重復(fù)序列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 美洲大蠊雌雄成蟲表達(dá)量前10的miRNATable 2 Top 10 highly expressed miRNAs in female and male Periplaneta americana

同時(shí)，在雄性成蟲和雌性成蟲中預(yù)測出152種和94種潛在的新miRNA，去除共有的miRNA，共發(fā)現(xiàn)152種新miRNA。這些新miRNA在miRBase 21數(shù)據(jù)庫中尚未報(bào)道，可能是美洲大蠊所特有的miRNA。

進(jìn)一步對預(yù)測的潛在的新miRNA的靶基因進(jìn)行預(yù)測，將miRanda和RNAhybrid預(yù)測結(jié)果整合后，共計(jì)在雄性成蟲和雌性成蟲中預(yù)測了5 389個(gè)和3 823個(gè)靶基因，去掉共有的，共5 822個(gè)靶基因。KEGG通路分析(圖2)顯示了這些靶基因的功能，其中，信號(hào)轉(zhuǎn)導(dǎo)通路注釋到的靶基因最多，且雄性成蟲比雌性成蟲多，說明美洲大蠊miRNA調(diào)控的靶基因主要是參與信號(hào)轉(zhuǎn)導(dǎo)通路的基因。

2.3 已知的miRNA在雌雄美洲大蠊中的差異表達(dá)

雌雄美洲大蠊成蟲差異表達(dá)的miRNA的散點(diǎn)圖顯示僅有1條顯著差異表達(dá)的miRNA：miR-750，其在雌性成蟲中的表達(dá)顯著高于雄性成蟲(圖3)。

3 討論

miRNA越來越受到研究者的關(guān)注，這是由于miRNA被發(fā)現(xiàn)在越來越多生物過程中發(fā)揮著重要作用，如發(fā)育(Reinhartetal.，2000)、細(xì)胞分化(Kawasaki & Taira，2003)、細(xì)胞增殖(Brenneckeetal.，2003)等。很多昆蟲的miRNA已被報(bào)道(Weietal.，2009；Jagadeeswaranetal.，2010；Skalskyetal.，2010；Zhangetal.，2012；Huangetal.，2014；Heetal.，2016)。本研究是首次對美洲大蠊的miRNA進(jìn)行報(bào)道，在雄性和雌性成蟲中分別發(fā)現(xiàn)了2 247 751條和2 051 813條miRNA序列。在miRNA序列長度分布上，顯示出2個(gè)峰值，在21～23 nt的峰代表的是miRNA，在28～30 nt的峰代表的可能是Piwi蛋白互作RNA(piRNA)，這是一種與miRNA完全不同的RNA沉默因子(Kawaokaetal.，2011)。這一現(xiàn)象在其他昆蟲，包括飛蛾Manducasexta(Zhangetal.，2012)、桔小實(shí)蠅Bactroceradorsalis(Huangetal.，2014)等都存在相似的結(jié)果。與德國小蠊的miRNA(Cristinoetal.，2011)相比，美洲大蠊表達(dá)量最高的前10個(gè)miRNA中，有4個(gè)與德國小蠊相同：miR-1、miR-8-5p、miR-10-5p和miR-31a。在德國小蠊中表達(dá)、在美洲大蠊中不表達(dá)的miRNA有miR-184、miR-276a、miR-279b、miR-375和miR-iab-4as-5p。然而，更多的miRNA在美洲大蠊中表達(dá)而未在德國小蠊中發(fā)現(xiàn)，如miR-133、miR-277-3p、miR-278、miR-2765-5p、miR-281-5p、miR-316和miR-9b-5p。其中，miR-1在心臟和骨骼肌祖細(xì)胞的發(fā)育上發(fā)揮重要作用(Zhao & Srivastava，2007)；miR-10-5p的靶基因Hox基因決定體軸的正確前后樣式(Lund，2010)；miR-8-5p參與眾多生物學(xué)過程，包括發(fā)育(Bellesetal.，2012)、神經(jīng)肌肉協(xié)調(diào)(Karresetal.，2007；Rubioetal.，2013)、卵及卵巢形成(Lucasetal.，2015)等；miR-31a可能在蛻皮過程中調(diào)控上皮代謝(Yuetal.，2008)；miR-276通過調(diào)控1個(gè)轉(zhuǎn)錄共活化基因調(diào)控孵化同步性(Heetal.，2016)；miR-133的其中1個(gè)靶基因是結(jié)締組織生長因子(CTGF)，其在纖維化過程中發(fā)揮重要作用(Duistersetal.，2009)；miR-281的靶基因是登革熱病毒血清2型(DENV-2)，增強(qiáng)在白紋伊蚊Aedesalbopictus中DENV-2的病毒復(fù)制(Zhouetal.，2014)。本研究結(jié)果顯示，美洲大蠊的miRNA廣泛參與了昆蟲組織和細(xì)胞的生長、發(fā)育、代謝等多種重要的調(diào)控機(jī)制。此外，美洲大蠊和德國小蠊都存在一些特有的miRNA調(diào)控的基因，顯示2種蜚蠊目昆蟲存在不同的調(diào)控機(jī)制，這或許是造成二者形態(tài)和生理差異的原因之一。

圖2 KEGG通路中美洲大蠊雌雄成蟲miRNA的靶基因數(shù)量
Fig. 2 Numbers of target genes of miRNA in female and malePeriplanetaamericanaannotated by KEGG

圖3 美洲大蠊雌雄成蟲已知miRNA的差異表達(dá)散點(diǎn)圖Fig. 3 The scatter figure of differentially expressed known miRNAs between female and male Periplaneta americana

差異表達(dá)分析顯示，美洲大蠊雌雄成蟲miRNA在組成和表達(dá)量上幾乎沒有差異，只有miR-750在雌性成蟲中高表達(dá)。miR-750在被白斑綜合癥病毒W(wǎng)SSV感染的斑節(jié)對蝦Penaeusmonodon中表達(dá)量顯著降低，其預(yù)測的靶基因是ATG14基因，但隨后的實(shí)驗(yàn)證實(shí)ATG14基因不是miR-750的靶基因(Kaewkascholkuletal.，2016)。Rebijith等(2016)通過RNAi實(shí)驗(yàn)在棕櫚薊馬Thripspalmi中證明JH候選受體——過剩氣門蛋白是miR-750的靶基因，表明miR-750可能通過調(diào)節(jié)卵黃蛋白原基因參與免疫、應(yīng)激反應(yīng)和激素信號(hào)傳導(dǎo)。卵黃蛋白原主要存在于非哺乳類性成熟的卵生雌性動(dòng)物的血淋巴、脂肪體和卵中，所以miR-750在雌性成蟲中高表達(dá)。總之，對雌雄美洲大蠊進(jìn)行了小RNA的測序，檢測出了57種已知的miRNA，并預(yù)測出152種潛在的新miRNA，但這些數(shù)據(jù)是用軟件預(yù)測出來的，這些潛在的新miRNA是否具有功能及其相應(yīng)的靶基因，還需要進(jìn)一步的驗(yàn)證。這些miRNA的發(fā)現(xiàn)為將來美洲大蠊中miRNA的研究提供了基礎(chǔ)，其表達(dá)量和功能的研究將推進(jìn)對于美洲大蠊調(diào)控網(wǎng)絡(luò)的理解，并為美洲大蠊的蟲害防治和在中藥中的應(yīng)用提供新的視角。