趙 芳，蔣 凱，蘇田娟，何 波，吳 清，林恭華，黃族豪

（1.井岡山大學生命科學學院，江西，吉安 343009；2.江西康之康中藥科技有限公司，江西，吉安 343009）

蛭類是環(huán)節(jié)動物門（Annelida）蛭綱（Hirudinea）動物的統(tǒng)稱，分布于除南極洲之外的各大洲[1]。全世界現存蛭類約680 種[2]，我國約有100 種[3]。一些蛭類物種以哺乳動物血液為食，能分泌抗凝血活性的物質，具有顯著的醫(yī)用價值和藥用價值，在蛭類研究中備受人們關注[4]。歐洲國家從17 世紀開始就廣泛使用蛭類治療各類炎癥和促進外科手術成功率，其中，最常用的物種是歐洲醫(yī)蛭（Hirudo medicinalis）[5]。在我國，蛭類是傳統(tǒng)醫(yī)藥中常用的破血、逐瘀、通經良藥[6]。據統(tǒng)計，我國以蛭類為原料藥材的中成藥有90 余種，其中通心絡膠囊、腦血康滴丸、腦心通膠囊等中成藥已被納入2019年的《國家基本醫(yī)療保險、工傷保險和生育保險藥品目錄》[7]。

現行《中國藥典》2020 年版收載的蛭類藥材種類為水蛭、螞蟥和柳葉螞蟥[8]，分別屬于日本醫(yī)蛭（Hirudo nipponica）、寬體金線蛭（Whitmania pigra）和尖細金線蛭（Whitmania acranulata）。不同蛭類物種之間習性不一，其抗凝血成分組成差異也很大[6]。為了保障臨床用藥安全，需要對蛭類物種進行有效鑒定[7]，而闡明不同物種之間的系統(tǒng)發(fā)育關系是達到這一目標的前提基礎。早期的研究以形態(tài)學、生理學和發(fā)育特征為主要分類依據[1]，隨著分子生物學的發(fā)展，以一個或少數幾個基因序列為標記的分類方法被廣泛應用到蛭類的分類工作中[9]。然而，由于蛭類的演化歷史極為復雜多變，現有的分類方法涉及的信息量較少，對一些類群難以有效劃分[10]。

醫(yī)蛭形亞目（Hirudiniformes）是蛭類中的代表類群，許多常見的物種如歐洲醫(yī)蛭、日本醫(yī)蛭、菲牛蛭（Hirudinaria manillensis）等，都隸屬于此亞目[1,3]。盡管醫(yī)蛭形亞目因其重要的經濟價值備受關注，然而，有關其分類學問題卻未得到有效解決。以醫(yī)蛭屬（Hirudo）和牛蛭屬（Hirudinaria）為例，有的學者認為兩者屬于親緣關系很近的姐妹群[11]，而另一些研究則發(fā)現兩者親緣關系很遠，甚至與黃蛭科（Haemopidae）或石蛭科（Erpobdellidae）的物種混雜，形成并系結構[12-13]。

近年來，高通量測序技術飛速發(fā)展，測序成本也直線下降。以基因組和轉錄組技術為代表的高通量測序分析技術，已廣泛應用于系統(tǒng)發(fā)育分析和分子進化分析研究[14-15]。目前，歐洲醫(yī)蛭[16]和菲牛蛭[17]基因組已被測序，也有一些物種如寬體金線蛭、棒紋牛蛭（Hirudinaria javanica）等的轉錄組相繼被測序[18-19]。盡管這些研究主要是分子進化分析，然而，所涉及的數據為系統(tǒng)發(fā)育分析提供了重要的基礎資料。本研究利用GenBank 上提交的醫(yī)蛭形亞目及其近緣類群的原始數據，開展系統(tǒng)發(fā)育基因組學研究，探討常見物種的進化地位，為這些物種的資源開發(fā)利用和保護提供理論依據。

1 材料與方法

1.1 數據來源

本研究涉及醫(yī)蛭形亞目的10 個蛭類物種，外加石蛭形亞目（Erpobdelliformes）和吻蛭目（Rhynchobdellida）各1 種作為外類群（表1），數據全部來源于GenBank。歐洲醫(yī)蛭、菲牛蛭和粗壯澤蛭（Helobdella robusta）的數據類型為已組裝的基因組，包含結構注釋信息和編碼區(qū)序列（coding sequence,CDS）。其中，歐洲醫(yī)蛭的CDS 序列從論文提供的鏈接中直接下載[16]；菲牛蛭基因組的發(fā)布者只提供基因組組裝序列（ GenBank No.GCA_015345955.1）和結構注釋（gff）文件[17]，需要利用gffread 軟件[20]進一步提取CDS 序列；粗壯澤蛭的CDS 序列直接從 GenBank 中下載（GCF_000326865.1）。其余物種GenBank 中只有轉錄組原始reads 數據，需要下載每個物種的Sequence Read Archive (SRA) 格 式 文 件， 用SRA-Toolkit 軟件包中的fastq-dump 軟件轉換成fastq 格式的reads 序列文件。用Trinity 軟件[21]對每個物種的reads 序列進行從頭組裝，然后用GeneMarkS-T 軟件[22]進行結構注釋，推測編碼區(qū)CDS 序列。

表1 本研究所用的蛭類基因組或轉錄組數據Table 1 Genomic or transcriptomic data of leeches used in this study

結合前人的研究，我們重點關注如下3 類基因：（1）COI（cytochrome c oxidase subunit I），即線粒體細胞色素c 氧化酶亞基I 基因。此基因變異速率適中，常被作為DNA 條碼序列鑒定物種[23]，頻繁應用于蛭類分子鑒定和系統(tǒng)發(fā)育研究。（2）18 S（18 S ribosomal DNA），即18 S 核糖體DNA。這個基因相對比較保守，但序列較長，常被應用于構建親緣關系較遠的系統(tǒng)發(fā)育關系[9]。（3）核基因編碼序列中的直系同源基因（orthologs），這類基因數量眾多，提供大量的變異信息位點，目前已成為系統(tǒng)發(fā)育基因組學的主要分析對象[14]。

1.2 分析方法

1.2.1 COI 和18 S 序列的提取和比對

COI 和18S 序列通過本地blast 的方式獲得。首先，從GenBank 中下載歐洲醫(yī)蛭的COI（AF003272）和18 S（AF116011）序列，并利用本地blast 軟件包NCBI-BLAST[24]中的makeblastdb 軟件構建每個基因組或轉錄組的本地blast 庫。以歐洲醫(yī)蛭序列為查詢序列，分別用NCBI-BLAST 軟件包中的tblastn和blastn 軟件，從本地blast 庫中檢索每個物種的COI 和18S 序列，同時確定匹配序列區(qū)段在基因組或轉錄組中的位置信息。用seqkit 軟件[25]的subseq工具，基于序列匹配信息，從基因組或轉錄組中提取各物種的目標序列及其上下游500 bp 區(qū)域。將每個基因的多個物種序列合并，用MEGA 軟件[26]分別進行比對，去除不匹配的區(qū)域。

1.2.2 核基因直系同源序列提取和比對

用cd-hit-est 軟件[27]去除每一個物種CDS 序列中的冗余序列，用seqkit 軟件中的fx2tab 工具統(tǒng)計所有非冗余CDS 序列（nonredundant CDS）的總體GC 含量。用OrthoFinder 軟件[28]提取所有物種的直系同源序列，用MACSE 軟件[29]對每個基因的直系同源序列進行比對，用Gblocks 軟件[30]去除比對質量較差的區(qū)段。最后，用seqkit 軟件中的concat 工具將所有比對好的基因按物種合并成多基因序列（命名為Orthologs 序列）。

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

用DAMBE 軟件[31]統(tǒng)計COI、18 S 和Orthologs序列的堿基組成、變異位點等信息，并進行飽和度檢驗。用IQ-TREE 軟件[32]對COI、18 S 和Orthologs序列分別進行最大似然法模型檢驗，基于貝葉斯信息量（Bayesian information criterion，BIC）選擇最適堿基替換模型?；谟米畲笏迫环嫿ㄏ到y(tǒng)發(fā)育樹，自舉檢驗次數為1000 次，支持率以百分數形式表示，最大值為100%。同時，用jModeltest 軟件[33]篩選可用于貝葉斯算法的最適堿基替代模型，用MrBayes 軟件[34]對COI、18 S 和Orthologs 序列構建貝葉斯樹，迭代次數為1000 萬次，支持率以小數形式表示，最大值為1。

2 結果

2.1 序列基本信息

從12 個物種的基因組或轉錄組數據中成功提取得到COI 和18 S 序列。比對后的COI 序列長度為1533 bp，變異位點數為747，占總位點數的48.7%；未發(fā)現插入/缺失位點。比對后的18 S 序列長度為1832 bp，變異位點數為385，占總位點數的21.0%；其中，插入/缺失位點數為119，占總位點數的6.5%。經提取和比對得到345 個直系同源序列，合并后的Orthologs 總長為220,095 bp，總變異位點數為112,240，占總長度的51%。

COI 序列A、C、G、T 四種堿基的含量(χ±SD)%分別為(29.0±1.3)%、(15.2±2.4)%、(15.8±1.5)%、(40.0%±3.0)%；G 和C 堿基含量占比（GC%）為(31.0±3.2)%，具有明顯的低GC 特征。18 S 序列A、C、G、T 四種堿基的含量分別為(24.8±0.3)% 、 (22.6±0.3)% 、 (28.0±0.4)% 、(24.6%±0.4)%；GC%為(50.6±0.6)%，堿基組成比較均勻。Orthologs 序列A、C、G、T 四種堿基的含量分別為（30.2±0.4）%、（20.7±0.4）%、(23.7±0.4)%、(25.3±0.4)%；GC%為(44.5±0.8)%，有少量的低GC 特征。此外，基于所有非冗余CDS序列的統(tǒng)計顯示，這些物種的GC 含量為(45.1±1.7)%，與Orthologs 序列的GC 含量最為接近。每個物種的GC 含量信息在表2 中列出。

表2 蛭類COI、18 S、Orthologs 及所有非冗余CDS 序列的GC 含量信息Table 2 GC contents of COI,18S,Orthologs and all nonredundant CDS sequences of leeches

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

飽和度檢驗分析顯示，COI 的Iss 和Iss.c 值分別為0.2753 和0.7898，18S 的Iss 和Iss.c 值分別為0.0971 和0.7986，Orthologs 的Iss 和Iss.c 值分別為0.2317 和0.8447。這3 類序列的Iss

IQ-TREE 軟件自帶的堿基替換模型分析顯示，COI、18S 和Orthologs 的最適堿基替換模型分別是TIM+F+I+G4、TIM3e+R2 和GTR+F+R3。jModeltest軟件的模型檢驗顯示，COI、18S 和Orthologs 的可用于 MrBayes 軟件的最適模型分別是TPM1uf+I+G、TrNef+G 和GTR+I+G。

基于COI 序列的系統(tǒng)發(fā)育分析顯示，IQ-TREE軟件構建的最大似然樹和MrBayes軟件構建的貝葉斯樹拓撲結構完全一致，但枝長和部分節(jié)點的支持率有明顯差別（圖1）?；?8S 序列的最大似然樹和貝葉斯樹之間在拓撲結構和支持率上有明顯差異，但枝長比較一致（圖2）?；贠rthologs 序列的最大似然樹和貝葉斯樹之間在拓撲結構、枝長和支持率上都高度一致，僅1 個節(jié)點（歐洲醫(yī)蛭和東方醫(yī)蛭交匯節(jié)點）上的支持率數值有少量差異（圖3）。

圖1 基于COI 序列的最大似然樹（A）和貝葉斯樹（B）Fig.1 Maximum likelihood tree(A)and Bayesian tree(B)based on COI sequences

圖2 基于18 S 序列的最大似然樹（A）和貝葉斯樹（B）Fig.2 Maximum likelihood tree(A)and Bayesian tree(B)based on 18 S sequences

基于Orthologs 序列的兩種系統(tǒng)發(fā)育樹在幾乎所有節(jié)點上支持率都高達100%，因此，本研究主要以此結果探討物種之間的系統(tǒng)發(fā)育關系。研究發(fā)現，歐洲醫(yī)蛭、東方醫(yī)蛭和側紋醫(yī)蛭組成一個單系群，寬體金線蛭、尖細金線蛭和日本醫(yī)蛭組成另一單系群，而馬蛭與上述6 個物種形成姐妹群。菲牛蛭和棒紋牛蛭組成一個單系群，與上述7 個物種形成姐妹群。山蛭科的洞穴山蛭在本研究所涉及的10個醫(yī)蛭形亞目種類中處于最基部位置。

3 討論

本研究利用GenBank 上公開發(fā)表的基因組或轉錄組數據，探討醫(yī)蛭形亞目代表類群的系統(tǒng)發(fā)育特征。與傳統(tǒng)的單個基因相比，本研究所用的Orthologs 序列在構建系統(tǒng)發(fā)育關系方面具有明顯優(yōu)勢。由于Orthologs 序列長度遠大于COI 或18 S，提供了大量的變異位點信息，使絕大部分節(jié)點的支持率達到100%，增加了結論的可靠性。同時，Orthologs 序列是從整個基因組或轉錄組中獲得，更能代表整個物種的遺傳特征。例如Orthologs 序列的GC 含量為(44.5±0.8)%，與蛭類所有CDS 的GC含量(45.1±1.7)%最為接近，可見Orthologs 序列充分體現了各物種的遺傳屬性，提高結論的可信度。

蛭類的演化歷史極其復雜，許多類群的分類地位存在爭議[1]。歐洲醫(yī)蛭和日本醫(yī)蛭是醫(yī)藥學領域最具代表性的物種。歐洲醫(yī)蛭主要分布在歐洲及北非地區(qū)，應用歷史長達4 個世紀[5]。盡管在醫(yī)學和生物學等領域受到廣泛關注，然而有關其分類地位的問題目前并未得到有效解決。最近的分子系統(tǒng)學研究發(fā)現，一些長期被認為是歐洲醫(yī)蛭變種的群體，其實是多個獨立的物種，包括側紋醫(yī)蛭、東方醫(yī)蛭和北非醫(yī)蛭（Hirudo troctina）[35]。這些物種之間在染色體數量、唾液成分和抗凝血基因組成等方面已出現明顯的分化[36]，盡管如此，一些物種之間（尤其是歐洲醫(yī)蛭、側紋醫(yī)蛭和東方醫(yī)蛭）的系統(tǒng)發(fā)育關系仍然存在爭議。例如，基于4 個基因（ITS2、5.8 S rRNA、12 S rRNA 和COI）的系統(tǒng)發(fā)育分析顯示，歐洲醫(yī)蛭和東方醫(yī)蛭的親緣關系最近[37]，而另一項基于4 個基因（18 S rDNA、28S rDNA、12 S rDNA 和COI）的研究則認為側紋醫(yī)蛭和東方醫(yī)蛭的親緣關系最近[38]。本研究中，COI 和18 S 序列支持側紋醫(yī)蛭和東方醫(yī)蛭親緣關系最近的觀點，然而Orthologs 的結果卻顯示，歐洲醫(yī)蛭和東方醫(yī)蛭的親緣關系最近。由于COI 和18 S 的系統(tǒng)發(fā)育樹在相關節(jié)點上的支持率都低于70%，而Orthologs 的系統(tǒng)發(fā)育樹在相關的支持率分別為88%和100%，因此本研究支持Trontelj and Utevsky 的研究[37]，即歐洲醫(yī)蛭和東方醫(yī)蛭的親緣關系最近。

日本醫(yī)蛭在我國廣泛分布，在日本、蒙古和俄羅斯遠東地區(qū)也有分布。日本醫(yī)蛭是傳統(tǒng)中藥水蛭的主要類群，作為活血化瘀類動物藥材的代表藥，具有破血、逐瘀、通經的功效[39]。寬體金線蛭和尖細金線蛭屬于黃蛭科下的物種[3]，盡管不吸食哺乳動物的血液，但具有一定的抗凝血活性和抗血小板聚集等作用[40-42]，因此也被《中國藥典》2020 年版收載[8]。前人的研究發(fā)現，醫(yī)蛭科和黃蛭科可能并非兩個獨立的單系群，而是具有相互交叉的并系關系[18,39]。我們的研究再次證實這個結論，即，日本醫(yī)蛭和這兩種金線蛭的親緣關系更近，而與其它幾種醫(yī)蛭科物種的親緣關系更遠。相反，同屬黃蛭科的馬蛭，卻與寬體金線蛭系統(tǒng)發(fā)育關系相去甚遠。我們認為，醫(yī)蛭科和黃蛭科的分類屬性需要重新修訂。這個結果也解釋了為什么不吸食血液的金線蛭也具有抗凝和抗栓活性——金線蛭屬的祖先也是吸血類型，現有的抗凝血和抗栓能力可能是其進化上的殘留性狀。

菲牛蛭隸屬醫(yī)蛭科牛蛭屬，因模式產地在菲律賓呂宋島而得名。菲牛蛭分布于東南亞國家包括菲律賓、馬來西亞、越南等，以及我國部分省份如廣西、海南、福建等[3]。菲牛蛭也是一味傳統(tǒng)中藥，其抗凝血活性明顯優(yōu)于《中國藥典》所收錄的3 個物種（日本醫(yī)蛭、寬體金線蛭和尖細金線蛭），被作為地方特色藥物收錄于《廣西壯族自治區(qū)壯藥質量標準》[43]和《云南省中藥材標準》[44]。當前的分類體系將菲牛蛭等牛蛭屬物種置于醫(yī)蛭科下，然而，本研究發(fā)現，牛蛭屬與醫(yī)蛭科醫(yī)蛭屬的分化程度甚至超過黃蛭科的馬蛭，這與前人的研究結論不一致[18,39]。因此，目前關于醫(yī)蛭科和黃蛭科的分類關系是混亂的，有必要對相關物種進行深入研究，進而對整個醫(yī)蛭形亞目的分類系統(tǒng)進行修訂。此外，山蛭科的洞穴山蛭在本研究所涉及的醫(yī)蛭形亞目種類中處于最基部位置，此結果與經典的Sawyer分類系統(tǒng)[1]一致。需要指出的是，作為中藥成分的動植物資源，其有效成分含量在不同物種甚至不同品種[45]之間都有可能出現明顯差別。本研究的結果不僅有利于闡明醫(yī)蛭形亞目的分類系統(tǒng)，同時，有助于蛭類中醫(yī)藥資源的開發(fā)、利用和保護。