3個(gè)蜈蚣科物種線粒體基因組測序及唇足綱的系統(tǒng)發(fā)育分析

潘陽洋，陳天韻，尤春雪，蔣 超

（1.天津農(nóng)學(xué)院動(dòng)物科學(xué)與動(dòng)物醫(yī)學(xué)學(xué)院天津市農(nóng)業(yè)動(dòng)物繁育與健康養(yǎng)殖重點(diǎn)實(shí)驗(yàn)室，天津 300384；2.中國中醫(yī)科學(xué)院中藥資源中心道地藥材品質(zhì)保障與資源持續(xù)利用全國重點(diǎn)實(shí)驗(yàn)室，北京 100700）

蜈蚣是唇足綱蜈蚣目動(dòng)物的統(tǒng)稱，目前已記錄5科27 屬.作為節(jié)肢動(dòng)物進(jìn)化的重要環(huán)節(jié)，蜈蚣的遺傳分類對物種進(jìn)化研究具有重要意義.目前已有許多基于單基因序列或多基因聯(lián)合分析蜈蚣進(jìn)化、系統(tǒng)發(fā)育和分類鑒定的報(bào)道，尤其是基于DNA 分子系統(tǒng)發(fā)育的研究[1].隨著二代DNA 測序技術(shù)的快速發(fā)展，越來越多的學(xué)者使用線粒體基因組（mitochondrial genome）研究物種分類、生物體的遺傳發(fā)育和生物系統(tǒng)發(fā)育[2].區(qū)別于核基因，線粒體基因組是母系遺傳，其結(jié)構(gòu)簡單，易于分析[3]. 線粒體基因組通常是1個(gè)環(huán)狀雙鏈DNA 分子，包括13個(gè)蛋白質(zhì)編碼基因（PCG）、22個(gè)轉(zhuǎn)移RNA 基因（tRNA）、2個(gè)核糖體RNA 基因（rRNA）和1個(gè)控制區(qū)（CR）[4-5].雖然線粒體基因組的結(jié)構(gòu)相對保守，近年來許多研究卻發(fā)現(xiàn)，在蜈蚣的線粒體基因組中存在基因重排現(xiàn)象，并且伴隨有tRNA 二級結(jié)構(gòu)的缺失等變化[6].目前已公布的唇足綱物種的線粒體基因組序列僅有14個(gè)，在蜈蚣目物種的線粒體基因組序列中，有4個(gè)屬于蜈蚣屬（Scolopendra），1個(gè)屬于尖盲蜈蚣屬（Scolopocryptops），缺乏科內(nèi)其他屬的數(shù)據(jù)，因此難以進(jìn)行跨種屬線粒體結(jié)構(gòu)分析.

本研究對蜈蚣科的3個(gè)物種，即衛(wèi)蜈蚣屬長足衛(wèi)蜈蚣（Rhysida longipes（Newport，1845））、篩孔蜈蚣屬篩孔蜈蚣sp.1（Ethmostigmussp.1）和蜈蚣屬日本蜈蚣（Scolopendra japonicaL.Koch，1878）進(jìn)行全長線粒體基因組測序，分析其線粒體基因組的基因重排模式，并依據(jù)13個(gè)PCGs 進(jìn)行唇足綱的系統(tǒng)發(fā)育分析，為蜈蚣的形態(tài)學(xué)分類提供數(shù)據(jù)支持，并有助于了解蜈蚣目獨(dú)有的基因重排模式.

1 材料與方法

1.1 樣品采集和DNA 提取

長足衛(wèi)蜈蚣（CMMI：20200511008），采集于云南省紅河哈尼彝族自治州；篩孔蜈蚣sp. 1（CMMI：20191212012），采集于云南省德宏傣族景頗族自治州盈江縣太平鎮(zhèn)；日本蜈蚣（CMMI：20190904004），采集于廣西壯族自治區(qū)來賓市金秀瑤族自治縣.樣本保存于70%～75%的乙醇中.

取每個(gè)樣本的3 條足，干燥后研磨成粉末放入2.0 mL 的離心管中，然后使用Promega WizardRSV Genomic DNA Purification kit 試劑盒（美國Promega 公司）提取蜈蚣DNA，DNA 樣品儲(chǔ)存于-20 ℃冰箱待用.

1.2 線粒體基因組的測序和組裝

用InvitrogenCollibriNGS建立基因組DNA文庫，用illuminaHiSeq2500 平臺(tái)，按照PE150bp 規(guī)程進(jìn)行測序，由北京諾禾致源科技股份有限公司完成.每個(gè)物種的線粒體基因組數(shù)據(jù)為5×109～6×109.用GetOrganelle[7]軟件進(jìn)行序列組裝，獲得線粒體基因組序列.

1.3 線粒體基因組的注釋

應(yīng)用Mitos2[8]注釋tRNAs、rRNA 及PCGs 的位置，對于未注釋上的tRNA，應(yīng)用Phylosuite 1.2.2[9]軟件并利用近緣物種比對的方式進(jìn)行預(yù)測.

通過rtools[10]網(wǎng)站（http：//rtools.cbrc.jp/cgi-bin/index.cgi）預(yù)測tRNA 的二級結(jié)構(gòu)，NCBI 網(wǎng)站Conserved Domains[11-14]（https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi）驗(yàn)證13個(gè)PCGs 的編碼區(qū)域并進(jìn)行調(diào)整，OGDRAW[15]網(wǎng)站（https：//chlorobox.mpimp-golm.mpg.de/OGDraw.html）繪制線粒體基因組圈圖.

應(yīng)用AT-Skew 和GC-Skew[16]評估所有線粒體基因組的鏈不對稱性，應(yīng)用CodonW[17]軟件計(jì)算密碼子偏好性（RSCU）.

式中：A、T、G、C 均為堿基.

1.4 系統(tǒng)發(fā)育分析

為判斷線粒體基因組數(shù)據(jù)與傳統(tǒng)的形態(tài)學(xué)系統(tǒng)發(fā)育關(guān)系之間的聯(lián)系，選取3種蜈蚣和已公布的唇足綱共14個(gè)物種（表1）的全長線粒體基因組序列建立系統(tǒng)發(fā)育樹，并且選取蚰蜒目的蚰蜒（Scutigera coleoptrata（Linnaeus，1758），AJ507061.2）作為外類群.利用CLUSTAL W[18]在BIOEDIT 7.1.3.0[19]中比對13個(gè)PCGs.用Phylosuite 1.2.2[9]軟件分別構(gòu)建貝葉斯推斷法（Bayesian inference，BI）和最大似然法（maximum likelihood，ML）的系統(tǒng)發(fā)育樹，并且利用ModelFinder[20]功能選擇合適的模型，最大似然法的模型為GTR+F+R3，貝葉斯推斷法的模型為GTR+F+I+G4. 使用IQTREE 構(gòu)建最大似然樹，采用超快自展（ultrafast bootstraping）法重復(fù)100 000 次計(jì)算分支置信度；應(yīng)用Mr-Bayes 構(gòu)建貝葉斯樹，利用4 條蒙特卡羅馬爾可夫鏈（Markov Chain Monte Carlo，MCMC）同時(shí)運(yùn)行1 000 000 代，抽樣頻率為1 000，摒棄25%的老化樣本.

表1 本研究中的3種蜈蚣和已發(fā)布的11個(gè)唇足綱物種的線粒體基因組數(shù)據(jù)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 結(jié)果與分析

通過高通量測序技術(shù)、序列組裝和基因組注釋，獲得長足衛(wèi)蜈蚣、篩孔蜈蚣sp.1 和日本蜈蚣的全長線粒體基因組，均為環(huán)狀DNA，其序列長度分別為14 992、15 002 和14 425 bp.由表1 可以看出，除石蜈蚣目的長角石蜈蚣（Cermatobius longicornis（Takakuwa，1939））線粒體基因組長度為16 833 bp 外，大部分的唇足綱物種線粒體基因組長度為14 000～16 000 bp.

2.1 線粒體基因組結(jié)構(gòu)及組成分析

3種蜈蚣的線粒體基因組均包含37個(gè)基因，包括13個(gè)PCGs、22個(gè)tRNA 基因和2個(gè)rRNA 基因，如圖1 所示.蜈蚣目物種線粒體基因組中，不同PCGs 的長度差異顯著，13個(gè)PCGs 長度最小的為ATP8，8個(gè)物種的平均值為158 bp，最大的為NAD5，平均值為1 622 bp，如圖2 所示.多數(shù)PCGs 在不同蜈蚣物種之間長度存在差異，如日本蜈蚣ATP6長度最大，達(dá)696 bp，而篩孔蜈蚣sp.1 的ATP6長度最小，僅361 bp.同其他PCGs 相比，ATP8的長度最穩(wěn)定.13個(gè)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 含量，表現(xiàn)為AT 富集，各物種的GC 含量在21.22%（少棘蜈蚣）～32.52%（篩孔蜈蚣sp.1）之間，其中，篩孔蜈蚣sp.1 的GC 含量與長足衛(wèi)蜈蚣（32.27%）的幾乎相等.

圖3 8個(gè)蜈蚣目物種線粒體基因組的核苷酸組成百分比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個(gè)蜈蚣目物種中，線粒體基因組RSCU>1 的密碼子有26～33個(gè)，而在長足衛(wèi)蜈蚣、篩孔蜈蚣sp.1 和日本蜈蚣中，有26～31個(gè)密碼子被高度使用.其中，篩孔蜈蚣sp.1 和長足衛(wèi)蜈蚣的重疊密碼子共有28個(gè)，而3種蜈蚣的重疊密碼子為18個(gè).除此之外，長足衛(wèi)蜈蚣和日本蜈蚣還重疊1個(gè)密碼子，即GUA.另外，蜈蚣的氨基酸RSCU 值均以2 作為分水嶺，因此將氨基酸RSCU>2 作為篩選條件，可以得到3種蜈蚣高度使用的氨基酸是丙氨酸（Ala）、精氨酸（Arg）、甘氨酸（Gly）、異亮氨酸（Ile）、亮氨酸（Leu（taa）和（tag））、脯氨酸（Pro）、絲氨酸（Ser（tag））、蘇氨酸（Thr）、纈氨酸（Val）.這10個(gè)氨基酸中，除了赤蜈蚣的亮氨酸Leu（taa）的RSCU 低于2 之外，其他氨基酸在8個(gè)蜈蚣目物種中也是被高度使用的.

圖4 8個(gè)蜈蚣目物種線粒體基因組的氨基酸同義密碼子相對使用度Fig.4 RSCU of mitochondrial genomes of eight species from Scolopendromorpha

在8個(gè)蜈蚣目物種中，13個(gè)PCGs 和2個(gè)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的轉(zhuǎn)錄方向與其他序列相反.少棘蜈蚣的排列順序與其他蜈蚣不一致，為COX1-COX2-ATP8-ATP6-COX3-NAD3-NAD1-rrnL-rrnS-NAD5-NAD4-NAD4L-NAD6-COB-NAD2，說明NAD1-rrnL-rrnS基因簇發(fā)生了一次整體易位.

在3種蜈蚣的線粒體基因組的基因排序中，PCGs的排列基本穩(wěn)定，重排現(xiàn)象主要出現(xiàn)在tRNA 的排列中，與其他序列的基因排序相比，trnQ、trnY、trnE的排序較為不固定，并且在rrnL和rrnS中穿插有trnV.

另外，長足衛(wèi)蜈蚣的14個(gè)tRNA 基因（trnI、trnM、trnW、trnK、trnD、trnG、trnA、trnR、trnN、trnS（gct）、trnE、trnT、trnS（tga）、trnL（tag））和8個(gè)PCGs（COX1、COX2、COX3、NAD2、NAD3、NAD6、ATP6、ATP8）都在正鏈轉(zhuǎn)錄，其余8個(gè)tRNA 基因、5個(gè)PCGs 和2個(gè)rRNA 基因（rrnS、rrnL）在負(fù)鏈轉(zhuǎn)錄. 相比之下，日本蜈蚣的線粒體基因組的基因大多從正鏈轉(zhuǎn)錄，為18個(gè)tRNA基因（trnI、trnM、trnW、trnL（tag）、trnL（taa）、trnK、trnD、trnG、trnA、trnR、trnN、trnS（tga）、trnS（gct）、trnH、trnT、trnY、trnV、trnE） 和8個(gè)PCGs（NAD2、NAD3、NAD6、COX1、COX2、COX3、ATP6、ATP8），其余4個(gè)tRNA 基因、5個(gè)PCGs 和2個(gè)rRNA 基因（rrnS、rrnL）在負(fù)鏈轉(zhuǎn)錄.

此外，在篩孔蜈蚣sp. 1 基因組中，trnR、trnP、trnS、trnV缺少TΨC 臂或者DHU 臂，并且日本蜈蚣同樣也有trnR、trnG、trnA、trnH、trnT、trnL（taa）缺少TΨC 臂或DHU 臂的現(xiàn)象.

2.3 唇足綱的系統(tǒng)發(fā)育關(guān)系

本研究構(gòu)建的ML 和BI 系統(tǒng)發(fā)育樹均高度支持（PP=1，BS=100）蜈蚣科分為蜈蚣亞科和耳孔蜈蚣亞科，如圖6 所示，由圖6 可以看出，在蜈蚣亞科蜈蚣屬中，哈氏蜈蚣、少棘蜈蚣和海南蜈蚣親緣關(guān)系很近，并且與日本蜈蚣互為姊妹群（PP=1，BS=100），與傳統(tǒng)的形態(tài)學(xué)分類結(jié)果相符.此外，系統(tǒng)發(fā)育樹顯示，地蜈蚣目和蜈蚣目互為姊妹群.地蜈蚣目和蜈蚣目均為整形亞綱分類群，符合傳統(tǒng)唇足綱系統(tǒng)發(fā)育的研究結(jié)果.

3 討論與結(jié)論

本研究首次對蜈蚣科長足衛(wèi)蜈蚣、篩孔蜈蚣sp.1和日本蜈蚣的線粒體基因組進(jìn)行測序、組裝和注釋，確定了22個(gè)tRNAs、13個(gè)PCGs 和2個(gè)rRNAs 的位置，并且發(fā)現(xiàn)大多數(shù)基因位于正鏈上. 在利用Mitos2網(wǎng)站對篩孔蜈蚣sp.1 進(jìn)行注釋時(shí)，發(fā)現(xiàn)PCGs 有大量的重疊區(qū)域，因此有必要采用近緣物種比對的方法進(jìn)行準(zhǔn)確注釋.與蜈蚣目其他物種相比，篩孔蜈蚣sp.1 中NAD5、NAD4、NAD4L、NAD1、rrnL、rrnS等多個(gè)PCGs轉(zhuǎn)錄方向相反，這種情況是否是篩孔蜈蚣屬的共有特征有待采集該屬其他物種進(jìn)行進(jìn)一步研究.

與PCGs 的排序相比，蜈蚣tRNA 的重排現(xiàn)象更頻繁，預(yù)測的tRNA 無法形成正確二級結(jié)構(gòu)的現(xiàn)象也很普遍. 在Park 等[21]發(fā)表的有關(guān)穴石蜈蚣屬（Bothropolyssp.）線粒體基因組的分析中，有13個(gè)tRNA 的TΨC 臂或者DHU 臂無法形成正確結(jié)構(gòu). 此外，在左寶平[22]推測的少棘蜈蚣線粒體基因組的tRNA 二級結(jié)構(gòu)中，存在堿基錯(cuò)配的現(xiàn)象，并且這種現(xiàn)象出現(xiàn)于多個(gè)tRNA 的氨基酸臂上，甚至哈氏蜈蚣的線粒體基因組還存在tRNA 缺失的現(xiàn)象[23]. 在節(jié)肢動(dòng)物中，由于tRNA 的次級莖和環(huán)結(jié)構(gòu)存在不對稱性，也會(huì)出現(xiàn)某些tRNA 結(jié)構(gòu)不符合典型的“三葉草”形狀[24-25]，本研究中篩孔蜈蚣sp. 1 和日本蜈蚣也出現(xiàn)了這一現(xiàn)象，某些tRNA 無法形成正確二級結(jié)構(gòu).

線粒體基因組是研究生物系統(tǒng)發(fā)育的重要工具.本研究發(fā)現(xiàn)，長足衛(wèi)蜈蚣、篩孔蜈蚣sp.1 和日本蜈蚣之間的親緣關(guān)系依次降低，并且驗(yàn)證了地蜈蚣目和蜈蚣目之間是姊妹群的關(guān)系，這符合傳統(tǒng)的形態(tài)分類學(xué)研究結(jié)果，同時(shí)還與Fernández 等[26]基于轉(zhuǎn)錄組中單拷貝基因數(shù)據(jù)建立的系統(tǒng)發(fā)育樹相符.基于轉(zhuǎn)錄組的系統(tǒng)發(fā)育基因組學(xué)研究需要使用活體標(biāo)本或液氮冷凍標(biāo)本，而線粒體基因組數(shù)據(jù)可以從陳舊或館藏標(biāo)本中獲得，極大降低了開展蜈蚣系統(tǒng)發(fā)育研究的門檻.

由于已公布的蜈蚣目物種線粒體基因組序列數(shù)量太少，利用不同屬、不同科甚至不同目的線粒體基因組數(shù)據(jù)進(jìn)行比較研究，建立完善的唇足綱動(dòng)物系統(tǒng)發(fā)育樹仍需要更多的基礎(chǔ)測序數(shù)據(jù).另外，少棘蜈蚣作為《中華人民共和國藥典》規(guī)定收載的傳統(tǒng)中藥[27]，有很高的研究價(jià)值，其他6個(gè)物種（如多棘蜈蚣、哈氏蜈蚣等）在我國不同地區(qū)也存在藥用習(xí)慣[28].通過補(bǔ)充與其親緣關(guān)系相近的物種各方面的數(shù)據(jù)，也有利于增強(qiáng)對我國藥用蜈蚣資源的掌握.