陳蔚濤，段辛斌，高 雷，李新輝，楊計平，汪登強

1. 中國水產(chǎn)科學(xué)研究院珠江水產(chǎn)研究所/農(nóng)業(yè)農(nóng)村部珠江中下游漁業(yè)資源環(huán)境科學(xué)觀測實驗站，廣東 廣州 510380

2. 中國水產(chǎn)科學(xué)研究院長江水產(chǎn)研究所，湖北 武漢 430223

鳤 (Ochetobius elongatus) 隸屬于鯉形目、鯉科、鳤屬，曾是廣泛分布于我國南方長江、珠江等諸多水系的重要經(jīng)濟魚類[1-3]。受到涉水工程、江河環(huán)境污染、棲息地破壞、過度捕撈等因素的影響，其資源量萎縮極其嚴重，在很多河流江段已多年未被監(jiān)測到[4-5]。2016年修訂的《中國脊椎動物紅色名錄》已將鳤的保護等級上升至極度瀕危[6]；近期頒布的《長江保護法》還特別針對鳤提出了重要指示，因此亟需給予重視和保護。然而，由于鳤樣本較難獲得，目前其相關(guān)研究鮮見報道，零星研究僅集中在大階元的系統(tǒng)發(fā)育進化關(guān)系[7-8]、小樣本的遺傳變異分析[5]和局域分布鳤種群的遺傳多樣性[9]方面，這極大制約了有效保護策略的制定。

珠江和長江是我國的兩大重要水系，在魚類生物多樣性和種質(zhì)資源孕育方面起著重要作用[3,10]。南嶺山脈作為珠江和長江兩大重要水系的分水嶺，是阻隔兩個水系魚類遷移擴散的天然屏障。已有大量遺傳進化研究表明，南嶺山脈導(dǎo)致珠江和長江兩個水系共同分布的魚類發(fā)生了遺傳分化甚至形成高度分化的遺傳譜系[11-13]。在此背景下，長江鳤種群與珠江鳤種群是否形成顯著的遺傳分化，遺傳分化的程度如何，若存在遺傳分化是在什么時間節(jié)點開始的，驅(qū)動分化的環(huán)境因素是什么，這些關(guān)鍵問題均尚未找到答案。

利用遺傳進化方法揭示物種的遺傳結(jié)構(gòu)對于物種的管理保護具有重要意義[14]。本文使用珠江和長江兩個水系鳤群體的2個線粒體基因 [細胞色素b (Cytb) 和酰胺腺嘌呤二核苷酸 (NADH) 脫氧酶亞單位2 (ND2)] 和2個核基因 [肌球蛋白重鏈6(myh6) 和重組激活基因2 (RAG2)] 的序列，利用系統(tǒng)發(fā)育分析、單倍型網(wǎng)狀圖、分化時間估算等方法分析了兩個水系鳤種群的遺傳結(jié)構(gòu)，估算了分化時間，探討了導(dǎo)致其分化的可能環(huán)境因素，以期為鳤的管理和保護提供科學(xué)依據(jù)。

1 材料與方法

1.1 樣品采集

長江的7尾鳤樣本 (體長21.0～24.2 cm，體質(zhì)量74.5～112.0 g) 由中國水產(chǎn)科學(xué)研究院長江水產(chǎn)研究所于2020年在長江湖北公安江段常規(guī)漁業(yè)調(diào)查采集獲得 (圖1和表1)。剪取少量胸鰭妥善保存于無水乙醇中，用以基因組DNA的提取。處理后的鳤樣本帶回中國水產(chǎn)科學(xué)研究院長江水產(chǎn)研究所，用甲醛固定制作成標本用以科研和科普。另外，在珠江廣東肇慶江段采集鳡 (Elopichthys bambusa)1尾，作為系統(tǒng)發(fā)育分析的外類群。

表1 采樣信息和GenBank序列號Table 1 Sampling information and GenBank No.

圖1 采樣示意圖Fig. 1 Map of sampling sites

1.2 基因組DNA提取、擴增與測序

采用高鹽法提取基因組DNA，并利用質(zhì)量分數(shù)為1.2%瓊脂糖凝膠電泳檢測DNA質(zhì)量。本文分別選取了2個線粒體基因 (Cytb和ND2) 和2個核基因 (myh6和RAG2)。PCR反應(yīng)體積為30 μL，其中 Taq PCR Master Mix (上海生工) 15 μL，包括Taq DNA 聚合酶、dNTP、PCR Buffer、PCR stabilizers、Gel loading和核酸染液，模板 DNA 1 μL，正反向引物各 1 μL (濃度為 10 μmol·L?1)，加滅菌超純水至總體積30 μL。鳤和鳡4個基因擴增與測序的引物信息和PCR反應(yīng)體系參考楊計平等[9]。PCR產(chǎn)物經(jīng)質(zhì)量分數(shù)為1.2%瓊脂糖凝膠電泳檢測，對于擴增效果良好的產(chǎn)物交由測序公司雙向測序。

1.3 數(shù)據(jù)處理與分析

測序所得序列使用Lasergene軟件包 (DNASTAR Inc., Madison, WI) 來檢測測序質(zhì)量并生成一致序列。所得的一致序列使用MUSCLE[14]進行比對，然后剪掉首尾兩端噪音序列至同樣長度。在Gen-Bank核苷酸序列數(shù)據(jù)庫下載了珠江52尾鳤樣本已公布的4個基因序列，具體采集信息參考楊計平等[9](圖1和表1)。這些公布序列將與新獲得序列比對后一起用于后續(xù)分析。為了獲得更多的變異信息，分別合并了2個線粒體基因和2個核基因用以后續(xù)分析。線粒體DNA單倍型和核基因等位基因均使用DnaSP 6.12[15]進行劃分。利用DnaSP 6.12統(tǒng)計單倍型數(shù)、計算單倍型多樣性和核苷酸多樣性，并估算了兩個水系鳤群體的遺傳分化狀況。使用Network 10.2 (http://www.fluxus-engineering.com)構(gòu)建了Cytb+ND2的單倍型中介網(wǎng)狀圖[16]和myh6+RAG2的等位基因中介網(wǎng)狀圖[16]。

為了揭示珠江和長江鳤種群的系統(tǒng)發(fā)育關(guān)系，本文通過MrBayes 3.2.7[17]構(gòu)建貝葉斯樹，并以鳡作為外類群。合并了4個基因進行系統(tǒng)發(fā)育樹的構(gòu)建。貝葉斯樹使用馬爾科夫鏈蒙特卡洛方法(Markov chain Monte Carlo, MCMC) 分析運行108萬代，每1 000代取樣1次，并舍棄初始的25×105代。另外，本文使用同樣的參數(shù)設(shè)置對Cytb+ND2和myh6+RAG2也分別進行了貝葉斯樹的構(gòu)建。不同基因組合的最優(yōu)核苷酸替代模型分別利用MrModeltest 2.3[18]計算獲得 (表2)。最后，基于Kimura雙參數(shù)模型 (Kimura-2-parameter,K2P)[19]，利用MEGA 7.0[20]計算兩個水系群體之間的遺傳距離。

表2 系統(tǒng)發(fā)育樹構(gòu)建和分化時間估算的最優(yōu)模型Table 2 Optimal model of phylogenetic tree construction and differentiation time estimation

本文聯(lián)合4個基因，使用BEAST v.1.8.1[21]軟件中的貝葉斯馬爾可夫鏈蒙特卡洛方法估算兩個水系鳤群體的分化時間，采用嚴格分子鐘模型、Yule先驗?zāi)Ｐ秃虶TR+I+G最優(yōu)核苷酸替代模型。由于缺乏有效的化石記錄，本文采用固定進化速率來估算分化時間。聯(lián)合基因的進化速率由Cytb的進化速率估算獲得。由于聯(lián)合基因的K2P平均遺傳距離 (0.29%) 與Cytb的K2P平均遺傳距離 (0.43%) 的比值為0.674，因此利用Cytb的進化速率乘以比值便可獲得聯(lián)合基因的進化速率。1%～2%是鯉科魚類Cytb基因普遍接受的一個進化速率范圍[22-24]，因此聯(lián)合基因的進化速率為0.674%～1.348%。分別利用0.674%和1.348%的進化速率進行運算，每次運算共運行1×109代，每1 000代進行1次抽樣，舍棄初始5×108代。運行結(jié)束后，利用Tracer 1.6[25]評價每個參數(shù)的穩(wěn)定性狀況。最后使用TreeAnnotator v.1.8.1[21]總結(jié)出一致樹，在Figtree v.1.3.1 (http://tree.bio.ed.ac.uk/software/figtree) 顯示節(jié)點的分化時間和后驗概率。

2 結(jié)果

2.1 序列信息

本文分析數(shù)據(jù)包括新獲得的長江7尾鳤的全部4個基因片段以及楊計平等[9]52尾珠江樣本的對應(yīng)4個基因片段，共計59尾樣本。序列比對后共獲得895 bp的Cytb序列，包含了30個變異位點；935 bp的ND2序列，包含了26個變異位點；763 bp的myh6序列，包含了3個變異位點；857 bp的RAG2序列，包含了5個變異位點。在4個基因片段中均未發(fā)現(xiàn)缺失和插入。此外，成功獲得了鳡的4個基因序列，GenBank序列號分別為MW657589、MW657597、MW657605和 MW657613。

2.2 遺傳多樣性、系統(tǒng)發(fā)育分析與遺傳結(jié)構(gòu)

基于線粒體聯(lián)合基因和核基因聯(lián)合基因估算的單倍型數(shù)目、單倍型多樣性和核苷酸多樣性見表3?；诰€粒體基因估算的兩個水系鳤的單倍型多樣性為0.950，核苷酸多樣性為0.494%?；诤嘶蛴嬎愕膯伪缎投鄻有院秃塑账岫鄻有悦黠@低于線粒體聯(lián)合基因，這主要是與核基因較慢的進化速率有關(guān)?；诰€粒體基因組合發(fā)現(xiàn)兩個水系之間的遺傳分化系數(shù) ФST為0.863 (P<0.001)，基于核基因組合獲得兩個水系之間的遺傳分化系數(shù)FST為0.518(P<0.001)。

表3 鳤種群的單倍型與遺傳多樣性指數(shù)Table 3 Haplotypes and genetic diversity indexes of O. elongatus populations

線粒體基因組合共界定了29個單倍型，其中珠江群體26個 (H1—H26)，長江群體3個 (H27—H29，圖2-a)。對于核基因組合，共鑒定了13個等位基因型 (A1—A13)，其中珠江和長江群體分別界定了5和8個等位基因型 (圖2-b)?；诰€粒體基因組合的單倍型中介網(wǎng)狀圖表明，珠江和長江的鳤群體形成顯著的遺傳分支結(jié)構(gòu)，并且兩個分支之間相差22個突變位點 (圖2-a)。基于核基因組合的等位基因中介網(wǎng)狀圖雖然未形成明顯的遺傳分支結(jié)構(gòu)，但兩個水系群體之間卻各自形成了特有的等位基因型 (圖 2-b)。

圖2 單倍型/等位基因中介網(wǎng)狀圖注：圓的大小代表單倍型頻率，顏色代表所屬群體，空心圓代表未檢測到的單倍型。Fig. 2 Haplotype/allele median joining networkNote: The size of the circle represents the haplotype frequency; the color represents the group, and the hollow circle represents the undetected haplotype.

基于4個基因組合和線粒體基因組合獲得的系統(tǒng)發(fā)育樹拓撲結(jié)構(gòu)一致，均支持珠江和長江形成了顯著的支系結(jié)構(gòu) (圖3-a和圖3-c)?；诤嘶蚪M合構(gòu)建的貝葉斯樹并未形成可信的支系結(jié)構(gòu) (圖3-b)?；?個基因組合、線粒體基因組合和核基因組合估算的兩個水系群體間的遺傳距離分別為0.87%、1.52%和0.14%。

圖3 基于不同基因組合的貝葉斯系統(tǒng)發(fā)育樹Fig. 3 Bayesian phylogenetic trees based on different gene combinations

2.3 分化時間估算

Beast運算的有效取樣大小 (Effect sample size,ESS) 參數(shù)均大于200，表明運算是收斂的?；诓煌M化速率估算的分化時間顯示，珠江和長江的鳤群體形成了2個進化支系，并獲得100%的后驗概率，分化時間大概在 0.38～0.76 Ma (百萬年前) (圖4)。

圖4 基于線粒體基因組合的分化時間結(jié)果Fig. 4 Estimation of divergence time based on combination of two mitochondrial genes

3 討論

3.1 譜系分化

基于線粒體基因組合的貝葉斯系統(tǒng)發(fā)育樹、單倍型網(wǎng)狀圖和基于4個聯(lián)合基因的系統(tǒng)發(fā)育樹均發(fā)現(xiàn)，兩個水系的鳤群體存在高度分化的遺傳譜系，且形成了嚴格的珠江-長江分化格局。此外，兩個水系群體之間的遺傳距離在線粒體基因?qū)用嫔弦渤尸F(xiàn)較高水平 (1.52%)。這與范鳳娟和章群[5]的研究結(jié)果相類似，其獲得了4尾長江樣本和2尾珠江樣本，利用線粒體Cytb基因發(fā)現(xiàn)兩個水系樣本各自形成了兩個獨立遺傳譜系，且譜系間的遺傳距離為1.5%。這些研究結(jié)果表明兩個水系鳤群體經(jīng)歷了長期的地理隔離。作為珠江和長江水系的分水嶺，南嶺山脈可能是導(dǎo)致其分化的重要推動力。類似的例子也發(fā)生在珠江和長江共存的其他鯉科魚類上，比如? (Hemiculter leucisculus)[26]、馬口魚(Opsariichthys bidens)[27]。基于核基因組合獲得的系統(tǒng)發(fā)育樹和等位基因網(wǎng)狀圖發(fā)現(xiàn)，兩個水系鳤群體雖未形成獨立的分支結(jié)構(gòu)，卻形成了特有的等位基因型，說明兩個水系群體之間不存在基因交流現(xiàn)象。然而，由于本文中涉及的樣本量較少，后期需要使用更多樣本量來驗證這個觀點。

草魚 (Ctenopharyngodon idella)、青魚 (Mylophayngodon piceus)、鰱 (Hypophthalmichthys molitrix) 和鳙 (H. nobilis) 是我國重要的經(jīng)濟魚類，在珠江和長江均有廣泛分布[1-3]。基于草魚[28]、青魚[29]、鰱[30]和鳙[30]的相關(guān)遺傳分化研究發(fā)現(xiàn)，雖然珠江和長江野生群體形成了一定程度的群體分化，但均未形成獨立的遺傳譜系。本文發(fā)現(xiàn)兩個水系的鳤群體形成了兩個深度分化譜系，表明鳤的分化程度顯著高于兩個水系共同分布的這幾種重要經(jīng)濟魚類。鳤有限的分布范圍和較小的種群數(shù)量可能是造成這種差異遺傳模式的原因[1,3]。早期的地質(zhì)學(xué)和魚類區(qū)系研究認為，珠江和長江水系的部分河流相隔距離很近，甚至有少數(shù)河流直接連通[31-33]。四大家魚分布廣泛、種群數(shù)量龐大，促使兩個水系群體之間可以通過一些相鄰的河流、江段進行遷移擴散，降低了群間的遺傳分化。另外，四大家魚人工繁育種群數(shù)量龐大，不可避免的養(yǎng)殖逃逸和人工增殖放流也會導(dǎo)致兩個水系四大家魚的遺傳種質(zhì)趨于同質(zhì)。相比之下，鳤主要分布在大型河流江段與湖泊，種群規(guī)模極其有限，加上人為干預(yù)極少，因此兩個水系群體之間進行相互擴散的概率較小。

已有大量研究表明古地質(zhì)歷史事件是促進物種分化的重要驅(qū)動因子[26-27,34-35]。本文發(fā)現(xiàn)珠江和長江鳤群體的譜系分化大概發(fā)生在更新世中期。這一時期剛好與上新世后期開始的青藏高原快速隆升階段第二幕的昆侖-黃河運動時間段吻合 (1.1～0.6 Ma)。研究表明，青藏高原在上新世后期進行快速抬升，共分為3個階段[36-37]：第一階段是青藏運動，時間介于3.6～1.7 Ma；第二階段是昆侖-黃河運動，時間介于1.1～0.6 Ma；第三階段是共和運動，時間大致在0.15 Ma。青藏高原的快速隆升重塑了東亞地區(qū)的景觀特征，包括河流和山脈。早期研究也證明了青藏高原的快速抬升導(dǎo)致了一些魚類的分化甚至成種[26,38-40]。因此，兩個水系鳤的譜系分化很可能與青藏高原的快速隆升階段相關(guān)。

3.2 保護建議

針對特定物種展開遺傳結(jié)構(gòu)解析能夠有效輔助物種的管理和保護[26,41-42]。鳤的資源極其缺乏，近幾年關(guān)于鳤的零星發(fā)現(xiàn)都能引起社會的廣泛關(guān)注。由于鳤的野外資源稀少以及人工繁殖尚未成功，細化保護單元并展開就地保護是一個重要舉措。因此，通過遺傳進化分析解析鳤在不同水系的遺傳結(jié)構(gòu)對于細化其保護單元非常必要。針對瀕危或稀有物種的保育，許多學(xué)者認為演化顯著單位 (Evolutionary Significant Unit, ESU) 為適當?shù)谋Ｗo管理單元[43]。ESU是指一個種群在線粒體基因或者核基因水平已經(jīng)形成獨立的單系群 (Monophyly)[44]。本文發(fā)現(xiàn)珠江和長江鳤群體形成了獨立的線粒體基因遺傳譜系和特有的等位基因型，因此可以看作兩個不同的ESU，在保護策略制定和實施上應(yīng)該區(qū)別對待。然而，由于本研究樣本量和分子標記有限，今后需要通過擴大研究范圍和樣本量，結(jié)合更多的標記 (如核基因組數(shù)據(jù)) 來全面評估鳤的種群遺傳結(jié)構(gòu)，以期為鳤的保育策略的制定提供更系統(tǒng)的支撐。