携带与未携带松材线虫的松墨天牛转录组差异表达基因分析

李鸣宇, 周 娇, 王海香, 赵莉蔺,*

(1. 山西农业大学林学院, 山西太谷 030801;2. 中国科学院动物研究所, 农业虫害鼠害综合治理研究国家重点实验室, 北京 100101)

松材线虫Bursaphelenchusxylophilus属于侧尾腺纲(Secernentea)滑刃目(Aphelenchida)滑刃科(Aphelenchoididae)伞滑刃属Bursaphelenchus(Nickleetal., 1981)，是一种林业重大外来入侵害虫，造成的松材线虫病，即松树的“癌症”，具有传染快致死快，防治困难，一经感染便无法治愈等特点(Mamiya, 1972, 1983; 宁眺等, 2005)。松材线虫原产于北美洲，于20世纪初传入亚洲，近些年来又扩散到了欧洲国家，除北美洲原发地以外，对许多国家的森林生态系统造成了严重的破坏。松材线虫于1982年首次在我国南京中山陵发现，而后不断扩散，快速蔓延至江苏、浙江、安徽、广州和山东等53个县市区，累计死亡松树3 500万株，近年来疫情逐步加剧，向中国北部蔓延(Yanoetal., 1913; 孙永春, 1982; Futai, 2008; Nakamura-Matorietal., 2008; Shin, 2008; Karnkowski and Sahajdak， 2010)。目前，松材线虫已成为我国第一大林业害虫，也是我国最为严重的林业外来有害生物，给病害发生区带来了严重的损失(张星耀和骆有庆, 2003)。随着国际贸易的全球化及我国加入世贸组织，松材线虫由其他疫区传入我国非疫区的可能性不断加大，风险不断加剧，而且由于国内各种基建项目的施行及国内经济贸易活动的活跃，松材线虫由我国国内疫区传入非疫区的可能性也在逐渐加大。因此，松材线虫在我国仍具有十分广泛的定殖及扩散的可能(谢丙炎等, 2009; 魏亚楠, 2015)。

在中国和日本，松材线虫的感染周期与松墨天牛Monchamusalternatus的羽化过程同步，并且没有松墨天牛的携带，松材线虫便无法被传播。松材线虫与松墨天牛间存在特殊的共生生活史。松材线虫生活史包含繁殖和扩散两个周期，繁殖周期主要扩大种群密度，致死松树；扩散周期则主要通过媒介昆虫转移至新寄主。当松木内松材线虫大量繁殖，其对松树的危害加剧，造成树势衰弱甚至枯死，从而吸引松墨天牛前来产卵。秋末冬初，此时树势已不利于繁殖型松材线虫生长，繁殖型2龄幼虫向扩散型3龄幼虫(LⅢ)转变。此时，松墨天牛孵化、幼虫蜕皮；翌年早春，树体内松墨天牛进入蛹期，松材线虫扩散型3龄幼虫(LⅢ)向松墨天牛蛹室聚集，松树蒸腾作用降低、针叶萎蔫黄化。松材线虫扩散型3龄幼虫(LⅢ)向扩散型4龄幼虫(LⅣ)的转变与松墨天牛蛹的发育期相一致，翌年春天，待松墨天牛从死树中羽化，可耐干旱及恶劣环境，也适合通过媒介昆虫传播的松材线虫扩散型4龄幼虫(LⅣ)同步形成，并通过松墨天牛气管转移至新的健康松树，每头松墨天牛平均携带量约18 000条，最多可达289 000条(Mamiya, 1983)。

松材线虫与松墨天牛的生活史显示，两者间特殊的共生关系构成了松材线虫病形成的生物学基础。松材线虫扩散型4龄幼虫(LⅣ)通过其关键媒介松墨天牛的气管被携带传播，而松墨天牛与松材线虫转型扩散发育进程的一致性是松材线虫扩散表型增强的关键。松墨天牛在羽化过程中释放大量C15-C18脂肪酸乙酯，启动松材线虫Insulin和DA/DAF-12信号通路，改变发育途径，进而诱导扩散型3龄幼虫(LⅢ)形成扩散型4龄幼虫(LⅣ)(Zhaoetal., 2013)。由LⅢ分泌的蛔甙能够诱导松墨天牛中的蜕皮激素生成并上调蜕皮激素相关基因表达，进而促进松墨天牛蛹化。一旦天牛进入成虫阶段，天牛分泌蛔甙并吸引LⅣ，从而促进松材线虫进入松墨天牛气管并被转运至新松树上 (Zhaoetal., 2016)。

总之，深入研究松墨天牛与松材线虫“协同互作”关系，对解释松材线虫扩散型4龄幼虫(LⅣ)被松墨天牛携带机制具有重要意义，对控制松墨天牛携带松材线虫的传播具有积极的理论和实际意义。因此，我们利用高通量Illumina HiSeqTM2000测序技术、通过组装、功能注释和分析，全面探讨松墨天牛携带松材线虫后各个组织的转录组信息，挖掘松墨天牛应对松材线虫的压力响应基因，从分子水平上提供松墨天牛携带松材线虫后维持稳态的基因信息，也为深入研究松墨天牛能够携带松材线虫的机制奠定基础。

1 材料与方法

1.1 供试昆虫

携带松材线虫的松墨天牛采集自浙江富阳(30°07′N, 119°95′E)死亡的松树上，对照组所用的未携带任何线虫的松墨天牛收集自实验室种群。实验室人工饲养条件：幼虫放在装有木屑饲料的10 cm的指型管中，每7 d换一次饲料；成虫置于塑料盒中，每7 d换一次新鲜松枝，以上培养过程均置于人工气候箱中，饲养温度为20～25℃，相对湿度60%～70%，光周期12L∶12D。所有天牛在实验前均解剖出气管并检测携带线虫情况及数量后，再收集组织用于文库构建。所选携带松材线虫的松墨天牛样本，体内携带松材线虫的数量为1 000～3 000头。

1.2 文库构建

收集未携带松材线虫的松墨天牛成虫表皮(EP-CK)、气管(Tr-CK)、脂肪体(Fb-CK)以及携带松材线虫的松墨天牛表皮(EP-Bx)、脂肪体(Fb-Bx)和气管(Tr-Bx)，用PBS洗涤，随后使用Trizol(Thermo Fisher Scientific)从组织中分离总RNA，并在使用前储存在-80℃冰箱。共构建6个基因文库，每个文库构建使用5头天牛的组织合并为一个样品。剩余的天牛尸体浸入PBS中统计线虫的数量。用ND-1000分光光度计(NanoDrop Technologies, Wilmington, DE, 美国)测定RNA样品浓度和纯度。用Dynabeads mRNA纯化试剂盒(Ambion; Thermo Fisher Scientific, 美国)纯化mRNA。RNA-Seq文库制备试剂盒(Compatible with Illumina, San Diego, CA, 美国，Illumina Whole-Genome Gene Expression BeadChip)用于转录组文库构建，使用寡核苷酸(dT)磁珠(Thermo Fisher Scientific, 美国)富集含Poly A的mRNA，然后对其进行片段化，反转录合成cDNA。使用MiniElute凝胶提取纯化试剂盒320-420 bp产物(Qiagen, Hilden, 德国)，测序是在中国科学院基因组研究所的Illumina HiSeq 2000测序平台上完成。

1.3 转录组拼接及注释

在组装之前，去除低质量的原始reads，过滤去除微生物(细菌/病毒/真菌)和线虫污染，获得的clean reads用于后续拼接。短reads在Trinity(v.2013)中进行Denovo组装 (Grabherretal., 2011)。组装的序列进一步通过CD-HIT软件(Li and Godzik, 2006) 进行聚类，以精简冗余序列并人工校对。为了使转录本尽可能全面并避免增加组装的复杂性，将6个库合并组装，使用Burrows-Wheeler Aligner(BWA)软件(Li and Durbin, 2009)将测序读数重新整合到unigenes上，评估组装序列的准确性。

使用BlastX软件，将拼接转录本与NCBI(http:∥www.ncbi.nlm.nih.gov/)非冗余蛋白数据库(NR)和Swiss-Prot数据做比对，设置E-value<10-5为阈值，以第一匹配序列的信息作为注释结果。使用Trinity软件中的RSEM评估基因丰度水平，同时使用具有P<0.001的DEG-Seq 软件鉴定差异表达的基因(DEG)，并使用以下规则标记。将FDR(false discovery rate)作为筛选的关键指标，将|log2Ratio|≥2且Q-value≤0.05的基因作为差异表达基因，计算出上调与下调的基因数量。用KEGG(http:∥www.kegg.jp)和GO(http:∥wego.genomics.org.cn/)数据库(Ogataetal., 1999)进行功能富集分析。使用R软件Pheatmap包(http:∥www.r-project.org/; R Foundation for Statistical Computing, Wien, 澳大利亚)绘制热图。本研究测得的转录组数据已上传美国国家生物技术信息中心(NCBI)SRA数据库，BioProject PRJNA374773。

1.4 实时荧光定量PCR验证差异表达基因的表达

采用实时荧光定量PCR的方法检测1.3节筛选的在携带与未携带松材线虫的松墨天牛间差异表达的基因在松墨天牛不同组织中的相对表达量。这些基因包括气管中表达上调的核酸内切酶V(endonuclease V)基因(ENV)、细胞周期蛋白依赖性激酶1(cyclin-dependent kinase 1)基因(CDK1)、热激蛋白70-C(heat shock protein 70-C)基因(HSP70-C)、热激蛋白70(heat shock protein 70)基因(HSP70)、双重氧化酶(dual oxidase)基因(DUO)、多腺苷酸结合蛋白(PABPN1 protein)基因(PABPN1)、胰岛素样生长因子结合蛋白复合物酸不稳定亚基(insulin-like growth factor-binding protein complex acid labile subunit)基因(IGFP)、热激蛋白75(heat shock protein 75)基因(HSP75)等。使用Primer Premier 6设计内参基因Actin引物序列和目的基因引物序列，以1.2节合成的cDNA为模板，引物信息详见表1。所需引物由北京六合华大基因科技有限公司合成。PCR反应体系: 2×Talent qPCR Premix 10 μL, 上下游引物(10 μmol/L)各0.8 μL, cDNA 2 μL, 超纯水补至20 μL。反应在ABI 7300Plus荧光定量PCR仪上进行。反应程序: 95℃预变性30 s; 95℃变性5 s, 55℃退火30 s, 72℃延伸30 s, 40次循环；接着进行溶解曲线分析。表达结果依据Ct值, 采用2-△△Ct法计算相对表达量。使用SPSS19.0软件对所得数据进行统计分析。

2 结果

2.1 RNA-Seq数据质控与评估

未携带松材线虫的松墨天牛气管、表皮及脂肪体的Illumina测序共得到60 734 045条raw reads，去除不可用数据509 836条，最后得到60 224 209条clean reads。携带松材线虫的松墨天牛气管、表皮及脂肪体的Illumina测序共得到66 296 151条raw reads，去除不可用数据715 587条，最后得到65 580 564条clean reads(表2)。上述结果说明可用reads在99%以上，本研究测得的转录组数据质量优，可用于进一步分析。

对转录本进行聚类和组装分析共得到55 059条unigenes，最大长度2 845 bp，最小长度为300 bp，平均长度为1 536 bp，中位数长度为906 bp，N50长度为2 516 bp。超过1 kb的unigenes有12 394条，占22.5%，超过2 kb的unigenes有13 219条，占24.0%。这些数据说明组装效果理想。

表1 引物信息Table 1 Primer information

表2 松墨天牛成虫各组织转录组数据评估Table 2 Analysis of the transcriptome data of adult tissues of Monochamus alternatus

2.2 携带和未携带松材线虫的松墨天牛成虫组织转录组中的差异表达基因

以未携带松材线虫的松墨天牛成虫各组织的转录组数据为对照，携带松材线虫天牛成虫的气管中有差异表达基因7 089个，其中上调基因4 526个，下调基因2 563个。表皮中有差异表达基因3 923个，其中上调基因2 575个，下调基因1 348个。脂肪体中有差异表达基因5 485个，其中上调基因1 748个，下调基因3 737个。气管与表皮中共有的差异表达基因1 871个；其中二者全部上调的基因849个，全部下调的基因425个；气管中上调而表皮中下调的基因265个，气管中下调而表皮中上调的基因332个。表皮和脂肪体中共有的差异表达基因1 587个；其中二者全部上调的基因480个，全部下调的基因676个；表皮中上调而脂肪体中下调的基因214个，表皮中下调而脂肪体中上调的基因217个。气管与脂肪体中共有的差异表达基因1 683个；其中二者全部上调的基因374个，全部下调的基因545个；气管中上调而脂肪体中下调的基因481个，气管中下调而脂肪体中上调的基因283个。气管、表皮与脂肪体中共有的差异表达基因1 446个；三者全部上调的基因为488个，全部下调的基因为367个；气管与表皮中上调而脂肪体中下调的基因为103个，气管与脂肪体中上调表达而表皮中下调表达的基因为78个；表皮与脂肪体中上调而气管中下调的基因为89个；气管中上调而表皮和脂肪体中下调的基因为163个，表皮中上调而气管与脂肪体中下调的基因为58个；脂肪体中上调而气管与表皮中下调的基因为100个(图1)。以上结果说明，松墨天牛气管相对于表皮及脂肪体是主要响应松材线虫的组织。

2.3 差异表达基因的功能注释

图1 携带和未携带松材线虫的松墨天牛成虫转录组差异表达基因韦恩图Fig. 1 Venn diagram of differentially expressed genes in Monochamus alternatus adults carrying and without carrying Bursaphelenchus xylophilusTr: 气管Trachea; Ep: 表皮Epidermis; Fb: 脂肪体Fat body. 向上和向下的箭头分别表示上调和下调基因。Up and down arrows indicate up-regulated and down-regulated genes, respectively.

利用BLASTX程序将55 059条unigenes序列与NR(14 234条unigenes), GO(4 022条unigenes)和KEGG(6 098条unigenes)数据库进行比对，成功注释到24 354条unigenes，占44.23%。

KEGG数据库注释到6 098条基因，占11.1%。属于气管与表皮中的共有2 696条基因；经过比对，相较于未携带松材线虫的松墨天牛，携带松材线虫的松墨天牛气管中上调表达基因与下调表达基因最多的均为复制和修复(replication and repair)相关的，分别为114和152条；表皮组织中上调表达基因最多的为折叠、分类和降解(folding, sorting and degradation)相关的，为54条，下调表达基因最多的为信号分子与相互作用(signaling molecules and interaction)相关的，为44条(图2)。

GO数据库注释到4 022条基因，占7.3%；相较于未携带松材线虫的松墨天牛，携带松材线虫的松墨天牛表皮与气管中上调的基因(图3)中属于细胞组分(cellular component)分类的主要有细胞(cell)、细胞外区域组分(cell part)和细胞内区域组分(intracellular part)等；属于分子功能(molecular function)分类的主要为催化活性(catalytic activity)、水解酶活性(hydrolase activity)和转移酶活性(transferase activity)等；属于生物学过程(biological process)分类的有代谢过程(metabolic process)、含氮化合物代谢过程(nitrogen compound metabolic process)和有机物代谢过程(organic substance metabolic process)等。下调的基因(图4)中属于细胞组分分类的主要有细胞(cell)、细胞外区域组分(cell part)和细胞内液(intracellular fluid)等；属于分子功能分类的主要为催化活性(catalytic activity)、结合(binding)和杂环基因的结合(heterocyclic compound binding)等；属于生物学过程分类的有代谢过程(metabolic process)、含氮化合物代谢过程(nitrogen compound metabolic process)和有机物代谢过程(organic substance metabolic process)等。

图3 在携带松材线虫松墨天牛成虫表皮与气管转录组中上调基因的GO注释分布Fig. 3 GO annotation distribution of up-regulated genes in epidermis and trachea transcriptomes of Monochamus alternatus adults carrying Bursaphelenchus xylophilusTr-u: 气管中上调表达的基因Up-regulated genes in trachea; Ep-u: 表皮中上调表达的基因Up-regulated genes in epidermis.

图4 在携带松材线虫松墨天牛成虫表皮与气管转录组中下调基因的GO注释分布Fig. 4 GO annotation distribution of down-regulated genes in epidermis and trachea transcriptomes of Monochamus alternatus adults carrying Bursaphelenchus xylophilusTr-d: 气管中下调表达的基因Down-regulated genes in trachea; Ep-d: 表皮中下调表达的基因Down-regulated genes in epidermis.

其中，GO功能富集分析结果表明，气管样品中上调的基因主要富集在分子功能中的催化活性、结合、生物学过程中的代谢过程和细胞过程(图3)；下调的基因则主要富集在细胞组分中的细胞、细胞外区域组分、细胞内区域组分以及分子功能分类中的结合杂环基因的结合、环类有机物的结合和离子作用(图4)。表皮样品中上调的基因同样主要富集在分子功能中的催化活性、结合以及生物学过程中的代谢过程、有机物代谢过程及初级代谢过程(图4)；下调的基因则富集在分子功能中的催化活性、结合以及生物过程中的代谢过程、含氮化合物代谢过程和有机物代谢过程(图4)。

2.4 松墨天牛胁迫响应相关基因的鉴定

通过GO, KEGG和NR注释分析显示，松墨天牛携带松材线虫后，复制和修复相关功能通路基因表达显著上调，热图分析进一步表明其在携带松材线虫的松墨天牛的表皮和气管中显著上调(图5)。其中，气管中上调表达的共济失调蛋白(ataxin-3-like)基因、锌指蛋白(zinc finger protein)基因、小核糖核蛋白颗粒蛋白SmE(small ribonucleoprotein particle protein SmE)基因，IGFP基因、ENV基因、CDK1基因、HSP70-C基因、HSP70基因、DUO基因、中链酰基辅酶A脱氢酶(medium-chain specific acyl-CoA dehydrogenase, mitochondri)基因、锚蛋白重复和FYVE域蛋白(ankyrin repeat and FYVE domain-containing protein 1)基因、CREB结合蛋白(CREB-binding protein)基因、参与囊泡运输的内质网(ER)蛋白 (vesicle-associated membrane protein-associated protein B)基因，HSP75基因，可能的RNA解旋酶(putative RNA helicase mov-10-B.1)基因和半胱氨酸蛋白酶抑制剂(Cystatin-C)基因等，与胁迫响应相关(表3)。经过定量PCR验证，其中CDK1基因、HSP70基因、HSP75基因、DUO基因、PABPN1基因和IGFP基因的表达量在携带松材线虫的松墨天牛成虫气管中相比未携带松材线虫的松墨天牛气管中的显著上调(图6)。

图5 胁迫响应相关基因在携带与未携带松材线虫的松墨天牛成虫组织中的表达热图Fig. 5 Heat map showing expression of genes related to stress response in the adult tissues of Monochamus alternatus carrying and without carrying Bursaphelenchus xylophilus

表3 携带和未携带松材线虫的松墨天牛成虫各组织中显著上调的压力调节相关基因Table 3 Significantly up-regulated genes related to pressure regulation in adult tissues of Monochamus alternatus carrying and without carrying Bursaphelenchus xylophilus

续表3 Table 3 continued

图6 胁迫响应相关基因在携带与未携带松材线虫的松墨天牛成虫气管中表达的实时荧光定量PCR验证Fig. 6 Expression verification of stress response-related genes in trachea of Monochamus alternatus adults carrying and without carrying Bursaphelenchus xylophilus by real-time quantitative PCRA: ENV; B: CDK1; C: HSP70-C; D: HSP70; E: DUO; F: PABPN1; G: IGFP; H: HSP75. 基因表达量以Tr-CK(对照)中的为基准，柱上星号和双星号分别表示两组间差异显著(P<0.05)和极显著(P<0.01)(独立样本t检验)。The gene expression level in Tr-CK(control) was normalized to 1. The asterisk and double asterisk above bars indicate significant difference (P<0.05) and extremely significant difference (P<0.01), respectively, between two groups (independent samples t-test).

3 讨论

目前，在全球范围内植物病原体已经带来了许多毁灭性的病害，而昆虫媒介负责传播一些世界上最致命的病原体，且某些病原需要特定的昆虫才能传播(Suzukietal., 2006)。媒介昆虫对植物病原体的携带主要集中在果蝇或半翅目昆虫的肠道和唾液腺中(Bassetetal., 2000)。昆虫对线虫的携带研究主要在昆虫对寄生性线虫和昆虫病原线虫之间展开(Castilloetal., 2011)。然而，媒介昆虫对植物病原线虫的携带，尤其是在气管系统中则鲜有报道。媒介昆虫松墨天牛能安全地携带成千上万头植物病原线虫——松材线虫而不影响其寿命，所以研究松墨天牛与松材线虫间的共生互作关系十分必要。本研究利用Illumina测序技术对携带与未携带松材线虫的松墨天牛进行转录组测序和分析，在无参考基因组转录组测序情况下，共获得55 059条unigenes。其中4 022条unigene在GO数据库中获得注释，6 098条unigene在KEGG数据库中获得注释。通过分析后发现松墨天牛携带松材线虫后，差异基因主要为松墨天牛气管中压力调节、组织与DNA修复相关通路，为进一步发掘重要的功能基因奠定基础。

通过对比携带及未携带松材线虫的松墨天牛的气管、表皮及脂肪体的转录组数据，发现除脂肪体外，携带松材线虫的松墨天牛气管与表皮组织中的差异表达基因均以上调为主，其中尤以气管中的基因表达量变化最大，高于其余两组织中的差异表达基因总数及上调基因数。对差异表达基因的GO和KEGG分析结果显示，携带松材线虫的松墨天牛气管及表皮中富集到的功能通路相一致，且表达模式也相一致。说明尽管松材线虫在天牛气管中被携带，昆虫气管系统是由外胚层细胞内陷形成，其内壁与表皮相连，与表皮构造相同，昆虫表皮及气管系统都具有一层上皮细胞且相连，作为昆虫的第一道防线抵挡外来物侵染。前期研究显示，松墨天牛对松材线虫的响应过程中，总体上松墨天牛的免疫处于压制状态，表皮和气管中一部分免疫基因上调表达，大部分免疫基因表达水平不变，松墨天牛通过活性氧调节维持二者间免疫稳态(Zhouetal., 2017, 2018)。在本研究中，我们通过进一步转录组分析，筛选到大量胁迫响应基因在携带松材线虫的松墨天牛中上调，说明松材线虫对于松墨天牛更像是一种压力因子的存在，而不是一种外来侵染物，可能是二者长期协同进化的过程中建立的特殊共生关系的体现。

根据KEGG分析及GO分析可看出，当松材线虫进入天牛气管后，表达量上调的共济失调蛋白、酪氨酸激酶受体等与DNA修复和细胞损伤修复有关，其中尤以PABPN1蛋白表达量变化最为明显。而CREB结合蛋白、锌指蛋白、PABPN1蛋白、RNA结合蛋白与真核基因表达调控以及RNA合成与代谢有关等方面的基因表达量增加；同时作为生殖系干细胞调节基因的核小核糖核蛋白颗粒SmE则与生殖活动相关(Ishigaketal., 2015; Dyson and Wright, 2016; Banerjeeetal., 2017)，而催化活性过程中的上调的FYVE域蛋白则与激素相应和应对胁迫压力方面有潜在作用(Wangetal., 2000; Yaoetal., 2016)。同时不难看出许多上调的基因均与DNA复制、损伤修复相关，如Endonuclease V与DNA的损伤修复相关；Cystatin-C基因的表达产物被认为对半胱氨酸蛋白酶的局部调节有着十分重要的作用；Putative helicase mov-10-B.1作为一种可能的RNA解旋酶，是miRNA介导的基因沉默过程所必需的；而Heat shock protein 75 kD, mitochondrial(precursor)表达形成的TNFR及其下游激酶则直接与细胞增殖调节与胁迫响应有关(Copeland, 2012)。所以我们推测，作为松材线虫的传播媒介，松墨天牛在携带松材线虫后，通过增强遗传物质与细胞组织修复及真核基因转录调控等方面的基因表达，来提高自身对于松材线虫的耐受性，抵消松材线虫寄生所带来的不利影响，重新维持了自身正常生命稳态，从而为松材线虫与松墨天牛之间形成共生关系提供了可能。

本研究从转录组测序方面探究了携带与不携带松材线虫状态下的松墨天牛气管、表皮及脂肪体3个组织间表达基因的变化，为探究二者之间的共生关系提供了基础数据支持，为松材线虫防治工作提供了新的参考思路。后续我们将对重点差异表达基因的功能进行分析，更全面深入地解析松墨天牛在携带松材线虫后所发生的生理变化的分子机制，为将来的实际应用打下良好基础。