携带与非携带松材线虫的松墨天牛miRNA表达谱比较分析

宁 静, 张 宾, 田浩楷,3, 柳小龙,4, 杨炳琰,3, 赵莉蔺*

1河北大学生命科学学院,河北 保定 071002； 2中国科学院动物研究所,农业虫害鼠害综合治理研究国家 重点实验室,北京 100101； 3中国科学院大学,北京 100049； 4山西农业大学林学院，山西 太谷 030801

松墨天牛MonochamusalternatusHope，属鞘翅目Coleoptera天牛科Cerambycidae墨天牛属Monochamus，是松树的主要蛀干害虫，主要为害衰弱木、濒死木和新伐木(王玲萍,2004； Tomiczek & Hoyer-Tomiczek,2008)。在亚洲，松墨天牛与林业重大外来入侵害虫松材线虫BursaphelenchusxylophilusSteiner et Bührer(pine wood nematode,PWN)建立了共生互作关系(Zhaoetal.,2014)。夏季，松材线虫连续经过卵、繁殖型1～4龄幼虫发育为成虫并快速繁殖。如遇到低温、干燥或种群密度过高等不利的环境条件，松材线虫就转型为扩散型3龄幼虫LⅢ；LⅢ在媒介天牛存在的情况下，进入天牛气管并转型为扩散型4龄幼虫(dispersal fourth stage larva,LⅣ)(Tomminenetal.,1991)。LⅣ随着媒介天牛从树势衰弱的病死木中转移到新的健康松树上，线虫在补充营养时离开天牛体内，进入新的寄主松树中，转型为繁殖型线虫在松树中不断传代，进而导致松材线虫病的不断扩散与蔓延(Dwinell,1997; Mamiya,1983; Tomiczek & Hoyer-Tomiczek,2008)。因此，从不同角度阐明松墨天牛和松材线虫的互作机制对于理解松材线虫病的传播扩散机理和发展新的防控技术至关重要。

从化学生态学的角度来说，松材线虫的传播主要靠松墨天牛产生的脂肪酸类物质、挥发性萜烯类物质、碳氢化合物和CO2进行调控(郑雅楠等,2014)。脂肪酸类物质可刺激线虫聚集，挥发性萜烯类物质、部分碳氢化合物以及CO2对线虫具有吸引作用(Futai,2008； Miyazakietal.,1978; Shuto & Watanabe,1987； Zhaoetal.,2007)。从分子生物学的角度来讲，天牛与线虫的协同进化使天牛产生了许多不同于其他昆虫的生物学特征。例如，松墨天牛为了适应松材线虫的进入，在JNK和STAT等免疫通路上产生了许多特异表达的基因(Zhouetal.,2017)。然而，表观遗传因子在松墨天牛—松材线虫互作中的作用鲜有报道。

MicroRNA (miRNA)是一种有重要功能的非编码RNA，成熟miRNA约22个核苷酸(nt)。它们通过与靶基因5′非翻译区(5′UTR)、3′非翻译区(3′UTR)及编码区(CDS)结合，对基因的转录后调控发挥非常重要的作用(Axtelletal.,2011； Bartel & Chen,2004)。在动物体内，miRNA主要通过“种子序列”即5′端的2～8个碱基与靶序列互补而实现对靶基因的调控(Brenneckeetal.,2005； Lewisetal.,2005)。在最新的miRBase数据库(http:∥www.mirbase.org/,release21,October 2017)中，已经收录了28645 种miRNA。对昆虫miRNA的研究表明，其参与调控细胞的增殖、分化、凋亡，以及机体的生理代谢、免疫防御等几乎所有的生物过程(刘永平等,2013)。由此可知，对携带和未携带松材线虫的松墨天牛miRNA进行鉴定将有助于研究天牛—线虫互作的表观遗传分子调控机理。本实验通过构建8个miRNA文库，并使用illumina 高通量测序平台进行测序，得到松墨天牛的保守miRNA和新miRNA,比较研究8个小RNA库中的miRNA表达谱，并预测其潜在的靶基因。由于miRNA在物种进化过程中较为保守，也可为其他鞘翅目昆虫miRNA的研究提供参考。

1 材料与方法

1.1 材料

实验所用松墨天牛成虫于2015年6月采自浙江富阳(N30°70′、E119°90′)野外。选取携带松材线虫的天牛(以PWN表示)和不携带松材线虫的天牛(以CK表示)雄虫，在无菌条件下解剖，在放有无菌水的培养皿中快速冲洗各组织，去掉附着在组织表面的线虫(对冲洗过组织的无菌水进行镜检，确定天牛有无携带线虫)。分别完成各个体的表皮(Epidermis,Ep)、脂肪体(Fat body,Fb)、中肠(Midgut,Mg)和气管(Trachea,Tr)的取样工作。依次编号，并放入 EP管中，冻存于液氮，置于-80 ℃保存备用。

1.2 样本RNA提取、建库及质量检测

RNA 提取使用RNeasy Micro Kit试剂盒(Qiagen，德国)。RNA 提取后，使用Nano Drop DN-1000(Nano Drop Technologies，美国)检测样品质量，选取质量好的样本进行建库。参考转录组数据上传至NCBIbioproject PRGNA374773。使用Agilent 2100 Bioanalyzer和ABI StepOnePlus Real-Time PCR System检测文库质量和产量，并用illuminaHiSeq 2000系统进行高通量测序，相关检测由深圳华大基因科技服务有限公司完成。

1.3 数据过滤

去除杂质，去除无插入片段序列、插入片段过长的序列、低质量序列、polyA序列和小片段序列，得到clean reads。统计小RNA(sRNA)的序列种类(用unique表示)及序列数量(用total表示)，并统计小RNA的序列长度分布情况。

1.4 小RNA分类注释

对所有小RNA与各类RNA进行Genbank(ftp:∥ftp.ncbi.nlm.nih.gov/genbank/)比对和注释，由于某些小RNA可能会比对上多个不同RNA的注释结果，按照rRNAetc>known miRNA>piRNA>repeat>exon>intron的优先级顺序对sRNA进行遍历，去除重复注释的小RNA。

1.5 已知miRNA比对

通过blast将sRNA和miRBase数据库(http:∥www.mirbase.org/)中所有昆虫(去除线虫)miRNA进行比对，鉴定出已知miRNA。

1.6 新miRNA预测

本实验主要是针对未注释上任何RNA且比对上基因组外显子反义链、内含子、基因间区的小RNA，通过选用软件mirdeep(Friedländeretal.,2012)筛选miRNA的生物特征而得到。miRNA初始转录位点多位于基因间隔区、内含子以及编码序列的反向重复序列上，其前体具有标志性的发夹结构，成熟体的形成是由Dicer酶的剪切而实现。

1.7 miRNA差异表达分析

利用ExpDiff方法对已知miRNA进行差异表达分析，将携带松材线虫的各组织miRNA库(Ep-PWN、Fb-PWN、Mg-PWN、Tr-PWN)和未携带松材线虫的miRNA库(Ep-CK、Fb-CK、Mg-CK、Tr-CK)进行比对。使用参数log2-ratio比较两者共同表达的miRNA表达量差异。具体方法：将两个样品(CK和PWN)归一化到同一个量级(Zhouetal.,2010)。

归一化表达量=(某一miRNA序列数×1000000)/总序列数。

对归一化后的数据进行fold change和P-value统计(Audic & Claverie,1997)。差异表达miRNA筛选条件：P-value≤0.05且|fold change|≥1。

1.8 已知miRNA的靶基因预测

使用miRanda (Enrightetal.,2003)和targetscan (Lewisetal.,2003)进行靶基因预测,取交集或并集作为预测结果。交集都只取相同miRNA预测的相同靶基因。

1.9 差异表达的已知miRNA的靶基因功能预测

对差异表达的已知miRNA靶基因进行GO分析，数据来自http:∥www.geneontology.org/。将靶基因与参考基因进行比较，选择靶基因中显著富集的几个GO功能条目，并筛选出与其显著相关的生物学功能(Sherlock,2009)。用KEGG数据提供代谢通路信息进行Pathway富集分析，当Qvalue≤0.05时，表示差异表达基因在该通路中显著富集(Kanehisaetal.,2008)。

2 结果与分析

2.1 sRNA测序结果

采用illuminaHiSeq 2000系统对4个未携带和4个携带松材线虫的松墨天牛成虫组织样品(Ep、Fb、Mg、Tr)进行高通量测序。结果(表1)表明，8个库中未经过滤的原始序列分别为22100500、25746674、18042623、23615366、19824853、20166204、22909595、28461495个。去除低质量reads后Ep-CK中包含22073317个高质量reads，过滤去除3′adapter、5′adapter和polyA reads后得到21624215个clean reads，为高质量reads的97.97%；在Ep-PWN样本中包含25620440个高质量reads、24777733个clean reads，clean reads为高质量reads的96.71%；在Fb-CK样本中包含18019573个高质量reads、16545231个clean reads，clean reads为高质量reads的91.82%；在Fb-PWN样本中包含23544786个高质量reads、23076468个clean reads，clean reads为高质量reads的98.01%；在Mg-CK样本中包含19800737个高质量reads、18946038个clean reads，clean reads为高质量reads的95.68%；在Mg-PWN样本中包含20138493个高质量reads、18791765个clean reads，clean reads为高质量reads的93.31%；在Tr-CK样本中包含22882734个高质量reads、20930087个clean reads，clean reads为高质量reads的91.47%；在Tr-PWN样本中包含28378622个高质量reads、27990568个clean reads，clean reads为高质量reads的98.63%。

表1 松墨天牛miRNA文库测序质量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等。在本实验中，8个样本的小RNA长度分布总体相似，主峰集中于22 nt(图1)，说明miRNA在样品小RNA中所占的比例较大，测序质量较好。

图1 松墨天牛小RNAs测序结果长度分布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文库的总数据和unique sequences。将可比对的唯一去重复的小RNA数据与Genebank和Rfam数据库中的非编码RNA(包括tRNA、rRNA、snoRNA和snRNA)进行比对。最后，剔除其他非编码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分析。结果表明，在5类小RNA库中，miRNA的丰度最高，约60%(图2)。

表2 比对到松墨天牛基因组的序列数量和种类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数据库进行比对，鉴定出已知miRNA，8个sRNA库中共鉴定出941个miRNA(扫描右侧二维码，查看详情)。其中，未携带松材线虫的天牛表皮、脂肪体、中肠、气管鉴定出的miRNA数分别为784、723、713、837个，携带松材线虫的天牛表皮、脂肪体、中肠、气管鉴定出的miRNA数分别为780、802、617、762个。未携带松材线虫的天牛表皮、脂肪体、中肠、气管的已知miRNA总表达量分别为792808、854688、832588、1787879，携带松材线虫的天牛表皮、脂肪体、中肠、气管的已知miRNA总表达量分别为409794、240149、55264、1224809。因此，从组织特异性角度来看，两组气管样本的总表达量远高于其他样本；从有无携带松材线虫角度来看，未携带松材线虫的样本miRNA总表达量远高于携带松材线虫的样本。

在鉴定到的已知miRNA中，表达量前四的miRNA分别为miR-281-5p、miR-281-2-5p、miR-10-5p和miR-31-5p(表3)。

表3 8个小RNA文库中预测出的前10位保守miRNA及其序列、表达量和同源性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预测结果

由于松墨天牛基因组测序尚未完成，新miRNA的预测就相对困难。表4中的新miRNA在8个新miRNA库中表达量最高，这些新miRNA的长度为21～23 nt，最小折叠自由能为-128.9～-189.5 kJ·mol-1。

表4 松墨天牛的8个小RNA文库中预测出的前10位新miRNA及其序列和表达量Table 4 Sequences and abundance of top ten predicted novel miRNA candidates in eight small RNA libraries of M. alternatus

2.4 差异表达miRNA统计

使用ExpDiff法对携带和未携带松材线虫的松墨天牛各组织进行差异表达分析。结果表明，携带松材线虫的松墨天牛表皮共同表达的miRNA中差异表达的已知miRNA数为65个，新miRNA数为34个；脂肪体中差异表达的已知miRNA数为77个，新miRNA数为42个；中肠中差异表达的已知miRNA数为85个，新miRNA数为22个；气管中差异表达的已知miRNA数为43个，新miRNA数为47个。松墨天牛携带松材线虫后，某些miRNA表达量会显著上调，如miR-14、miR-279和miR-312等。

对以上结果进行聚类分析(图3)，发现Tr-CK与Tr-PWN聚成一类，Mg-PWN、Fb-PWN和Ep-PWN聚成一类，Mg-CK、Fb-CK和Ep-CK聚成一类。此外，对携带与未携带松材线虫的松墨天牛miRNA差异表达进行比较分析发现，携带松材线虫的天牛miRNA表达量总体下调。

图3 松墨天牛中差异表达的miRNA聚类分析Fig.3 Heat map of the expression profiles of verified miRNAs in different small RNA libraries of M. alternatus

2.5 已知miRNA的靶基因预测

从表5可看出，携带松材线虫的松墨天牛表皮(Ep-PWN)的780个保守miRNA对应81778个靶基因，脂肪体(Fb-PWN)的802个保守miRNA对应80868个靶基因，中肠(Mg-PWN)的617个保守miRNA对应27467个靶基因，气管(Tr-PWN)的762个保守miRNA对应80453个靶基因；未携带松材线虫的松墨天牛表皮(Ep-CK)的784个保守miRNA对应79578个靶基因，脂肪体(Fb-CK)的723个保守miRNA对应79993个靶基因，中肠(Mg-CK)的713个保守miRNA对应75788个靶基因，气管(Tr-CK)的837个保守miRNA对应79001个靶基因。

2.6 差异表达的已知miRNA的靶基因功能预测

从图4中可发现，这些差异表达miRNA靶向的基因功能大都相近，主要在合成、代谢和免疫方面起作用。但其在参与生物过程(biological process)中的细胞过程(cellular process)的基因数量有12856个，参与代谢过程(metabolic process)的基因数量为9751个;此外参与催化活性(catalytic activity)、免疫系统过程(immune system process)的基因也相对较多。该结果说明松墨天牛在携带线虫过程中可能影响了自身的代谢和免疫反应。

KEGG通路富集分析结果(图5)显示，在表皮富集差异表达miRNA较高的通路中，与微生物感染有关的通路有5个，如致病性大肠杆菌感染(pathogenicEscherichiacoliinfection)、志贺杆菌病(shigellosis)等；脂肪体富集差异表达miRNA较高的通路与表皮类似；在中肠富集差异表达miRNA较高的通路中，可找到免疫相关通路NF-κB信号通路(NF-kappa B signaling pathway)和JAK/STAT信号通路(JAK/STAT signaling pathway)；气管富集表达miRNA较高的通路中主要与代谢相关，包括丙酮酸代谢(pyruvate metabolism)、糖酵解/糖异生(glycolysis/gluconeogenesis)、果糖和甘露糖代谢(fructose and mannose metabolism)以及次生代谢产物的生物合成(biosynthesis of secondary metabolites)。

表5 松墨天牛保守miRNAs的靶基因预测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 松墨天牛差异表达miRNA靶基因的GO功能分类结果Fig.4 Result of GO function classification of target genes of differential expression miRNA in M. alternatus

3 讨论

本研究在对松墨天牛的sRNA进行测序和分析后发现，8个miRNA文库所有总表达量结果中，两组气管样本(Tr-CK、Tr-PWN)的总表达量远高于其他组织样本；未携带松材线虫的样本miRNA总表达量远高于携带松材线虫的样本。这组结果与后文的差异表达miRNA聚类分析结果相符，说明：(1)气管相较于其他组织可能有较为特异的表达模式，这可能与气管是携带松材线虫的主要器官有关；(2)携带松材线虫的天牛中下调表达的miRNA可能在携带松材线虫过程中起关键作用；(3)在差异表达聚类热图中，携带松材线虫的3个组织(Ep、Tr、Fb)聚成一类，未携带松材线虫的3个组织(Ep、Tr、Fb)聚成一类，说明有无携带松材线虫在表观遗传角度对松墨天牛的影响较大。

携带松材线虫的松墨天牛中，某些已知miRNA表达量会显著升高，如miR-14、miR279和miR-312等。在果蝇Drosophila中，miR-14的靶基因是drice和sugarbabe，参与细胞凋亡及物质代谢(Vargheseetal.,2010; Xuetal.,2003)。miR-279的靶基因是nerfin-1，可调控嗅觉神经元(Cayirliogluetal.,2008)，天牛—线虫互作过程中对某些化学信息物质的感受可能与此有关；miR-279还可调节果蝇对CO2的感觉回路(Hartletal.,2010)，天牛释放CO2吸引线虫的过程也可能与此有关。miR-312的靶基因是crebA，能影响昆虫表皮几丁质的发育(Burgler & Macdonald,2005)。

在生物过程GO分类中，细胞过程、代谢过程和单组织过程(single-organism process)3个分类单元富集水平最高。以往的研究表明，松墨天牛携带松材线虫过程中具有重要功能的几类物质，如脂肪酸类物质、碳氢化合物的合成及分解，均属于代谢过程。本研究结果表明，携带松材线虫的松墨天牛与未携带松材线虫的松墨天牛具有不同的代谢特点。在分子功能GO分类中，结合(binding)、催化活性和转运活性(transporter activity)的富集水平最高；这3个分类单元均在代谢过程中起重要作用，进一步证实了GO富集的准确性。在KEGG富集分析中，气管的次级代谢产物的生物合成、糖酵解/糖异生、丙酮酸代谢及果糖和甘露糖代谢通路的富集水平均较高，这可能与松材线虫通过进入松墨天牛的气管来完成扩散传播有关。在中肠富集差异表达miRNA较高的通路中，可找到免疫相关通路NF-κB和JAK/STAT信号通路，可能与松材线虫进入松墨天牛时导致的松墨天牛免疫反应有关。总体来说，松材线虫影响了松墨天牛的新陈代谢和免疫反应，其机理有待日后进一步研究。

致谢: 山西大学附属中学校谢武韬同学协助采样和整理数据，特此表示感谢！

刘永平, 杨静, 刘蕴， 2013. 昆虫microrna的研究进展. 昆虫学报, 56(9): 1026-1037.

王玲萍， 2004. 松墨天牛生物学特性的研究. 福建林业科技, 31(3): 23-26.

郑雅楠, 杨忠岐, 王小艺， 2014. 松墨天牛携带松材线虫传播的化学生态学机制. 植物保护, 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.