美洲大蠊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的研究提供了基礎，其表達量和功能的研究將推進對于美洲大蠊調控網絡的理解，并為美洲大蠊的蟲害防治和在中藥中的應用提供新的視角。