3个蜈蚣科物种线粒体基因组测序及唇足纲的系统发育分析

潘阳洋，陈天韵，尤春雪，蒋 超

（1.天津农学院动物科学与动物医学学院天津市农业动物繁育与健康养殖重点实验室，天津 300384；2.中国中医科学院中药资源中心道地药材品质保障与资源持续利用全国重点实验室，北京 100700）

蜈蚣是唇足纲蜈蚣目动物的统称，目前已记录5科27 属.作为节肢动物进化的重要环节，蜈蚣的遗传分类对物种进化研究具有重要意义.目前已有许多基于单基因序列或多基因联合分析蜈蚣进化、系统发育和分类鉴定的报道，尤其是基于DNA 分子系统发育的研究[1].随着二代DNA 测序技术的快速发展，越来越多的学者使用线粒体基因组（mitochondrial genome）研究物种分类、生物体的遗传发育和生物系统发育[2].区别于核基因，线粒体基因组是母系遗传，其结构简单，易于分析[3]. 线粒体基因组通常是1个环状双链DNA 分子，包括13个蛋白质编码基因（PCG）、22个转移RNA 基因（tRNA）、2个核糖体RNA 基因（rRNA）和1个控制区（CR）[4-5].虽然线粒体基因组的结构相对保守，近年来许多研究却发现，在蜈蚣的线粒体基因组中存在基因重排现象，并且伴随有tRNA 二级结构的缺失等变化[6].目前已公布的唇足纲物种的线粒体基因组序列仅有14个，在蜈蚣目物种的线粒体基因组序列中，有4个属于蜈蚣属（Scolopendra），1个属于尖盲蜈蚣属（Scolopocryptops），缺乏科内其他属的数据，因此难以进行跨种属线粒体结构分析.

本研究对蜈蚣科的3个物种，即卫蜈蚣属长足卫蜈蚣（Rhysida longipes（Newport，1845））、筛孔蜈蚣属筛孔蜈蚣sp.1（Ethmostigmussp.1）和蜈蚣属日本蜈蚣（Scolopendra japonicaL.Koch，1878）进行全长线粒体基因组测序，分析其线粒体基因组的基因重排模式，并依据13个PCGs 进行唇足纲的系统发育分析，为蜈蚣的形态学分类提供数据支持，并有助于了解蜈蚣目独有的基因重排模式.

1 材料与方法

1.1 样品采集和DNA 提取

长足卫蜈蚣（CMMI：20200511008），采集于云南省红河哈尼彝族自治州；筛孔蜈蚣sp. 1（CMMI：20191212012），采集于云南省德宏傣族景颇族自治州盈江县太平镇；日本蜈蚣（CMMI：20190904004），采集于广西壮族自治区来宾市金秀瑶族自治县.样本保存于70%～75%的乙醇中.

取每个样本的3 条足，干燥后研磨成粉末放入2.0 mL 的离心管中，然后使用Promega WizardRSV Genomic DNA Purification kit 试剂盒（美国Promega 公司）提取蜈蚣DNA，DNA 样品储存于-20 ℃冰箱待用.

1.2 线粒体基因组的测序和组装

用InvitrogenCollibriNGS建立基因组DNA文库，用illuminaHiSeq2500 平台，按照PE150bp 规程进行测序，由北京诺禾致源科技股份有限公司完成.每个物种的线粒体基因组数据为5×109～6×109.用GetOrganelle[7]软件进行序列组装，获得线粒体基因组序列.

1.3 线粒体基因组的注释

应用Mitos2[8]注释tRNAs、rRNA 及PCGs 的位置，对于未注释上的tRNA，应用Phylosuite 1.2.2[9]软件并利用近缘物种比对的方式进行预测.

通过rtools[10]网站（http：//rtools.cbrc.jp/cgi-bin/index.cgi）预测tRNA 的二级结构，NCBI 网站Conserved Domains[11-14]（https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi）验证13个PCGs 的编码区域并进行调整，OGDRAW[15]网站（https：//chlorobox.mpimp-golm.mpg.de/OGDraw.html）绘制线粒体基因组圈图.

应用AT-Skew 和GC-Skew[16]评估所有线粒体基因组的链不对称性，应用CodonW[17]软件计算密码子偏好性（RSCU）.

式中：A、T、G、C 均为碱基.

1.4 系统发育分析

为判断线粒体基因组数据与传统的形态学系统发育关系之间的联系，选取3种蜈蚣和已公布的唇足纲共14个物种（表1）的全长线粒体基因组序列建立系统发育树，并且选取蚰蜒目的蚰蜒（Scutigera coleoptrata（Linnaeus，1758），AJ507061.2）作为外类群.利用CLUSTAL W[18]在BIOEDIT 7.1.3.0[19]中比对13个PCGs.用Phylosuite 1.2.2[9]软件分别构建贝叶斯推断法（Bayesian inference，BI）和最大似然法（maximum likelihood，ML）的系统发育树，并且利用ModelFinder[20]功能选择合适的模型，最大似然法的模型为GTR+F+R3，贝叶斯推断法的模型为GTR+F+I+G4. 使用IQTREE 构建最大似然树，采用超快自展（ultrafast bootstraping）法重复100 000 次计算分支置信度；应用Mr-Bayes 构建贝叶斯树，利用4 条蒙特卡罗马尔可夫链（Markov Chain Monte Carlo，MCMC）同时运行1 000 000 代，抽样频率为1 000，摒弃25%的老化样本.

表1 本研究中的3种蜈蚣和已发布的11个唇足纲物种的线粒体基因组数据Tab.1 Mitochondrial genome data of the three species in this study and the other 11 species from Chilopoda data of which have been published

2 结果与分析

通过高通量测序技术、序列组装和基因组注释，获得长足卫蜈蚣、筛孔蜈蚣sp.1 和日本蜈蚣的全长线粒体基因组，均为环状DNA，其序列长度分别为14 992、15 002 和14 425 bp.由表1 可以看出，除石蜈蚣目的长角石蜈蚣（Cermatobius longicornis（Takakuwa，1939））线粒体基因组长度为16 833 bp 外，大部分的唇足纲物种线粒体基因组长度为14 000～16 000 bp.

2.1 线粒体基因组结构及组成分析

3种蜈蚣的线粒体基因组均包含37个基因，包括13个PCGs、22个tRNA 基因和2个rRNA 基因，如图1 所示.蜈蚣目物种线粒体基因组中，不同PCGs 的长度差异显著，13个PCGs 长度最小的为ATP8，8个物种的平均值为158 bp，最大的为NAD5，平均值为1 622 bp，如图2 所示.多数PCGs 在不同蜈蚣物种之间长度存在差异，如日本蜈蚣ATP6长度最大，达696 bp，而筛孔蜈蚣sp.1 的ATP6长度最小，仅361 bp.同其他PCGs 相比，ATP8的长度最稳定.13个PCGs 在不同蜈蚣物种间的长度差异分别为：ATP8，12 bp；ATP6，335 bp；COB，31 bp；COX1，45 bp；COX2，23 bp；COX3，53 bp；NAD1，96 bp；NAD2，169 bp；NAD3，122 bp；NAD4，189 bp；NAD4L，81 bp；NAD5，527 bp；NAD6，331 bp.

图1 3种蜈蚣的线粒体基因组圈图Fig.1 Mitochondrial genome composition of three species of Scolopendridae

蜈蚣目物种的全线粒体基因组的核苷酸组成如图3 所示.由图3 可以看出，AT 含量均大于GC 含量，表现为AT 富集，各物种的GC 含量在21.22%（少棘蜈蚣）～32.52%（筛孔蜈蚣sp.1）之间，其中，筛孔蜈蚣sp.1 的GC 含量与长足卫蜈蚣（32.27%）的几乎相等.

图3 8个蜈蚣目物种线粒体基因组的核苷酸组成百分比Fig.3 Nucleotide composition percentage of mitochondrial genomes of eight species from Scolopendromorpha

2.2 同义密码子相对使用度（RSCU）及基因重排分析

蜈蚣目8种蜈蚣线粒体基因组编码氨基酸的RSCU 在0～3 之间（图4），通常小于1，因此，以密码子的RSCU>1 筛选使用频率较高的密码子.8个蜈蚣目物种中，线粒体基因组RSCU>1 的密码子有26～33个，而在长足卫蜈蚣、筛孔蜈蚣sp.1 和日本蜈蚣中，有26～31个密码子被高度使用.其中，筛孔蜈蚣sp.1 和长足卫蜈蚣的重叠密码子共有28个，而3种蜈蚣的重叠密码子为18个.除此之外，长足卫蜈蚣和日本蜈蚣还重叠1个密码子，即GUA.另外，蜈蚣的氨基酸RSCU 值均以2 作为分水岭，因此将氨基酸RSCU>2 作为筛选条件，可以得到3种蜈蚣高度使用的氨基酸是丙氨酸（Ala）、精氨酸（Arg）、甘氨酸（Gly）、异亮氨酸（Ile）、亮氨酸（Leu（taa）和（tag））、脯氨酸（Pro）、丝氨酸（Ser（tag））、苏氨酸（Thr）、缬氨酸（Val）.这10个氨基酸中，除了赤蜈蚣的亮氨酸Leu（taa）的RSCU 低于2 之外，其他氨基酸在8个蜈蚣目物种中也是被高度使用的.

图4 8个蜈蚣目物种线粒体基因组的氨基酸同义密码子相对使用度Fig.4 RSCU of mitochondrial genomes of eight species from Scolopendromorpha

在8个蜈蚣目物种中，13个PCGs 和2个rRNAs的线粒体基因组排列顺序类似，除少棘蜈蚣之外，以COX1为起始，各物种的排列顺序均为COX1-COX2-ATP8-ATP6-COX3-NAD3-NAD5-NAD4-NAD4L-NAD6-COBNAD1-rrnL-rrnS-NAD2（图5），筛孔蜈蚣sp.1 线粒体基因组中NAD5、NAD4、NAD4L、NAD1、rrnL、rrnS的转录方向与其他序列相反.少棘蜈蚣的排列顺序与其他蜈蚣不一致，为COX1-COX2-ATP8-ATP6-COX3-NAD3-NAD1-rrnL-rrnS-NAD5-NAD4-NAD4L-NAD6-COB-NAD2，说明NAD1-rrnL-rrnS基因簇发生了一次整体易位.

在3种蜈蚣的线粒体基因组的基因排序中，PCGs的排列基本稳定，重排现象主要出现在tRNA 的排列中，与其他序列的基因排序相比，trnQ、trnY、trnE的排序较为不固定，并且在rrnL和rrnS中穿插有trnV.

另外，长足卫蜈蚣的14个tRNA 基因（trnI、trnM、trnW、trnK、trnD、trnG、trnA、trnR、trnN、trnS（gct）、trnE、trnT、trnS（tga）、trnL（tag））和8个PCGs（COX1、COX2、COX3、NAD2、NAD3、NAD6、ATP6、ATP8）都在正链转录，其余8个tRNA 基因、5个PCGs 和2个rRNA 基因（rrnS、rrnL）在负链转录. 相比之下，日本蜈蚣的线粒体基因组的基因大多从正链转录，为18个tRNA基因（trnI、trnM、trnW、trnL（tag）、trnL（taa）、trnK、trnD、trnG、trnA、trnR、trnN、trnS（tga）、trnS（gct）、trnH、trnT、trnY、trnV、trnE） 和8个PCGs（NAD2、NAD3、NAD6、COX1、COX2、COX3、ATP6、ATP8），其余4个tRNA 基因、5个PCGs 和2个rRNA 基因（rrnS、rrnL）在负链转录.

此外，在筛孔蜈蚣sp. 1 基因组中，trnR、trnP、trnS、trnV缺少TΨC 臂或者DHU 臂，并且日本蜈蚣同样也有trnR、trnG、trnA、trnH、trnT、trnL（taa）缺少TΨC 臂或DHU 臂的现象.

2.3 唇足纲的系统发育关系

本研究构建的ML 和BI 系统发育树均高度支持（PP=1，BS=100）蜈蚣科分为蜈蚣亚科和耳孔蜈蚣亚科，如图6 所示，由图6 可以看出，在蜈蚣亚科蜈蚣属中，哈氏蜈蚣、少棘蜈蚣和海南蜈蚣亲缘关系很近，并且与日本蜈蚣互为姊妹群（PP=1，BS=100），与传统的形态学分类结果相符.此外，系统发育树显示，地蜈蚣目和蜈蚣目互为姊妹群.地蜈蚣目和蜈蚣目均为整形亚纲分类群，符合传统唇足纲系统发育的研究结果.

3 讨论与结论

本研究首次对蜈蚣科长足卫蜈蚣、筛孔蜈蚣sp.1和日本蜈蚣的线粒体基因组进行测序、组装和注释，确定了22个tRNAs、13个PCGs 和2个rRNAs 的位置，并且发现大多数基因位于正链上. 在利用Mitos2网站对筛孔蜈蚣sp.1 进行注释时，发现PCGs 有大量的重叠区域，因此有必要采用近缘物种比对的方法进行准确注释.与蜈蚣目其他物种相比，筛孔蜈蚣sp.1 中NAD5、NAD4、NAD4L、NAD1、rrnL、rrnS等多个PCGs转录方向相反，这种情况是否是筛孔蜈蚣属的共有特征有待采集该属其他物种进行进一步研究.

与PCGs 的排序相比，蜈蚣tRNA 的重排现象更频繁，预测的tRNA 无法形成正确二级结构的现象也很普遍. 在Park 等[21]发表的有关穴石蜈蚣属（Bothropolyssp.）线粒体基因组的分析中，有13个tRNA 的TΨC 臂或者DHU 臂无法形成正确结构. 此外，在左宝平[22]推测的少棘蜈蚣线粒体基因组的tRNA 二级结构中，存在碱基错配的现象，并且这种现象出现于多个tRNA 的氨基酸臂上，甚至哈氏蜈蚣的线粒体基因组还存在tRNA 缺失的现象[23]. 在节肢动物中，由于tRNA 的次级茎和环结构存在不对称性，也会出现某些tRNA 结构不符合典型的“三叶草”形状[24-25]，本研究中筛孔蜈蚣sp. 1 和日本蜈蚣也出现了这一现象，某些tRNA 无法形成正确二级结构.

线粒体基因组是研究生物系统发育的重要工具.本研究发现，长足卫蜈蚣、筛孔蜈蚣sp.1 和日本蜈蚣之间的亲缘关系依次降低，并且验证了地蜈蚣目和蜈蚣目之间是姊妹群的关系，这符合传统的形态分类学研究结果，同时还与Fernández 等[26]基于转录组中单拷贝基因数据建立的系统发育树相符.基于转录组的系统发育基因组学研究需要使用活体标本或液氮冷冻标本，而线粒体基因组数据可以从陈旧或馆藏标本中获得，极大降低了开展蜈蚣系统发育研究的门槛.

由于已公布的蜈蚣目物种线粒体基因组序列数量太少，利用不同属、不同科甚至不同目的线粒体基因组数据进行比较研究，建立完善的唇足纲动物系统发育树仍需要更多的基础测序数据.另外，少棘蜈蚣作为《中华人民共和国药典》规定收载的传统中药[27]，有很高的研究价值，其他6个物种（如多棘蜈蚣、哈氏蜈蚣等）在我国不同地区也存在药用习惯[28].通过补充与其亲缘关系相近的物种各方面的数据，也有利于增强对我国药用蜈蚣资源的掌握.