東帕米爾高原地區(qū)藏兔群體遺傳進(jìn)化分析

周世玉，米熱姑麗·麥麥提，劉 鵬，孟 陽，聶文悅，滕培晨，單文娟

(新疆大學(xué)生命科學(xué)與技術(shù)學(xué)院 新疆生物資源基因工程重點實驗室，烏魯木齊 830046)

藏兔(LepustibetanusWaterhouse, 1841)隸屬于兔形目(Lagomorpha)兔科(Leporidae)兔屬(Lepus)。藏兔身型纖細(xì)，身體毛色主要呈沙黃色或淺褐色。分布于阿富汗、印度、吉爾吉斯斯坦、巴基斯坦、塔吉克斯坦、中國以及蒙古[1]。在中國境內(nèi)主要分布于新疆西南部帕米爾高原地區(qū)，與新疆特有兔種塔里木兔(Lepusyarkandensis)的分布區(qū)有一定程度的重疊，易產(chǎn)生雜交和基因滲透現(xiàn)象[2]。藏兔主要以荒漠植物為食，但由于其分布區(qū)存在大量農(nóng)田，也會對農(nóng)作物產(chǎn)生危害，且藏兔皮毛及兔肉具有一定的經(jīng)濟價值，對當(dāng)?shù)匦竽赁r(nóng)業(yè)發(fā)展有一定影響。

由于藏兔的分類地位一直存在爭議，因此曾將藏兔和托氏兔(Lepustolai)作為非洲草兔(Lepuscapensis)的亞種進(jìn)行研究[3]。但是，羅澤珣[4]和程承等[5]的研究表明，非洲草兔與新疆“草兔”間的形態(tài)學(xué)指標(biāo)存在差異。Wang和Yang[6]基于線粒體全基因組水平的研究指出，中國“草兔”與南非草兔親緣關(guān)系較遠(yuǎn)。Lado等[7]基于核基因和線粒體基因?qū)Ψ侵薏萃玫倪M(jìn)化歷史進(jìn)行了研究，結(jié)果表明非洲草兔與中國兔屬物種間的親緣關(guān)系較遠(yuǎn)。單文娟等[8-11]基于線粒體基因及線粒體全序列探討了新疆“草兔”的亞種分化及新疆兔屬物種的分類，支持將新疆帕米爾高原的野兔劃分為藏兔帕米爾亞種。以上研究表明，無論是基于形態(tài)學(xué)的頭骨測量數(shù)據(jù)還是基于分子生物學(xué)的系統(tǒng)發(fā)育關(guān)系，均支持將生活在新疆帕米爾高原地區(qū)的野兔歸為藏兔。但目前對藏兔的研究主要集中于分類地位的探討，對該地區(qū)藏兔的群體遺傳學(xué)研究鮮有報道。

帕米爾高原是由印度板塊與歐亞板塊碰撞形成的一處山結(jié)，平均海拔4 500 m，主要山峰均在6 000 m 以上，屬于顯著的大陸性高山氣候[12]，其中約三分之一坐落于中國新疆西南部，即東帕米爾高原地區(qū)，該地區(qū)包括塔什庫爾干縣全部，阿克陶縣大部及烏恰縣南部3個行政區(qū)。東帕米爾高原地區(qū)環(huán)境的主要特點為：海拔高、氣壓低、含氧量低、降水少、輻射多、氣候寒冷，這種極端生境對當(dāng)?shù)貏游锏纳媸且环N挑戰(zhàn)。近年來，高原極端環(huán)境對動物遺傳多樣性與遺傳結(jié)構(gòu)的影響成為廣受關(guān)注的研究熱點，如有學(xué)者基于線粒體基因以及全基因組SNP等分子遺傳工具對藏豬、藏雞、藏狐以及藏羚羊等高原物種的遺傳多樣性與遺傳結(jié)構(gòu)進(jìn)行了探討，其結(jié)果表明高原環(huán)境是影響群體遺傳多樣性與遺傳結(jié)構(gòu)的一項重要因素[13-16]，也是進(jìn)一步研究極端環(huán)境適應(yīng)性的基礎(chǔ)。但高原極端環(huán)境會對藏兔的遺傳多樣性、遺傳結(jié)構(gòu)及遺傳分化產(chǎn)生何種影響尚未見報道。

線粒體分子標(biāo)記具有母系遺傳、變異性強、成本低等特點，張玉琮等[17]在新疆兔屬物種DNA條形碼篩選的研究中發(fā)現(xiàn)，CO1與ND4基因表現(xiàn)出更強的變異水平和更高的鑒定準(zhǔn)確率。因此，本研究選用CO1與ND4基因作為分子標(biāo)記對采集自東帕米爾高原地區(qū)的藏兔進(jìn)行研究，旨在探討高原極端環(huán)境對藏兔的遺傳多樣性、遺傳結(jié)構(gòu)及遺傳分化產(chǎn)生的影響。

1 材料與方法

1.1 樣本信息

本研究在東帕米爾高原及周邊的4個行政區(qū)共收集了基于64項頭骨指標(biāo)[18]和線粒體基因[8]初步鑒定為藏兔的37例野兔整只樣本，以下簡稱藏兔。將全部藏兔樣本按照其地理位置及海拔高度分為4個地理種群，分別為位于東帕米爾高原南部的塔什庫爾干縣紅其拉甫口岸種群(TXH)共13只樣本，當(dāng)?shù)仄骄０? 700 m，是典型高原地區(qū)；阿克陶縣及與其接壤的塔什庫爾干縣大同鄉(xiāng)種群(AKT)共8只 樣本，當(dāng)?shù)仄骄０? 800 m；位于東帕米爾高原北部烏恰縣種群(WQ)共8只樣本，當(dāng)?shù)仄骄０? 100 m；位于東帕米爾高原周邊的喀什市疏附縣、疏勒縣及澤普縣種群(KS)共8只樣本，當(dāng)?shù)仄骄０? 300 m，海拔高度相對較低(表1)。

表1 東帕米爾高原地區(qū)藏兔樣本采集信息表

1.2 基因組DNA提取

本試驗選擇藏兔的耳尖組織通過試劑盒法提取全基因組DNA，所用試劑盒為天根生化科技(北京)有限公司生產(chǎn)的血液/細(xì)胞/組織基因組DNA提取試劑盒(離心柱型)，提取后的全基因組DNA通過0.8%的瓊脂糖凝膠進(jìn)行電泳檢測。

1.3 目的基因PCR擴增

選擇電泳條帶清晰明亮的藏兔基因組總DNA樣本作為模板，使用表2中已發(fā)表的引物[17]和寶日醫(yī)生物技術(shù)(北京)有限公司生產(chǎn)的Premix Taq，用25 μL標(biāo)準(zhǔn)PCR反應(yīng)體系(其中模板DNA 1 μL，終濃度為1.0 μmol·L-1的上、下游引物各1 μL，Premix Taq酶13 μL，去離子水9 μL)擴增CO1與ND4兩個基因。反應(yīng)條件：98 ℃預(yù)變性5 min；95 ℃變性10 s，兩種引物退火條件相同均為51 ℃ 30 s，72 ℃延伸60 s，進(jìn)行30個循環(huán)；72 ℃終末延伸5 min。擴增產(chǎn)物經(jīng)1.5%的瓊脂糖凝膠電泳檢測合格后，送至生工生物工程(上海)股份有限公司進(jìn)行測序。

表2 藏兔線粒體基因PCR擴增反應(yīng)所用引物

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

表3 本研究所用線粒體基因GenBank登錄號

2 結(jié) 果

2.1 全基因組DNA提取與目的基因PCR擴增

以藏兔耳尖組織提取全基因組DNA作為模板，分別加入CO1、ND4特異性引物后進(jìn)行PCR反應(yīng)，擴增產(chǎn)物經(jīng)電泳檢測與預(yù)期產(chǎn)物大小一致，如圖1所示，選擇條帶清晰明亮的擴增產(chǎn)物進(jìn)行測序分析。

M.DNA相對分子質(zhì)量標(biāo)準(zhǔn)；a.藏兔CO1基因擴增產(chǎn)物；b.藏兔ND4基因擴增產(chǎn)物

2.2 序列特征分析

將藏兔兩個線粒體基因CO1與ND4拼接，得到總長為1 145 bp的序列。堿基平均含量分別為A：27.4%，T：33.4%，C：25.2%，G：14.0%，T+A：60.8%，C+G：39.2%。共發(fā)現(xiàn)多態(tài)性位點(variable sites)77 個，突變位點(mutations sites)78個，簡約信息位點(parsimony informative sites)59個，轉(zhuǎn)換(transition)12個、顛換(transversion)1個，轉(zhuǎn)換顛換比(Ti/Tv)12，插入和缺失位點(insertion and deletion)共9個。

2.3 遺傳多樣性參數(shù)

基于37條藏兔CO1、ND4基因串聯(lián)序列共得到17種單倍型，分別按照4個不同地理種群計算各項遺傳多樣性參數(shù)，結(jié)果見表4?？偤塑账岫鄻有?π)為(0.012±0.006)，其中WQ種群核苷酸多樣性相對最低，為(0.003±0.002)；TXH與AKT種群核苷酸多樣性相對較高，分別為(0.009±0.005)、(0.011±0.006)。KS種群核苷酸多樣性相對最高，為(0.021±0.012)?？倖伪缎投鄻有?h)為(0.935±0.021)，其中TXH種群的相對較低為(0.833±0.082)，AKT、WQ與KS 3個種群的單倍型多樣性相對較高，分別為(0.893±0.086)、(0.893±0.111)、(0.893±0.111)。中性檢驗Tajima`s D值中，TXH、KS種群為正，AKT、WQ種群及總體結(jié)果為負(fù)；Fu`s Fs值中，除WQ種群為負(fù)值以外，其它均為正值，所有結(jié)果P值均未達(dá)到顯著水平。

表4 基于藏兔CO1、ND4基因計算的遺傳多樣性參數(shù)

2.4 遺傳結(jié)構(gòu)

2.4.1 系統(tǒng)發(fā)育分析 以家兔為外群，基于藏兔17種串聯(lián)序列單倍型，通過ML法構(gòu)建的系統(tǒng)發(fā)育樹顯示(圖2)，全部單倍型聚為明顯的兩枝，其中一個大枝又以較高的置信度分出兩個小枝，共3個小枝，分別用Clade A-C表示。Clade A枝中的單倍型囊括了全部4個采樣地的樣本，且不同地理種群之間有共享單倍型；Clade B中僅包含來自TXH種群的4個單倍型；Clade C則由兩個來自KS地區(qū)的單倍型和一個來自AKT地區(qū)的單倍型組成。

2.4.2 遺傳距離 本研究中藏兔4個地理種群之間的遺傳距離見表5，結(jié)果顯示CO1基因計算出的遺傳距離相對ND4較小，但變化趨勢基本一致。KS種群與另外3個地理種群間的遺傳距離均相對較大，分別為KS-TXH=1.63%/1.96%，KS-AKT=1.35%/1.67%，KS-WQ=1.21%/1.44%；WQ-AKT群體遺傳距離最小為0.57%/0.82%。種群內(nèi)遺傳距離計算結(jié)果中，KS種群相對最大為1.93%/2.23%，WQ種群相對最小為0.21%/0.38%。

表5 基于藏兔CO1、ND4基因計算的遺傳距離

2.4.3 單倍型中介網(wǎng)絡(luò)圖 基于藏兔17種串聯(lián)序列單倍型構(gòu)建的中介網(wǎng)絡(luò)圖如圖3所示，結(jié)果與藏兔系統(tǒng)發(fā)育樹一致，全部單倍型共聚為左右兩枝，左邊又可以細(xì)分為兩小枝，左下角的小枝由來自4個地理種群的混合單倍型組成；左上角的小枝由來自TXH種群的4種單倍型組成；右邊小枝由來自KS種群的2種單倍型與來自AKT種群的1種單倍型組成，且與另外兩枝間距離較遠(yuǎn)。

圖2 基于藏兔17種串聯(lián)序列單倍型構(gòu)建的ML樹

圖3 基于藏兔17種串聯(lián)序列單倍型構(gòu)建的中介網(wǎng)絡(luò)圖

2.5 遺傳分化與基因交流

2.5.1 兔屬物種系統(tǒng)發(fā)育樹 選取圖2藏兔系統(tǒng)發(fā)育樹CladeA、B兩枝中的部分樣本以及Clade C中全部樣本的CO1基因，與GenBank中發(fā)表的其他兔屬物種的CO1基因一起構(gòu)建了ML樹(圖4)。結(jié)果顯示，Clade C中全部3個樣本與塔里木兔聚成一枝，顯示出與塔里木兔更近的親緣關(guān)系。CladeA、Clade B的藏兔聚成一枝，加上Clade C以及塔里木兔一起聚為一個的大枝，該枝與雪兔及托氏兔親緣關(guān)系稍近，與非洲草兔親緣關(guān)系最遠(yuǎn)。

圖4 基于兔屬物種CO1基因構(gòu)建的ML樹

2.5.2 AMOVA及分化指數(shù)FST分子方差分析(Analysis of Molecular Variance, AMOVA)中，將樣本分為4個基因池[TXH][AKT][WQ][KS]和2個基因池[TXH][AKT WQ KS]，[AKT][TXH WQ KS]、[WQ][TXH AKT KS]、[KS][TXH AKT WQ]分別計算分子變異與FST值，結(jié)果如表6所示。按照不同地理種群分為4個基因池時，種群間變異小于種群內(nèi)變異，F(xiàn)ST值為0.168；將每個地理種群與剩余地理種群分為2個基因池時，種群間變異均小于種群內(nèi)變異，TXH種群FST值相對最高為0.228(P<0.01)，AKT、WQ和KS種群FST值相對較小，分別為0.002(P<0.01)、0.063(P<0.05)、0.097(P<0.01)。

表6 藏兔分子方差分析

2.5.3 遺傳分化指數(shù)與基因流 藏兔不同地理種群間的遺傳分化指數(shù)FST與基因流Nm結(jié)果顯示(表7)，TXH與WQ種群間FST值最大為0.344，AKT與WQ種群間FST值最小為0.015，各種群間FST的P值均小于0.01。根據(jù)FST計算出的基因流中，AKT-WQ組基因流最大為16.495，其次為AKT-KS組為4.500，以及WQ-KS組為2.077，TXH-AKT、TXH-WQ、TXH-KS 3個組的基因流相對較小，均低于1。

表7 藏兔遺傳分化指數(shù)FST與基因流

3 討 論

3.1 藏兔的群體遺傳多樣性

單倍型多樣性(h)與核苷酸多樣性(π)是反映遺傳多樣性的重要指標(biāo)，遺傳多樣性高代表物種基因庫積累的變異豐富，應(yīng)對環(huán)境變化能力強；反之則說明該物種受到了選擇壓力的影響或經(jīng)歷了種群收縮。本研究中，藏兔樣本總核苷酸多樣性及海拔較高的TXH、AKT、WQ 3個種群的核苷酸多樣性均小于0.012，低于已經(jīng)報道的其他兔屬物種如塔里木兔(0.045±0.022)、托氏兔(0.033±0.018)、歐兔(0.027±0.003)、雪兔(0.057±0.029)等，僅海拔較低的KS種群核苷酸多樣性較高為(0.021±0.012)與其他兔屬物種接近。相應(yīng)的，海拔最高的TXH種群的單倍型多樣度也較低為(0.833±0.082)，AKT、WQ和KS 3個種群單倍型多樣度相對較高為(0.893±0.111)，但仍低于塔里木兔(0.979±0.014)和歐兔(0.988±0.004)[25-27]。中性檢驗結(jié)果的兩個指標(biāo)均未檢測到近期種群擴張的顯著證據(jù)。因此認(rèn)為在線粒體基因水平上，東帕米爾高原地區(qū)的藏兔隨著海拔的升高，核苷酸多樣性與單倍型多樣性均有降低，且當(dāng)?shù)夭赝眠z傳多樣性水平總體偏低，應(yīng)對環(huán)境變化能力弱，可能需要補充制定相關(guān)保護(hù)措施。

3.2 藏兔的群體遺傳結(jié)構(gòu)

地理隔離影響物種遺傳結(jié)構(gòu)、基因交流與物種分化，帕米爾高原山地中溝壑縱橫且分布有大量現(xiàn)代冰川，是形成地理隔離的因素之一[28]。系統(tǒng)發(fā)育分析表明(圖2)，東帕米爾高原地區(qū)的藏兔樣本共分成3枝Clade A-C，除Clade C枝中的樣本被懷疑受到了雜交的影響外，其余兩枝由海拔相異的4個地理種群的藏兔相互混雜在一起，并未呈現(xiàn)出較強的地理分布模式。中介網(wǎng)絡(luò)圖(圖3)也呈現(xiàn)出與此一致的結(jié)果。值得注意的是，海拔最高的TXH種群在系統(tǒng)發(fā)育樹中被分為兩部分，一部分與其他地理種群一起聚在Clade A中，另一部分單獨聚類在Clade B中，同一個地區(qū)的樣本出現(xiàn)了相互獨立的兩個進(jìn)化枝。王玉濤等[29]在對東帕米爾高原雪雞的系統(tǒng)發(fā)育研究中，也發(fā)現(xiàn)了類似的聚類情況。Sarabia等[30]在關(guān)于非洲金狼的研究中同樣發(fā)現(xiàn)了類似的聚類情況，并認(rèn)為這種現(xiàn)象是由于物種占據(jù)了不同的生態(tài)型而產(chǎn)生的，雖然為同一個采樣地的同一個物種，但其種群內(nèi)遺傳結(jié)構(gòu)仍存在差異，物種在不斷擴張和收縮的過程中，不同種群由于食物和棲息地競爭而形成了食性不同或棲息地碎片化的不同生態(tài)型，由于高原環(huán)境容易形成地理屏障，阻止了不同生態(tài)型群體之間的基因交流，進(jìn)而使同一地理區(qū)域分布的同物種產(chǎn)生了一定的分化。

在此基礎(chǔ)上，本研究分析了4個藏兔地理種群兩兩群體間的遺傳距離，遺傳距離反映不同群體間的差異大小和親緣關(guān)系，遺傳距離越近，兩個群體相似程度越高，親緣關(guān)系也越近；反之則表示兩群體間差異越大，親緣關(guān)系越遠(yuǎn)[31-32]。從結(jié)果可以看出(表5)，WQ-AKT組遺傳距離最小。KS種群與其他地理種群間遺傳距離最大，且是唯一一個種群內(nèi)遺傳距離大于種群間遺傳距離的群體，這可能是受到了Clade C中兩個來自KS地區(qū)雜交樣本的影響。

3.3 藏兔的群體遺傳分化

為了進(jìn)一步檢測高原環(huán)境影響下不同群體間的差異與基因交流情況，本研究進(jìn)行了AMOVA分子方差分析，并計算了固定系數(shù)FST和基因流。固定系數(shù)取值范圍從1到0，值越高，則說明群體間差異越大，分化程度也越大，一般認(rèn)為，F(xiàn)ST值大于0.25就說明群體間已經(jīng)產(chǎn)生了明顯的分化[33]。本研究分析結(jié)果顯示(表6)，4個基因池分析時，4個地理種群之間FST值為0.168，2個基因池分析時，TXH種群FST值最大為0.228，P<0.01，其他種群FST相對較小(均低于0.1，P<0.05)，說明不同種群間分子差異主要來源于TXH種群。各種群間分化指數(shù)與基因流的分析結(jié)果也顯示(表7)，TXH與AKT、WQ、KS種群間的FST較高(均大于0.4)，且基因流較小(均小于1)，表明TXH種群分化程度最大。從地理上看，由于TXH種群采樣地的海拔為所有采樣點中最高(4 700 m)，當(dāng)?shù)夭赝眯杳媾R更為嚴(yán)苛的高原環(huán)境的挑戰(zhàn)，且Naldi等[34]在其研究中指出，高原陡坡會影響歐兔的遷移和分布，紅其拉甫口岸的極端高海拔環(huán)境形成的地理屏障加大了當(dāng)?shù)赝脤傥锓N與周邊地區(qū)兔屬物種進(jìn)行基因交流的難度，使遺傳分化指數(shù)升高。

而AKT-WQ組的FST值相對最小，基因流最大(FST=0.015，Nm=16.495)，表明AKT與WQ種群間分化程度不明顯，且保持著密切的基因交流，這與表5中計算的AKT-WQ組遺傳距離結(jié)果一致。AKT-KS、WQ-KS兩組的基因流也相對較高，分別為4.500和2.007。且兩地之間多有綠洲、耕地相連。AKT種群采樣地阿克陶縣雖然海拔相對較高(平均海拔3 800 m)，但北部與烏恰縣(平均海拔2 100 m)相連，南部有葉兒羌河將其與東部喀什地區(qū)的澤普縣(平均海拔1 300 m)相連。說明AKT、WQ和KS 3個種群的藏兔雖然海拔高度有差異，但依托周圍地理條件，仍保持著一定的基因交流。

3.4 藏兔線粒體基因滲透現(xiàn)象

雖然本研究中的全部樣本基于形態(tài)學(xué)及分子生物學(xué)數(shù)據(jù)被初步劃分為藏兔，但基于藏兔線粒體基因串聯(lián)序列單倍型構(gòu)建的系統(tǒng)發(fā)育樹(圖2)及中介網(wǎng)絡(luò)圖(圖3)顯示，Clade C枝中2個KS地區(qū)的藏兔樣本與1個AKT地區(qū)的藏兔樣本均以單系分布，與其他地區(qū)藏兔樣本形成的另外兩小枝組成的一大支呈并系關(guān)系。由于全世界其他兔屬物種CO1基因數(shù)據(jù)較為全面，因此選用CO1基因構(gòu)建世界兔屬物種系統(tǒng)發(fā)育樹(圖4)，結(jié)果顯示Clade C枝中的3個樣本與數(shù)據(jù)庫中下載的塔里木兔親緣關(guān)系最近?；诓赝脙蓚€線粒體基因(CO1、ND4)分別計算的遺傳距離結(jié)果也顯示，KS種群與其他種群遺傳距離較大，且種群內(nèi)遺傳距離大于種群間遺傳距離，說明KS種群內(nèi)的樣本可能發(fā)生過雜交。從地理上看，該枝中的3個樣本所處的地理位置(喀什市與阿克陶縣)位于塔里木盆地邊緣，與塔里木兔的分布范圍重疊[8]，是兩個物種的分布接觸區(qū)。

Alves等[35]早在2003年研究兔屬物種遺傳結(jié)構(gòu)時就發(fā)現(xiàn)該屬存在由于雜交導(dǎo)致的線粒體基因滲透(introgression of mitochondria DNA)現(xiàn)象，F(xiàn)erreira等[36]在其關(guān)于歐兔與伊比利亞野兔的研究中指出，在兩種野兔分布范圍的接觸區(qū)，由于物種入侵和收縮，導(dǎo)致了反復(fù)雜交引起的線粒體基因滲透現(xiàn)象，Wu等[2]及Liu等[37]曾報道，中國兔屬不同種之間存在雜交和線粒體基因滲透現(xiàn)象，可以看出，在兔屬物種接觸區(qū)發(fā)生雜交及線粒體基因滲透是頻繁且常見的現(xiàn)象。綜上所述，認(rèn)為Clade C中的3個藏兔樣本可能受到了與塔里木兔雜交而產(chǎn)生的線粒體基因滲透的影響。

4 結(jié) 論

本研究基于東帕米爾高原地區(qū)藏兔線粒體CO1及ND4兩個基因的分析表明，TXH、AKT、WQ種群的藏兔受高原極端環(huán)境影響，遺傳多樣性水平較低，應(yīng)對環(huán)境變化能力較弱，應(yīng)加強當(dāng)?shù)氐姆N質(zhì)資源保護(hù)意識，制定相關(guān)保護(hù)政策。在遺傳結(jié)構(gòu)上不同采樣地種群沒有呈現(xiàn)出較強的系統(tǒng)地理分布模式。紅其拉甫口岸的藏兔種群由于高原地理屏障作用，可能形成了不同生態(tài)型且與周邊地區(qū)藏兔存在一定程度的分化。此外，本研究還發(fā)現(xiàn)，藏兔與塔里木兔分布范圍接觸區(qū)可能存在線粒體基因滲透現(xiàn)象。