基于高通量测序技术的松材线虫研究进展

丁晓磊，张 悦，林司曦，叶建仁

(南京林业大学，南方现代林业协同创新中心，南京林业大学林学院，江苏 南京 210037)

林业有害生物危害已成为当前制约我国林业建设和经济发展的重要障碍[1]。松材线虫病(pine wilt disease)作为林业毁灭性病害之一，具有传播迅速、防治困难，以及一旦染病全株迅速枯死且不可逆转等特点。松树感病后最快40 d即可枯死，如不及时除治，从单棵松树发病到数百公顷的松林毁灭仅需3～5 a的时间。正是由于该病的严重性，其病原松材线虫(Bursaphelenchusxylophilus)已成为国际上公认的重要检疫性有害生物。松材线虫原产于北美，对乡土松种并未造成明显的损失[2]。但从20世纪初起，松材线虫先后传入日本、中国和韩国等东亚地区，给该地区造成了巨大的经济损失和生态破坏[3]，并于1999年传入葡萄牙，随后扩散至西班牙，危害欧洲大片松林[4-5]。1982年我国首次在南京发现松材线虫病，并在随后的40年内迅速扩展蔓延，至2021年底，我国有19个省(自治区、直辖市)的728个县级行政区被划分为疫区，防控形势异常严峻[6]。

各国学者一直试图通过阐明松材线虫的致病机制来推进病害防控，研究者们围绕该病害已经开展了许多研究也取得了诸多成果，但其致病机制目前尚无定论[7-16]。随着分子生物学的不断发展，松材线虫致病的分子机制研究不断取得突破。近20年来，纤维素酶、果胶裂解酶、细胞色素酶等松材线虫致病关键基因及其潜在调控功能研究已见报道[14-17]，而高通量测序技术的发展和普及为研究该病提供了新的方法。首代短读长测序组装的松材线虫基因组草图在2011年发布[18]，为后续分子致病机理的研究提供了重要的基础数据支持。此后关于该病害的致病机理研究主要聚焦在致病基因、非编码RNA和种群变异等方面[8, 19-22]，但大多以解析具有潜在调控作用的基因为主[20-24]。近年来，三代长读长测序和染色体空间结构测序技术将分子生物学研究带入了新的阶段，使高质量染色体级别基因组的组装变得切实可行，越来越多物种的基因组也同步得到了更新[25-27]。日本学者和笔者所在研究团队分别在2021年、2022年公布了2个版本的松材线虫染色体级别基因组，为开展松材线虫的分子生物学研究奠定了重要基础[28-29]。在此，笔者对近年来松材线虫高通量测序相关研究进行综述，全面概括利用该技术在松材线虫分子致病研究中取得的重要进展，以期为后续开展松材线虫致病的相关分子生物学研究提供借鉴。

1 早期松材线虫病相关致病基因研究

在高通量测序技术问世之前，松材线虫致病相关分子机制主要采用基因同源克隆、抑制消减杂交文库和基因芯片等手段进行相关基因功能研究[14-17, 30]。如2006年Kikuchi等[14]对松材线虫2个果胶裂解酶基因进行了同源克隆，发现这2个基因主要在松材线虫的食道中特异表达，证明了它们在线虫侵入寄主过程中发挥着重要作用。Kang等[15]和Hirao等[17]分别于2009年和2012年对不同龄期的松材线虫进行了差减文库分析，发现了19个特异表达的基因，其中半胱氨酸蛋白酶、毒液过敏原相关基因等与致病性密切相关；Qiu等[30]通过基因芯片技术对接种前后的松材线虫进行基因表达分析，发现1 310个基因中，633个基因上调、677个基因下调。这些差异基因中包含大量致病相关基因(如果胶裂解酶基因、细胞色素P450、UGTs、ABC转运蛋白基因等)，注释结果表明这些致病相关基因主要参与了细胞壁降解和生殖过程。以上研究发现了多个与松材线虫致病相关的关键基因,为阐明松材线虫致病机理奠定了重要基础，但由于当时缺少基因组数据，学者们难以深入挖掘松材线虫特有或未知基因，阻碍了对该病致病机理的进一步认识。

2 基于初代基因组的相关致病基因研究

日本学者在2011年公布了首个基于短读长测序的松材线虫基因组草图，为广大研究者提供了约75 Mb、5 527个scaffolds和12 568个注释基因的基础数据[18]。自此，松材线虫的分子生物学相关研究进入了有参基因组时代，极大地方便了研究人员的实验设计和对关键致病基因的挖掘。

2011年Li等[31]在松材线虫基因组中挖掘到了1个新的2-半胱氨酸过氧化物酶BRP-RX，进一步的研究发现在BRP-RX的一级结构中有1个由23个连续疏水氨基酸组成的跨膜α-螺旋，推测其可能在线虫抵御寄主植物 ROS 防御系统中起着关键作用，是一个潜在的促进松材线虫侵染松树的遗传因子。2012年Wang等[32]对松材线虫精氨酸激酶BxAK1进行了克隆和RNA干扰验证，结果表明沉默该基因可以显著提升松材线虫的死亡率并降低其繁殖率，揭示了其作为防控靶标的潜在可能性。2015年Xu等[33]对3个松材线虫致病相关细胞色素酶P450基因进行了克隆，相关功能分析表明上述基因与松材线虫传播、繁殖和致病力相关。2020年Liu等[34]开展了松材线虫自噬基因相关研究，发现松材线虫自噬基因与宿主体内活性氧代谢密切相关，且在氧胁迫下能够诱导自噬基因的表达，明确了自噬基因在松材线虫致病过程中可能发挥的作用。除了以上研究成果，越来越多关于松材线虫分子生物学探索中开始使用高通量测序技术进行转录组水平的差异基因挖掘和基因组水平的遗传变异研究。基于参考基因组和海量测序数据的帮助，松材线虫分子生物学相关研究在近10年来取得了长足进展。

3 松材线虫高通量转录组测序相关研究

3.1 基于转录组测序的松材线虫功能基因挖掘

转录组测序能够在海量数据的支持下揭示松材线虫致病相关基因的表达模式，以及预测基因功能和发现未知基因，达到对转录水平变化的高覆盖度分析。2012年Yan等[35]对松材线虫与拟松材线虫的转录组(RNA-seq)进行de-novo组装的尝试，基于松材线虫与拟松材线虫分别获得了23 765和21 782个转录本片段，并发现了6个尚未报道的类毒液过敏原蛋白(VAPs)基因，认为松材线虫与拟松材线虫可能在进化过程中形成了相似的分子机制以适应在宿主体内的生活。2016年Tsai等[8]利用RNA-seq分析了松材线虫在接种48 h后基因转录的变化，发现有60个与解毒物质相关的基因在接种后明显上调表达，而下调表达的基因则主要属于胶原蛋白家族，并且雌性松材线虫的尾尖突在接种后产生了较为明显的表型变化。因此，作者认为松材线虫在侵染松树的过程中可能会通过自身表型的变化来适应外界环境的变化。2016年，Ding等[36]对不同毒力的4个松材线虫虫株进行RNA-seq分析，发现了差异表达的238个转录本和84个外显子。差异基因功能分析发现，不同毒力虫株表现出不同的生长、繁殖特征和氧化还原酶活性，明确了以上基因与独立分化之间的关系。Espada等[9]在2016年利用RNA-seq对松材线虫早期侵染阶段进行效应子的筛选，发现多个在感染早期高表达的效应子并进行了初步的功能研究，这也是松材线虫效应子相关研究的首次报道。Hu等[37]和Zhao等[38]分别于2019年和2020年通过RNA-seq成功鉴定了松材线虫中的BxSapB1和BxSapB2等2个效应子。其中BxSapB1可以导致本氏烟草细胞死亡，并对松材线虫毒力有一定的影响；BxSapB2是一种分泌蛋白，在松材线虫的食管腺细胞和化感器中高表达，且在松材线虫与寄主的互作中扮演重要角色。2019年日本研究者Tanaka等[7]对松材线虫不同发育阶段进行RNA-seq测序，发现超过9 000个基因在不同发育阶段有差异表达，其中与线虫生物学有关的基因多涉及生殖和蜕皮方面。2019年Wang等[39]通过RNA-seq测序来研究松材线虫在低温下的分子反应模式，分析了6个低温相关BxGPCR基因的转录丰度后发现，BxGPCR基因在低温下转录水平显著升高，而BxGPCR17454基因沉默后的松材线虫其低温存活率显著下降。这一研究结果部分揭示了松材线虫的低温适应分子机制。2021年，Chen等[40]对扩散性阶段的松材线虫进行转录组测序分析，发现其主要功能富集在细胞自噬和寿命相关通路。后续实验发现，松材线虫中海藻糖相关基因的协同表达促进了3龄幼虫的形成，并帮助其在-20 ℃下存活。其他转录组相关研究还包括石榴苷、甲维盐和蒎烯处理对松材线虫相关基因表达量的影响[41-43]。为了阐明细菌在松材线虫病中可能发挥的作用，陈阳雪等[44]采用对野生型和含有马尾松内生细菌GD2的松材线虫分别进行接种实验，对接种后松材线虫进行转录组测序，比较了不同处理后筛选出的143个差异表达基因，其后续功能分析表明菌株GD2可以通过提高松材线虫的免疫调节能力来影响松材线虫病的发生。

3.2 松材线虫非编码RNA鉴定

目前关于松材线虫非编码RNA的研究较少，且多数都围绕在基于高通量测序的miRNA鉴定和表达水平分析等方面，很少涉及某个miRNA的调控机制研究。如2010年Huang等[45]首次发现松材线虫体内存在非编码RNA，有49个miRNA在种间高度保守，可能调控979个靶标基因。其中14个miRNA参与调控了松材线虫的耐热性，初步揭示了非编码RNA与松材线虫环境适应性之间的关系。2015年，Ding等[22]构建了松材线虫侵染松树后4个不同阶段的miRNA文库，获得了数百个进化上相对保守的miRNA和未知候选miRNA数据集，表达丰度分析发现大部分miRNA在松材线虫病发病中期呈现高水平表达，揭示了miRNA在松材线虫致病过程中的潜在调控作用。Modesto等[46]研究了miRNAs在松材线虫侵染Pinuspinaster早期发挥的调节作用，鉴定出105个响应松材线虫侵染的miRNAs，并根据预测的靶标，认为这些miRNAs中的一部分与茉莉酸反应途径、ROS解毒和萜类生物合成的作用有关。此外，通过对抗病和感病植物的比较，还鉴定出了8个在植物防御和抗病方面具有潜在功能的miRNAs。2017年，Xie等[47]以马尾松为供体植物，利用高通量测序技术研究了松材线虫侵染后前3 d针叶中miRNA的表达。结果表明，松材线虫侵染后3 d内，差异表达的miRNAs在样品中呈现先增加后减少的趋势，其中以感染后第2天采集的样品表达最高。对感染后不同时间针叶中miRNA靶基因的分析表明，植物激素信号转导和RNA转运两条途径富含miRNAs。Cai等[48]于2022年利用miRDeep2、mireap和miR分析器等miRNA预测工具或通过同源作图，已在至少36种寄生蠕虫(包括松材线虫)中鉴定出了不同的miRNA种群。

4 松材线虫伴生微生物群落高通量测序相关研究

松材线虫伴生微生物群落一直以来被一些学者认为可能与松材线虫病毒力和繁殖力等密切相关。基于高通量测序的宏基因组技术可以全面分析松材线虫体表和体内携带的微生物群落情况，明确线虫和微生物之间的互作机制。2012年，Cheng等[49]采用核糖体RNA通用引物扩增后对扩增产物进行RNA-seq测序的方法分析了松材线虫伴生微生物群落,结果发现松材线虫微生物群落的解毒机制以苯甲酸降解代谢为主，不同于松材线虫基于细胞色素CYP450代谢的解毒通路。实验分析还表明，大多数可培养的细菌不仅可以耐受0.4%质量分数的α-pinene，还能作为其生长的碳源，故而认为松材线虫和其共生微生物之间已经建立了互惠互利的解毒机制。2014年，Xiang等[50]通过细菌16S rDNA高通量测序研究了4株不同毒力的松材线虫和4株拟松材线虫的内生细菌结构差异，基于运算分类单元(operational taxonomic unit，OTU)的结果表明，松材线虫强毒虫株AMA3和ZL1的菌群丰富度均高于中毒虫株AA3、HE2和所有拟松材线虫虫株，不同种群线虫所携带的内生细菌差异显著。2016年，Wu等[51]利用Biolog结合高通量测序分析了松材线虫繁殖型和扩散型之间伴生细菌群落差异，高通量分析表明寡养单胞菌(Stenotrophomonassp.)、无色杆菌(Achromobactersp.)和鞘氨醇杆菌(Sphingobacteriumsp.)在繁殖型松材线虫中更丰富，而柠檬酸杆菌(Citrobacteramalonaticus)、弧菌(Vibriosp.)和肠杆菌科(Enterobacteriaceae)在扩散型松材线虫中含量更高，这种变化可能有助于松材线虫适应环境的改变。Xue等[42]为探讨松材线虫内生细菌嗜麦芽窄食单胞菌(Stenotrophomonasmaltophilia)对松材线虫致病力的影响，对无菌线虫和携带嗜麦芽窄食单胞菌的松材线虫进行了转录组测序，结果表明嗜麦芽窄食单胞菌可通过调节松材线虫寄生、免疫和致病相关基因的表达来影响松材线虫的毒力和繁殖，此结果为进一步探讨松材线虫病的发病机制提供了依据。为了弄清松材线虫侵染寄主对寄主内生细菌的影响，芦伟等[52]以 2年生马尾松苗为试验材料，采用 Illumina MiSeq 高通量测序技术分析受松材线虫侵染后的马尾松茎干内生细菌群落结构组成与多样性的变化情况，发现在属水平上马尾松茎干内生细菌组成种类与丰度均发生明显改变，表明松材线虫侵染会对马尾松茎干内生细菌产生显著影响。尹诗恒等[53]为了探究松材线虫侵染下松干、松针与根系部位内生细菌菌群的结构变化，充分挖掘马尾松内生细菌资源，采用变性梯度凝胶电泳(DGGE)及高通量测序技术，对松材线虫侵染下3年生马尾松不同部位内生细菌的结构多样性进行分析，结果表明受侵染后的马尾松根部基本不受影响，但树干部位内生细菌菌群结构变化极显著，其次是松针。Deng等[54]通过从微孔膜和0.5 g土壤中提取基因组DNA,对细菌16S rDNA基因的V3—V4区进行扩增并进行高通量测序，对健康红松和红松早期及后期感染松材线虫的叶圈和根际细菌、真菌群落进行分析，发现随着松材线虫病感染程度的加重，大多数细菌群趋向于共排除而不是共生。

5 松材线虫高通量基因组测序相关研究

5.1 松材线虫短读长基因组重测序

在基因组水平上，高通量测序能够检测松材线虫不同种群间的核苷酸位点差异和基因家族扩张现象，为解析该病害的进化过程、明确关键性状的生物靶标及分析种间差异提供有效的技术手段。2015年Palomares-Rius等[20]将日本境内不同地理来源的松材线虫进行了基于单核苷酸多态性(SNPs)的差异分析。该研究首次报道了松材线虫体内含有约3 Mb的SNPs位点，占到了基因组的4.1%。这些高水平的遗传变异也表明松材线虫在毒力与生态适应性方面具有较高的多样性。2020年Zhang等[55]通过开展松材线虫、动物寄生线虫和植物寄生线虫等线虫种的比较基因组分析，发现松材线虫的51个基因家族都存在基因扩张情况，基中参与外源解毒途径的黄素单加氧酶(FMO)、细胞色素P450(CYP450)、醇脱氢酶(ADH)、葡萄糖醛基转移酶(UGT)和谷胱甘肽转移酶(GST)等基因的扩张现象尤为显著。虽然大多数的扩张可能是由DNA片段重复产生，但其中9个ADH基因可能是通过水平基因转移现象从真菌中获得。研究认为宿主防御化学物质导致了松材线虫外源解毒途径的基因家族扩展，促进了其在树脂道中的生存，提高了线虫在宿主体内的适应性。2021年Filipiak等[56]为了对4株具有不同毒力的虫株进行重测序分析，发现1 576个与毒力差异有潜在关联的位点，为在分子水平上区分不同毒力松材线虫虫株提供了分子靶标。2019年，黄金思等[57]通过Illumina Hi-seq 4000平台对采集自广东省不同区域的30株松材线虫虫株进行基因组重测序，利用生物信息学软件对其SNPs位点信息及基因型进行分析，最后通过主成分和聚类分析将这30株虫株分为3类，完成了基于SNP标记的广东省松材线虫种群分化研究。2021年Ding等[28]在松材线虫新的染色体级别基因组基础上对来自中国和美国的181个菌株进行了重测序分析，发现了约780万个SNPs位点，并将中国松材线虫种群细化为4个类群。后续分析发现，不同类群线虫与地理来源关系密切，不同地区间存在严重的交叉感染现象。此外，该研究还建立了一种基于机器学习的松材线虫疫源追溯体系，为控制病害蔓延提供了有力的防控手段。

5.2 基于长读长测序技术的松材线虫基因组组装

第2代测序技术的出现，极大地推进了基因组层次研究的进展，使得基因组测序工作变得比较简单。但其缺点也日益凸显，最主要的问题是第2代测序读长短，测通困难[58]，这给序列拼接、组装及注释等分析带来巨大的困难。而第3代测序技术是在单一的、由 DNA 分子组成的阵列上进行的合成测序。这种技术可以大大增加独立分析的 DNA 片段的数目，增强产出序列的长度，极大地简化了基因组组装过程[59]。目前，在第3代测序技术的帮助下，已陆续完成了各物种基因组版本更新，提高了测序结果的完整度和精细度，且很多物种已有端粒至端粒的完整精细图谱。高质量的基因组无疑为后续所有分子生物学相关研究提供了重要的基础数据支持，Dayi等[29]在2020年将其发布的初代基因组使用Nanopore长读长序列和HI-C技术进行了更新，获得了松材线虫染色体级别基因组。几乎同时期，Ding等[28]利用PacBio+Illumina +Bionano+Hi-C多项技术对松材线虫进行了基因组组装，两项研究均获得了6条近完整的染色体序列，将基因组大小提高至78 Mb。以上新基因组的scaffold长度均达到了12 Mb，是早期基因组(0.2 Mb)的60倍，平均scaffold长度为其24倍。此外，功能注释和RNA-seq数据分析证明，借助新版本基因组可以鉴定出更多的新基因和亚型，这将有助于研究人员获取更完整的松材线虫转录组信息。以上两个版本的基因组数据可以完成对松材线虫复杂重复区域的高度覆盖，弥补了上一版本松材线虫基因组草图的不足，为后续相关研究奠定了重要的基础。

6 结 语

随着高通量测序技术的快速发展，越来越多的松材线虫关键致病基因、非编码RNA，以及种群遗传差异位点被报道，加深了对松材线虫分子致病机理的理解。尤其是随着松材线虫最新版本基因组数据的公布，研究者们将向着更为精细、前沿的研究领域探索。在转录组相关研究方面，可尝试开展空间转录组、单细胞测序等在松材线虫中尚未使用的新技术[60-62]，在染色体基因组的支持下寻找松材线虫中具有潜在调控功能的长链非编码RNA、融合基因和环状RNA等罕见核酸分子；挖掘可能存在的甲基化、RNA编辑、结构变异等未知现象。在基因组学研究方面，可以开展基因家族扩增、大样本全基因组关联分析和数量性状定位等研究，以挖掘与致病力和繁殖力等重要性状紧密关联的基因或位点，阐明松材线虫不同种群在进化过程中发生的适应性变化。以上新的研究方向将为松材线虫病相关分子研究提供全新的思路和切入点，有助于深入揭示其潜在的致病机理。