松树内生细菌GD2对松材线虫入侵寄主时转录组的影响

陈阳雪，赵晓佳，谈家金

南京林业大学林学院，南京 210018

松材线虫是一种危险性极高的外来有害生物，由它引起的松材线虫病是一种毁灭性的森林病害[1]。松材线虫病可使易感病松树在感染松材线虫后2～3个月内发生枯萎死亡[2]。自1905年在日本大暴发以来，该病害在世界各地被发现，截至目前该病已在欧亚和北美等多个国家或地区被发现和报道[3-4]。亚洲是松材线虫病流行危害最严重的地区，其中中国是目前遭受松材线虫病威胁最严重的国家，病害已经在中国10多个省造成了大面积松树死亡，现有6 000万hm2松林正面临着松材线虫病大流行的威胁[5]。经统计，该病害目前已在我国危害马尾松(Pimusmassoniana)、油松(P.tabulaeformis)、赤松(P.densiflora)、华山松(P.armandii)、黑松(P.thunbergii)等50多种松科松属树种和落叶松(Larixolgensis)、云杉(Piceaasperata)、冷杉(Abiesfabri)等10多种非松属树种[6]，对森林生态系统造成严重的破坏，也造成了严重的经济损失。

松材线虫病由松材线虫、寄主、传播媒介、细菌真菌和环境等多个因素构成了其复杂的病害系统，因此松材线虫的致病机制尚未完全弄清[7]。目前，国内外学者对松材线虫病的研究工作主要以松材线虫为病原展开[1，8-10]。近年来，已有研究显示松材线虫致病作用是松材线虫与其共生细菌共同致病[11]。Kawazu等[12-15]、洪英娣等[16-17]从松材线虫体表分离出细菌，并通过接种试验表明松材线虫与其伴生细菌共同导致松材线虫病。但是细菌在这种复杂病害中的特定参与特性仍然存在争议。Wang等[18]、He等[19]、Xue等[20]通过转录组测序研究了细菌对松材线虫基因表达的影响，发现细菌对松材线虫的致病相关基因有一定的影响。这些对于深入研究细菌在松材线虫病中的作用提供了重要的实验依据。菌株GD2是课题组前期从马尾松体内分离的内生细菌，初步研究显示其在松材线虫致病过程中发挥着重要作用[21]。本研究拟进一步探讨在松材线虫侵入寄主建立适应性的过程中，菌株GD2对线虫转录组的影响，这对全面认识松材线虫的致病机理和松材线虫-细菌微生态系统具有重要的理论意义。

1 材料与方法

1.1 材料与处理

供试菌株为马尾松内生细菌GD2。供试松材线虫AA3，分离自安徽的感病黄山松(P.taiwanensis)。虫株由南京林业大学森林病理实验室提供。供试松苗为生长较一致的 2 年生马尾松(中国，广西)，由江苏沂北园林苗木有限公司提供。

将松材线虫 AA3 接种到长满新鲜灰葡萄孢的 PDA 平板上，置于 25 ℃ 培养6 d 后，采用贝尔曼漏斗法分离。将线虫液收集到无菌的 15 mL 离心管中，3 000 r/min 离心 3 min,弃上清，用 0.1%硫酸链霉素处理 15 min，用无菌水冲洗离心3次。将松材线虫调成浓度为 15 000 条/mL 的虫悬液。

用接种环挑取1环菌株GD2的新鲜单菌落，接入含有20 mL NB培养液的锥形瓶中，27 ℃、200 r/min振荡摇培12 h，置入无菌的 50 mL离心管中，6 500 r/min 离心 3 min，弃上清，无菌去离子水洗 3 次。将菌体用无菌水重悬制成菌悬液(浓度为5×107CFU/mL)。

线虫 AA3 与菌株 GD2 预处理后，采用皮接法[22]接种到2年生马尾松。设3个处理，1个对照，每处理接种13棵马尾松：(1)AA3+GD2(T)：接种 1 mL 虫悬液和 1 mL 菌悬液；(2)AA3(CK)：接种 1 mL 虫悬液和 1 mL无菌去离子水；(3)菌株 GD2：接种1 mL 菌悬液和 1 mL无菌去离子水；(4)AG：接无菌去离子水共 2 mL。放于温室 25 ℃培养，24 h后将接种无菌水(AG)、AA3(CK)、AA3+GD2(T)的马尾松接种点部位剪下，-80 ℃保存备用。每个处理各3个生物学重复。剩余的马尾松每 7 d 观察接种寄主的感病症状，统计感病率和感病指数。根据病害分级标准和参照Yu等[23]方法计算感病指数。马尾松的感病严重程度分为5个级别，Ⅰ级：植物针叶绿色，健康生长，代表值0；Ⅱ级：0～25%针叶褪绿变黄，代表值1；Ⅲ级：25%～50% 针叶褪绿变黄，枝头弯曲，代表值2；Ⅳ级：50%～70% 针叶褪绿变黄，枝头下垂，代表值3；Ⅴ级：75%～100% 针叶失绿褐变，植株枯萎死亡，代表值4。计算公式如下：

感病率=(感病株数/接种株树) × 100% 感病指数= [∑(各级病树株数 × 代表值)/ (总株数 × 发病最高级代表值)] × 100

1.2 总RNA提取、文库构建及测序

将前面-80 ℃保存备用的样品采用OminPlant RNA kit(DNase I)(CWBIO)法提取各处理组的总RNA，用 1% 的琼脂糖凝胶电泳检测，2100 Bioanalyser (Agilent) 、ND-2000 (NanoDrop Technologies)方法检测RNA样品的质量。RNA文库的建立采用TruSeqTM RNA sample preparation Kit( Illumina,San Diego,CA)试剂盒。使用IlluminaNovaSeq 6000测序平台进行高通量测序(上海美吉生物医药科技有限公司)。数据已提交NCBI，登录号为PRJNA690180。

1.3 数据处理

首先将测得马尾松(AG组)的Raw reads先进行过滤处理得到 Clean reads，再用CK组和T组的Raw reads直接比对马尾松的Clean reads，去除比对上的reads，将剩下的reads经过滤后比对到参考基因组松材线虫基因组(参考基因组版本：GCA_000231135.1;参考基因组来源：https://parasite.wormbase.org/Bursaphelenchus_xylophilus_prjea64437/Info/Index)进行序列比对。

1.4 生物信息学分析

使用FPKM(Fragments Per Kilobases per Millionreads)软件进行定量计算获得基因的表达量，再用软件DESeq2进行差异基因分析，并用BH(fdr correction with Benjamini/Hochberg)进行多重检验。将P<0.05 且 |log2FC|≥0.585设为条件筛选差异表达基因并进行GO (Gene Ontology,http://www.geneontology.org/)、COG (clusters of orthologous groups,http://www.ncbi.nlm.nih.gov/research/cog-project/)以及KEGG (Kyoto Encyclopedia of Genes and Genomes，http://www.genome.jp/kegg/)分析。

1.5 实时荧光定量 PCR 验证差异基因

选取8个差异表达基因(5个上调基因、3个下调基因)，利用SIGMA在线设计软件(http://www.oligoarchitect.com/LoginServlet)为 8个差异表达基因设计引物(表1)，以松材线虫肌动蛋白基因Actin(GenBank EU100952)为内参基因，每个样品3个技术学重复，3个生物学重复，并采用 2-ΔΔCt算法分析数据。

2 结果与分析

2.1 菌株GD2影响松材线虫病的发生

温室接种试验表明，菌株GD2增强了松材线虫在马尾松中的致病力。AA3 处理组、GD2+AA3混合处理组都使 2 年生马尾松盆栽苗出现感病症状，并且所致萎蔫症状一致，均为针叶先褪色、黄化，进而枯萎变成红褐色；而 GD2 处理组和对照无菌水处理组均保持健康，马尾松松针绿色。接种后21 d，接种 AA3 的马尾松感病率为20%，感病指数为5，而混合接种AA3和菌株GD2的马尾松感病率达到40%，感病指数为10；接种后 42 d， AA3处理组的感病率为 60%，感病指数 52.5，而GD2+AA3 处理组的感病率为 70%，感病指数为55(表2)。

表1 松材线虫 8 个差异表达基因的 qRT-PCR 引物 Table 1 The qRT-PCR primers for 8 differentially expressed genes of Bursaphelenchus xylophilus

表2 接种不同处理的松材线虫的马尾松的感病率和感病指数 Table 2 Disease incidence and disease index (DI) of Pinus massoniana inoculated with different treatments of Bursaphelenchus xylophilus

2.2 差异表达基因的筛选

对含松材线虫的RNA样本进行转录组测序，去除比对到的马尾松reads,剩下的 reads经过滤后得到的clean reads质量达到 Q20 以上的碱基比例大于 80%，达到 Q30 以上的碱基比例大于85%，每个样本碱基 GC 含量较为一致，且测序碱基平均错误率均在0.1%以下，说明测序数据良好(表 3)。将过滤后的数据与rRNA数据库进行比对(表4)，再利用HISAT2 软件将质控后的原始数据，即Clean reads，与参考基因组(https://parasite.wormbase.org/Bursaphelenchus_xylophilus_prjea64437/Info/Index)比对(表5)。将P<0.05 且|log2FC|≥0.585设为条件共筛选到143个差异表达基因，其中 63 个基因表达上调，80个基因表达下调(图1)。这些差异基因与内生细菌GD2协同松材线虫侵染马尾松相关。

表3 转录组序列质量统计表 Table 3 The data statistics for RNA-sequencing

表4 过滤后数据与rRNA数据库比对统计表 Table 4 Results of comparison with rRNA database

表5 序列与松材线虫参考基因组比对结果统计 Table 5 Results of comparison with the reference genomeof Bursaphelenchus xylophilus

图1 差异基因火山图Fig.1 Differential gene volcano

2.3 差异基因 GO 分析

对差异表达基因进行 GO 功能注释与分类统计。差异基因被归类到GO 3个组分的38个功能分类中(图2)。细胞组成(cellular component)中细胞组分(cell part)和膜组分(membrane part)最多，分别有52个基因和44个基因；分子功能(molecular functions)中结合(binding)和催化活性(catalytic activity)最多，分别有60个基因和51个基因，生物代谢过程(biological processes)中细胞进程(cellular process)和代谢过程(metabolic process)最多，分别有48个基因和41个基因。此外，还检测到一些涉及先天免疫的生物过程，如免疫系统过程、行为、对刺激的反应等。

2.4 差异表达基因 COG 数据库注释分析

将差异表达基因比对到 COG 数据库当中，预测基因的功能并进行分类和统计，可分为16类(图3)。除了未知功能，所占数量最多的功能是蛋白质翻译后修饰与转运，分子伴侣(posttranslational modification,protein turnover,chaperones)，注释到此类的差异基因有11个。其次是信号转导机制(signal transduction mechanisms)、碳水化合物运输与代谢(carbohydrate transport and metabolism)和细胞内运输、分泌和囊泡运输(intracellular trafficking,secretion,and vesicular transport)。

图2 差异基因GO功能分析Fig.2 GO functional analysis of differential gene

2.5 差异基因 KEGG 功能分析

将差异表达基因进行KEGG生物通路富集分析，为挖掘细菌GD2增强松材线虫寄主适应性相关基因提供依据。全部差异表达基因共映射到106个代谢通路，图4显示了富集程度在前50的通路，其中新陈代谢(metabolism)、有机系统(organismal systems)、人类疾病(human diseases)涉及的通路最多。ECM-receptor互作(ECM-receptor interaction)、细胞凋亡(cellular senescence)、肾素-血管紧张素系统(renin-angiotensin system)和PPAR信号通路(PPAR signaling pathway)、AMPK 信号通路(AMPK signaling pathway)、谷胱甘肽代谢(glutathione metabolism)等通路主要涉及机体的免疫调节、延长寿命和生长代谢。

B:染色质结构和变化 Chromatin structure and dynamics； C:能量生产和转化 Energy production and conversion； E:氨基酸运输和代谢 Amino acid transport and metabolism； F:核苷酸运输和代谢 Nucleotide transport and metabolism； G:碳水化合物运输和代谢 Carbohydrate transport and metabolism； I:脂质运输和代谢 Lipid transport and metabolism； K:转录 Transcription； M:细胞壁/膜/包膜生物发生 Cell wall/membrane/envelope biogenesis； N:细胞运动 Cell motility； O:蛋白质翻译后修饰与转运,分子伴侣 Posttranslational modification,protein turnover,chaperones； P:无机离子转运与代谢 Inorganic ion transport and metabolism； Q:次级代谢物的生物合成、运输和代谢 Secondary metabolites biosynthesis,transport and catabolism； S:功能未知 Function unknown； T:信号转导机制 Signal transduction mechanisms； U:细胞内运输、分泌和囊泡运输 Intracellular trafficking,secretion,and vesicular transport； Z:细胞骨架 Cytoskeleton.

图4 差异基因 KEGG 功能分类Fig.4 KEGG functional classification of differential gene

2.6 qPCR验证

通过RT-qPCR对8个表达丰度不同的转录本(5个上调和3个下调)进行验证，发现基因的qPCR 验证结果与转录组中的趋势一致，说明转录组测序数据较为准确，可用于下一步研究(图5)。

A.转录组所测8 个差异表达基因 Read count 值； B. RT-qPCR 检测所得8 个差异表达基因的相对表达量。*P<0.05，**P<0.01，***P<0.001，****P<0.000 1。 A. Read count of the 8 differentially expressed genes measured in the transcriptome;B. The relative expression levels of the 8 differentially expressed genes detected by RT-qPCR. *P<0.05,**P<0.01,*** P<0.001,**** P<0.000 1.

3 讨 论

RNA测序(RNA-Seq)是一种用高通量测序技术进行测序分析的转录组测序技术，目前已成为所有物种全基因组分析的第一方法，而鉴定比较组间的差异表达基因是RNA-Seq 最典型的应用[24-25]。尽管RNA-Seq是一项正在积极开发的技术，但它提供的基因关键功能仍使之成为目前最为广泛使用的测序技术[26]。本研究通过转录组测序分析了细菌GD2对松材线虫入侵寄主马尾松24 h后的影响，但是我们发现松材线虫转录组序列与参考基因组的比对率比较低。原因可能是①没有确保每棵马尾松都准确接种到15 000条松材线虫；②在松材线虫接种到马尾松24 h 后，没有将松材线虫用贝尔曼漏斗法漏下线虫进行RNA 提取，而是直接将接种部位取下，这时候已经有部分线虫已经扩散到松树其他地方，再加上马尾松和松材线虫混样提取RNA，所以导致提取的RNA 里马尾松的RNA 占比很大，只有较少的为松材线虫的 RNA；③由于混样提取RNA，为了剔除马尾松的同源序列，我们将测得的松材线虫的reads 与马尾松的Clean reads 进行了比对，去除了比对上的reads，将剩下的reads 进行过滤后得到的松材线虫的Clean reads 再与参考基因组比对。所以可能是松材线虫的Clean reads 的基因覆盖度比较低，导致了松材线虫转录组序列与参考基因组的比对率比较低。

通过本次转录测序，从中筛选出143个差异表达显著的基因，并对这些差异基因进行 GO、COG以及 KEGG 分析。

从GO功能注释可以看出，差异表达基因主要集中在代谢过程、结合、催化活性、膜等与药物的吸收、运输、代谢有关的功能，这些基因可能参与了外源有害物质(马尾松抗病反应中产生的植保素)的解毒作用。此外，涉及的GO功能免疫和行为是动物保护自己免受伤害的2种机制[27]。这可能是松材线虫应对寄主的防御反应所产生的。

COG功能注释结果表明，细菌GD2与松材线虫混合接种马尾松后，GD2帮助线虫增强信号转导机制来感受寄主的防御反应，然后线虫通过改变自身的能量运输与代谢以及蛋白质翻译后修饰与转运等方面来加强对寄主的适应性。

对差异基因富集的KEGG 通路进行分析，发现差异基因集中于视黄醇新陈代谢、酪氨酸代谢、果糖和甘露糖的代谢、鞘脂类代谢等代谢通路。可见细菌GD2对松材线虫的生长代谢产生了影响。此外，一些差异基因如BXY_0706100涉及的PI3K-Akt信号通路在与细菌互作及对抗逆境胁迫中发挥着重要作用[28]，涉及的胰岛素信号通路调节固有免疫[29]；BXY_1558900涉及的ECM-receptor互作通路被认为是免疫相关通路[30]；BXY_1767700涉及的细胞色素P450对外源物质的代谢通路可催化杀虫剂、植物次生代谢物等外源物质的代谢，使生物体适应其生存环境或产生抗逆性[31]。

本研究通过 RNA-Seq，初步揭示了接种马尾松后，马尾松内生细菌GD2对松材线虫转录组的影响，发现菌株GD2影响着松材线虫的143个基因，这些DEGs显著富集于一些与代谢、免疫相关的GO、COG类别和KEGG途径，这些与松材线虫的生命活动密切相关。推测菌株GD2在松材线虫抵御寄主的防御反应时增强了线虫的信号传导机制，刺激了线虫的免疫相关反应，提高了线虫的能量运输和对寄主所产生的次生代谢物代谢的能力，从而加重了松材线虫病的发生。这为进一步从分子层面弄清细菌增强松材线虫对寄主的适应性机制奠定了基础。