趙文浩 易少奎 蘇君曉 周 瓊 沈建忠 李大鵬 周小云

(1. 華中農業(yè)大學水產學院, 農業(yè)農村部淡水生物繁育重點實驗室, 武漢 430070; 2. 湖州師范學院生命科學學院, 湖州 313000)

葉爾羌高原鰍(Triplophysa yarkandensis), 俗稱“狗頭魚”, 是新疆南部塔里木河水系的特有魚類,隸屬于鯉形目鰍科(Cobitidae), 條鰍亞科(Noemacheilinae), 是高原鰍屬(Triplophysa)魚類中個體最大、生長速度最快的經(jīng)濟魚種, 迄今發(fā)現(xiàn)的最大個體長30 cm、重305 g[1]。近幾十年來, 塔里木河干流水流量日益減少, 鹽堿化程度加劇, 加之大量濫捕、外來入侵魚類捕食等, 導致葉爾羌高原鰍種群數(shù)量急劇減少, 即將成為繼扁吻魚和塔里木裂腹魚之后的又一瀕危物種[2]。遺傳多樣性是物種繁衍、抵抗疾病和適應環(huán)境變化的物質基礎, 因此, 開展葉爾羌高原鰍野生群體遺傳多樣性研究, 從分子水平分析其種群遺傳結構特征和種群歷史動態(tài), 進而制定合理有效的保護方案, 對葉爾羌高原鰍的種質資源保護和科學開發(fā)利用有重要的意義。

線粒體DNA(mtDNA)是一種核外遺傳物質, 因具有基因組小、結構簡單、進化速度快、幾乎不發(fā)生重組以及不同的區(qū)域進化速率存在差異等特點, 成為物種鑒定和群體遺傳學研究的理想分子標記[3]。在mtDNA中, 細胞色素b(Cytochromeb,Cytb)基因的進化速度適中, 許多學者選擇其作為分子標記, 用于分析種群母系起源、系統(tǒng)進化、遺傳結構和遺傳多樣性等[4—6]; 控制區(qū)(D-loop)是線粒體上的一段非編碼區(qū), 由于受選擇壓力較小, 其進化速率是線粒體其他區(qū)域的5—10倍, 研究發(fā)現(xiàn), 用D-loop序列進行群體遺傳多樣性和系統(tǒng)進化研究,獲得的信息更全面、結果更準確[7,8]。近年來, 通過分析Cytb和D-loop序列的堿基信息研究水產動物遺傳多樣性已被廣泛應用, 如Guo等[9]用Cytb和D-loop序列對雅魯藏布江6個地理群體的拉薩裸裂尻魚(Schizopygopsis younghusbandi)進行遺傳多樣性分析, 發(fā)現(xiàn)各群體的遺傳多樣性水平均較低, 并且遺傳變異主要來自群體內部; 杜景豪等[8]用Cytb和D-loop序列對5個日本囊對蝦(Penaeus japonicus)群體分析結果顯示, 群體間遺傳分化較明顯, 近期未曾經(jīng)歷瓶頸效應和種群擴張。上述研究都為種質資源的保護與開發(fā)利用提供了重要的數(shù)據(jù)基礎。

目前, 關于葉爾羌高原鰍的研究主要集中在基礎生物學[11,12]、繁殖習性[1,10]、肌肉營養(yǎng)成分[13]和耐鹽堿性[14,16]等方面, 而在遺傳多樣性方面, 僅見王錦秀等[15]用SSR標記對塔里木河5個群體的研究報道。為了更多地了解塔里木河葉爾羌高原鰍的遺傳多樣性現(xiàn)狀, 本研究采用線粒體Cytb基因和Dloop序列對塔里木河8個地理群體的遺傳多樣性進行分析, 研究其種群結構和歷史動態(tài)信息等, 旨在為葉爾羌高原鰍的資源保護與科學開發(fā)利用提供理論依據(jù)。

1 材料與方法

1.1 樣品采集

本研究所用葉爾羌高原鰍為2019年5至10月從塔里木河支流渭干河和車爾臣河采集, 其中拜城大宛其(DWQ)22尾、托克遜(TKX)20尾、且末第二分水樞紐(JDT)22尾、阿克提坎村(AKT)11尾、蘇外斯埂管護站(SWS)23尾、且末五葦場(WWC)28尾、瓦石峽鄉(xiāng)(WSX)21尾和臺特瑪湖(TTM)27尾, 共174尾(圖 1)。根據(jù)形態(tài)學特征確認物種后,取鰭條于100%乙醇中暫存, 運回實驗室后置于–80℃冰箱中保存?zhèn)溆谩?/p>

圖1 葉爾羌高原鰍采樣點Fig. 1 The sampling locations of T. yarkandensis

1.2 DNA提取、序列擴增及測序

取鰭條約30 mg, 用醋酸銨/異戊醇法提取基因組DNA[17], 分別用1%的瓊脂糖凝膠電泳和Nano-Drop 2000紫外分光光度計(Thermo Scientific,USA)檢測DNA的質量和濃度。將合格的DNA稀釋至50 ng/μL用于后續(xù)實驗。

分別采用Xiao等[18]和Liu等[19]的引物擴增Cytb和D-loop序列, 具體為: CytbF: 5′-GACTTG AAAAACCACCGTTG-3′, R: 5′-CTCCGATCTCCG GATTACAAGAC-3′; D-loop F: 5′-ACCACTGGC TCCCAAAGC-3′, R: 5′-ATCTTAGCATCTTCAGTG-3′。PCR反應體系30 μL, 包括: DNA模板3 μL、上下游引物各0.75 μL(10 μmol/L), 2×PCR Master Mix 15 μL, ddH2O 10.5 μL。PCR程序為: 94℃預變性5min, 35個循環(huán)[94℃變性30s, 退火30s(Cytb為56℃,D-loop為54℃), 72℃延伸50s], 72℃延伸10min。PCR產物經(jīng)1%瓊脂糖凝膠電泳檢測為單一明亮的條帶后送公司測序。

1.3 數(shù)據(jù)分析

測序結果用Bio Edit 7.0.5軟件[20]進行編輯并人工校對, 然后用MEGA 7.0[21]將其與NCBI中葉爾羌高原鰍mtDNA(KT224367)中的Cytb和D-loop序列同源比對, 確認測序結果的準確性, 并統(tǒng)計序列的堿基組成、變異和群體間遺傳距離; 用鄰接法(Neighbor-joining, NJ)基于Kimura2-parameter(K2-p)模型構建樣本系統(tǒng)發(fā)育樹, 用Bootstrap(1000次)檢驗聚類樹各分支的置信度。用DnaSP 5.1[24]統(tǒng)計多態(tài)位點數(shù)目(Number of polymorphic sites, S)、核苷酸多樣性指數(shù)(Nucleotide diversity,Pi)、單倍型數(shù)目(Number of haplotype,H)、單倍型多樣性指數(shù)(Haplotype diversity,Hd)等。用PopART[22]構建單倍型網(wǎng)絡進化圖。用Barrier v.2.2[23]分析群體間的遺傳障礙。用Arlequin 3.5[25]進行AMOVA分析(Analysis of molecular variance)計算群體內和群體間的遺傳變異, 并進行Tajima’sD[26]和Fu’sFs[27]中性檢驗以及核苷酸不配對分布分析(Mismatch distribution)研究群體是否發(fā)生擴張。

2 結果

2.1 序列組成及變異

測序結果經(jīng)校正和比對后, 分別選取1079 bp的Cytb和887 bp的D-loop序列用于后續(xù)分析。結果顯示, 兩序列的A、T、C、G平均含量分別為24.0%、33.1%、25.7%、17.2%和32.5%、31.1%、21.3%、15.1%、63.6%, A+T含量(Cytb, 57.1%; D-loop,63.6%)均顯著高于G+C含量(Cytb, 42.9%; D-loop,36.4%)。Cytb的多態(tài)性位點數(shù)為50個, 包括40個簡約信息位點和10個單一突變位點; D-loop的多態(tài)位點數(shù)為57個, 包括47個簡約信息位點和10個單一突變位點; 兩序列的堿基轉換與顛換比分別為9.8﹕1和2.9﹕1。堿基突變位點如圖 2所示。

2.2 群體遺傳多樣性

基于Cytb序列共檢測到41個單倍型(圖 2), 其中, Hap-41為8個群體共有單倍型。單倍型多樣性指數(shù)(Hd)為0.9247±0.0124, 核苷酸多樣性指數(shù)(Pi)為0.0065±0.0034; SWS群體的Hd和Pi指數(shù)最高,JDT群體最低?；贒-loop序列共檢測到51個單倍型, 沒有發(fā)現(xiàn)所有群體的共享單倍型,Hd為0.9776±0.0032,Pi為0.0115±0.0058; AKT群體的Hd指數(shù)最高, 而TKX群體的Pi指數(shù)最高。此外, 基于D-loop分析得到的單倍型數(shù)量、Hd和Pi均高于基于Cytb序列的分析結果(表 1)。

表1 葉爾羌高原鰍各群體的遺傳多樣性參數(shù)Tab. 1 Genetic diversity parameters in each population of T. yarkandensis

圖2 葉爾羌高原鰍Cytb (a)與D-loop (b)序列的變異位點圖Fig. 2 The genetic variation loci of T. yarkandensis in Cytb (a) and D-loop (b) sequences

2.3 群體遺傳結構與系統(tǒng)發(fā)育

遺傳分化指數(shù)(Fst)表示群體間的遺傳分化程度, 通常分為3個等級: 弱分化(Fst<0.05)、中度分化(0.05≤Fst≤0.25)和高度分化(Fst>0.25)[28]。在本研究中, 除TTM和TKX群體與其他群體間的Fst在0.05—0.25, 屬于中度分化外, 其他群體間的Fst均小于0.05, 屬于弱分化水平(圖 3)。同時, 各群體間的遺傳距離小, 其范圍分別為0.0055—0.0071(Cytb)和0.0092—0.0128(D-loop; 表 2), 表明群體間具有較高的遺傳同質性。

表2 基于Cytb(對角線下)和D-loop(對角線上)序列的葉爾羌高原鰍群體間遺傳距離Tab. 2 Genetic distance between populations of T. yarkandensis based on Cytb (under the diagonal) and D-loop (on the diagonal)sequences analysis

圖3 基于Cytb(a)和D-loop(b)序列的群體間遺傳分化指數(shù)(Fst)Fig. 3 Pairwise Fst between pairs of populations based on Cytb (a) and D-loop (b) sequences

AMOVA分析結果顯示, 基于Cytb和D-loop序列獲得的群體內變異貢獻率分別為96.57%和95.25%, 而群體間的變異貢獻率僅為3.43%和4.75%(表 3)。為了解支流間的遺傳變異, 本研究將8個群體按支流分為兩個亞群(渭干河亞群和車爾臣河亞群), AMOVA分析發(fā)現(xiàn), 遺傳變異依然主要來自群體內部(Cytb, 93.55%; D-loop, 92.65%), 而支流間(Cytb, 5.40%; D-loop, 4.71%)和支流內群體間(Cytb,1.05%; D-loop, 2.64%)遺傳變異較小。但是, 遺傳分化指數(shù)分析結果顯示, 兩支流之間出現(xiàn)了一定程度的遺傳分化(Cytb,Fct=0.0540,P<0.05; D-loop,Fct=0.0471,P<0.05; 表 3)。

表3 基于Cytb和D-loop序列的葉爾羌高原鰍群體間AMOVA分析結果Tab. 3 AMOVA analysis results of inter-population based on Cytb and D-loop sequences

基于各樣本的Cytb和D-loop序列構建的NJ系統(tǒng)發(fā)育樹顯示, 8個群體的樣本交錯分布, 沒能形成明顯的地理聚類(圖 4)。為了驗證結果的準確性,本研究進一步構建了單倍型網(wǎng)絡進化圖, 結果顯示,不同群體間共享單倍型較多, 也沒能形成明顯的地理聚群和譜系結構(圖 5)。這與上述分析中各群體間遺傳分化指數(shù)Fst值小、遺傳距離值低的結果相一致。用BARRIER軟件分析兩支流的群體間是否存在遺傳障礙, 基于Cytb和D-loop序列的分析結果一致, 均顯示渭干河與車爾臣河群體之間存在潛在的遺傳障礙(圖 6)。

圖4 基于Cytb (a)和D-loop (b)序列構建的葉爾羌高原鰍各樣本NJ系統(tǒng)發(fā)育樹Fig. 4 NJ phylogenetic trees of T. yarkandensis samples based on Cytb (a) and D-loop (b) sequences

圖5 基于Cytb (a)和D-loop (b)序列構建的葉爾羌高原鰍單倍型網(wǎng)絡圖Fig. 5 Median-Joining network of Cytb (a) and D-loop (b) haplotypes from eight populations of T. yarkandensis

圖6 基于Cytb (a)和D-loop (b)遺傳距離的遺傳障礙圖Fig. 6 Map of genetic barriers based on Cytb (a) and D-loop (b)genetic distances圓點表失8個葉爾羌高原鰍所在地理位置, 紅色箭頭a, b代表遺傳障礙Circle points represent the locations of the eight T. yarkandensis labeled by abbreviations of their populations, the red arrows a and b represent genetic barriers

2.4 群體歷史動態(tài)信息

Tajima’sD和Fu’sFs中性檢驗結果如表 4所示,除SWS和WWC兩個群體的Fs值為負值且統(tǒng)計學結果顯著(P<0.05)外, 其他6個群體的D值和Fs值的統(tǒng)計學結果都不顯著(P>0.05)。隨后, 我們將所有個體作為一個整體進行Tajima’sD和Fu’sFs檢驗, 得到的D值和Fs值均為負值, 但D值的統(tǒng)計學結果不顯著(P>0.05), 結合大部分群體D值與Fs值都為正值且統(tǒng)計學不顯著的結果, 說明沒有明顯偏離中性突變。此外, 雖然錯配分析擬合度檢驗結果支持擴張假設(P>0.05; 表 4), 但是基于Cytb與D-loop所得核苷酸錯配分布曲線都為多峰形(圖 7a和7b)。綜合數(shù)據(jù)結果表明, 整個葉爾羌高原鰍種群相對穩(wěn)定,近期沒有經(jīng)歷過瓶頸效應和明顯擴張。

圖7 基于Cytb (a)和D-loop (b)序列的葉爾羌高原鰍8個群體核苷酸錯配分布圖Fig. 7 Mismatch distribution of eight T. yarkandensis populations based on Cytb (a) and D-loop (b) sequences

表4 葉爾羌高原鰍8個群體的中性檢驗及擬合度檢驗Tab. 4 Neutrality test and goodness-of-test for eight T. yarkandensis populations

3 討論

3.1 序列組成及突變

本研究基于線粒體Cytb和D-loop序列對塔里木河支流渭干河和車爾臣河8個葉爾羌高原鰍群體的遺傳多樣性進行了分析。兩序列的堿基組成分析發(fā)現(xiàn), A+T含量均顯著高于G+C含量, 呈現(xiàn)明顯的偏隨機分布, 這與多數(shù)脊椎動物線粒體基因的編碼模式一致[29]。Meyer[30]提出, 在線粒體DNA進化過程中, 發(fā)生轉換的頻率遠高于發(fā)生顛換的頻率。葉爾羌高原鰍Cytb基因和D-loop序列的堿基轉換與顛換比值分別為9.8﹕1和2.9﹕1, 表明兩序列的堿基突變類型均以轉換為主, 這與Meyer的結論一致[30]。

3.2 種群遺傳多樣性

遺傳多樣性是物種進化和適應環(huán)境變化的物質基礎, 也是評價種群資源狀況的重要依據(jù), 一個物種的遺傳多樣性水平越高, 其對環(huán)境的適應能力就越強, 就能獲得更多的生存空間使得種群得到更好的發(fā)展, 相反, 遺傳多樣性缺乏則會導致物種資源衰退甚至滅絕[31]。線粒體DNA的核苷酸多樣性(Pi)和單倍型多樣性(Hd)是評價物種群體遺傳變異程度的重要指標, 也是衡量一個群體遺傳多樣性和遺傳分化的重要參數(shù)[32]。Grant和Bowen[33]認為, 如果群體的核苷酸多樣性指數(shù)Pi>0.05且單倍型多樣性指數(shù)Hd>0.5, 則表明該群體的遺傳變異較豐富,遺傳多樣性也較高。在本研究中, 8個葉爾羌高原鰍群體的Pi指數(shù)均大于0.005、Hd指數(shù)均大于0.5,因此均屬于高核苷酸多樣性、高單倍型多樣性群體, 這與王錦秀等[15]用SSR標記對塔里木河5個群體的葉爾羌高原鰍的研究結果類似。塔里木河葉爾羌高原鰍較高的遺傳多樣性水平, 可能是由于其對環(huán)境的適應能力較強, 在水體發(fā)生匯流時能很好地進行遷移, 使得群體間能夠進行基因交流, 進而促進了群體的遺傳變異[14—16]。此外, 基于D-loop序列獲得的各群體單倍型數(shù)、單倍型多態(tài)性指數(shù)和核苷酸多態(tài)性指數(shù)均大于基于Cytb序列分析的結果, 這可能與D-loop區(qū)不編碼蛋白質, 所受的選擇壓力小, 進化速率會更快等因素有關[34,35], 這也與杜景豪等[8]在日本囊對蝦和Guo等[9]在拉薩裸裂尻魚中的研究結果一致。

3.3 群體遺傳結構

遺傳距離大小是衡量物種間或同一物種不同群體間親緣關系遠近的常用指標, 種群內的遺傳距離通常小于10%, 大于6%的群體間會有明顯的亞種或種的分化[36]。在本研究中, 8個群體間的遺傳距離值均較低(Cytb, 0.0055—0.0071; D-loop, 0.0092—0.0128), 表明相互間的遺傳距離較近。根據(jù)Wright[28]的遺傳分化理論, 遺傳分化指數(shù)Fst<0.05代表遺傳分化水平較低, 在本研究中, 大部分群體間的遺傳分化指數(shù)Fst都小于0.05, 因此為低分化水平。王錦秀等[15]用SSR分析塔里木河5個葉爾羌高原鰍群體的遺傳結構發(fā)現(xiàn), 其遺傳距離在0.606—1.901, 明顯高于本研究結果, 并且多數(shù)群體間的遺傳分化指數(shù)在0.05—0.15, 達到中等程度的遺傳分化。產生上述結果差異的原因, 一方面可能是由于線粒體DNA標記(Cytb和D-loop)與核DNA標記(SSR)的進化速率不同[37], 因此, 用兩種類型的遺傳標記對群體遺傳結構解析時可能會出現(xiàn)一定的差異; 另一方面,由于種群的遺傳結構受地理隔離、生存環(huán)境、種群瓶頸和基因流等多種因素的影響, 在王錦秀等[15]的研究中, 其5個采樣點分別位于不同的支流(阿爾干、臺特瑪湖、車爾臣河、阿克蘇河和臺南河),相互之間的地理位置距離較遠, 群體間的基因交流相對較少。此外, AMOVA分析發(fā)現(xiàn), 遺傳變異主要來源于群體內部, 而群體間遺傳變異很小, 僅3.43%(Cytb)和4.75%(D-loop); 基于樣本構建的系統(tǒng)發(fā)育樹和單倍型構建的網(wǎng)絡進化圖均顯示, 各群體的樣本沒能形成明顯的地理聚群和譜系結構。綜合上述結果表明, 本研究的8個葉爾羌高原鰍群體間存在基因交流。雖然渭干河和車爾臣河屬于塔里木河中游和下游的不同支流, 但自2000年起, 塔里木河下游啟動了生態(tài)輸水工程[38], 這可能促進了渭干河和車爾臣河間的葉爾羌高原鰍的基因交流。王錦秀等[15]在其研究中也發(fā)現(xiàn), 不同支流的葉爾羌高原鰍群體間存在基因交流, 并指出, 這種交流增加了群體的等位基因數(shù), 是其群體遺傳多樣性高的原因之一。

3.4 群體歷史動態(tài)

Tajima’sD與Fu’sFs檢驗在檢測群體是否發(fā)生擴張方面比較敏感, 當明顯偏離中性突變時,D值與Fs值為顯著負值, 表明發(fā)生過群體擴張事件[39]; 核苷酸錯配分布曲線呈單峰時表明群體發(fā)生過擴張,呈雙峰或多峰形時表明群體未經(jīng)歷擴張, 而是保持相對穩(wěn)定的狀態(tài)[40]。本研究基于Cytb和D-Loop序列的分析結果顯示, 大部分群體的D值與Fs值為正值且沒有顯著性差異(P>0.05), 因此不能拒絕中性突變, 即沒有發(fā)生過群體擴張事件; 同時, 核苷酸錯配分布曲線為雙峰形。這些結果表明, 本研究中的8個葉爾羌高原鰍群體沒有經(jīng)歷擴張, 而是保持著一個相對穩(wěn)定的狀態(tài)。

綜上所述, 塔里木河葉爾羌高原鰍各群體的遺傳多樣性水平較高, 但遺傳變異主要來源于群體內部, 群體間的遺傳分化程度較小, 各樣本沒能形成明顯的遺傳分支與地理群體聚類。雖然兩支流群體間存在潛在的遺傳障礙, 并發(fā)生了一定的遺傳分化, 但其遺傳距離只在0.05附近(Cytb,Fct=0.0540;D-loop,Fct=0.0471)。因此, 根據(jù)進化上具有關鍵意義的單元ESU(Evolutionarily Significant Unit)設置原則[41], 本研究建議將塔里木河葉爾羌高原鰍種群作為一個保護管理單元進行保護。