陳蔚濤，段辛斌，高 雷，李新輝，楊計(jì)平，汪登強(qiáng)

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

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

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

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

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

1 材料與方法

1.1 樣品采集

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

表1 采樣信息和GenBank序列號(hào)Table 1 Sampling information and GenBank No.

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

1.2 基因組DNA提取、擴(kuò)增與測(cè)序

采用高鹽法提取基因組DNA，并利用質(zhì)量分?jǐn)?shù)為1.2%瓊脂糖凝膠電泳檢測(cè)DNA質(zhì)量。本文分別選取了2個(gè)線粒體基因 (Cytb和ND2) 和2個(gè)核基因 (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個(gè)基因擴(kuò)增與測(cè)序的引物信息和PCR反應(yīng)體系參考楊計(jì)平等[9]。PCR產(chǎn)物經(jīng)質(zhì)量分?jǐn)?shù)為1.2%瓊脂糖凝膠電泳檢測(cè)，對(duì)于擴(kuò)增效果良好的產(chǎn)物交由測(cè)序公司雙向測(cè)序。

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

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

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

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

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

2 結(jié)果

2.1 序列信息

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

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

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

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

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

圖2 單倍型/等位基因中介網(wǎng)狀圖注：圓的大小代表單倍型頻率，顏色代表所屬群體，空心圓代表未檢測(cè)到的單倍型。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個(gè)基因組合和線粒體基因組合獲得的系統(tǒng)發(fā)育樹(shù)拓?fù)浣Y(jié)構(gòu)一致，均支持珠江和長(zhǎng)江形成了顯著的支系結(jié)構(gòu) (圖3-a和圖3-c)。基于核基因組合構(gòu)建的貝葉斯樹(shù)并未形成可信的支系結(jié)構(gòu) (圖3-b)。基于4個(gè)基因組合、線粒體基因組合和核基因組合估算的兩個(gè)水系群體間的遺傳距離分別為0.87%、1.52%和0.14%。

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

2.3 分化時(shí)間估算

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

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

3 討論

3.1 譜系分化

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

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

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

3.2 保護(hù)建議

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