曾媛玲，吴文涛，周绍芳，闫曦蕊，陈荣春，王雪兰，王扬

（云南农业大学 植物保护学院，云南 昆明 650201）

根结线虫寄主范围广泛，对蔬菜，薯类，豆类，棉麻类，香料和咖啡等作物都构成了严重威胁［1］。根结线虫种类超过90种，种间致病性差异显著，在生产上引起的危害程度也有所差异，不同地区间的优势种群也有所不同［2-3］。目前我国记载中值得关注的根结线虫种类—象耳豆根结线虫（Meloido⁃gyne enterolobii）具有致病力强、寄主范围广、并且能克服线虫抗性基因的特性［4］，逐渐变为热带和亚热带地区最重要的根结线虫之一［5］。目前，象耳豆根结线虫在亚洲、非洲、欧洲和美洲皆有分布［6］，仅象耳豆根结线虫就可导致蔬菜高达65%的损失，高于迄今为止研究的任何其它根结线虫［7-8］。

象耳豆根结线虫的生活史与其它种类的根结线虫相似［9］，属于土壤习居性植物寄生线虫，生活周期可能因环境条件的不同而不同，大概是3~6周。土壤中的卵在适宜的环境条件下孵化成二龄幼虫（J2）在土壤中移动［10］。在土壤中的J2通过头部化感器可以感受到土壤中寄主根部分泌出的化学物质，向寄主植物根部趋向性移动，并从根尖部位侵入根组织，经过平行于根体轴的方向移动后，最终到达中柱鞘开始寄生。在寄生启始阶段，线虫将口针刺入中柱鞘分生细胞注入其食道腺分泌物，这些分泌物不仅诱导线虫头部周围的细胞膨大形成为其终生提供营养的取食位点—巨型细胞（Giant Cells，GC），并且还会使线虫周围的寄主细胞分裂加快，导致细胞壁的分解、细胞核的扩大以及细胞质组成的变化，从而使线虫侵染部分的根组织膨大，形成明显的根结［10-11］。研究表明，根结线虫食道腺细胞分泌物是决定其寄生性和致病性的关键因素，而其中的分泌蛋白则是作用主体成分［12］。随着根结线虫致病性研究的深入，研究者开始接受效应子（Effectors）的概念，Hogenhout等［13］将效应子定义为“植物病原物分泌的可以改变寄主细胞结构和功能的蛋白质或小分子物质”。线虫分泌的效应子可以促进取食位点的建立和帮助线虫在植物组织中迁移。植物线虫的头感器、表皮、排泄系统和食道腺都能产生效应子，其中食道腺是植物寄生线虫主要的合成效应子的组织，食道腺细胞由2个亚腹食道腺细胞和1个背食道腺细胞组成［14］。有研究认为食道腺细胞的活性在不同发育阶段有一定差异，例如亚腹食道腺在线虫早期侵染过程中活性最高，分泌与根的侵入、线虫迁移相关基因；而背食道腺在线虫寄生阶段的活性更高，分泌与取食位点形成和维持相关的基因［15］。根据目前已有研究可以初步判断，根结线虫分泌效应子具有降解或修饰细胞壁、参与调控寄主细胞代谢和信号途径和对抗或抑制寄主防卫反应等方面的功能［9］。

番茄是受根结线虫病危害最严重的作物之一，利用抗性品种是防控番茄根结线虫病最环保和经济的方法［15］，但目前番茄上可供应用的抗线虫品种匮乏。Mi基因家族作为番茄上目前已报道发现的唯一抗线虫来源，其中来源于秘鲁番茄（Solanum peruvianum）的Mi⁃1基因已通过胚胎拯救的方法转入栽培番茄，该基因可介导对南方根结线虫（M. incognita）、花生根结线虫（M. arenar⁃ia）、爪哇根结线虫（M. javanica）等常见根结线虫的抗性，目前生产中使用的抗根结线虫番茄品种均来源于此，相关研究表明象耳豆根结线虫能够克服Mi⁃1基因的抗性［16］，这不仅增加了番茄根结线虫病的防治难度，也使Mi基因的应用前景受到挑战。本团队前期发现Mi⁃1纯合基因型番茄品种VFNT在接种象耳豆根结线虫出现亲和性反应，而接种南方根结线虫则产生的非亲和性反应，进行组织病理学观察发现接种线虫后2~6 d是亲和性与非亲和反应分化的关键阶段［17］。可以推测象耳豆根结线虫侵染后的2~6 d是其克服Mi⁃1基因抗性相关效应子基因增强表达阶段。鉴于此，本研究通过转录组测序技术比较分析象耳豆根结线虫侵染番茄根系2~6 d关键时期基因表达情况，筛选差异表达基因并利用生物信息学软件筛选出其中具有信号肽且无跨膜结构域的象耳豆根结线虫候选效应子，旨在为解析象耳豆根结线虫的致病机理以及靶向防治根结线虫病害提供理论依据。

1 材料和方法

1.1 植物材料及线虫的处理

供试番茄材料为Mi⁃1纯合基因型栽培品种VFNT（Lycopersicon esculentum）［18］，由 云 南 农 业大学园艺学院提供。供试线虫为象耳豆根结线虫（Meloidogyne enterolobii），由本实验室采集纯化，活体保存于易感番茄（Rutgers）根系的种群。参照Abdel-Momen［19］的 方 法 进 行 植 物 培育与线虫接种。主要步骤简述如下：将番茄种子用次氯酸钠消毒洗净后放置于垫有灭菌滤纸的培养皿内置于28 ℃度温箱催芽。2~5 d后将露白的种子移栽到育苗盘于日光温室内进行育苗，待幼苗出芽长至3~4片真叶换大的灭菌花盆进行移栽培育，育苗基质为腐殖土∶壤土∶石英砂（3∶1∶1）比例混合，使用前高温灭菌（121 ℃，60 min），保证无其它微生物和线虫。番茄苗长出4~5片真叶时便可用于接种根结线虫。

次氯酸钠分散法［20］分离从用于扩繁的番茄根系上的根结线虫卵，置于25 ℃度温孵化，过筛收集二龄幼虫，经定量后以每盆植株的接种1000条虫的虫量定量接种番茄，按常规方法进行水肥管理。

1.2 总RNA的提取与数据组装

测序样品分别在番茄接种象耳豆根结线虫的2 dpi（days post-inoculation，dpi）、6 dpi取样，每处理3次生物学重复取样，选用全式金公司（北京）的试剂TransZol Up提取根系组织总RNA。提取的RNA需进行浓度和纯度检测，用质检合格的RNA样品用于mRNA文库制备。样品交由武汉康测科技有限公司采用基于UID建库并基于Illumina高通量测序平台构建的Illumina PE文库进行测序。使用Trimmomatic软件对raw data进行Clean，去掉接头、低质量的reads后，需要对clean reads进行新一轮质控，得到clean reads的数据信息，然后再将对比去除来自番茄植株reads（Accession：PRJ⁃NA119，2018），用去除后的clean data与象耳豆根结线虫基因组（Accession：PRJEB36431，2020）进行序列比对，获得全面的转录本信息以及相应的基因信息。

1.3 基因表达定量和差异表达分析

针对基因组注释的蛋白编码基因的mRNA的出现次数（用count表示）进行分析，对基因比对结果在不同区域富集到的片段数进行统计、并以此评估组内及组间样品基因表达特征的相关性，通过edgeR软件分析基因在各样本中的差异表达情况，以以log2FC的绝对值>1且FDR<0.05作为筛选差异基因的标准，并对筛选到的差异基因进行GO功能注释以及KEGG富集分析。

1.4 利用生物信息学软件对差异基因进行致病相关基因的筛选与分析

通过转录组比较分析，筛选出象耳豆根结线虫特异性上调的差异基因，并通过signalP5.0预测具有信号肽的基因，随后利用TMHMM2.0在线软件预测无跨膜结构域的差异基因，筛选出具有信号肽且无跨膜结构域的候选效应子，再利用In⁃terpro及NCBI在线网站对筛选出的基因进行特殊结构域分析及序列比对分析，最后通过String在线软件预测候选基因的蛋白互作网络。

1.5 qPCR验证

利用CTAB方法提取植株2 dpi和6 dpi根部总RNA，利用iScript™ gDNA Clear cDNA Synthe⁃sis Kit（Biorad）反转录试剂盒将提取的RNA逆转录为cDNA，再使用SYBR Premix EX Taq（Taka⁃ra）试剂盒进行荧光定量试验，检测仪器为Quant⁃Studio 6 Flex（ABI）荧光定量PCR仪，选用primer5软件对从预测的效应子中选取的5个基因及内参基因Me-actin进行引物设计（表1），检测所选基因在象耳豆根结线虫侵染番茄植株2 dpi、6 dpi的表达量，以验证高通量测序结果的准确性。

表1 引物序列Table 1 Primer sequence

2 结果与分析

2.1 转录组测序质量分析

通过质检发现样品RNA量都大于1 μg；样品OD260/280值大于1.8，说明样品无污染，电泳结果显示样品RNA条带清晰、没有拖带，说明RNA完整性很好，符合建库要求。对6个样品测序数据进行结果统计（表2），有效数据占原始数据的90%以上，为后期转录组拼接提供了良好条件，除此外碱基质量达到Q20、Q30的百分比均高于98%，质控后Cleandate的GC值含量均高于45%。此外样本的主成分分析发现（图1），同一处理的不同样本之间空间分布差异小，表明具有较好的重复性，2 dpi与6 dpi两个处理之间样本相互分离表明其基因表达相差大，有较好的特异性。综合以上数据说明转录组测序质量较高，可用于深入分析。

图1 样品基因表达水平的主成分分析（PCA）的聚类图Fig.1 Cluster diagram of principal component analysis （PCA） of gene expression levels of samples

表2 转录组测序结果统计表Table 2 Statistical table of transcriptome sequencing results

2.2 差异表达基因的筛选

差异表达基因是指样品间或者同一样品经过不同处理后，上调表达的基因和下调表达基因的汇总。通常从差异倍数和显著水平2方面进行评估，对差异表达基因进行筛选，通过比较2 dpi（为象耳豆根结线虫侵染2 d后的样品）与6 dpi（为象耳豆根结线虫侵染6 d后的样品）之间的基因表达水平，通过鉴定共得到了967个差异表达基因（图2）。其中有647个基因表达量上调，320个基因表达量下调，这些差异基因可能在线虫寄生过程中发挥重要功能。

图2 差异基因表达火山图Fig.2 Volcano of differentially expressed genes

2.3 差异表达基因GO注释分析

对差异表达的967个差异基因进行GO功能注释，共注释到656个差异表达基因，共划分为162个功能组，这些功能组都分布在3个主要类别中，且注释到的上调基因数目多于下调基因数目。将差异表达基因按照参与的生物学过程、构成细胞组分、分子功能进行分类统计（图3）。图中显示2 dpi与6 dpi相比，上调基因被注释到细胞组分分类中的差异表达基因在胞内细胞器（intracellular organelle）、细 胞 器（organelle）和 蛋 白 质 复 合 体（protein complex）这3个层级中的比例较大；在分子功能分类中，比例较为突出的3个层级分别是：catalytic activity（催化活性），氧化还原酶活性（oxi⁃doreductase activity）和脂质结合（lipid binding）；在生物过程分类中，比例较为突出的3个层级分别是：单有机物过程（single-organism process），细胞组分组织或生物合成（cellular component organiza⁃tion or biogenesis）、碳水化合物代谢过程（carbohy⁃drate metabolic process）。从对差异上调基因的GO功能注释结果看出，象耳豆根结线虫在侵染携带Mi-1抗性基因的番茄后极有可能通过上调与氧化还原酶活性等相关基因来提高对活性氧耐受力，从而加强其对寄主的侵染性。

图3 差异基因GO功能注释图Fig.3 Statistics of GO annotation classification of differentially ex⁃pressed genes

2.4 差异表达基因KEGG注释分析

对差异表达的967个差异基因进行KEGG富集分析，共有532个基因得到富集，下调数目多于上调，统计到在上调、下调差异基因中分别显著富集到51和59条通路，富集基因数分别为364和168，共涉及到102个信号通路。由图4可见，上调基因主要富集在碳代谢（Carbon metabolism）、核糖 体（Ribosome）、糖 酵 解/糖 异 生（Glycolysis /Gluconeogenesis）、氨基酸的生物合成（Biosynthe⁃sis of amino acids）、剪切体（Spliceosome）、丙酮酸代谢（Pyruvate metabolism）、蛋白质在内质网上加工（Protein processing in endoplasmic reticulum）等，其中注释到最多基因数的通路为碳代谢和核糖体，而富集程度最大的是谷胱甘肽代谢。

图4 差异基因KEGG富集分析气泡图Fig.4 Bubble diagram of KEGG enrichment analysis of differential⁃ly expressed genes

2.5 致病相关基因的筛选与分析

通过KEGG、GO功能分析结果看出，上调基因中极有可能存在潜在的致病相关基因，通过sig⁃nalP5.0在线软件预测所有上调基因，发现具有信号肽的差异基因共88个，利用TMHMM2.0在线软件预测具有信号肽且无跨膜结构域的差异基因共55个，其中上调基因中具信号肽且无跨膜结构域的差异基因共29个。对预测到的29个致病相关基因进行了基因序列比对分析及保守结构域预测（表3）。除此外，通过String预测候选基因的蛋白互作网络（图5）发现，CAD218859.1与CAD2177646.1之间可能存在基因共现性与共表达的互作类型，CAD2179971.1与基因CAD2177646.1之间也可能存在共表达的互作类型。

图5 候选基因的蛋白互作网络预测图Fig.5 Protein interaction network prediction map of candidate genes

表3 筛选到的致病相关基因信息表Table 3 Information table of pathogenicity⁃related genes screene

由于通过在NCBI中进行基因序列比对分析及结构域分析发现有3个基因具有类RxLR结构域，分 别 是CAD2177575.1、CAD2175494.1、CAD2161020.1，且筛选到的基因中只有2个基因与已报道的南方根结线虫候选效应子类似，分别是CAD2205632.1和CAD2171197.1，且这2个基因差异倍数较高，综合分析这5个基因具有一定的后续研究价值，因此选取这5个基因进行后续qP⁃CR验证。

2.6 实时荧光定量PCR验证

通过对转录组数据的分析及基因序列比对的从预测到的29个致病相关基因中挑选了5个基因，使用qPCR的方法验证转录组数据的准确性，并对验证的5个候选致病相关基因表达量与转录组测序数据进行了相关性分析（图6），CAD2175632.1与CAD2175494基因表达量验证结果与转录组数据之间为负相关关系，但差异不显著，qPCR和RNA-seq在基因表达定量上存在一定差异可能是qPCR的定量只在基因的局部区域 进 行 测 量 所 致，CAD2205632、CAD2177575、CAD2161020基因表达量验证结果与转录组数据之间为正的强相关关系，差异不显著。此外，从图7中可以看出5个候选致病相关基因表达量都有所上调，其中6 dpi与2 dpi相比，CAD2177575、CAD2205632的基因表达量差异极显著，CAD2171197.1与CAD2175494基因表达量差异显著，CAD2161020的基因表达量差异不显著。综上所诉，5个候选致病相关基因经qPCR验证均为象耳豆根结线虫侵染番茄VFNT后在侵染早期表达量上调的基因，且与转录组测序数据具有一定相关性，验证了转录组数据的可靠性。

图6 5个候选基因RNA-seq与qPCR验证结果的相关性热图Fig.6 Correlation heatmap between RNA-seq and qPCR valida⁃tion results of five candidate genes

图7 转录组中预测的5个候选致病相关基因的qPCR验证Fig.7 qPCR verification of five pathogenicity-related gene predict⁃ed in transcriptome

3 讨论

由于植物寄生线虫效应子蛋白在线虫与寄主互作中的重要作用，挖掘线虫基因组中的效应蛋白编码基因成为了研究线虫与寄主互作的重要途径。随着近年来分子生物学技术的发展，也开发了各种研究技术来筛选与鉴定线虫效应蛋白，包括生物信息学、转录组学、基因组学、蛋白质组学等。例如Lee等［21］利用对南方根结线虫分别接种甘薯感病和抗病品种的根系进行转录组分析，确定了可能触发特定信号通路变化的候选致病基因；Pogorelko等［22］使用了感染丁香假单胞菌和荧光假单胞菌的本氏烟草植株以确定大豆孢囊线虫中干扰植物免疫反应的甘氨酸氨基转移酶类效应蛋白，最终筛选鉴定出3个抑制ETI反应的效应蛋白，7个抑制PTI反应的效应蛋白；史倩倩基于南方根结线虫基因组信息，结合生物信息学软件预测了110个具有信号肽，无跨膜结构的候选效应子［23］。此外，一些具有特殊结构域的效应蛋白也常可能为关键致病基因，例如具有RxLR结构域、SPRY结构域、纤维素结合结构域、RD激酶结构域、FLAK结构域等［24-26］。其中具有RxLR结构域的效应子，可以抑制由PTI/ETI引发的HR反应，实现免疫抑制作用，张文玥等［25］通过基于南方根结线虫野生型和Me3毒性群体基因组重测序及转录组分析，初步预测其存在24个候选的类RxLR效应子基因并进行了初步功能验证。所以本试验基于转录组测序数据分析筛选致病相关基因时除信号肽、跨膜结构域预测以外还进行了其它结构域的预测及基因序列比对，从筛选的29个候选致病相关基因中，发现有3个基因也具有RxLR结构域，其它候选基因还具有泛素折叠类似结构域、纤维素酶结构域等与致病相关的结构域，通过NCBI序列比对分析还发现有两个基因与已经报道的南方根结线虫候选效应子具有一定的相似度。

前期研究进行了针对亲和线虫（象耳豆根结线虫）和非亲和线虫（南方根结线虫）接种含有Mi⁃1抗性基因的番茄植株VFNT的组织病理学观察，发现在2 dpi时两者几乎无差别，而非亲和互作6 d时不能像亲和互作一样在根组织内形成巨型细胞，根外形成明显膨大的根结［17］，根据Jones等［27］提出的植物免疫反应zigzag模型，可以推测此时间段内象耳豆根结线虫分泌的某些特有的效应子能够使线虫成功避开寄主植物的抗性反应，实现对寄主植物的寄生。本研究针对象耳豆根结线虫侵染番茄VFNT关键时期进行转录组分析，对效应蛋白相关基因进行初步预测筛选，结果发现在侵染关键阶段差异上调的基因共有656个，GO注释和KEGG注释都说明上调基因功能与象耳豆根结线虫致病力具有一定关系，例如GO功能注释中氧化还原酶活性相关基因的显著上调，同先前的研究显示象耳豆根结线虫对过氧化氢的耐受力强于南方根结线虫的研究结果一致［17］。此外，在KEGG代谢通路分析中发现，富集程度最大为谷胱甘肽代谢，有研究发现谷胱甘肽（GSH）氧化和还原的形式对线虫侵染植物必不可少，因为它调节巨型细胞的代谢活动［18］。其它注释和富集到的功能未见特殊功能富集，但上调的候选致病相关基因可挖掘性强，极有可能存在影响象耳豆根结线虫发育和侵染的关键基因，因此研究针对侵染关键时期所有上调的差异表达基因进行了筛选。

当然在转录组水平对植物寄生线虫效应蛋白编码基因进行初步预测的方法也有其局限性［28］，因为有研究发现有的效应子不具有典型的信号肽，例如Huang等［29］就发现50%的褐飞虱唾液腺蛋白基因不具有信号肽序列。除此以外，这些转录本的基因序列并不完整，还需通过RACE等技术进行全长扩增进行进一步的研究，并且后续还需对候选致病相关基因进行原位杂交、基因沉默等试验对其功能进行进一步验证。

综上，随着响应国家减肥减药双减政策要求和粮食生产需求，番茄生产系统中防治对经济具强破坏性影响的象耳豆根结线虫，了解其寄生致病分子机制，开发有效防治新手段极具研究价值。

4 结论

研究通过象耳豆根结线虫侵染携带Mi⁃1抗性基因的番茄植株关键时期进行转录组分析，并结合生物信息学软件对上调的差异表达基因进行初步筛选，筛选出29个具有信号肽且无跨膜结构域的候选致病相关基因，从中挑选出保守结构域与致病相关的基因以及与已报道的南方根结线虫效应子具有一定同源性的基因进行qPCR试验，经验证5个候选基因表达趋势与转录组一致，皆为上调，且与转录组数据也具有一定相关性，验证了转录组数据的具有可信度的同时，对后续试验挑选候选致病相关基因进行功能验证提供了研究思路。本研究基于转录组筛选出的候选致病相关基因可能在象耳豆根结线虫侵染抗病植株过程中发挥正向调节作用，为研究象耳豆根结线虫致病基因提供了新思路及新可能。