扬子鳃蛭线粒体基因组全序列测定及其系统学地位——基于Pacbio和Illumina技术*

谷生丽，刘谏君，孙恩涛，湛孝东，聂刘旺

(1.安徽师范大学 生命科学院，安徽 芜湖 241000；2.皖南医学院 医学寄生虫学教研室，安徽 芜湖 241002)

扬子鳃蛭(OzobranchusjantseanusOka，1912)是一种淡水吸血蛭类，隶属蛭纲(Hirudinea)吻蛭目(Rhynchobdellida)鳃蛭科(Ozobranchidae)鳃蛭属(Ozobranchus)[1]，常寄生于乌龟(Mauremysreevesii)的颈部及四肢皱褶区.目前仅日本和我国有关于扬子鳃蛭分布的报道[2-3].2014年，日本学者Suzuki[4]首次报道未经任何预处理的扬子鳃蛭经液氮冷冻24 h后，可100%复苏存活，且在-90 ℃的温度下可长期存活32个月，引起了研究者极大的兴趣.迄今，有关扬子鳃蛭线粒体基因组的研究甚少，Liu等[5]对扬子鳃蛭线粒体全序列数据进行了简短报道，但没有对其结构进行分析；Nakano等[6]基于线粒体COX1基因分析了扬子鳃蛭的遗传多样性.本研究采用Pacbio三代单分子测序技术以及Illumina二代测序技术对扬子鳃蛭线粒体基因组进行测序和注释，并与目前NCBI nucleotide数据库中已完成测序的其他10种蛭类的线粒体基因组全序列进行比较、分析，以期为深入了解扬子鳃蛭线粒体基因组的结构及其系统学地位提供基础数据.

1 材料与方法

1.1 材料收集及DNA提取

本实验用扬子鳃蛭采自安徽省无为市某龟类养殖场.取吸附于乌龟体表的扬子鳃蛭成虫1只，用去离子水冲洗，经SDS/蛋白酶K裂解，再使用酚/氯仿提取总DNA，1%琼脂糖凝胶电泳检测DNA完整性.

1.2 Illumina和Pacbio建库及测序

取1 μg DNA，用Covaris M220将DNA超声打断至300～500 bp的短片段，用TruSeqTMNano DNA Sample Prep Kit补平并在3’端加A，然后连接index接头，富集DNA后进行PCR扩增.利用2%琼脂糖回收目的条带，TBS380定量后按比例混合上机，使用Illumina Hiseq平台对文库进行测序.

参照文献[7],取1 μg DNA，用G-tubes方法将基因组处理成15～20 kb的片段，末端补平后，两端分别加上环状接头，最后使用Pacbio Sequel平台进行DNA测序.

1.3 数据质控及统计

参照文献[8]，Illumina测序的原始数据使用Trimmomatic去除接头及低质量序列.Pacbio测序的原始数据使用SMART Analysis将bam文件格式转化成fasta格式文件，并过滤到小于100 bp和质量小于0.8的Polymerase reads，去除adapter序列.

1.4 线粒体基因组组装及注释

参照文献[9]，利用ABySSV2.0组装Illumina测序数据，然后利用blasR比对Pacbio三代数据，根据比对结果对单分子测序数据进行一次矫正与纠错，最后利用纠正过的单分子测序数据与二代数据进行混合组装.使用SPAdes挑选覆盖深度足够高且组装长度较长的序列作为候选序列[10]，比对NT库进行确认.再次利用Illumina数据进行校验.得到完整的基因组数据之后，使用MITOS进行基因组的初步注释[11]，再参考GenBank数据库中相似物种的线粒体基因组数据进行最终注释，最后使用GCview绘制基因组图谱[12].

1.5 线粒体基因组数据比较分析

参照文献[13]，使用MAFFT version 7在线比对扬子鳃蛭的线粒体基因组序列与GenBank数据库中扬子鳃蛭线粒体基因组全序列(登录号KY861060.1)之间的差异.

登录NCBI nucleotide数据库，输入关键词“Hirudinea mitochondrion complete genome”，以检索蛭纲的线粒体基因组全序列.结果共检索到11个物种的线粒体基因组全序列数据，其中Poecilobdellamanillensis和Hirudinariamanillensis为同种异名，且基因序列完全相同，故有效种为10个.下载10个线粒体基因组序列，并加入本研究所获得扬子鳃蛭的线粒体基因组数据，进行蛭类线粒体基因组排序的比较分析.

为进一步确定扬子鳃蛭的分类地位，除下载以上全序列的10种蛭的线粒体基因组外，另加入了蛭纲中包含所有13个蛋白编码基因(PCGs)的不完整线粒体基因组序列，它们分别为吻蛭目舌蛭科(Glossiphoniidae)的Haementeriaofficinalis和Placobdellaparasitica，无吻蛭目医蛭科(Hirudinidae)的Hirudoverbana和Hirudomedicinalis，以及同属于环节动物门的其他外类群——寡毛纲的Lumbricusterrestris和Perionyxexcavatus及多毛纲的Orbinialatreillii的相应序列.参照文献[14]，利用MEGA6.0软件对18个物种的13个蛋白编码基因DNA序列构建Neighbor-Joining(NJ)树并进行分析.

2 结果

2.1 扬子鳃蛭线粒体基因组的组装

本研究采用Illumina和Pacbio平台对扬子鳃蛭的DNA进行测序，其中：Illumina测序得到原始数据4 432 Mb，质控之后剩余3 733 Mb，片段长度为450 bp；Pacbio得到过滤后Subreads数目为9 752，Subreads最大的长度为20 926 bp，Subreads的平均长度为3 834 bp.Pacbio的长度分布如图1所示.

图1 PacBio Subreads的长度分布Fig. 1 Length Distribution of PacBio Subreads

2.2 扬子鳃蛭线粒体基因组注释

经二代和三代数据组装获得扬子鳃蛭线粒体基因组序列长度为14 868 bp，A+T含量(72.4%)显著高于G+C含量(27.6%).所有基因都在同一条链，即J链(Majoritystrand)编码全部基因，包含13个蛋白编码基因、22个tRNA基因、2个rRNA基因及1个长非编码区(Long Non-coding Region， LNR)(图2).经过统计，扬子鳃蛭线粒体基因组中蛋白质编码基因占整个线粒体基因组长度的74.60%.

图2 扬子鳃蛭线粒体基因组结构Fig. 2 Schematic Map of the Mitochondrial Genome of O. jantseanus

本研究中扬子鳃蛭的编码基因共存在2种起始密码子，分别为ATG(ND5，ND4L，ND4，ND1，ND3，ND2，COX1，COX2，ATP8，ND6，CYTB，ATP6)和ATT(COX3)；终止密码子共有4种，分别为TAA(ND5，ND4L，ND1，COX2，ATP8，CYTB，ATP6)，TAG(ND3)，TA(ND4，ND6)和T(ND2，COX1，COX3).除最长非编码区为301 bp外，其他基因间隔共19处，最长为43 bp，其他间隔区总长度为162 bp.重叠区共6处，最大为-12 bp，总重叠区长度为-28 bp.具体结果详见表1.

表1 扬子鳃蛭线粒体基因组注释

2.3 基于线粒体mtDNA结构和13个蛋白编码基因DNA序列的蛭类系统发育分析

将NCBI核酸数据库中检索到的10种蛭纲物种的线粒体基因组全序列数据信息汇总，发现10种蛭纲物种均隶属于真蛭亚纲(Euhirudinea)的无吻蛭目(Arhynchobdellida)和吻蛭目(Rhynchobdellida)(表2).

表2 10种蛭纲物种的线粒体全序列基本数据

对目前已完成的10种蛭纲线粒体基因组的基因分布情况进行初步分析，结果见图3.

图3 蛭纲线粒体基因组基因序列差异Fig. 3 Differences in Gene Order in the Mitochondrial Genomes from the Hirudinea

如3所示，所有蛭的蛋白质编码基因在线粒体上的排列顺序均相同，但4个tRNAs有基因位置互换现象.蛭纲物种线粒体基因按照排序方式的不同可分为a和b 2种类型.a型：无吻蛭目中除日本石蛭(Erpobdellajaponica)外，均呈此排列方式，ATP8和COX3中间为tRNA(G)和tRNA(Y)，16S和ND1之间为tRNA(L1)，tRNA(S2)，tRNA(A)，tRNA(L2)，最长的非编码区(LNR)均小于100 bp.b型：所有吻蛭目物种如扬子鳃蛭以及无吻蛭目中的日本石蛭呈现此排列方式，即ATP8和COX3中间为tRNA(Y)和tRNA(G)，16S和ND1之间为tRNA(L1)，tRNA(A)，tRNA(S2)，tRNA(L2)，且LNR均大于299 bp.其中鱼蛭属(Zeylanicobdella)的Zeylanicobdellaarugamensis[15]的LNR达1 670 bp.

基于线粒体的13个蛋白编码基因DNA序列的蛭类系统发育分析结果如图4所示.

图4 基于线粒体13个蛋白编码基因DNA序列构建的蛭纲物种的NJ树Fig. 4 NJ Phylogenetic Tree of Hirudinea Species Based on 13 PCGs

由图4可知：本研究得到扬子鳃蛭与数据库中的扬子鳃蛭聚为一支，之后再与其他吻蛭目物种聚为一支；无吻蛭目中日本石蛭虽与其他无吻蛭目聚为一支，但其位于无吻蛭目最外侧，与亲缘关系较远且同为石蛭科的八目石蛭(Erpobdellaoctoculata)之间遗传距离也较远.

3 讨论

线粒体基因组测序对于物种的鉴定、系统发育以及生物的地理学研究具有重要的意义[16].目前发现环节动物的不同物种中蛋白编码基因相对保守.在无脊椎动物中进化较为高等的环节动物的基因顺序和方向具有高度保守性，这被认为是它们的原始特性[17-18].

本研究经Illumina二代测序技术和Pacbio三代单分子测序技术获得的扬子鳃蛭线粒体基因组全序列与GenBank中的Ozobranchusjantseanus(GenBank登录号KY861060.1)Blast比对相似度达到99.48%，全序列长度分别为14 868 bp和14 864 bp.通过MAFFT version 7软件比较2个序列发现，扬子鳃蛭线粒体基因组序列存在有8段gap，共14个位点，均位于非蛋白编码区；共有57个突变位点，其中11个位点涉及氨基酸编码的改变，ND4和ND6分别有2个，ND4L，ND1，COX1，ATP8无氨基酸变化，其余7个蛋白编码基因仅1个氨基酸的改变.另外，本研究扬子鳃蛭线粒体基因组全序列与GenBank中的Ozobranchusjantseanus在基因位点的划分方面也存在较多差异.如COX2的长度分别是684 bp和676 bp，后者缺少2个氨基酸和终止密码子；COX3分别是787 bp和781 bp，后者缺少2个氨基酸且无起始密码子.

本研究所采用的测序技术为Pacbio三代单分子测序技术和Illumina二代测序技术，该技术具有读取片段长和准确率高的特点[19-20].本研究表明，Subreads最大可达20 926 bp，所有候选序列长度均大于10 kb，可纠正步移法和Illumina二代测序可能产生的线粒体假基因[21]的偏差.通过以上可靠、高效的测序手段，笔者对扬子鳃蛭的线粒体基因组进行测序并组装，获得了高度精确的扬子鳃蛭mtDNA，并且对测序获得的扬子鳃蛭的线粒体基因组进行了较为详细的注释，通过软件分析结合近缘物种的比较，纠正了一些基因位点的不正确划分[5]，尤其是蛋白基因，更加准确地确定了扬子鳃蛭功能基因的位点.

笔者对10种蛭纲物种的线粒体基因排列进行比较分析并归类，首次发现蛭纲内物种的线粒体tRNA的排列方式有2种.有研究指出线粒体基因顺序在研究系统发生方面有着重要的作用，越相近的种相应的线粒体基因的顺序也越趋于一致[22-23].笔者归类的2种排列方式也与非编码区的长度以及基于线粒体的13个蛋白编码基因DNA序列所构建的NJ树的聚类较为一致，但是聚类结果与最近报道的通过进化分析发现的扬子鳃蛭聚类于无吻蛭目的研究结果[5]有根本的差别.本研究再次证明扬子鳃蛭无论在形态学分类还是在线粒体基因排列方式和NJ树进化分析上，都与吻蛭目物种聚为一类.另外，无吻蛭目中石蛭科的日本石蛭线粒体基因的排列方式为b型，在基于13个PCGs基因构建的NJ树分析中，日本石蛭与其他无吻蛭目以及同为石蛭科的八目石蛭的亲缘关系较远，说明石蛭科物种可能非单系起源，是否为有效种值得进一步分析研究.

目前，扬子鳃蛭作为蛭纲中报道的唯一能够耐受低温冷冻的蛭[4]，是蛭纲中的一个特例，但是关于扬子鳃蛭分子研究的报道仅有2篇[5-6].而非洲揺蚊幼虫和缓步虫作为耐低温和天体生物学研究中的模型生物[24-25]，相关研究报道涉及形态、生理、转录组和蛋白组等多方面[26-27].扬子鳃蛭与非洲揺蚊幼虫、缓步虫相比体积更大[1]，在今后的耐超低温以及冻后修复机制研究中具有更好的应用前景.在生物冷冻后的解冻过程中，组织会出现缺血缺氧和再灌注等损伤[28]，由此引起的一系列功能改变均与线粒体代谢密切相关[29].因此，对于扬子鳃蛭线粒体的后续研究以及蛋白组学和代谢组学的研究，可能会为人类探究低温休眠的生物钟调节机制，并最终实现整个生物体的冷冻保存及复苏提供帮助.