室內(nèi)飼養(yǎng)與野外采集的松墨天牛幼蟲microRNA表達譜比較分析

楊炳琰， 趙莉藺*

1中國科學(xué)院動物研究所，農(nóng)業(yè)蟲害鼠害綜合治理研究國家重點實驗室，北京100101；2中國科學(xué)院生物互作卓越創(chuàng)新中心，中國科學(xué)院大學(xué)，北京100049

松墨天牛Monochamus alternatusHope又名松褐天牛，屬鞘翅目Coleoptera天牛科Cerambycidae溝脛天牛亞科Lamiinae墨天牛屬Monochamus(宋士涵和崔錫明，1991)。寄主有馬尾松Pinus massonianaLamb.、黑松P.thunbergiiParl.、濕地松P.taedaL.、云南松P.yun-nanensisFranch.等，也會危害雪松Cedrus deodara(Roxb.)G.Don、冷杉Abies fabri(Mast.)Craib和落葉松Larix gmelinii(Rupr.)Kuzen(呂傳海等，2000)，是林業(yè)上重要的蛀干害蟲。松墨天牛幼蟲時期鉆蛀樹干和大枝條，取食韌皮部與木質(zhì)部(Wuet al.，2016)，蟲口密度較大時可導(dǎo)致松樹死亡。此外，松墨天牛還是國際公認的重要檢疫性有害生物松材線蟲Bursaphelenchus xylophilus(Steiner et Bührer)Nickel的主要傳播媒介，威脅全球森林生態(tài)系統(tǒng)(徐華潮等，2012)。自20世紀80年代初松材線蟲侵入我國境內(nèi)，松墨天牛協(xié)助松材線蟲擴散，最終暴發(fā)成災(zāi)，使松材線蟲病成為我國最危險的林業(yè)病害(駱有慶，2001)。

松墨天牛的老熟幼蟲化蛹時，會吸引大量松材線蟲向其聚集，聚集在蛹室的4齡線蟲被松墨天牛攜帶，隨著天牛補充營養(yǎng)或產(chǎn)卵侵染到健康的松樹中(楊寶君等，2003； Zhaoet al.，2014)。 此外，調(diào)節(jié)線蟲和松墨天牛發(fā)育的信息素ascarosides在野外松墨天牛老熟幼蟲體內(nèi)起重要作用(Zhaoet al.，2016)。因此，對松墨天牛幼蟲期的研究將有助于研究松材線蟲和松墨天牛的種間互作及ascarosides的合成機理，并有助于松材線蟲病的早期監(jiān)測及傳播阻斷(陳井榮和彭小忠，2014)。然而，室內(nèi)飼養(yǎng)與野外捕獲的松墨天牛在發(fā)育速率、產(chǎn)卵、壽命、體型大小等方面有較明顯的差異(吳桂康等，2019；徐金華等，2009)。發(fā)生這一現(xiàn)象的表觀遺傳響應(yīng)機制尚不明確。因此，本實驗通過研究microRNA的表達譜，以期為室內(nèi)飼養(yǎng)和野外松墨天牛幼蟲之間的代謝差異以及對ascarosides合成差異的研究提供新的思路與方法。

MicroRNAs(miRNA)由約22個核苷酸(nt)組成，是一類內(nèi)源性非編碼的單鏈小分子RNA，在動植物中廣泛表達，且結(jié)構(gòu)保守(Ambros，2004；Bartel，2004； Carrington&Ambros，2003； Lauet al.，2001)。研究發(fā)現(xiàn)，它們在多種生物過程中發(fā)揮關(guān)鍵作用，包括細胞分化、凋亡、增殖和脂肪儲存，并且可能在昆蟲對病毒抗性的產(chǎn)生中起作用(陳研等，2019；魯莎等，2019； 史宏利等，2019；Chenet al.，2010； Xuet al.，2003)。 因此，松墨天牛 miRNA的測序鑒定，對于室內(nèi)飼養(yǎng)與野外天牛發(fā)育與代謝的表觀遺傳調(diào)控研究將有很大幫助。

本實驗通過構(gòu)建4個miRNA文庫，2個南京采集的野外松墨天牛幼蟲組織樣品和2個實驗室飼養(yǎng)的松墨天牛幼蟲組織樣品，使用illuminaHiseq高通量測序，通過目標序列分類注釋，獲得樣品中包含的各組分及表達量信息。對未注釋的小RNA片段進行新miRNA的預(yù)測；對保守miRNA和新miRNA進行差異分析、靶基因預(yù)測和靶基因的GO功能注釋以及KEGG pathway注釋，旨在為松墨天牛幼蟲發(fā)育與代謝的表觀調(diào)控機制研究提供一定的參考依據(jù)。

1 材料與方法

1.1 實驗材料

野外松墨天牛是2016年10月于南京(江蘇省林科院提供)林間捕獲的5齡老熟幼蟲(NJ)。實驗室飼養(yǎng)的松墨天牛是在南京采回來的松墨天牛進行實驗室(溫度25℃，濕度35%，無光照)內(nèi)飼養(yǎng)10代的5齡老熟幼蟲(LKY)。將2種天牛分別在無菌條件下解剖，獲得表皮(epidermis，Ep)和中腸(mgmidgut， Mg)的樣本。

1.2 文庫構(gòu)建及測序

使用 RNeasy Micro Kit試劑盒(Qiagen，German)提取總RNA，并使用DNaseI消化DNA后，用Urea-PAGE純化總RNA中的sRNA，在純化后的sRNA的3′末端和5′末端分別連接特異性的接頭，再以已連接接頭的sRNA為模板，逆轉(zhuǎn)錄合成cDNA，再通過PCR擴增，從而完成整個文庫制備工作。構(gòu)建文庫，使用Agilent 2100 Bioanalyzer和ABI StepOnePlus Real-Time PCR System對其進行質(zhì)量和產(chǎn)量的檢測，因為松墨天牛的基因組數(shù)據(jù)尚未完全公布，所以本文參考光肩星天牛Anoplophora glabripennisMotsch的基因組數(shù)據(jù)(https:∥www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id＝217634)。最后用HiSeq 2000(illumina?，San Diego，CA，USA)平臺對cDNA文庫測序。一系列檢測及測序由深圳華大基因科技有限公司完成。

1.3 數(shù)據(jù)過濾及分類注釋

去掉接頭序列、污染、低拷貝、低質(zhì)量序列，得到可信的clean reads，統(tǒng)計小RNA(sRNA)的序列種類(unique)及序列數(shù)量(total)，并對sRNA做序列的長度和數(shù)量統(tǒng)計。將處理后的clean reads與各類RNA進行Genbank(http:∥ftp.ncbi.nlm.nih.gov/genbank/)和Rfam11.0(http:∥rfam.janelia.org/)比對和注釋，由于某些小RNA可能會比對上多個不同RNA的注釋結(jié)果，按照rRNAetc＞known miRNA＞piRNA＞repeat＞exon＞intron的優(yōu)先級順序?qū)RNA進行遍歷，沒有比上任何注釋信息的用unann表示(Burgeet al.，2013)。 由于 rRNAetc是由NCBI Genbank和Rfam 2個數(shù)據(jù)庫比對所得，規(guī)定這2個數(shù)據(jù)庫間的優(yōu)先級為Genbank＞Rfam11.0。

1.4 已知miRNA比對

將sRNA和miRBase數(shù)據(jù)庫(http:∥www.mirbase.org/)中所有昆蟲miRNAs進行比對，鑒定樣本中保守miRNA，并預(yù)測保守miRNA的靶基因(Ambros， 2003； Ana&Sam， 2010，2013； Griffithsjoneset al.，2006，2008； Meyerset al.，2009)。

1.5 新miRNA預(yù)測

新miRNA是指在與miRBase數(shù)據(jù)庫中所有昆蟲miRNAs進行比對時，未注釋上任何RNA且比對上基因組外顯子反義鏈、內(nèi)含子、基因間區(qū)的小RNA，通過選用軟件mitp(miRNA identification and target prediction pipeline)篩選miRNA的生物特征得到的(Friedl?nderet al.， 2008)。

1.6 miRNA差異表達分析

使用ExpDiff方法對miRNA進行差異表達分析，將野外采集的松墨天牛的miRNA庫(NJ)和室內(nèi)飼養(yǎng)的松墨天牛的miRNA庫(LKY)的各組織進行比對。對所有數(shù)據(jù)進行歸一化并對數(shù)轉(zhuǎn)化。對歸一化后的數(shù)據(jù)進行組間t檢驗，差異表達miRNA篩選條件:P-value≤0.05且fold change≥2(魏明輝等，2014)。

1.7 miRNA的靶基因預(yù)測及功能預(yù)測

使用miRanda軟件對差異表達的miRNA進行生物信息學(xué)分析，預(yù)測靶基因(Hsuet al.，2008)。將靶基因投射到 Gene Ontology數(shù)據(jù)庫(http:∥www.geneontology.org/)，進行基因GO分析(Harriset al.，2004)。將靶基因與參考基因進行比較，選擇靶基因中顯著富集的幾個GO功能條目，并篩選出與其顯著相關(guān)的生物學(xué)功能。用KEGG(Kyoto Encyclopedia Genes Genomes)數(shù)據(jù)庫(http:∥www.genome.ad.jp/kegg/)提供代謝通路信息進行Pathway富集分析，當Q-value≤0.05時，表示差異表達基因在該通路中顯著富集(Kanehisaet al.，2000，2008)。

2 結(jié)果與分析

2.1 miRNA測序的結(jié)果統(tǒng)計

應(yīng)用illumina測序平臺完成野外松墨天牛的表皮(NJ-Ep)和中腸(NJ-Mg)與室內(nèi)飼養(yǎng)的松墨天牛的表皮(LKY-Ep)和中腸(LKY-Mg)4個樣品的miRNA測序。經(jīng)分析，4個庫的 clean reads中，LKY-Ep、NJ-Ep、LKY-Mg、NJ-Mg 的 GC 含量分別為40.53593%、41.26336%、41.62256%、41.69211%，GC含量位于35%～65%之間，說明測序質(zhì)量好。測序長度集中分布在21～22 nt范圍內(nèi)，長度為22 nt的reads占總reads的比例最高。野外采集與室內(nèi)飼養(yǎng)的松墨天牛的測序片段分布具有明顯的相似性。

通過bowtie將sRNA定位到基因組上，4個庫中得到的原始序列分別為 11842208、11528297、11852104、12118488個，過濾低質(zhì)量(Q-value＜30)和接頭序列后，NJ-Ep樣品中，高質(zhì)量 reads為11692536個，clean reads為10976643個，占高質(zhì)量reads的93.88%；NJ-Mg的樣本中，高質(zhì)量 reads為11383034個，clean reads為 10750830個，占高質(zhì)量reads的94.45%；LKY-Ep的樣本中，高質(zhì)量的reads為11704237個，clean reads為11165380個，占高質(zhì)量reads的95.40%；LKY-Mg的樣本中，高質(zhì)量的reads為 11965143個，clean reads為 11579370個，占高質(zhì)量reads的96.78%。說明測序所獲得的數(shù)據(jù)量大，測序讀長質(zhì)量高。將得到的所有有效數(shù)據(jù)比對光肩星天牛的基因組，發(fā)現(xiàn)分別有25.62%、25.99%、37.06%、33.45%的 unique sRNA 在 NJ-Ep、NJ-Mg、LKY-Ep、LKY-Mg中能定位到基因組上。

通過blast將sRNA與Genbank和Rfam數(shù)據(jù)庫中的非編碼 RNA進行比對，包括核糖體 RNA(rRNA)、核仁小RNA(snoRNA)、核內(nèi)小RNA(snRNA)、轉(zhuǎn)運RNA(tRNA)。篩選和去除其中的其他非編碼RNA的干擾。其豐度分布如圖1所示，可以看出實驗室飼養(yǎng)的松墨天牛老熟幼蟲的miRNA豐度大于野外松墨天牛。

2.2 保守miRNA的鑒定

將sRNA和miRBase數(shù)據(jù)庫(http:∥www.mirbase.org/)中所有昆蟲miRNAs進行比對，構(gòu)建保守miRNA表達譜。

在 LKY-Ep、LKY-Mg、NJ-Ep、NJ-Mg 的樣本中，鑒定到保守 miRNA 分別為 16、14、13、13個，表達量分別為 1509、2296、1301、1621。 從組織特異性來看，中腸比表皮的miRNA表達量高；從室內(nèi)飼養(yǎng)與野外捕獲的角度來看，室內(nèi)飼養(yǎng)的松墨天牛miRNA表達量更高一些。

通過miRDeep2根據(jù)比對到基因組上sRNA的序列，識別保守miRNA的豐度和表達量情況。在表1中列出了在4個庫中表達量最高的10個保守miRNA。

圖1 室內(nèi)飼養(yǎng)與野外采集的松墨天牛幼蟲的小RNAs的分類圖Fig.1 Distribution of different small RNA classes in indoor rearing and wild M.alternatus larvae

表1 松墨天牛幼蟲的4個庫中保守miRNA中表達量最高10個miRNA的序列和表達量Table 1 Sequences and abundance of top ten conserved miRNA candidates in four small RNA libraries of M.alternatus larvae

2.3 新miRNA的預(yù)測及表達量結(jié)果

在 LKY-Ep、LKY-Mg、NJ-Ep、NJ-Mg 的樣本中，分別鑒定到 431、514、344、414 個新 miRNA，表達量分別為 3581、5050、2487、3217。 表 2中是從 4 個文庫中篩選出表達量前10位的新miRNA。

2.4 保守miRNA的靶基因預(yù)測

用miRanda對保守miRNA進行靶基因預(yù)測，選取能量值低于-84 kj·mol-1，且得分值大于140的進行統(tǒng)計，可得出室內(nèi)飼養(yǎng)的松墨天牛表皮(LKY-Ep)的16個保守miRNA對應(yīng)16673個靶基因，中腸(LKY-Mg)的14個保守miRNA對應(yīng)14451個靶基因。野外松墨天牛表皮(NJ-Ep)的13個保守miRNA對應(yīng)13763個靶基因，中腸(NJ-Mg)的13個保守miRNA對應(yīng)13763個靶基因。

表2 松墨天牛幼蟲的4個小RNA文庫中預(yù)測出的前10位新miRNA及其序列和表達量Table 2 Sequences and abundance of top ten predicted novel miRNA candidates in four small RNA libraries of M.alternatus larvae

2.5 差異表達的miRNA的靶基因功能預(yù)測

保守miRNA差異分析發(fā)現(xiàn)，分析到的保守miRNA無顯著差異。故對所有的miRNA進行分析，在室內(nèi)飼養(yǎng)后，中腸的miRNA表達量總體高于表皮，室內(nèi)飼養(yǎng)的松墨天牛的表達量較高于野外松墨天牛。某些microRNA表達量與野外捕獲的天牛有明顯變化，如 novel-mir-62127、novel-mir-184731、novel-mir-290819明顯上調(diào)，novel-mir-251851明顯下調(diào)(圖2)。

圖2 松墨天牛幼蟲4個庫中差異表達miRNA聚類分析Fig.2 Cluster analysis of the expression profiles of verified miRNAs in 4 small RNA libraries of M.alternatus larvae

從圖3中可以看出，這些差異表達miRNA靶向的基因功能大都相近，主要集中在氧化還原過程、細胞外組分、氧化還原酶活性3個分類單元。并且富集到很多代謝途徑，如氨基糖、幾丁質(zhì)、含氨基葡萄糖的化合物、碳水化合物代謝途徑等。說明室內(nèi)飼養(yǎng)可能影響了天牛的代謝與形態(tài)結(jié)構(gòu)，并且酶活性受到很大影響，如氧化還原酶、幾丁質(zhì)酶、水解酶、?；o酶A脫氫酶的活性等，可能是由于室內(nèi)飼養(yǎng)的恒溫條件與野外環(huán)境變化的溫度導(dǎo)致。

KEGG通路富集分析發(fā)現(xiàn)(圖4)，室內(nèi)飼養(yǎng)與野外天牛的表皮和中腸中，代謝水平發(fā)生了巨大的變化，如丙酸酯、甘油磷脂、脂肪酸、色氨酸代謝等，說明不同的生活環(huán)境對于天牛的代謝有重要影響。

3 討論

本研究通過高通量測序分析鑒定室內(nèi)飼養(yǎng)與野外捕獲松墨天牛的miRNA，得到基因序列的注釋信息。結(jié)果發(fā)現(xiàn)，從長度分布上可以看到，主要分布在21～23 nt，該長度符合miRNA長度分布模式(Lagos-Quintanaet al.， 2001； Lauet al.，2001)。 并且在4個miRNA文庫所有表達量結(jié)果中，室內(nèi)飼養(yǎng)后的松墨天牛miRNA的表達量高于野外采集天牛，2組中腸的表達量高于表皮，此結(jié)果與寧靜等(2018)的結(jié)果一致。這可能說明:從表觀遺傳的角度上，中腸相對表皮組織在冷適應(yīng)方面響應(yīng)更強烈，比表皮發(fā)揮更多的功能(Hunget al.，2000；Cristofolettiet al.，2003； Matthewset al.，2010)。 室內(nèi)飼養(yǎng)的松墨天牛miRNA上調(diào)表達可能會通過抑制一些生理過程的基因，從而影響其代謝、發(fā)育過程。此外，由于本研究用光肩星天牛的基因組作為miRNA的鑒定依據(jù)，限制了比對的范圍，導(dǎo)致鑒定到的總miRNA的數(shù)量較少。隨著高通量測序以及生物信息學(xué)的迅猛發(fā)展，待松墨天牛基因組測序完成，對于松墨天牛miRNA組學(xué)的研究將會更加簡便、全面。

圖3 松墨天牛幼蟲差異表達miRNA靶基因的GO功能分類結(jié)果Fig.3 Result of GO function classification of target genes of differential expression miRNA in M.alternatus larvae

圖4 差異表達miRNA靶基因的KEGG通路富集分析結(jié)果Fig.4 Results of the KEGG pathway enrichment analysis of target genes of differential expression miRNA

在室內(nèi)飼養(yǎng)后，有些miRNA表達量有明顯的變化，比如 novel-mir-62127、novel-mir-184731、novelmir-290819明顯上調(diào)。由于分析得到的差異表達miRNA未在昆蟲中注釋，故比對到其他物種鑒定同源miRNA。novel-mir-62127比對到的同源miRNA為osa-miR166i-5p，與耐熱性相關(guān)，可能由于實驗室與野外環(huán)境的溫度不同(Liuet al.，2017)。mir-62127等一系列的具有顯著差異的miRNA，目前在昆蟲中未有注釋，這些松墨天牛的新miRNA可能調(diào)節(jié)新的靶標基因，具有新的功能。本研究只是一個初步的數(shù)據(jù)分析，有關(guān)基因的挖掘及其表達和功能驗證還需進一步研究。

以往研究發(fā)現(xiàn)，室內(nèi)飼養(yǎng)會很大程度影響松墨天牛的代謝以及形態(tài)結(jié)構(gòu)的變化，會使其發(fā)育歷期縮短，體型變小，產(chǎn)卵量減少(徐金華等，2009)。GO功能預(yù)測發(fā)現(xiàn)，差異基因很多與糖代謝途徑相關(guān)，比如:氨基糖代謝途徑、氨基糖分解代謝過程、幾丁質(zhì)代謝途徑、幾丁質(zhì)分解代謝過程、含氨基葡萄糖的化合物代謝途徑、碳水化合物代謝途徑等代謝途徑中有較多的富集，表明實驗室飼養(yǎng)對松墨天牛的糖代謝有重要影響，這可能與信息素ascarosides的產(chǎn)生或者代謝相關(guān)。幾丁質(zhì)代謝途徑、幾丁質(zhì)分解代謝過程的富集意味著在蛻皮時間上有所調(diào)控，如野外的松墨天牛幼蟲是一年一代，從12月到來年的4月一直處于老熟幼蟲期，5月才開始變態(tài)發(fā)育進行蛻皮，但是在室內(nèi)恒溫培養(yǎng)下發(fā)育蛻皮速度加快，一年為3～4代，因此，這2個信號途徑有助于進一步探究其在蛻皮速度調(diào)控中的功能。在KEGG富集分析中也發(fā)現(xiàn)，室內(nèi)飼養(yǎng)對脂肪酸代謝、甘油磷脂代謝等脂代謝相關(guān)過程有影響，這為后續(xù)實驗提供了思路。利用室內(nèi)miRNA抑制的代謝靶標基因進一步研究及開發(fā)代謝抑制劑或干擾劑，可擾亂松墨天牛代謝過程以達到防治目標。所以本研究對野外干擾劑的開發(fā)具有一定的指導(dǎo)意義。并且在KEGG通路中發(fā)現(xiàn)染色體的富集水平較高，室內(nèi)外的溫度差異可能是對染色體產(chǎn)生影響的主要原因(Bora&Soper，1971)?？傮w來說，室內(nèi)飼養(yǎng)與野外捕獲的松墨天牛存在一定的差異，希望本研究結(jié)果為后續(xù)挖掘兩者不同表型的表觀遺傳調(diào)控機制及防治害蟲研究提供一定的幫助。

致謝:江蘇省林業(yè)科學(xué)院幫助采樣，特此表示感謝!