線粒體Cyt b基因對九種穿山甲鑒別力的分析

關鍵詞：穿山甲；細胞色素b基因（ Cyt b）；遺傳距離；物種鑒定

穿山甲主要分布于亞洲、非洲的熱帶及亞熱帶地區，屬鱗甲目（Pholidota）哺乳動物，含穿山甲屬（Manis）、地穿山甲屬（Smutsia）和長尾穿山甲屬（Phataginus）3屬。目前穿山甲公認有8個種，分別是印度穿山甲（Manis crassicaudata）、菲律賓穿山甲（M. culionensis）、馬來穿山甲（M. javanica）、中華穿山甲（M. pentadactyla）、長尾穿山甲（Phataginus tetra?dactyla）、樹穿山甲（Ph. tricuspis）、大穿山甲（Smutsiagigantea）和南非穿山甲（S. temminckii）［1?2］。2023年，Gu et al.［3］通過對海關查獲的穿山甲進行群體基因組學分析，發現了一個新的物種，命名為神秘穿山甲（Manis mysteria），使其成為第9 種穿山甲或第5 種亞洲穿山甲，但其分布地區目前仍不清楚。

已知的8種穿山甲在國際自然保護聯盟（IUCN）紅色名錄中分別被定義為易危（VU）、瀕危（EN）或極危（CR）［4］。近年來，大量穿山甲制品自非洲、東南亞通過走私流入亞洲和歐美市場，其中穿山甲甲片是海關查獲最多的產品［5?6］。在打擊走私過程中，對查獲的甲片進行物種鑒定是其關鍵環節。但是，穿山甲屬內的物種分化程度較低，如菲律賓穿山甲曾一度被認為是馬來穿山甲，而最新發現的神秘穿山甲，過去也始終認為是馬來穿山甲。主要原因是這兩種穿山甲的形態與馬來穿山甲非常相近，通過甲片形態更是難以鑒定到種。因此，進行DNA鑒定就成為當前的最佳選擇。

甲片雖然是完全角質化的材料，但也含有少量的DNA。其中mtDNA占有優勢，能夠通過PCR技術擴增特定的片段［7?8］，進而可以測定其序列。mtDNA中的細胞色素b（cytochrome b，Cyt b）基因進化速率與物種分化速率大體相當，含有豐富的物種識別信息，被廣泛應用于包括穿山甲在內的物種鑒定［9?11］。迄今為止，穿山甲物種鑒定的方法仍存在局限性，如常用的PCR擴增片段都在300 bp以上，這對于高度片段化的甲片DNA 來說過長，經常得不到擴增產物。縮短目的片段雖然可以提高擴增和測序的成功率，但減少了物種信息含量，提高了誤判的風險，需要妥善解決該問題，才能建立起適合各類穿山甲樣本的物種鑒定方法。

物種鑒別主要基于這樣一個假設：同一物種的個體之間親緣關系比不同物種的個體之間更近。親緣關系可以用遺傳距離來衡量，兩兩個體之間的遺傳距離可以通過mtDNA序列計算。在衡量一段DNA序列能否可靠地判定物種時，一般通過比較同種個體間遺傳距離（intra-specific distances，intra-ds）的分布與異種個體間遺傳距離（inter-specific distances，inter-ds）的分布，如果二者的范圍有明顯的差異，重疊很少，則表明這一序列具有良好的鑒別力，反之鑒別力較低［12］。通常而言，各部分基因序列的進化速度并不均勻，選擇目的片段不當會增加物種判別的不確定性，而加大擴增片段的大小又影響其對高度降解DNA 的適用性。另一方面，遺傳距離也受到計算方法的影響，對物種信息解析較好的方法才能得到更有利于物種判別的遺傳距離［13?14］。

本研究以線粒體Cyt b 基因序列為對象，探究適合穿山甲物種鑒定的最佳遺傳距離計算方法，評估并篩選鑒別力較優的區段，為建立可靠的穿山甲物種鑒定方法提供參考。

1 材料與方法

1. 1 穿山甲Cyt b 基因序列

從NCBI的Nuleotide Database（https：//www. ncbi.nlm. nih. gov/nuccore）下載9 種穿山甲54 條線粒體Cyt b 基因的完整序列（1 140 bp），詳見表1。根據王辰等［15］的研究結果對其中2個錯誤標注的序列重新歸類：JN411577歸入馬來穿山甲，NC004027歸入樹穿山甲。

1. 2 遺傳距離計算模型有效性的測試

將54 條穿山甲Cyt b 序列導入MEGA 11［16］軟件，使用23種方式計算兩兩序列間的遺傳距離：在估計方差（estimate variance）欄中選擇自助抽樣法（bootstrap method）并迭代1 000次；在替代模型（substitutionmodel）欄中分別選擇使用Nucleotide和Syn-Nonsynonymous 兩種替代類型進行計算，其中在Nucleotide 選項下使用p-distance、Kimura 2-parametermodel（K2P）［17］、maximum composite likelihood（MCL）［18］和LogDet（Tamura-Kumar）［19］4 種模型及其下替代選項中Transitions（TS）、Transversions（TV）、TS + TV、TS / TV 的4種替換類型和算法進行計算；在Syn-Nonsynonymous 選項下使用標準遺傳密碼表，Nei-Gojobori method（NG）［20］、modified Nei-Gojoborimethod（Jukes-Cantor，JC）［21］和Kumar method（Kimura2-para）［22］ 3 種模型及其下替代選項中Synonymousonly（dS）、Nonsynonymous only（dN）、Difference（dSdN）、At 0-fold degenerate sites only（0-fold deg. sites）和4-fold degenerate sites only（4-fold deg. sites）5種算法計算；其他選項和參數均使用默認值。

1. 3 不同區段有效性的滑窗測試

使用Python腳本將54條Cyt b 序列切分為步長均為50，滑窗長度為100 bp 的序列22 組及長度為200 bp的序列20組（表2），將42組序列數據按物種分組后導入MEGA 11軟件。綜合先期測試數據和其他文獻研究［23］結果，選擇分辨率最優的K2P模型和TV值作為遺傳距離的計算方法。

1. 4 數據分析

將遺傳距離計算結果導出至Excel，使用Excel函數分別統計每個遺傳距離矩陣中種內遺傳距離（intra-ds）和種間遺傳距離（inter-ds）的平均值、標準差和95% 置信區間。同時計算inter-ds 的下限與intra-ds的上限間的差值。該差值可以定義為物種分辨率，當其大于0時，表示鑒別力較好；小于0時表示鑒別力不足。由于本研究涉及不同模型間跨坐標軸比較遺傳距離，且在遺傳距離計算中不總是呈現intra-ds在左（較小）、inter-ds在右（較大）的分布格局（如計算TS / TV時），故采用置信區間比率（R）反映物種分辨率。

R 值越小對應相對物種分辨率越大，鑒別力越好，當R 小于100%時，說明鑒別力極佳；等于100%時，說明鑒別力較好；大于100% 時則鑒別力不足。使用Origin 2021計算intra-ds和inter-ds頻率分布并繪圖。

使用MEGA 11 軟件基于鄰接法（neighborjoiningmethod）［24］繪制聚類樹，選擇K2P方法及TV堿基替換類型進行計算，自助法迭代1 000次以檢驗各分支置信度，其他選項和參數使用默認值。

2 結果

2. 1 不同算法模型對遺傳距離計算結果的影響

使用核苷酸類模型計算基于9種穿山甲Cyt b 全長序列的intra-ds和inter-ds，結果見表3及圖1。總體來說，各種方法計算的intra-ds和inter-ds置信區間基本無交集，但是二者分布都有所重疊，說明總體上的物種分辨率不夠高。在K2P、MCL和p-distance三種模型中，TS / TV 作為堿基替換模式時intra-ds 和inter-ds幾乎完全重疊，沒有物種分辨力。用K2P和p-distance兩種模型計算TV時，其intra-ds置信上限與inter-ds置信下限比率為8%，物種分辨率最好。其他方法置信比率在17% ～ 27%，相對差距不顯著，其中MCL 模型計算的物種分辨率較為穩定（R = 17%），LogDet算法分辨率一般（R = 20%）。

使用Kumar、JC 和NG 三種Syn-Nonsynonymous遺傳距離模型計算9種穿山甲的intra-ds和inter-ds結果如表4和圖2所示。各方法計算的遺傳距離均呈現出intra-ds在左（較小）、inter-ds在右（較大）的分布格局。Kumar模型的5種不同計算方法間的置信區間比率為14% ～ 17%，是3種模型中整體表現最好的，其中dN 的比率最低（R = 14%），物種分辨率最高。Kumar與NG模型的dN置信比率（14%、15%）均低于dS置信比率（17%、23%），證明本數據集中dN算法能更好地對種內種間進行區分。與上述2種模型相比，JC 模型的遺傳距離置信比率明顯偏大（35% ～ 39%）。

綜合來看，采用Cyt b 全序列鑒別9種穿山甲時，JC模型得到的物種鑒別力較低，另外所有TS / TV替換模式物種分辨力均比較有限。相比之下，K2P、pdistance兩種模型效果在采用TV替換類型時最好，可得到較高的物種判別力。

2. 2 Cyt b 基因不同區段對穿山甲鑒定的有效性

使用滑窗對54條1 140 bp的穿山甲Cyt b 基因序列進行切分，使用K2P + TV的方式獲取遺傳距離計算結果，研究不同區段對intra-ds和inter-ds差值的影響。滑窗長度分別為100、200 bp，步長均為50 bp的遺傳距離置信區間分布情況見圖3，其中橫坐標的起始位點是指片段首個位點在Cyt b 基因片段上的位置，各起始位點對應的實際序列片段見表2。整體遺傳距離置信差值呈現為200 bp長度的片段更為穩定、浮動較小；更短的100 bp片段浮動較大且出現更高遺傳距離差值，其4個波峰分別出現在堿基序位101 ～ 200號、301 ～ 400號、551 ～ 650號和851 ～ 950號，其物種分辨率分別是Cyt b 全長序列差值（0. 05）的1. 00、1. 14、1. 20、1. 22倍。100 bp片段遺傳距離置信區間上下限差值的最大值（851 ～ 950號堿基序位，0. 061）高于200 bp 片段的最大值（851 ～ 1050 號，0. 052），且100 bp 片段差值最小值（751 ～ 850 號，0. 025）低于200 bp 片段最小值（651 ～ 850 號，0. 032）。綜上證明遺傳距離差值數據與序列片段長度無直接關系，且序列片段的選擇將在一定程度上影響遺傳距離計算結果。

由圖3可見，沿Cyt b 基因，大部分片段的物種分辨率都低于全長Cyt b 基因的分辨率（0. 05），如在100 bp的滑窗中檢測到751 ～ 850 bp序列最為保守，其物種分辨率僅為全長序列差值的50%。100 bp片段鑒別力高但穩定性較低，200 bp片段與之相反，因此可推測使用部分100 bp片段組成不連續的200 bp序列能夠提供一種更理想的方法。從100 bp片段中提取4個置信差值位于波峰的片段，即堿基序位為101 ～ 200、301 ～ 400、551 ～ 650和851 ～ 950號堿基的4個片段。統計不同片段組合對遺傳距離的影響（圖4），片段組合的鑒別力隨原片段鑒別力增加而遞增，相比Cyt b 全長序列，片段組合的鑒別力提升至1. 08 ～ 1. 24倍，其中（551 ～ 650）+（851 ～ 950）號堿基片段的組合獲得了最高置信區間差值0. 062，1. 24倍于全長序列物種鑒別力，且高于單個100 bp片段的鑒別力最大值（0. 061），說明相較于使用連續的200 bp片段，由高鑒別力離散片段組合而成的序列能夠最大程度提升鑒別能力且提供穩定結果。

為檢驗上述得到的物種分辨率較優的片段在實際鑒別中的作用，采用鄰接法及K2P + TV算法繪制9種穿山甲系統樹，代表性結果見圖5。堿基序位在第101 ～ 200號堿基和301 ～ 400號堿基的2個小片段都可以單獨把各個物種區分開，而（101 ～ 200）+（301 ～ 400）號堿基及（301 ～ 400）+（851 ～ 950）號堿基2 個片段組合的表現最好。其他結果表明在自401號堿基開始至1 050號堿基的區域內，馬來穿山甲與菲律賓穿山甲間的遺傳距離均為0，在551 ～650號堿基片段中神秘穿山甲與前二者的遺傳距離也為0，因此如需鑒定相關穿山甲應盡可能避免僅使用上述序列。另外，原標注為中華穿山甲的NC016008應為馬來穿山甲。值得說明的是，不同區段上突變的特點和頻率有差異，故此其系統樹的拓撲形狀與實際的物種樹存在差異。這些結果主要反映各個片段的物種鑒別力，即可否把近緣物種有效分開，而不探討是否符合實際的系統發育關系。

3 討論

雖然法醫動物學已經進入基因組時代，但線粒體Cyt b 基因片段分析簡單，結果可靠，耗費計算資源少，對樣本質量要求寬泛，尤其是分析成本低，適合構建物種大數據，是鑒定實踐中開展物種鑒定的主要選擇［25］。通過Cyt b 序列解析物種間的差異有很多方法，如在GenBank數據庫中通過BLAST搜索尋找序列最為相似的物種記錄［26］，通過宏觀的遺傳距離估算，基于同種個體間遺傳距離小于異種間遺傳距離的基本前提來確定給定樣品的來源物種，也有通過序列的相似性或遺傳距離與多個可疑物種構建系統發育樹來觀察所屬分支等［13］。遺傳距離法需要滿足種內遺傳距離最大值小于種間遺傳距離最小值的方法判別物種［27］。因此基于遺傳距離的物種鑒定方法有兩個問題必須探討：一是何種方法計算的遺傳距離最能展現出intra-ds和inter-ds的這種差異，以盡可能降低錯誤判斷物種的概率；二是使用Cyt b哪個（些）區段才能獲得最佳的物種分辨率。

針對這兩個問題，本研究通過9種穿山甲Cyt b基因序列開展了系統分析和評估。首先，關于遺傳距離的估算方法，測試了23 種估算方法，包括p-distance、Kimura 2-parameter model（K2P）、maximumcomposite likelihood（MCL）、LogDet（Tamura-Kumar）、Nei-Gojobori method（NG） 、modified Nei-Gojoborimethod（Jukes-Cantor，JC）和Kumar method（Kimura 2-para）等距離估算模型，以及Transitions（TS）、Transversions（TV）、TS + TV、TS / TV、Synonymous only（dS）、Nonsynonymous only（dN）、Difference（dS-dN）、At 0-fold degenerate sites only（0-fold deg. sites）和4-folddegenerate sites only（4-fold deg. sites）九種堿基替代模式和類型。相比傳統觀念中置信差值評估物種間遺傳距離的方式，采用置信比率方式更能直接反映不同模型生成的遺傳距離矩陣中intra-ds和inter-ds置信區間在橫坐標軸上的相對距離。在可能的組合（表3，表4）中，Nucleotide類模型的TS / TV 替換模式估算的距離完全沒有物種分辨力，相比之下，K2P、p-distance兩種模型效果在采用TV 替換類型時置信區間比率為8%，其物種鑒別力最高。在Syn-Nonsynonymous 類中Kumar模型整體的鑒別力較好（14% ～ 17%），而JC 模型得到的物種鑒別力較低（35% ～ 39%）。因此，建議在采用穿山甲全長Cyt b序列時，盡量避免使用JC 模型和TS / TV 堿基替換模式。

前述研究結果顯示，在穿山甲的物種鑒定中，Cyt b 全長序列分辨率比較有限。因此，采用距離估算效果較好的K2P + TV 方法評估Cyt b 各個區段的物種分辨率。沿Cyt b 基因，大部分片段的物種分辨率都低于全長Cyt b基因的分辨率（0. 05），如在100 bp的滑窗中檢測到751 ～ 850 號堿基間序列最為保守，其物種分辨率僅為全長序列差值的50%；但另一些100 bp的短片段也具有較高的鑒別力，如堿基序位在851 ～ 950號的片段，其物種分辨率是Cyt b 全序列的1. 22倍；類似的片段還有堿基序位為551 ～ 650號、301 ～ 400 號，物種分辨率分別是Cyt b 全序列的1. 20 倍和1. 14 倍。因此，如果采用多片段復合擴增［28］和擴增子高通量測序技術，把多個這樣的片段同時進行分析，能夠去掉很多噪聲，提高鑒別力。由于自401號堿基序位后無法鑒別部分物種，必須避免使用只包含該部分堿基的片段，故而物種分辨率達全長序列的1. 2倍的（301 ～ 400）+（851 ～ 950）號堿基序位的片段組合為最佳理論選擇。基于遺傳距離的系統發育樹（圖5）也直觀地支持了這一結論，這為設計更有效的穿山甲鑒定方法提供了依據。已發表的幾種使用短片段的穿山甲鑒定方法中Kocher et al.和Yeo et al. 設計的引物［29?30］得到的長度為307 bp（96 ～ 402號堿基序位）的序列范圍包含兩段高分辨率區域，推測鑒定效果近似或略優于Cyt b 全長序列；L14724 / H15149引物［31?32］得到的長度為402 bp（1 ～ 402 號堿基序位）的序列包含0 ～ 100 及200 ～300號堿基兩段低鑒別力區域，推測分辨率不及上述方法；王辰等［15］設計的一種產物為131 bp（70 ～ 200號堿基序位）長度的引物基本位于低分辨率區域，鑒別力有限。