魚類環(huán)境DNA metabarcoding片段的近緣物種識別差異

陳治 ，馬春來 ，葉樂 ，楊超杰 ，王海山 *

( 1.海南熱帶海洋學院 熱帶海洋生物資源利用與保護教育部重點實驗室，海南 三亞 572022；2.海南熱帶海洋學院 海南省熱帶海洋漁業(yè)資源保護與利用重點實驗室，海南 三亞 572022)

1 引言

近年來，環(huán)境 DNA（Environmental DNA，eDNA）metabarcoding技術開始廣泛應用于漁業(yè)生態(tài)學領域，在物種多樣性調查、外來魚種監(jiān)測、瀕危魚類保護方面展現了極大的應用價值和發(fā)展?jié)摿1-2]。其中，metabarcoding片段的選擇對該技術的有效使用至關重要[2-3]。經過眾多研究者的不懈努力，目前已篩選出一些物種識別率較高、對應引物擴增能力較強的優(yōu)質片段[3]。如Miya等[4]基于880種硬骨魚類線粒體12S序列設計的MiFish-U，擴增片段長度控制在163～185 bp范圍內，對沖繩水族館的魚類檢出數目高達230種；Vences等[5]針對線粒體16S設計的脊椎動物Vert-16S，不僅擴增片段長度更長（約260 bp），而且對魚類、兩棲類均具有較好的適用性—以德國境內的兩條河流為驗證水域時，物種檢出率高達82%～93%；Balasingham等[6]設計的PS1，將擴增片段選定于線粒體COI基因上，在獲得較多魚類可操作分類單元（Operational Taxonomic Units，OTUs）的同時，可高效利用公共數據庫內的參考序列進行物種注釋。截至2020年，被報道的魚類eDNA metabarcoding片段已多達數十個，其中表現優(yōu)異的片段約10余個[3]。各片段的篩選背景不同，實際表現也各具特色，這為豐富魚類eDNA metabarcoding庫、推動魚類多樣性研究奠定了堅實基礎。

然而，迄今為止的魚類eDNA研究尚不健全，如何根據研究水域和目標類群合理選擇metabarcoding片段仍處于探索階段[3]。其中，metabarcoding片段對近緣物種的適用性就是一個亟需解決的問題。物種識別率和引物通用性是選擇魚類eDNA metabarcoding片段的兩大核心指標[3-4]。為了同時兼顧上述兩方面需求，主流的eDNA metabarcoding片段（AcMDB07、Ac12S、MiFish-U、Tele02等）主要位于線粒體核糖體（12S、16S）高變區(qū)[3]，雖然能基本滿足魚類多樣性調查需求，但卻存在著少數近緣魚類因條形碼片段變異速率較低而難以識別的風險?；谝延形锓N監(jiān)測結果可知，采用核糖體區(qū)片段進行魚類eDNA調查可能會造成“物種丟失”。如MiFish-U對平鲉屬（Sebastes）、東方鲀屬（Takifugu）、蝴蝶魚屬（Chaetodon）[7]，Tele01對玉筋魚屬（Ammodytes）、杜父魚屬（Cottus）、雅羅魚屬（Leuciscus）[8]等類群的OTUs注釋結果均出現了此問題。表明“物種丟失”可能是一種比較常見的現象。然而，由于本底資料缺乏、物種數量不足等原因，已知的魚類eDNA調查只是簡單提及了此現象，未能充分揭露所用片段的近緣物種識別缺陷；不同片段間的近緣物種識別差異更缺乏相互比較。究竟哪些片段的“物種丟失”風險更低，更適合近緣魚類多樣性調查？目前尚不得而知。

隨著水體污染、過度捕撈、棲息地喪失等威脅的加劇，世界各地土著魚類的物種多樣性、遺傳多樣性正急劇下降[9-11]；而外來魚種卻因水產養(yǎng)殖、觀賞活動釋放等原因大幅增加[11]。這就對當前的魚類多樣性監(jiān)測提出了更高的要求。而全面、準確的物種調查則是開展各類生物多樣性保護的基礎和前提[1-3]。在近緣魚類識別是魚類eDNA調查無法回避的大背景下，本研究將系統探究主流metabarcoding片段的“物種丟失”風險、比較不同片段的近緣物種識別差異，以期能夠從中篩選出物種識別率最高的片段，最終為完善魚類eDNA metabarcoding技術、提高魚類多樣性調查結果的準確性提供一定技術支撐。

2 材料與方法

2.1 序列下載及篩選

自然水域難以出現近緣魚類大規(guī)模出現的情況，結果具有較大的偶然性。因此，本研究采用MitoFish公共數據庫內的線粒體全序比較各片段的近緣物種識別差異。具體操作如下：登錄MitoFish數據庫（http://mitofish.aori.u-tokyo.ac.jp/download.html），下 載“All Mitogenomes”選項內的3 081種魚類的線粒體全序fasta壓縮文件（截至2022年1月23日）；同步下載txt文檔，以方便核對序列來源、GenBank號、上傳者單位等信息。使用MegAlign軟件將下載的fasta文件轉換為seq文件。基于Fishbase數據庫（https://www.fishbase.de/）核對物種分類有效性。剔除雜交種、爭議種及含有簡并堿基的序列。為降低序列截齊的工作量，本研究僅將物種較多（5種及以上）的屬挑出用于后續(xù)分析。物種分類階元的確認依據Fishbase數據庫及中國動物主題數據庫（China Animal Scientific Database，http://www.zoology.csdb.cn/）。

2.2 metabarcoding片段選擇

查閱2008?2021年期間涉及魚類eDNA metabarcoding技術的有關研究，特別是參考Zhang等[3]的研究結果對其中表現較好的metabarcoding片段及對應引物進行匯總。以花鰻鱺（Anguilla marmorata）線粒體全序（序列號：NC006540）為參照，核對正、反向引物方向標注是否正確。按所在位置對metabarcoding片段進行編號、排序。

2.3 序列截齊

使用DNAstar軟件包中的Seqman程序，對挑選出的序列進行截齊：（1）擴增片段序列截齊后用于近緣魚類條形碼相似度分析；（2）正、反向引物序列截齊后用于引物變異度分析。

2.4 遺傳距離分析

采用PAUP 4.0軟件構建最大簡約法（Maximum Parsimony，MP）系統發(fā)育樹，系統樹可信度均采用Bootstrap檢驗，經1 000次重復抽樣檢驗得到分支樹節(jié)點支持率。物種的聚類情況即物種識別率[12]：若不同物種聚為一支（取遺傳距離等于0），認為該物種識別失敗；若不同物種聚為不同支（遺傳距離大于0），認為該物種識別成功（由于不同metabarcoding片段的種間差異標準不同，難以給出準確可靠的判定閾值，因此本研究只探究遺傳距離是否為0的情況）；同時，基于正、反向引物聯合序列計算總平均遺傳距離[7]，從而判斷引物的通用程度。

運用SPSS25.0 軟件基于單因素方差分析（Oneway ANOVA）對不同引物的物種識別率等進行差異顯著性分析。顯著差異閾值為p=0.05，試驗結果以“平均值±標準差”表示。

2.5 非度量多維尺度分析

使用非度量多維尺度（Non-Metric Multidimensional Scaling，NMDS）[13]方法對物種識別結果 [0,1]矩陣進行排序分析，該分析在Canoco 5.0軟件中完成。選用Bray-Curtis距離，因為其在大范圍和小范圍的坐標軸上都具有穩(wěn)健性[12,14]。分析結果以脅強系數（stress）作為評判標準：當stress＜0.2時，認為可以用NMDS的二維點圖表示，該圖形有一定的解釋意義；當stress＜0.1時，認為該排序是一個好的排序；而當stress＜0.05時，則認為該排序結果具有很好的代表性[15]。NMDS運算的步驟如下：（1）將物種識別率處理為[0,1]矩陣，鑒定成功為1，鑒定失敗為0；以物種名為首行、metabarcoding片段名為首列儲存.xlsx格式的[0,1]矩陣；（2）打開Canoco 5.0軟件，加載.xlsx文檔，給定table name，選擇“import all species as factors”，其余選擇默認參數，完成矩陣導入；（3）“analysis”選項下，選擇“Canoco Adviser”，調出 NMDS 分析程序；（4）“treatment of ties in distance”選擇“secondary”， “stress formula”選擇“type 1”，其余選擇默認參數。

3 結果

3.1 用于分析的近緣魚類簡介

本研究用于分析的魚類共計2綱、20目、52科、106屬、935種。其中鯉形目（Cypriniformes）種類最多，共計3科、37屬、337種，鱸形目（Perciformes）次之，共計 18科、22屬、190種。扁鯊目（Squatiniformes）、脂鯉目（Characiformes）、胡瓜魚目（Osmeriformes）的數量最少，各只有1科、1屬、5種。物種名稱及對應Genbank序列號見表S1、表S2。

3.2 引物簡介

經文獻查閱及匯總，特別是參照Zhang等[3]對22個主流片段（引物）的全面比較結果，本研究用于分析的魚類eDNA metabarcoding片段共計15個。涵蓋線粒體12S、16S、COI 3個目標基因。其中以12S為目標基因的片段最多（8個），16S次之（6個），COI最少（僅1個）。各片段、對應引物位置及信息見圖1及表1。

表1 本研究中15個魚類eDNA metabarcoding片段簡介Table 1 Summary of 15 fish eDNA metabarcoding fragments analyzed in this study

圖1 本研究中15個魚類eDNA metabarcoding片段和對應引物在目標基因上的位置Fig.1 Locations of the 15 fish eDNA metabarcoding fragments and primer pairs on the target mitochondrial genes

3.3 物種識別率及引物變異情況

不同片段的物種識別率存在極顯著差異（F=32.39，p=1.5×10?76），蛋白質編碼基因（片段 15）的物種識別率最高，達到（82.76±24.66）%。核糖體基因中，片段03、片段05的表現最為優(yōu)異，識別率分別為（68.9±30.81）%和（71.62±27.3）%；片段 04、片段 08和片段11的識別率較低，平均值未超過42.92%；其余片段的物種識別率較為中等，平均值在54.58%～64.6%間（圖2，表2）。

基于15個片段對應正、反向引物聯合序列構建的MP系統發(fā)育樹總平均遺傳距離依次為1.52%、1.84%、4.58%、0.88%、0.97%、2.46%、5.78%、3.32%、9.32%、2.29%、6.21%、4.14%、1.71%、3.17%和15.70%（圖2）。表明蛋白質編碼基因（片段15）的對應引物通用性最差。片段09、片段11、片段07、片段03和片段12的引物變異程度也較高（總平均遺傳距離≥4.14%），實際使用過程中應考慮其eDNA擴增效率。

圖2 本研究中15個魚類eDNA metabarcoding片段的物種識別率及對應引物序列總平均遺傳距離Fig.2 Fish species resolution rates of 15 eDNA metabarcoding fragments and overall mean distances of primer pairs in this study

不同類群的識別率也存在極顯著差異（F=14.40,p=7×10?157）：旗鳉屬（Aphyosemion）、舌鰨屬（Cynoglossus）、鏢鱸屬（Etheostoma）等屬內近緣物種較容易區(qū)分，全部15個片段的平均識別率均在90%以上；而金槍魚屬（Thunnus） 、白鮭屬（Coregonus）等類群的目標基因同源性較高，大部分片段對其識別能力較差（圖3，表2）。

圖3 本研究中106屬魚類的物種識別率Fig.3 Fish species resolution rates of 106 genera in this study

表2 metabarcoding片段對各類群的識別率Table 2 The species resolution rates of metabarcoding fragments for different groups

3.4 NMDS結果

根據935種魚類的識別結果對15個metabarcoding片段做NMDS排序，片段04、片段08和片段11明顯遠離其余片段，說明此3個片段物種識別結果較為不同（圖4）。其余11個核糖體片段根據所在目標基因聚類成A、B兩組，表明不同目標基因影響物種識別。NMDS分析結果的脅強系數為0.13，說明將15個metabarcoding片段劃歸不同組具有一定的解釋意義。

4 討論

本研究用于分析的魚類數量龐大、階元眾多（見3.1節(jié)及表S1、表S2），涉及軟骨魚類及硬骨魚類、淡水魚類及海水魚類、定居性魚類及洄游性魚類等各類群。樣本量具有較好的統計學意義，結果能夠比較充分地反映15個主流eDNA metabarcoding片段的優(yōu)劣。本研究是目前國內第一篇針對性探究魚類eDNA metabarcoding近緣物種適用性的報告。

4.1 COI片段的引物缺陷

片段15（PS1，COI）的物種識別率顯著高于剩余14個12S/16S片段，是唯一物種識別率超過80%的metabarcoding片段。這與Balasingham等[6]認為的魚類線粒體蛋白質編碼基因比核糖體基因能夠更有效避免北美五大湖內部分入侵物種eDNA注釋結果假陰性的結論一致，也與經典的DNA條形碼普遍采用線粒體細胞色素C氧化酶亞基I基因的事實相符[24]。此外，COI基因在魚類eDNA metabarcoding參考數據庫的全面性方面也具有12S、16S基因難以比擬的優(yōu)勢[3]。然而，引物序列的高變性可能會阻礙片段15的廣泛應用。本研究中該片段正、反向引物序列的總平均遺傳距離高達15.70%。中國團扇鰩（Platyrhina sinen-sis）、細條銀口天竺鯛（Jaydia lineata）、土佐鰧（Uranoscopus tosae）、橫帶髭鯛（Hapalogenys analis）的 COI特異性引物、探針設計結果表明：當總平均遺傳距離大于8.21%～9.65%時，高純度DNA模板會出現明顯的彌散擴增；總平均遺傳距離大于12.47%～13.77%時，引物會對部分近緣物種產生絕對特異[7]。因此，片段15存在eDNA低效擴增/偏倚擴增的缺陷。Collins等[25]、Zhang等[3]分別對英吉利海、北京水系的魚類多樣性調查顯示，COI片段的物種檢出數目不及所用的12S片段，其結果支持本研究觀點。

此外，其他的COI或Cytb metabarcoding片段可能也難以解決引物通用性的問題。魚類eDNA metabarcoding技術需要依托 2×150 bp、2×250 bp兩種高通量測序（Next-generation Sequencing，NGS）平臺實現[3]，而Miya等[4]對880種硬骨魚類線粒體基因組全序的篩查結果顯示，COI、Cytb等蛋白質編碼基因在該讀長限定范圍內不存在成對的側翼保守區(qū)，理論上難以設計出與經典的魚類DNA條形碼（基于Sanger測序，讀長可達2×600 bp）媲美的成對引物。Menning等[26]、Jennings等[27]、Sultana等[28]的序列比對結果也表明，整個Cytb基因基本為高變序列，COI基因則只在近5′端前350 bp范圍內變異較小，相對適合設計短片段通用引物。結合本研究（圖1）和Zhang等[3]對22對引物的匯總結果可知，片段15等少數幾對魚類COI metabarcoding優(yōu)質片段恰位于此相對保守區(qū)，其引物通用性較差，而其他區(qū)域metabarcoding片段的引物序列總平均遺傳距離可能更大。因此，本研究認為后續(xù)研究沒有必要進一步比較其他COI、Cytb片段的近緣物種識別差異；同時認為以片段15為代表的蛋白質編碼片段更適合作為eDNA研究的輔助metabarcoding，并且需要根據研究水域、目標類群降低引物調查范圍，從而對引物序列進行針對性設計或優(yōu)化。

4.2 核糖體片段的比較

片段長度影響物種識別（表1，圖2，圖4）。本研究中識別率最低的3個metabarcoding片段（片段08、片段11和片段04）長度均不超過110 bp，而片段05、片段03等高識別率片段則在200 bp以上。Balasingham等[6]、Gantner等[29]的研究也認為12S和16S擴增片段越短，物種鑒定的準確性越低。雖然由于水體中的痕量eDNA通常存在嚴重降解，主流的觀點認為200 bp以內的微條形碼可能具有更高的PCR成功率[4,30]，但Bylemans等[19]、Zhang等[3]多個研究卻發(fā)現：片段05、片段03等的高通量測序數據量和物種檢出數目并不低于片段08、片段04或片段11，只是在定量分析等方面相互存在較大差異。魚類eDNA的產生、降解動力學過程極為復雜，目前沒有研究能夠準確闡明eDNA的微觀存在、分布和變化規(guī)律，Deiner等[31]、Bylemans等[19]推測長片段eDNA可能反而可以更長久地存在水體中。本研究不推薦使用metabarcoding短片段（片段08、片段11和片段04）進行近緣魚類多樣性調查。

圖4 不同metabarcoding片段的非度量多維尺度（NMDS）分析Fig.4 Analysis of non-metric multidimensional scaling (NMDS) for different metabarcoding fragments

續(xù)表2

續(xù)表2

續(xù)表2

續(xù)表2

除去COI和3個短片段，片段05的物種識別率最高，其引物通用性（總遺傳距離=0.97%）也僅次于片段04。片段05是近緣魚類多樣性調查的第一選擇，該結果與Zhang等[3]的研究結論一致。片段09、片段07、片段03、片段12的物種識別率雖然也較高，但在本研究中，其引物序列總平均遺傳距離較大（≥4.14%），存在類似于片段15的“物種丟失”風險；其他metabarcoding引物比較研究也表明其物種檢出數目不及片段06、片段13[3]。片段14、片段10雖然長度大于200 bp，但物種識別率和引物通用性卻并不出眾；并且可能因為引物序列3′位置存在變異[7]，這兩個片段在個別研究中的實際表現甚至不如片段15[3]。而片段01和片段02則在引物序列、條形碼性能方面基本一致，二者存在明顯的相互替代性，本研究更傾向二者中選用片段01。綜上所述，片段05、片段06、片段01、片段13是本研究篩選出的4個最優(yōu)metabarcoding片段。然而，受高通量測序成本和參考數據庫的影響，在單一研究中使用的metabarcoding片段通常不超過3個[3-4]。雖然本研究及Zhang等[3]、Bylemans等[19]的研究均表明片段06的物種識別率和引物通用性優(yōu)異，但NMDS分析（圖4）顯示片段06在散點圖上與片段05位置較近，因此同時使用這兩個片段進行魚類多樣性研究可能并不能顯著提高物種檢出數目。NMDS分析也顯示，不同基因、同一基因不同片段間均存在較大差異（圖4）。因此不僅需要多片段聯合應用，而且所用片段要有足夠的多基因代表性。本研究傾向于以片段05、片段01為主，片段13等為輔，進行近緣魚類多樣性調查。

4.3 本研究所用片段的不足及新片段開發(fā)的可能性

本研究對106屬魚類的識別結果差異極顯著（見3.3節(jié)），表明物種類群直接影響eDNA metabarcoding調查效果。從15個片段的總識別率角度比較（表2），僅15屬魚類物種識別率超過90%，多達38屬魚類識別率不足50%。由于935種魚類的線粒體全序列皆來自公共數據庫，其中一些序列可能存在同物異名錯誤—如鯧屬魚類的銀鯧（Pampus argenteus）與鐮鯧（Pampus echinogaster）[32]，因此該屬魚類識別率偏低；但更可能的解釋是低識別類群物種間遺傳差異較小，導致metabarcoding片段從源頭上就難以完全區(qū)分—如白鮭屬（Coregonus）[33]、紅點鮭屬（Salvelinus）[34]的眾多物種集中分布于50°N以北環(huán)北極圈水域，物種間可能存在較為頻繁的基因交流；平鲉屬（Sebastes）則廣泛分布于北太平洋且仔魚營漂浮生活，不同群體間也不存在明顯的地理隔離[35]。Miya等[4]在篩選Mi-Fish-U（片段01）引物時已發(fā)現，金槍魚屬等大洋性魚類的大部分線粒體基因高度保守，只有借助變異速率極高的線粒體NADH脫氫酶亞基5（NADH Dehydrogenase Subunit 5，ND5）條形碼才能進行有效區(qū)分。因此，受制于物種自身的遺傳背景，本研究所用片段尚無法對全部魚類進行100%區(qū)分。后續(xù)仍需不斷嘗試篩選新片段從而盡量提高eDNA metabarcoding技術的魚類多樣性調查能力。

metabarcoding片段的開發(fā)與目標基因內莖環(huán)結構的數量和分布情況密切相關。斑鰶（Konosirus punctatus）[36]、大口鰜（Psettodes erumei）[37]、北極茴魚（Thymallus arcticus）[38]等線粒體基因二級結構示意圖顯示，16S基因僅在中后段存在1～2個大型莖環(huán)結構，相比之下，12S基因存在2～4個。由于12S基因的總長度僅為16S的60%（約900 bp : 1 600 bp），因此該基因內的莖環(huán)結構分布更為密集和均勻。莖環(huán)結構有利于尋找穩(wěn)定的“保守?高變?保守”序列，Evans和Lamberti[39]、H?nfling 等[40]、Bylemans等[19]據此認為，線粒體12S基因比16S基因更適合開發(fā)魚類eDNA metabarcoding標記，本研究支持此結論。然而，基于14個核糖體片段的位置示意圖（圖1）和簡約信息（表1）可知，魚類16S基因的后半段大型莖環(huán)及12S基因的主要莖環(huán)皆已經被篩選出metabarcoding片段。圍繞單一莖環(huán)結構開發(fā)新片段不僅可能性極小，而且片段識別率可能也不理想（如片段07、片段10）。相比片段01、片段04和片段08等，片段05和片段06等則是由多個中小型莖環(huán)結構組合形成，其物種識別率反而更高。以此類推，聯合多個相鄰莖環(huán)區(qū)域可能是形成新的高識別率片段的手段之一。本研究中，12S片段對應引物存在明顯的位置重合及序列共用（圖1，表1），說明這些區(qū)域的引物保守性得到了較大認可。因此，可以嘗試通過這些引物的組合應用開發(fā)更長、更高識別率的metabarcoding片段。如片段01、片段02的正向引物與片段06、片段07正向引物的反向互補序列共用，PCR產物約420 bp（圖1）。該片段由1個大型莖環(huán)加2個小型莖環(huán)結構組成，其物種識別率可能高于片段05。

綜合物種識別率、引物通用性等多方面因素，本研究推薦2×150 bp測序平臺使用片段01（MiFish-U）、2×250 bp 測序平臺使用片段 05（Ac12S），輔以片段13（Vert-16S-eDNA）進行近緣魚類多樣性調查；新片段的開發(fā)也是不斷完善eDNA metabarcoding技術的重要工作。然而，魚類eDNA研究受多方面因素的影響，基于935種近緣魚類線粒體全序的比較結果可能會與實際表現有差異。后期還需尋找近緣魚類廣泛分布的水域，在建立龐大的本底資料數據庫的基礎上驗證本研究的結論和推測。

補充材料

表S1 本研究中106屬魚類簡介

表S2 本研究中935種魚類學名及對應序列號

補充材料可通過https://www.hyxbocean.cn/獲取。補充材料未進行排版和編輯，內容的準確性和科學性由作者承擔。