攜帶與非攜帶松材線蟲的松墨天牛miRNA表達(dá)譜比較分析

寧 靜, 張 賓, 田浩楷,3, 柳小龍,4, 楊炳琰,3, 趙莉藺*

1河北大學(xué)生命科學(xué)學(xué)院,河北 保定 071002； 2中國科學(xué)院動(dòng)物研究所,農(nóng)業(yè)蟲害鼠害綜合治理研究國家 重點(diǎn)實(shí)驗(yàn)室,北京 100101； 3中國科學(xué)院大學(xué),北京 100049； 4山西農(nóng)業(yè)大學(xué)林學(xué)院，山西 太谷 030801

松墨天牛MonochamusalternatusHope，屬鞘翅目Coleoptera天?？艭erambycidae墨天牛屬M(fèi)onochamus，是松樹的主要蛀干害蟲，主要為害衰弱木、瀕死木和新伐木(王玲萍,2004； Tomiczek & Hoyer-Tomiczek,2008)。在亞洲，松墨天牛與林業(yè)重大外來入侵害蟲松材線蟲BursaphelenchusxylophilusSteiner et Bührer(pine wood nematode,PWN)建立了共生互作關(guān)系(Zhaoetal.,2014)。夏季，松材線蟲連續(xù)經(jīng)過卵、繁殖型1～4齡幼蟲發(fā)育為成蟲并快速繁殖。如遇到低溫、干燥或種群密度過高等不利的環(huán)境條件，松材線蟲就轉(zhuǎn)型為擴(kuò)散型3齡幼蟲LⅢ；LⅢ在媒介天牛存在的情況下，進(jìn)入天牛氣管并轉(zhuǎn)型為擴(kuò)散型4齡幼蟲(dispersal fourth stage larva,LⅣ)(Tomminenetal.,1991)。LⅣ隨著媒介天牛從樹勢(shì)衰弱的病死木中轉(zhuǎn)移到新的健康松樹上，線蟲在補(bǔ)充營養(yǎng)時(shí)離開天牛體內(nèi)，進(jìn)入新的寄主松樹中，轉(zhuǎn)型為繁殖型線蟲在松樹中不斷傳代，進(jìn)而導(dǎo)致松材線蟲病的不斷擴(kuò)散與蔓延(Dwinell,1997; Mamiya,1983; Tomiczek & Hoyer-Tomiczek,2008)。因此，從不同角度闡明松墨天牛和松材線蟲的互作機(jī)制對(duì)于理解松材線蟲病的傳播擴(kuò)散機(jī)理和發(fā)展新的防控技術(shù)至關(guān)重要。

從化學(xué)生態(tài)學(xué)的角度來說，松材線蟲的傳播主要靠松墨天牛產(chǎn)生的脂肪酸類物質(zhì)、揮發(fā)性萜烯類物質(zhì)、碳?xì)浠衔锖虲O2進(jìn)行調(diào)控(鄭雅楠等,2014)。脂肪酸類物質(zhì)可刺激線蟲聚集，揮發(fā)性萜烯類物質(zhì)、部分碳?xì)浠衔镆约癈O2對(duì)線蟲具有吸引作用(Futai,2008； Miyazakietal.,1978; Shuto & Watanabe,1987； Zhaoetal.,2007)。從分子生物學(xué)的角度來講，天牛與線蟲的協(xié)同進(jìn)化使天牛產(chǎn)生了許多不同于其他昆蟲的生物學(xué)特征。例如，松墨天牛為了適應(yīng)松材線蟲的進(jìn)入，在JNK和STAT等免疫通路上產(chǎn)生了許多特異表達(dá)的基因(Zhouetal.,2017)。然而，表觀遺傳因子在松墨天?！刹木€蟲互作中的作用鮮有報(bào)道。

MicroRNA (miRNA)是一種有重要功能的非編碼RNA，成熟miRNA約22個(gè)核苷酸(nt)。它們通過與靶基因5′非翻譯區(qū)(5′UTR)、3′非翻譯區(qū)(3′UTR)及編碼區(qū)(CDS)結(jié)合，對(duì)基因的轉(zhuǎn)錄后調(diào)控發(fā)揮非常重要的作用(Axtelletal.,2011； Bartel & Chen,2004)。在動(dòng)物體內(nèi)，miRNA主要通過“種子序列”即5′端的2～8個(gè)堿基與靶序列互補(bǔ)而實(shí)現(xiàn)對(duì)靶基因的調(diào)控(Brenneckeetal.,2005； Lewisetal.,2005)。在最新的miRBase數(shù)據(jù)庫(http:∥www.mirbase.org/,release21,October 2017)中，已經(jīng)收錄了28645 種miRNA。對(duì)昆蟲miRNA的研究表明，其參與調(diào)控細(xì)胞的增殖、分化、凋亡，以及機(jī)體的生理代謝、免疫防御等幾乎所有的生物過程(劉永平等,2013)。由此可知，對(duì)攜帶和未攜帶松材線蟲的松墨天牛miRNA進(jìn)行鑒定將有助于研究天?！€蟲互作的表觀遺傳分子調(diào)控機(jī)理。本實(shí)驗(yàn)通過構(gòu)建8個(gè)miRNA文庫，并使用illumina 高通量測(cè)序平臺(tái)進(jìn)行測(cè)序，得到松墨天牛的保守miRNA和新miRNA,比較研究8個(gè)小RNA庫中的miRNA表達(dá)譜，并預(yù)測(cè)其潛在的靶基因。由于miRNA在物種進(jìn)化過程中較為保守，也可為其他鞘翅目昆蟲miRNA的研究提供參考。

1 材料與方法

1.1 材料

實(shí)驗(yàn)所用松墨天牛成蟲于2015年6月采自浙江富陽(N30°70′、E119°90′)野外。選取攜帶松材線蟲的天牛(以PWN表示)和不攜帶松材線蟲的天牛(以CK表示)雄蟲，在無菌條件下解剖，在放有無菌水的培養(yǎng)皿中快速?zèng)_洗各組織，去掉附著在組織表面的線蟲(對(duì)沖洗過組織的無菌水進(jìn)行鏡檢，確定天牛有無攜帶線蟲)。分別完成各個(gè)體的表皮(Epidermis,Ep)、脂肪體(Fat body,Fb)、中腸(Midgut,Mg)和氣管(Trachea,Tr)的取樣工作。依次編號(hào)，并放入 EP管中，凍存于液氮，置于-80 ℃保存?zhèn)溆谩?/p>

1.2 樣本RNA提取、建庫及質(zhì)量檢測(cè)

RNA 提取使用RNeasy Micro Kit試劑盒(Qiagen，德國)。RNA 提取后，使用Nano Drop DN-1000(Nano Drop Technologies，美國)檢測(cè)樣品質(zhì)量，選取質(zhì)量好的樣本進(jìn)行建庫。參考轉(zhuǎn)錄組數(shù)據(jù)上傳至NCBIbioproject PRGNA374773。使用Agilent 2100 Bioanalyzer和ABI StepOnePlus Real-Time PCR System檢測(cè)文庫質(zhì)量和產(chǎn)量，并用illuminaHiSeq 2000系統(tǒng)進(jìn)行高通量測(cè)序，相關(guān)檢測(cè)由深圳華大基因科技服務(wù)有限公司完成。

1.3 數(shù)據(jù)過濾

去除雜質(zhì)，去除無插入片段序列、插入片段過長的序列、低質(zhì)量序列、polyA序列和小片段序列，得到clean reads。統(tǒng)計(jì)小RNA(sRNA)的序列種類(用unique表示)及序列數(shù)量(用total表示)，并統(tǒng)計(jì)小RNA的序列長度分布情況。

1.4 小RNA分類注釋

對(duì)所有小RNA與各類RNA進(jìn)行Genbank(ftp:∥ftp.ncbi.nlm.nih.gov/genbank/)比對(duì)和注釋，由于某些小RNA可能會(huì)比對(duì)上多個(gè)不同RNA的注釋結(jié)果，按照rRNAetc>known miRNA>piRNA>repeat>exon>intron的優(yōu)先級(jí)順序?qū)RNA進(jìn)行遍歷，去除重復(fù)注釋的小RNA。

1.5 已知miRNA比對(duì)

通過blast將sRNA和miRBase數(shù)據(jù)庫(http:∥www.mirbase.org/)中所有昆蟲(去除線蟲)miRNA進(jìn)行比對(duì)，鑒定出已知miRNA。

1.6 新miRNA預(yù)測(cè)

本實(shí)驗(yàn)主要是針對(duì)未注釋上任何RNA且比對(duì)上基因組外顯子反義鏈、內(nèi)含子、基因間區(qū)的小RNA，通過選用軟件mirdeep(Friedl?nderetal.,2012)篩選miRNA的生物特征而得到。miRNA初始轉(zhuǎn)錄位點(diǎn)多位于基因間隔區(qū)、內(nèi)含子以及編碼序列的反向重復(fù)序列上，其前體具有標(biāo)志性的發(fā)夾結(jié)構(gòu)，成熟體的形成是由Dicer酶的剪切而實(shí)現(xiàn)。

1.7 miRNA差異表達(dá)分析

利用ExpDiff方法對(duì)已知miRNA進(jìn)行差異表達(dá)分析，將攜帶松材線蟲的各組織miRNA庫(Ep-PWN、Fb-PWN、Mg-PWN、Tr-PWN)和未攜帶松材線蟲的miRNA庫(Ep-CK、Fb-CK、Mg-CK、Tr-CK)進(jìn)行比對(duì)。使用參數(shù)log2-ratio比較兩者共同表達(dá)的miRNA表達(dá)量差異。具體方法：將兩個(gè)樣品(CK和PWN)歸一化到同一個(gè)量級(jí)(Zhouetal.,2010)。

歸一化表達(dá)量=(某一miRNA序列數(shù)×1000000)/總序列數(shù)。

對(duì)歸一化后的數(shù)據(jù)進(jìn)行fold change和P-value統(tǒng)計(jì)(Audic & Claverie,1997)。差異表達(dá)miRNA篩選條件：P-value≤0.05且|fold change|≥1。

1.8 已知miRNA的靶基因預(yù)測(cè)

使用miRanda (Enrightetal.,2003)和targetscan (Lewisetal.,2003)進(jìn)行靶基因預(yù)測(cè),取交集或并集作為預(yù)測(cè)結(jié)果。交集都只取相同miRNA預(yù)測(cè)的相同靶基因。

1.9 差異表達(dá)的已知miRNA的靶基因功能預(yù)測(cè)

對(duì)差異表達(dá)的已知miRNA靶基因進(jìn)行GO分析，數(shù)據(jù)來自http:∥www.geneontology.org/。將靶基因與參考基因進(jìn)行比較，選擇靶基因中顯著富集的幾個(gè)GO功能條目，并篩選出與其顯著相關(guān)的生物學(xué)功能(Sherlock,2009)。用KEGG數(shù)據(jù)提供代謝通路信息進(jìn)行Pathway富集分析，當(dāng)Qvalue≤0.05時(shí)，表示差異表達(dá)基因在該通路中顯著富集(Kanehisaetal.,2008)。

2 結(jié)果與分析

2.1 sRNA測(cè)序結(jié)果

采用illuminaHiSeq 2000系統(tǒng)對(duì)4個(gè)未攜帶和4個(gè)攜帶松材線蟲的松墨天牛成蟲組織樣品(Ep、Fb、Mg、Tr)進(jìn)行高通量測(cè)序。結(jié)果(表1)表明，8個(gè)庫中未經(jīng)過濾的原始序列分別為22100500、25746674、18042623、23615366、19824853、20166204、22909595、28461495個(gè)。去除低質(zhì)量reads后Ep-CK中包含22073317個(gè)高質(zhì)量reads，過濾去除3′adapter、5′adapter和polyA reads后得到21624215個(gè)clean reads，為高質(zhì)量reads的97.97%；在Ep-PWN樣本中包含25620440個(gè)高質(zhì)量reads、24777733個(gè)clean reads，clean reads為高質(zhì)量reads的96.71%；在Fb-CK樣本中包含18019573個(gè)高質(zhì)量reads、16545231個(gè)clean reads，clean reads為高質(zhì)量reads的91.82%；在Fb-PWN樣本中包含23544786個(gè)高質(zhì)量reads、23076468個(gè)clean reads，clean reads為高質(zhì)量reads的98.01%；在Mg-CK樣本中包含19800737個(gè)高質(zhì)量reads、18946038個(gè)clean reads，clean reads為高質(zhì)量reads的95.68%；在Mg-PWN樣本中包含20138493個(gè)高質(zhì)量reads、18791765個(gè)clean reads，clean reads為高質(zhì)量reads的93.31%；在Tr-CK樣本中包含22882734個(gè)高質(zhì)量reads、20930087個(gè)clean reads，clean reads為高質(zhì)量reads的91.47%；在Tr-PWN樣本中包含28378622個(gè)高質(zhì)量reads、27990568個(gè)clean reads，clean reads為高質(zhì)量reads的98.63%。

表1 松墨天牛miRNA文庫測(cè)序質(zhì)量Table 1 Read quality in the miRNA libraries of M. alternatus

Ep-CK、Fb-CK、Mg-CK、Tr-CK分別表示未攜帶松材線蟲的松墨天牛表皮、脂肪體、中腸、氣管；Ep-PWN、Fb-PWN、Mg-PWN、Tr-PWN分別表示攜帶松材線蟲的松墨天牛表皮、脂肪體、中腸、氣管。

Ep-CK， Fb-CK， Mg-CK， Tr-CK indicate epidermis, fat body, midgut and trachea ofM.alternatusin the absence of its symbiont,B.xylophilus; Ep-PWN, Fb-PWN, Mg-PWN, Tr-PWN indicate epidermis, fat body, midgut and trachea ofM.alternatusin the presence of its symbiont,B.xylophilus.

一般而言，小RNA的長度為18～30 nt，因此可通過sRNA的長度分布圖初步判斷sRNA的類型，如miRNA一般集中在21或22 nt，siRNA集中在24 nt，piRNA集中在30 nt等。在本實(shí)驗(yàn)中，8個(gè)樣本的小RNA長度分布總體相似，主峰集中于22 nt(圖1)，說明miRNA在樣品小RNA中所占的比例較大，測(cè)序質(zhì)量較好。

圖1 松墨天牛小RNAs測(cè)序結(jié)果長度分布Fig.1 Length distribution of sequenced small RNAs from M. alternatusEp-CK、Fb-CK、Mg-CK、Tr-CK分別表示未攜帶松材線蟲的松墨天牛表皮、脂肪體、中腸、氣管；Ep-PWN、Fb-PWN、Mg-PWN、Tr-PWN分別表示攜帶松材線蟲的松墨天牛表皮、脂肪體、中腸、氣管。Ep-CK， Fb-CK， Mg-CK， Tr-CK indicate epidermis, fat body, midgut and trachea of M. alternatus in the absence of its symbiont, B. xylophilus; Ep-PWN, Fb-PWN, Mg-PWN, Tr-PWN indicate epidermis, fat body, midgut and trachea of M. alternatus in the presence of its symbiont, B. xylophilus.

如表2所示，通過SOAP和bowtie將sRNA定位到基因組上，得到各樣本小RNA文庫的總數(shù)據(jù)和unique sequences。將可比對(duì)的唯一去重復(fù)的小RNA數(shù)據(jù)與Genebank和Rfam數(shù)據(jù)庫中的非編碼RNA(包括tRNA、rRNA、snoRNA和snRNA)進(jìn)行比對(duì)。最后，剔除其他非編碼RNA的干擾，分別將Ep-CK、Ep-PWN、Fb-CK、Fb-PWN、Mg-CK、Mg-PWN、Tr-CK和Tr-PWN中的2386101、2936382、1766002、2771045、893938、1090703、1647849、2441619種unique reads用于接下來miRNA分析。結(jié)果表明，在5類小RNA庫中，miRNA的豐度最高，約60%(圖2)。

表2 比對(duì)到松墨天?；蚪M的序列數(shù)量和種類Table 2 Unique reads and total reads mapped to the genome of M. alternatus

Ep-CK、Fb-CK、Mg-CK、Tr-CK分別表示未攜帶松材線蟲的松墨天牛表皮、脂肪體、中腸、氣管；Ep-PWN、Fb-PWN、Mg-PWN、Tr-PWN分別表示攜帶松材線蟲的松墨天牛表皮、脂肪體、中腸、氣管。

Ep-CK， Fb-CK， Mg-CK， Tr-CK indicate epidermis, fat body, midgut and trachea ofM.alternatusin the absence of its symbiont,B.xylophilus; Ep-PWN, Fb-PWN, Mg-PWN, Tr-PWN indicate epidermis, fat body, midgut and trachea ofM.alternatusin the presence of its symbiont,B.xylophilus.

圖2 攜帶和未攜帶松材線蟲的松墨天牛中小RNAs的分類Fig.2 Distribution of different small RNA classes in M. alternatus in the absence and presence of its symbiont, B. xylophilusEp-CK、Fb-CK、Mg-CK、Tr-CK分別表示未攜帶松材線蟲的松墨天牛表皮、脂肪體、中腸、氣管；Ep-PWN、Fb-PWN、Mg-PWN、Tr-PWN分別表示攜帶松材線蟲的松墨天牛表皮、脂肪體、中腸、氣管。Ep-CK， Fb-CK， Mg-CK， Tr-CK indicate epidermis, fat body, midgut and trachea of M. alternatus in the absence of its symbiont, B. xylophilus; Ep-PWN, Fb-PWN, Mg-PWN, Tr-PWN indicate epidermis, fat body, midgut and trachea of M. alternatus in the presence of its symbiont, B. xylophilus.

2.2 保守miRNA鑒定

通過blast將sRNA和miRBase數(shù)據(jù)庫進(jìn)行比對(duì)，鑒定出已知miRNA，8個(gè)sRNA庫中共鑒定出941個(gè)miRNA(掃描右側(cè)二維碼，查看詳情)。其中，未攜帶松材線蟲的天牛表皮、脂肪體、中腸、氣管鑒定出的miRNA數(shù)分別為784、723、713、837個(gè)，攜帶松材線蟲的天牛表皮、脂肪體、中腸、氣管鑒定出的miRNA數(shù)分別為780、802、617、762個(gè)。未攜帶松材線蟲的天牛表皮、脂肪體、中腸、氣管的已知miRNA總表達(dá)量分別為792808、854688、832588、1787879，攜帶松材線蟲的天牛表皮、脂肪體、中腸、氣管的已知miRNA總表達(dá)量分別為409794、240149、55264、1224809。因此，從組織特異性角度來看，兩組氣管樣本的總表達(dá)量遠(yuǎn)高于其他樣本；從有無攜帶松材線蟲角度來看，未攜帶松材線蟲的樣本miRNA總表達(dá)量遠(yuǎn)高于攜帶松材線蟲的樣本。

在鑒定到的已知miRNA中，表達(dá)量前四的miRNA分別為miR-281-5p、miR-281-2-5p、miR-10-5p和miR-31-5p(表3)。

表3 8個(gè)小RNA文庫中預(yù)測(cè)出的前10位保守miRNA及其序列、表達(dá)量和同源性Table 3 Sequences, abundance and homologues of top ten predicted conserved miRNA candidates in eight small RNA libraries of M. alternatus

Ep-CK、Fb-CK、Mg-CK、Tr-CK分別表示未攜帶松材線蟲的松墨天牛表皮、脂肪體、中腸、氣管；Ep-PWN、Fb-PWN、Mg-PWN、Tr-PWN分別表示攜帶松材線蟲的松墨天牛表皮、脂肪體、中腸、氣管。

Ep-CK， Fb-CK， Mg-CK， Tr-CK indicate epidermis, fat body, midgut and trachea ofM.alternatusin the absence of its symbiont,B.xylophilus; Ep-PWN, Fb-PWN, Mg-PWN, Tr-PWN indicate epidermis, fat body, midgut and trachea ofM.alternatusin the presence of its symbiont,B.xylophilus.

2.3 新miRNA預(yù)測(cè)結(jié)果

由于松墨天?；蚪M測(cè)序尚未完成，新miRNA的預(yù)測(cè)就相對(duì)困難。表4中的新miRNA在8個(gè)新miRNA庫中表達(dá)量最高，這些新miRNA的長度為21～23 nt，最小折疊自由能為-128.9～-189.5 kJ·mol-1。

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

2.4 差異表達(dá)miRNA統(tǒng)計(jì)

使用ExpDiff法對(duì)攜帶和未攜帶松材線蟲的松墨天牛各組織進(jìn)行差異表達(dá)分析。結(jié)果表明，攜帶松材線蟲的松墨天牛表皮共同表達(dá)的miRNA中差異表達(dá)的已知miRNA數(shù)為65個(gè)，新miRNA數(shù)為34個(gè)；脂肪體中差異表達(dá)的已知miRNA數(shù)為77個(gè)，新miRNA數(shù)為42個(gè)；中腸中差異表達(dá)的已知miRNA數(shù)為85個(gè)，新miRNA數(shù)為22個(gè)；氣管中差異表達(dá)的已知miRNA數(shù)為43個(gè)，新miRNA數(shù)為47個(gè)。松墨天牛攜帶松材線蟲后，某些miRNA表達(dá)量會(huì)顯著上調(diào)，如miR-14、miR-279和miR-312等。

對(duì)以上結(jié)果進(jìn)行聚類分析(圖3)，發(fā)現(xiàn)Tr-CK與Tr-PWN聚成一類，Mg-PWN、Fb-PWN和Ep-PWN聚成一類，Mg-CK、Fb-CK和Ep-CK聚成一類。此外，對(duì)攜帶與未攜帶松材線蟲的松墨天牛miRNA差異表達(dá)進(jìn)行比較分析發(fā)現(xiàn)，攜帶松材線蟲的天牛miRNA表達(dá)量總體下調(diào)。

圖3 松墨天牛中差異表達(dá)的miRNA聚類分析Fig.3 Heat map of the expression profiles of verified miRNAs in different small RNA libraries of M. alternatus

2.5 已知miRNA的靶基因預(yù)測(cè)

從表5可看出，攜帶松材線蟲的松墨天牛表皮(Ep-PWN)的780個(gè)保守miRNA對(duì)應(yīng)81778個(gè)靶基因，脂肪體(Fb-PWN)的802個(gè)保守miRNA對(duì)應(yīng)80868個(gè)靶基因，中腸(Mg-PWN)的617個(gè)保守miRNA對(duì)應(yīng)27467個(gè)靶基因，氣管(Tr-PWN)的762個(gè)保守miRNA對(duì)應(yīng)80453個(gè)靶基因；未攜帶松材線蟲的松墨天牛表皮(Ep-CK)的784個(gè)保守miRNA對(duì)應(yīng)79578個(gè)靶基因，脂肪體(Fb-CK)的723個(gè)保守miRNA對(duì)應(yīng)79993個(gè)靶基因，中腸(Mg-CK)的713個(gè)保守miRNA對(duì)應(yīng)75788個(gè)靶基因，氣管(Tr-CK)的837個(gè)保守miRNA對(duì)應(yīng)79001個(gè)靶基因。

2.6 差異表達(dá)的已知miRNA的靶基因功能預(yù)測(cè)

從圖4中可發(fā)現(xiàn)，這些差異表達(dá)miRNA靶向的基因功能大都相近，主要在合成、代謝和免疫方面起作用。但其在參與生物過程(biological process)中的細(xì)胞過程(cellular process)的基因數(shù)量有12856個(gè)，參與代謝過程(metabolic process)的基因數(shù)量為9751個(gè);此外參與催化活性(catalytic activity)、免疫系統(tǒng)過程(immune system process)的基因也相對(duì)較多。該結(jié)果說明松墨天牛在攜帶線蟲過程中可能影響了自身的代謝和免疫反應(yīng)。

KEGG通路富集分析結(jié)果(圖5)顯示，在表皮富集差異表達(dá)miRNA較高的通路中，與微生物感染有關(guān)的通路有5個(gè)，如致病性大腸桿菌感染(pathogenicEscherichiacoliinfection)、志賀桿菌病(shigellosis)等；脂肪體富集差異表達(dá)miRNA較高的通路與表皮類似；在中腸富集差異表達(dá)miRNA較高的通路中，可找到免疫相關(guān)通路NF-κB信號(hào)通路(NF-kappa B signaling pathway)和JAK/STAT信號(hào)通路(JAK/STAT signaling pathway)；氣管富集表達(dá)miRNA較高的通路中主要與代謝相關(guān)，包括丙酮酸代謝(pyruvate metabolism)、糖酵解/糖異生(glycolysis/gluconeogenesis)、果糖和甘露糖代謝(fructose and mannose metabolism)以及次生代謝產(chǎn)物的生物合成(biosynthesis of secondary metabolites)。

表5 松墨天牛保守miRNAs的靶基因預(yù)測(cè)Table 5 Predicted target genes of the identified conserved miRNAs in the M. alternatus genome

Ep-CK、Fb-CK、Mg-CK、Tr-CK分別表示未攜帶松材線蟲的松墨天牛表皮、脂肪體、中腸、氣管；Ep-PWN、Fb-PWN、Mg-PWN、Tr-PWN分別表示攜帶松材線蟲的松墨天牛表皮、脂肪體、中腸、氣管。

Ep-CK， Fb-CK， Mg-CK， Tr-CK indicate epidermis, fat body, midgut and trachea ofM.alternatusin the absence of its symbiont,B.xylophilus; Ep-PWN, Fb-PWN, Mg-PWN, Tr-PWN indicate epidermis, fat body, midgut and trachea ofM.alternatusin the presence of its symbiont,B.xylophilus.

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

3 討論

本研究在對(duì)松墨天牛的sRNA進(jìn)行測(cè)序和分析后發(fā)現(xiàn)，8個(gè)miRNA文庫所有總表達(dá)量結(jié)果中，兩組氣管樣本(Tr-CK、Tr-PWN)的總表達(dá)量遠(yuǎn)高于其他組織樣本；未攜帶松材線蟲的樣本miRNA總表達(dá)量遠(yuǎn)高于攜帶松材線蟲的樣本。這組結(jié)果與后文的差異表達(dá)miRNA聚類分析結(jié)果相符，說明：(1)氣管相較于其他組織可能有較為特異的表達(dá)模式，這可能與氣管是攜帶松材線蟲的主要器官有關(guān)；(2)攜帶松材線蟲的天牛中下調(diào)表達(dá)的miRNA可能在攜帶松材線蟲過程中起關(guān)鍵作用；(3)在差異表達(dá)聚類熱圖中，攜帶松材線蟲的3個(gè)組織(Ep、Tr、Fb)聚成一類，未攜帶松材線蟲的3個(gè)組織(Ep、Tr、Fb)聚成一類，說明有無攜帶松材線蟲在表觀遺傳角度對(duì)松墨天牛的影響較大。

攜帶松材線蟲的松墨天牛中，某些已知miRNA表達(dá)量會(huì)顯著升高，如miR-14、miR279和miR-312等。在果蠅Drosophila中，miR-14的靶基因是drice和sugarbabe，參與細(xì)胞凋亡及物質(zhì)代謝(Vargheseetal.,2010; Xuetal.,2003)。miR-279的靶基因是nerfin-1，可調(diào)控嗅覺神經(jīng)元(Cayirliogluetal.,2008)，天?！€蟲互作過程中對(duì)某些化學(xué)信息物質(zhì)的感受可能與此有關(guān)；miR-279還可調(diào)節(jié)果蠅對(duì)CO2的感覺回路(Hartletal.,2010)，天牛釋放CO2吸引線蟲的過程也可能與此有關(guān)。miR-312的靶基因是crebA，能影響昆蟲表皮幾丁質(zhì)的發(fā)育(Burgler & Macdonald,2005)。

在生物過程GO分類中，細(xì)胞過程、代謝過程和單組織過程(single-organism process)3個(gè)分類單元富集水平最高。以往的研究表明，松墨天牛攜帶松材線蟲過程中具有重要功能的幾類物質(zhì)，如脂肪酸類物質(zhì)、碳?xì)浠衔锏暮铣杉胺纸猓鶎儆诖x過程。本研究結(jié)果表明，攜帶松材線蟲的松墨天牛與未攜帶松材線蟲的松墨天牛具有不同的代謝特點(diǎn)。在分子功能GO分類中，結(jié)合(binding)、催化活性和轉(zhuǎn)運(yùn)活性(transporter activity)的富集水平最高；這3個(gè)分類單元均在代謝過程中起重要作用，進(jìn)一步證實(shí)了GO富集的準(zhǔn)確性。在KEGG富集分析中，氣管的次級(jí)代謝產(chǎn)物的生物合成、糖酵解/糖異生、丙酮酸代謝及果糖和甘露糖代謝通路的富集水平均較高，這可能與松材線蟲通過進(jìn)入松墨天牛的氣管來完成擴(kuò)散傳播有關(guān)。在中腸富集差異表達(dá)miRNA較高的通路中，可找到免疫相關(guān)通路NF-κB和JAK/STAT信號(hào)通路，可能與松材線蟲進(jìn)入松墨天牛時(shí)導(dǎo)致的松墨天牛免疫反應(yīng)有關(guān)。總體來說，松材線蟲影響了松墨天牛的新陳代謝和免疫反應(yīng)，其機(jī)理有待日后進(jìn)一步研究。

致謝: 山西大學(xué)附屬中學(xué)校謝武韜同學(xué)協(xié)助采樣和整理數(shù)據(jù)，特此表示感謝！

劉永平, 楊靜, 劉蘊(yùn)， 2013. 昆蟲microrna的研究進(jìn)展. 昆蟲學(xué)報(bào), 56(9): 1026-1037.

王玲萍， 2004. 松墨天牛生物學(xué)特性的研究. 福建林業(yè)科技, 31(3): 23-26.

鄭雅楠, 楊忠岐, 王小藝， 2014. 松墨天牛攜帶松材線蟲傳播的化學(xué)生態(tài)學(xué)機(jī)制. 植物保護(hù), 40(1): 12-15.

AUDIC S, CLAVERIE J M, 1997. The significance of digital gene expression profiles.GenomeResearch, 7(10): 986-995.

AXTELL M J, WESTHOLM J O, LAI E C, 2011. Vive la différence: biogenesis and evolution of microRNAs in plants and animals.GenomeBiology, 12(4): 1-13.

BARTEL D P, CHEN C Z, 2004. Micromanagers of gene expression: the potentially widespread influence of metazoan microRNAs.NatureReviewsGenetics, 5(5): 396-400.

BRENNECKE J, STARK A, RUSSELL R B, COHEN S M, 2005. Principles of microRNA-target recognition.PLoSBiology, 3(3): e85.

BURGLER C, MACDONALD P M, 2005. Prediction and verification of microRNA targets by moving targets, a highly adaptable prediction method.BMCGenomics, 6(1): 88.

CAYIRLIOGLU P, KADOW I G, ZHAN X, OKAMURA K, SUH G S, GUNNING D, 2008. Hybrid neurons in a microRNA mutant are putative evolutionary intermediates in insect CO2sensory systems.Science, 319: 1256-1260.

DWINELL L D, 1997. The pinewood nematode: regulation and mitigation.AnnualReviewofPhytopathology, 35(1): 153-166.

ENRIGHT A J, JOHN B, GAUL U, TUSCHL T, SANDER C, MARKS D S, 2003. MicroRNA targets inDrosophila.GenomeBiology, 5(1): R1.

FUTAI K, 2008. Host preference ofBursaphelenchuslignicolus(Nematoda: Aphelenchoididae) andB.mucronatusshown by their aggregation to pine saps.AppliedEntomology&Zoology, 15(3): 193-197.

HARTL M, LOSCHEK L, KADOW I G, 2010. Prospero and miR-279 regulate CO2sensory circuit formation inDrosophila.JournalofNeurogenetics, 24(Suppl.1): 9-10.

KANEHISA M, ARAKI M, GOTO S, HATTORI M, HIRAKAWA M, ITOH M, 2008. Kegg for linking genomes to life and the environment.NucleicAcidsResearch, 36: 480-484.

LEWIS B P, BURGE C B, BARTEL D P, 2005. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets.Cell, 120(1): 15-20.

LEWIS B P, SHIH I H, JONES-RHOADES M W, BARTEL D P, BURGE C B, 2003. Prediction of mammalian microRNA targets.Cell, 115(7): 787-798.

MAMIYA Y, 1983. Pathology of the pine wilt disease caused byBursaphelenchusxylophilus.AnnualReviewofPhytopathology, 21(21): 201-220.

MIYAZAKI M, YAMAGUCHI A, ODA K, 1978. Behaviour ofBursaphelenchuslignicolusin response to carbon dioxide released by respiration ofMonochamusalternatuspupa.NipponRingakuKaishiJournaloftheJapaneseForestrySociety, 60: 203-208.

SHERLOCK G, 2009. Gene ontology: tool for the unification of biology.CanadianInstituteofFoodScience&TechnologyJournal, 22(4): 415.

SHUTO Y, WATANABE H, 1987. Attractants from a vector,Monochamusalternatus, for the pine wood nematode.JournaloftheAgriculturalChemicalSocietyofJapan, 51(5): 1457-1458.

TOMICZEK C, HOYER-TOMICZEK U, 2008. Biology studies relevant to the vector role ofMonochamusspecies for pine wood nematode∥MOTA M, VIEIRA P.Pinewiltdisease:Aworldwidethreattoforestecosystems. Dordrecht, The Netherlands: Springer: 215-221.

TOMMINEN J, HALIK S, BERGDAHL D R, 1991. Incubation temperature and time effects on life stages ofBursaphelenchusxylophilusin wood chips.JournalofNematology, 23(4): 477-484.

VARGHESE J, LIM S F, COHEN S M, 2010.DrosophilamiR-14 regulates insulin production and metabolism through its target, sugarbabe.Genes&Development, 24(24): 2748-2753.

XU P, VERNOOY S Y, GUO M, HAY B A, 2003. TheDrosophilamicroRNA miR-14 suppresses cell death and is required for normal fat metabolism.CurrentBiologyCb, 13(9): 790-795.

ZHAO L, MOTA M, VIEIRA P, BUTCHER R A, SUN J, 2014. Interspecific communication between pinewood nematode, its insect vector, and associated microbes.TrendsinParasitology, 30(6)： 299-308.

ZHAO L, WEI W, KANG L, SUN J, 2007. Chemotaxis of the pinewood nematode,Bursaphelenchusxylophilus, to volatiles associated with host pine,Pinusmassoniana, and its vectorMonochamusalternatus.JournalofChemicalEcology, 33(6): 1207-1216.

ZHOU L, CHEN J H, LI Z Z, LI X X, HU X D, HUANG Y, ZHAO X K, LIANG C Z, WANG Y, SUN L, SHI M, XU X H, SHEN F, CHEN M S, HAN Z J, PENG Z Y, ZHAI Q N, CHEN J, ZHANG Z F, YANG R L, YE J X, GUAN Z C, YANG H M, GUI Y T, WANG J, CAI Z M, ZHANG X Q, 2010. Integrated profiling of microRNAs and mRNAs: microRNAs located on xq27.3 associate with clear cell renal cell carcinoma.PLoSONE, 5(12): e15224.

ZHOU J, ZHAO L L, YU H Y, ZHANG W, AHMAD F, HU S N, ZHAO L L, ZOU Z, SUN J H, 2017. Comparative analysis of theMonochamusalternatusimmune system.InsectScience. DOI:10.1111/1744-7917.12453.