苜蓿滑刃线虫线粒体基因组及其系统发育研究

2021-08-22薛清杜虹锐薛会英王译浩王暄李红梅

生物技术通报 2021年7期

薛清杜虹锐薛会英王译浩王暄李红梅

（1. 南京农业大学植物保护学院农作物生物灾害综合治理教育部重点实验室，南京 210095；2. 西藏农牧学院资源与环境学院，林芝 860000）

滑刃线虫属（AphelenchoidesFischer，1894）已报道种类约200种，大多数种类取食真菌，有13个种类为兼性植物寄生线虫，具有广泛的寄主范围［1］。贝西滑刃线虫（A. besseyi）、草莓滑刃线虫（A.fragariae）和菊花滑刃线虫（A. ritzemabosi）因其对作物产量和品质造成严重影响而受到广泛的关注和研究，它们的典型寄主分别为水稻、草莓和菊花，并引起水稻干尖病、草莓夏矮病和菊叶线虫病［2-4］。针对其他植物寄生性滑刃线虫种类的研究则相对缺乏。

苜蓿滑刃线虫A. medicagus是新近描述的一种滑刃线虫［5］。该线虫分离自宁波口岸进境的美国干苜蓿草，可在链格孢菌（Alternariasp.）、灰葡萄孢菌（Botrytis cinerea）等多种真菌上完成生活史，同时对大豆和苜蓿具有弱寄生性［6］。系统发育研究显示，苜蓿滑刃线虫与贝西滑刃线虫、福建滑刃线虫（A. fujianensis）亲缘关系较近［5］，但其隶属的滑刃线虫属在垫刃亚目（Tylenchina）中的分类地位尚未完全解析，特别是与真滑刃线虫属（Aphelenchus）的高阶元系统发育关系仍然存在较多的争议［7-8］。现有的线虫系统发育研究大多基于核糖体RNA与线粒体COI基因，这些基因的片段一般较短（常用片段约为450-1 600 bp），所含信息量有限，难以完全解析各物种间的进化关系。动物线粒体基因组一般在10 kb以上，兼具编码与非编码区域，富含丰富的遗传信息，其中基因结构排列、核苷酸和氨基酸序列，都可成为研究进化关系的良好靶标。因此，基于线粒体基因组的系统发育已成为近年来分类进化研究方向的热点之一［9］。

线虫线粒体基因组的研究相比其他动物类群起步较晚，目前已完成的线粒体基因组测序绝大多数集中于重要的动物寄生线虫，而植物寄生类和土壤腐生类的线虫鲜有报道［10］。滑刃线虫属种类众多，然而目前仅有贝西滑刃线虫完成了线粒体全基因组测序［11］。造成这种现象的重要原因之一是现有的线虫线粒体基因组测序大多采用操作难度较大的长PCR扩增法，该方法是根据已有的线粒体基因组信息，设计多对特异性引物进行长片段的扩增与测序［11-13］。然而，线虫种间的线粒体基因差异较大，且可供参考的线虫线粒体基因组较少，因此可能会导致引物设计存在较大盲目性，使实践操作中成功扩增长片段变得困难。近年来，高通量测序技术的发展为线粒体基因组测序提供了新的策略，已经在多种生物中获得了较好的结果［14］，然而该方法尚未在线虫的线粒体基因组研究中得到广泛应用。目前采用高通量测序法仅完成了艾灵顿球孢囊线虫（Globodera ellingtonae）和哥伦比亚纽带线虫（Hoplolaimus columbus）线粒体基因组测序与组装［15-16］，尚未见其他线虫类群的国内外相关报道。

本研究利用二代测序技术Illumina平台对苜蓿滑刃线虫进行低覆盖全基因组测序，从获得的原始数据中提取线粒体信息，通过多种方法组装线粒体基因组后，对线粒体基因组的基因结构和排列、碱基组成，蛋白编码基因和密码子使用情况，以及RNA基因和非编码区的结构等进行预测分析，结合已发表的线虫线粒体基因组序列对苜蓿滑刃线虫的系统发育关系进行探讨，旨在为滑刃线虫的分类与线粒体基因组研究提供新的科学方法和依据。

1 材料与方法

1.1 材料

供试苜蓿滑刃线虫来自本实验室在灰葡萄孢平板上保存的纯培养，为该种描述时使用的模式群体［5］。挑取300条线虫用无菌蒸馏水清洗后，放入装有8 μL蒸馏水的200 μL PCR管中。采用液氮冻融法进行线虫DNA的粗提，将PCR管在液氮中冷冻1 min，然后迅速放入65℃水浴锅中溶解2-3 min，重复冻融步骤3次后，向PCR管中加入100 mg/mL蛋白酶K 1 μL和10×PCR缓冲液1 μL，混匀后于65℃保温1.5 h，然后于95℃保持15 min。使用Qubit 1x ds DNA HS检测试剂盒及 Qubit®3.0型荧光定量仪测定 DNA溶液的浓度。

1.2 方法

1.2.1 组装种子序列的获取本研究采用线粒体COI基因部分序列作为种子，用于线粒体基因组组装的起始序列。采用引物对COI-F1（5'-CCT ACT ATG ATT GGT GGT TTT GGT AAT TG-3'）与COI-V2（5'-GTA GCA GCA GTA AAA TAA GCA CG-3'）扩增COI基因片段，预测片段大小约为720 bp［17］。扩增体系为25 μL，含有DNA模板1 μL、浓度10 μmol/L的上游和下游引物各2 μL、100 μmol/L dNTPs 2 μL、10×ExTaq PCR缓冲液2.4 μL、25 mmol/L MgCl22 μL以及5 U/μL ExTaq酶0.1 μL，加ddH2O补足至25 μL。PCR反应条件为：94℃ 5 min；94℃ 30 s，51℃ 30 s，72℃ 2 min，循环42次；72℃延长10 min。PCR产物经凝胶电泳检测后，由北京擎科生物科技有限公司采用Sanger法进行测序。

1.2.2 高通量测序与线粒体基因组组装采用标准的 Illumina TruSeq Nano DNA LT文库制备实验流程构建苜蓿滑刃线虫DNA样本所需的基因组上机文库，由南京派森诺基因科技有限公司使用 Hiseq X ten PE150平台进行高通量测序，测序产生2 Gb数据量，得到低覆盖全基因组测序数据。采用FSATP［18］对下机数据进行质控，去除接头污染，过滤掉长度低于50 bp、Q值低于20或n数量大于3的reads。

从GenBank中下载贝西滑刃线虫（A. besseyi）、伤残短体线虫（Pratylenchus vulnus）、香蕉穿孔线虫（Radopholus similis）、燕麦真滑刃线虫（Aphelenchus avenae）、秀丽隐杆线虫（Caenorhabditis elegans）、大豆孢囊线虫（Heterodera glycines）以及花生根结线虫（Meloidogyne arenaria）等7种线虫的线粒体基因组序列作为参考基因组。在Linux系统下使用NextGenMap［19］根据线粒体基因组参考序列快速比对全基因组测序数据，比对阈值设定为0.3。通过SAMtools［20］提取比对结果中至少在单向符合要求的候选线粒体序列。利用上述获得的COI基因扩增序列作为“种子”，结合自定义循环脚本与NOVOPlasty的种子序列扩展（seed-extend）算法，对已提取的线粒体基因组进行组装。使用MitoZ［21］软件并结合已发表的贝西滑刃线虫［11］线粒体基因组信息，对获得的苜蓿滑刃线虫线粒体基因组进行注释。根据基因组初步组装结果，在nad4与cox1两个基因上分别设计正向引物AM_F（5'- TTA TAC TCT ATT CGT TTT AGT GTC AGG TTT-3'）与反向引物AM_R（5'- TCA AAC AGC AAA ACC ACC TTG ATA CT-3'），以扩增非编码片段。PCR采用 Prime STAR MAX酶，反应条件为：94℃ 5 min；94℃ 30 s，48℃ 30 s，62℃ 2 min，循环38次；62℃延长10 min。最后，将基因组组装结果与扩增获得的非编码区进行拼接，将拼接后的基因组序列上传至GenBank，登录号为MZ424209。

1.2.3 基因注释与系统发育研究采用MITOS在线服务（http：//mitos.bioinf.uni-leipzig.de/index.py）对组装的线粒体基因组进行注释［22］，并结合贝西滑刃线虫线粒体基因组已有的注释对其进行修正。利用MITOS软件的MiTFi 程序对tRNA二级结构进行预测［23］，并用FORNA（http：//rna.tbi.univie.ac.at/forna）［24］软件绘制tRNA二级结构。

为研究苜蓿滑刃线虫在线虫门的系统发育地位，使用线虫门31个不同种类的线粒体基因组作为参考序列，以科斯格罗夫异索虫（Thaumamermis cosgrovei）作为外群构建基于极大似然法的系统发育树。使用Geneious7.1.3软件将12个蛋白编码基因（PCGs）按照无脊椎动物线粒体密码子分别转换为氨基酸序列，并在MAFFT软件中进行多重比对［25］。采用SMS软件［26］对12个基因进行评估，选择的模型为MtZoa+G+I+F。按照nad6-nad5-nad4-nad4Lnad3-nad2-nad1-cytb-cox3-cox2-cox1-atp6的顺序构建序列矩阵，利用CIPRES运算平台［27］的RAxML v8.2.12［28］软件进行1 000次bootstrap运算，获得极大似然自举值（BS），构建极大似然系统发育树。

2 结果

2.1 苜蓿滑刃线虫线粒体基因组的基因结构与排列

苜蓿滑刃线虫DNA样本经Hiseq X ten PE150平台测序共获得213 159 776条原始序列，通过质控处理后得到211 293 116条高质量序列，占总原始序列的99.12%，测序总体质量较高。与7种线虫的线粒体基因组序列进行比对和过滤后，共获得线粒体序列1 430 492条，占获得高质量序列的0.68%。拼接共获得13 890 bp序列，为典型的闭合双链DNA分子（图1）。

图1 苜蓿滑刃线虫线粒体基因组结构与基因排列顺序Fig.1 Genome structure and gene arrangement in A. medicagus

通过线粒体基因组基因排列顺序研究发现，苜蓿滑刃线虫与贝西滑刃线虫、松材线虫和拟松材线虫的基因排列完全一致，而与燕麦真滑刃线虫及根结线虫、孢囊线虫、短体线虫等垫刃类线虫存在较大的差异（图2）。经过注释分别属于36个基因，其中PCGs 12个，分别是atp6、cob、cox1-3、nad1-6和nad4L，无atp8基因；编码线粒体核糖体RNA的2个rRNA（rrnS和rrnL）和22个tRNA基因。非编码区位于nad4与cox1基因间，该部分未能通过测序片段完成组装。采用结合非编码区两侧区域的引物AM_F与AM_R进行扩增，并对获得的片段进行测序，结果显示非编码区长度约1 600 bp，显著少于贝西滑刃线虫的3 142 bp。通过Sanger测序获得其中1 393 bp序列，其余部分由于存在大量AT串联重复，测序未能成功。对获得的非编码区进行分析显示，除AT短重复外，该区域至少存在9个103 bp的串联重复。将Sanger测序获得的非编码区与高通量测序的组装结果进行拼接，共获得14 411 bp的苜蓿滑刃线虫基因组，碱基A、T、C、G的含量分别为23.7%、55.9%、10.2%、10.2%，存在明显的AT 偏好性。

图2 垫刃亚目不同线虫种群的线粒体基因排列顺序Fig.2 Mitochondrial gene arrangement among different species in suborder Tylenchina

2.3 苜蓿滑刃线虫线粒体基因组的rRNA基因和tRNA基因结构

苜蓿滑刃线虫线粒体基因组中rrnS基因长度为910 bp，位于trnH和nad3基因间；rrnL基因长度为667 bp，位于trnE与trnS2之间，两者AT含量分别为84.9%和83.2%。MUSCEL序列比对显示，苜蓿滑刃线虫的rRNA与贝西滑刃线虫的rrnL和rrnS的相似度最高，分别为72.1%和77.2%。对苜蓿滑刃线虫基因组进行同源性比对分析，预测得到了22个tRNA基因的序列，有关tRNA基因的位置见表1。22个tRNA基因的长度在53-57 bp 范围之间，其中最长的为trnK，最短的为trnS1。

表1 苜蓿滑刃线虫的线粒体基因位置分布与起始终止密码子Table 1 Placements and start/stop codons of mitochondrial genes in A. medicagus

对tRNA基因二级结构的预测显示，得到的tRNA基因大部分均无法折叠形成典型的后生动物的三叶草结构，只有trnI基因序列的二级结构能折叠成典型的三叶草结构，其余21个tRNA 的二级结构均为非典型（图3）。这些非典型的二级结构均属于可变茎（variable arm）及配对臂缺失，呈TV-loop 结构。

图3 苜蓿滑刃线虫线粒体基因组tRNA二级结构预测Fig.3 Predicted secondary structure of tRNA in the mitochondrial genome of A. medicagus

2.4 苜蓿滑刃线虫系统发育分析

将12个PCGs（atp6、cob、cox1-3、nad1-6和nad4）基因分别转换为氨基酸序列，联合构建了包含3 948个信息位点的序列矩阵。采用极大似然法构建了包括线虫门两个亚目20个科共30种线虫的系统发育树（图4），从图中可以看出，苜蓿滑刃线虫与贝西滑刃线虫互为姐妹群，且获得高度支持（BS=100），且与松材线虫、拟松材线虫处在高度支持（BS=100）的滑刃线虫科（Aphelenchoididae）单系中。燕麦真滑刃线虫与滑刃科线虫的亲缘关系虽然较远，但还是以较高的支持率与小杆亚目（Rhabditina）、旋尾亚目（Spirurina）、滑刃线虫科、斯氏线虫科（Steinernematidae）以及全凹线虫科（Panagrolaimidae）的部分种类聚类在一个大分支中（BS=97）。以短体线虫科（Pratylenchidae）、异皮线虫科（Heteroderidae）与根结线虫科（Meloidogynidae）为代表的垫刃类线虫与寄生脊椎动物的盘尾丝虫科（Onchocercidae）互为姐妹群，但拓扑结构支持率较低（BS=53），且两者构成的分支在系统发育树的地位尚无法解析（BS=29）。

图4 基于12个蛋白编码基因的氨基酸序列构建的极大似然系统发育树Fig. 4 Phylogenetic tree based on amino acid sequences from concatenated 12 protein coding genes by maximum likelihood method

3 讨论

线粒体基因组信息在物种多样性与系统发育研究上有着广泛的应用前景，但因传统PCR扩增测序方法效率低，成本较高，目前仅有105个物种完成了线粒体基因组测序。相比超过25 000种已描述的线虫［29］，大多数线虫的线粒体基因信息仍然未知。已有研究显示，线虫的线粒体基因结构为一个双链闭合的单环或双环［10］，大小在12-22 kb，一般包括36个基因，即编码氧化磷酸化所需酶的12个蛋白编码基因（PCGs），分别是atp6、cob、cox1-3、nad1-6和nad4L，编码线粒体核糖体RNA组分的2个rRNA（rrnS和rrnL），以及翻译不同线粒体蛋白的22个tRNA。目前已知只有寄生脊椎动物的旋毛属（Trichinellaspp.）和鞭虫属（Trichurisspp.）包含atp8基因，而大多数线虫缺少atp8基因［12］。本研究共获得14 411 bp的苜蓿滑刃线粒体基因组序列，另有约1 600 bp序列因存在大量AT串联重复，测序未能成功。对基因组成分析显示，苜蓿滑刃线虫的线粒体基因组与其他大多数线虫较为相似，具有除atp8外的36个基因；在基因排列与组成上，与贝西滑刃线虫、松材线虫和拟松材线虫这3种滑刃科线虫的基因排序结构完全相同且序列相似度最高，而与燕麦真滑刃线虫以及伤残短体线虫、香蕉穿孔线虫、大豆孢囊线虫和花生根结线虫等垫刃类线虫存在较大差异，显示滑刃线虫科与真滑刃线虫科以及其他垫刃类线虫均存在较远的亲缘关系，这与前期线粒体基因组系统发育研究［8，13］与核糖体RNA系统发育结果一致，支持滑刃科线虫与垫刃类线虫分别有着独立的进化起源［30］。

对tRNA的比较分析显示，苜蓿滑刃线虫与贝西滑刃线虫［11］最为相似，两者具有相似的基因长度（分别为53-57 bp与51-59 bp），比哥伦比亚纽带线虫的基因长度变化范围（52-72 bp）更窄。苜蓿滑刃线虫tRNA二级结构中有一个为典型的三叶草结构（rrnI），其余均为可变茎及配对臂缺失，贝西滑刃线虫tRNA也均为可变茎及配对臂全部缺失，两者较为相近，而哥伦比亚纽带线虫中有7个tRNA具三叶草结构［16］，明显多于苜蓿滑刃与贝西滑刃线虫线虫，因此二级结构同样支持滑刃科线虫具有更为相近的亲缘关系，而与垫刃类线虫存在一定差异。

本研究中测序获得的原始数据大部分为细胞核基因，仅有0.68%属于线粒体基因，这与昆虫中0.5%-1.4%的比例类似［31-32］，但是明显高于哥伦比亚纽带线虫的0.2%线粒体序列含量［16］。截至到 2021年初，部分测序公司的建库费用已降低至200-300元/文库，单个样品二代测序费用亦低至50元/Gb。因此，随着测序成本的降低，通过增加测序通量弥补线粒体数据含量低，进而获得充足的测序深度和组装质量的策略在实践中成为可能。

本研究首次完成了苜蓿滑刃线虫的线粒体基因测序，结果显示，低覆盖全基因组测序法可以成功组装出线粒体基因组中大部分序列，证明了利用该方法获取线虫线粒体全基因组序列的可行性。同时，线粒体基因组与传统的核糖体RNA相比拥有更多的信息位点，构建的系统发育树也具有更高的支持率，因此低覆盖全基因组测序法在线虫线粒体基因组研究上将有广阔的应用前景。