基于簡(jiǎn)化基因組測(cè)序技術(shù)的喜馬拉雅旱獺遺傳多樣性分析

李 優(yōu)，南新?tīng)I(yíng)，趙坤鈺，李耀東

（青海大學(xué)，西寧，810016）

喜馬拉雅旱獺（Marmota himalayana）屬?lài)X目（Rodentia），松鼠科（Sciuridae），旱獺屬（Marmota），是高原環(huán)境中的世居動(dòng)物，也是高原環(huán)境中的特有物種，主要分布于我國(guó)青海、西藏、云南和青藏高原邊緣山地［1］。由于喜馬拉雅旱獺長(zhǎng)期生活在高原地區(qū)，對(duì)低氧環(huán)境有很好的適應(yīng)性，在生理生化上獲得了適應(yīng)高原地區(qū)低溫低氧的穩(wěn)定遺傳學(xué)特性，是研究高原低氧適應(yīng)的天然理想動(dòng)物模型［2］。目前，對(duì)于喜馬拉雅旱獺的研究有起源與分化［3］、轉(zhuǎn)錄組學(xué)［4］、實(shí)驗(yàn)動(dòng)物學(xué)［5］和線粒體基因組多態(tài)性及群體結(jié)構(gòu)［6］等，對(duì)于喜馬拉雅旱獺的研究較為成熟。近年來(lái)隨著分子生物學(xué)技術(shù)高速發(fā)展，基于第二代測(cè)序技術(shù)的簡(jiǎn)化基因組測(cè)序隨之產(chǎn)生并得到廣泛應(yīng)用，其中雙酶切RAD（dd-RAD）測(cè)序技術(shù)是在傳統(tǒng)二代測(cè)序RAD-Seq 的基礎(chǔ)上研發(fā)的雙酶切RAD-Seq［7］，被廣泛應(yīng)用于群體遺傳結(jié)構(gòu)和多樣性分析［8］以及群體遺傳進(jìn)化［9］等領(lǐng)域，能快速鑒別高密度單核苷酸多態(tài)性（SNP）位點(diǎn)［10］。SNP 位點(diǎn)作為基因組中最容易發(fā)生變異的類(lèi)型，可作為群體遺傳學(xué)和定量遺傳學(xué)的理想工具［11］，但迄今為止，尚未有利用dd-RAD測(cè)序技術(shù)對(duì)喜馬拉雅旱獺開(kāi)展遺傳多樣性研究的相關(guān)報(bào)道。本研究基于dd-RAD 測(cè)序技術(shù)對(duì)喜馬拉雅旱獺進(jìn)行群體遺傳分析，通過(guò)SNP 位點(diǎn)刻畫(huà)其遺傳特性，構(gòu)建系統(tǒng)進(jìn)化樹(shù)，分析種群結(jié)構(gòu)特征，從基因水平分析不同地區(qū)個(gè)體之間的遺傳進(jìn)化關(guān)系，旨在了解青藏高原本地野生動(dòng)物資源的遺傳潛力，為喜馬拉雅旱獺比較基因組學(xué)和旱獺屬動(dòng)物的資源合理利用提供參考。

1 研究方法

1.1 樣本采集

共使用44 例喜馬拉雅旱獺樣本，按地點(diǎn)分為3個(gè)地理種群，其中14例采自青海省玉樹(shù)州曲麻萊縣麻多鄉(xiāng)（簡(jiǎn)稱(chēng)玉樹(shù)，編號(hào)為YS1～YS14；海拔4 452 m；35°2'52″ N，96°37'39″ E），18例采自黃南州同仁縣瓜什則鄉(xiāng)（簡(jiǎn)稱(chēng)黃南，編號(hào)HN1～HN18；海拔3 273 m；35°49'35″ N，102°28'28″ E），12例采自青海省海東市樂(lè)都區(qū)引勝鄉(xiāng)（簡(jiǎn)稱(chēng)樂(lè)都，編號(hào)LD1～LD12；海拔 2 279 m；36°60'19″ N，102°40'37″ E）。取樣時(shí)喜馬拉雅旱獺均已用乙醚深度麻醉處死，符合動(dòng)物倫理要求，將得到的喜馬拉雅旱獺肝臟組織用液氮快速冷凍，放入實(shí)驗(yàn)室冰箱-80 ℃貯存?zhèn)溆谩?/p>

1.2 基因組DNA提取與質(zhì)量控制

用EasyPure 基因組提取試劑盒（北京全式金生物技術(shù)股份有限公司）提取喜馬拉雅旱獺基因組DNA，用1.5%瓊脂糖凝膠電泳和紫外分光光度計(jì)NanoDrop 2000C 儀器檢測(cè)提取的基因組DNA 質(zhì)量，基因組DNA 的A260/280值為1.7～2.0符合質(zhì)量要求，樣本基因組DNA于-20 ℃保存?zhèn)溆谩?/p>

1.3 dd-RAD文庫(kù)構(gòu)建

通過(guò)參考基因組序列模擬限制性?xún)?nèi)切酶的酶切位點(diǎn)的數(shù)量與分布，選擇EcoRⅠ限制性?xún)?nèi)切酶和NlaⅢ限制性?xún)?nèi)切酶進(jìn)行雙酶切試驗(yàn)并建庫(kù)。質(zhì)檢合格的基因組DNA 樣本采用dd-RAD 構(gòu)建文庫(kù)，具體方法：取喜馬拉雅旱獺基因組DNA 1 μg 作為起始量建庫(kù)，加入10×NE Buffer 5 μL，加入1 μL限制性?xún)?nèi)切酶EcoRⅠ（NEB）和1 μL 限制性?xún)?nèi)切酶NlaⅢ，最后加入無(wú)酶水使反應(yīng)體系達(dá)到50 μL。反應(yīng)體系在37 ℃下溫育15 min，在65 ℃溫育20 min，熱失活后進(jìn)行末端修復(fù)，加入接頭。使用2%的瓊脂糖凝膠電泳對(duì)連接產(chǎn)物進(jìn)行片段選擇，選擇400～600 bp 的條 帶，切膠回收酶切產(chǎn)物。文庫(kù)質(zhì)量控制合格后，用 Illumina TruSeq 測(cè)序平臺(tái)對(duì)構(gòu)建的DNA 文庫(kù)進(jìn)行簡(jiǎn)化基因組測(cè)序。

1.4 全基因組SNP篩查、分型與質(zhì)控

用GATK 軟件工具包檢測(cè)SNP［12］。以高質(zhì)量測(cè)序數(shù)據(jù)的結(jié)果為依據(jù)，為保證檢測(cè)得到的SNP 的準(zhǔn)確性，用Samtools 去重復(fù)，GATK 局部重比對(duì)、校正堿基質(zhì)量值，再使用GATK 進(jìn)行核苷酸多態(tài)性的檢測(cè)、過(guò)濾，得到最終的SNP位點(diǎn)集。質(zhì)控條件：丟棄低質(zhì)量（Q＜20）和缺失預(yù)期限制的位點(diǎn)，同時(shí)丟棄具有較小等位基因頻率（MAF＜0.05）和最大雜合度（He＞0.7）的SNP。此外，對(duì)每個(gè)個(gè)體獲得的基因型進(jìn)行支持每個(gè)等位基因的reads 數(shù)量的詢(xún)問(wèn)，丟棄少于20 個(gè)reads 支持的基因型和其中一個(gè)等位基因的覆蓋率比另一個(gè)等位基因高3 倍以上的基因型。最后，丟棄所有樣本中低于75%個(gè)體可以分型的位點(diǎn)。

1.5 數(shù)據(jù)統(tǒng)計(jì)與分析

利用BWA 軟件［13］、GATK 軟件、ADMIXTURE 軟件［14］、GCTA 軟件［15］和Fast Tree 軟件［16］分別對(duì)SNP注釋、遺傳分化系數(shù)統(tǒng)計(jì)，并對(duì)SNP位點(diǎn)觀測(cè)雜合度（Ho）、期望雜合度（He）、近交系數(shù)（Fis）、多樣性指數(shù)（Shannon）、多態(tài)信息含量（PIC）、有效等位基因數(shù)（Ne）和最小等位基因頻率（MAF）7 個(gè)遺傳多樣性參數(shù)進(jìn)行分析?；诘玫降腟NP，使用ADMIXTURE軟件計(jì)算群體結(jié)構(gòu)，以分群數(shù)（K值）為2～10 進(jìn)行聚類(lèi)，做出群體結(jié)構(gòu)聚類(lèi)分析。使用GCTA 軟件，通過(guò)主成分分析（PCA）獲得樣本遺傳距離遠(yuǎn)近關(guān)系，輔助進(jìn)化分析。

2 結(jié)果

2.1 測(cè)序數(shù)據(jù)統(tǒng)計(jì)

基于dd-RAD 測(cè)序技術(shù)構(gòu)建44 例喜馬拉雅旱獺的測(cè)序文庫(kù)，共獲得106 699 757條原始序列，平均每個(gè)樣本的測(cè)序序列為2 424 994 條，獲得總堿基數(shù)為30 729 530 016 個(gè)，平均值為698 398 409 個(gè)，喜馬拉雅旱獺樣本測(cè)序獲得的GC比例為39.53%～41.09%，平均GC 比例為40.39%，所有樣本的Q20是100%，所有的Q30在82.51%以上，表明測(cè)序結(jié)果較準(zhǔn)確，堿基錯(cuò)誤率低，數(shù)據(jù)質(zhì)量合格，可用于進(jìn)一步分析（表1）。

表1 喜馬拉雅旱獺dd-RAD測(cè)序數(shù)據(jù)Tab.1 dd-RAD sequencing data of Marmota himalayana

2.2 SNP檢測(cè)與位點(diǎn)開(kāi)發(fā)

44 例喜馬拉雅旱獺樣本共檢測(cè)出19 279 845 個(gè)SNP 位點(diǎn)。轉(zhuǎn)換SNP 為147 893～332 375，顛換SNP為81 923～186 554，發(fā)生轉(zhuǎn)換與顛換的比值為1.759～1.809，轉(zhuǎn)換/顛換均大于1.5，表明喜馬拉 雅旱獺與大多數(shù)脊椎動(dòng)物相同，存在轉(zhuǎn)換型顛倒 偏差現(xiàn)象［17］。SNP 遺傳多樣性分析結(jié)果顯示，雜合 突變數(shù)為21 598～59 268，純合突變數(shù)為208 218～463 280（表2）。

表2 喜馬拉雅旱獺SNP統(tǒng)計(jì)結(jié)果Tab.2 Statistical results of SNP of Marmota himalayana

2.3 遺傳進(jìn)化分析

2.3.1 遺傳統(tǒng)計(jì)量分析

根據(jù)得到的SNP結(jié)果，分析遺傳統(tǒng)計(jì)量，結(jié)果顯示：樂(lè)都（0.236）和玉樹(shù)（0.269）喜馬拉雅旱獺觀測(cè)雜合度均大于期望雜合度（0.232、0.235），黃南則是期望雜合度（0.254）大于觀測(cè)雜合度（0.240）。近交系數(shù)為0.014～0.145。3 個(gè)地區(qū)的喜馬拉雅旱獺的多態(tài)信息含量均小于0.250，為低度多態(tài)性。黃南有效等位基因數(shù)為1.435，最接近等位基因數(shù)2，3 個(gè)地區(qū)喜馬拉雅旱獺有效等位基因數(shù)在種群中分布均勻程度從大到小依次為黃南、玉樹(shù)、樂(lè)都。最小等位基因頻率為0.183～0.192。3個(gè)地區(qū)喜馬拉雅旱獺Shannon多樣性指數(shù)為2.855～3.322，從大到小依次為黃南、玉樹(shù)、樂(lè)都（表3）。

表3 遺傳統(tǒng)計(jì)量信息Tab.3 The results of genetic statistics

2.3.2 群體結(jié)構(gòu)聚類(lèi)

群體結(jié)構(gòu)分析結(jié)果表明，44 例喜馬拉雅旱獺遺傳背景一致的個(gè)體都聚類(lèi)在一起，說(shuō)明聚類(lèi)結(jié)果準(zhǔn)確。擁有最低交叉驗(yàn)證錯(cuò)誤率（CV error）的分群數(shù)被認(rèn)為是最優(yōu)分群，由圖1 可知，交叉驗(yàn)證錯(cuò)誤率曲線呈現(xiàn)上升趨勢(shì)，當(dāng)分群數(shù)（K）為2 時(shí)，最低交叉驗(yàn)證錯(cuò)誤率值最小，即K=2 為最佳分群數(shù)，3 個(gè)地區(qū)的喜馬拉雅旱獺被劃分為2 個(gè)不同的類(lèi)群，說(shuō)明樣本之間共享2個(gè)基因庫(kù)。

圖1 群體結(jié)構(gòu)（A）及不同K值對(duì)應(yīng)的交叉驗(yàn)證錯(cuò)誤率（B）Fig.1 Group structure（A）and cross validation error rate corresponding to different K values（B）

2.3.3 PCA分析

主成分分析結(jié)果顯示，選取的喜馬拉雅旱獺樣本存在3 個(gè)明顯的分組，樂(lè)都、玉樹(shù)和黃南的喜馬拉雅旱獺各自聚在一起（圖2），部分旱獺親緣關(guān)系較為緊密，而部分旱獺親緣關(guān)系較遠(yuǎn)。

圖2 44例喜馬拉雅旱獺PCA分析Fig.2 PCA analysis of 44 Marmota himalayana

2.3.4 系統(tǒng)進(jìn)化分析

系統(tǒng)進(jìn)化樹(shù)表明，44 例喜馬拉雅旱獺樣本匯聚成3 個(gè)較大的遺傳分支。第1 類(lèi)、第2 類(lèi)是黃南的喜馬拉雅旱獺，第3類(lèi)包括樂(lè)都和玉樹(shù)2個(gè)地區(qū)的喜馬拉雅旱獺（圖3）。

圖3 基于鄰接法的44例喜馬拉雅旱獺的系統(tǒng)發(fā)育樹(shù)Fig.3 The phylogenetic tree of 44 Marmota himalayana based on adjacency method

3 討論

青藏高原特有的地理環(huán)境使之成為高原自然界的原始“本底”，保存了許多珍稀瀕危生物，為開(kāi)展有關(guān)青藏高原生物學(xué)、高原醫(yī)學(xué)、動(dòng)植物遺傳多樣性及保護(hù)等學(xué)科的研究，提供了理想的基地和天然實(shí)驗(yàn)室［18］。喜馬拉雅旱獺作為青藏高原特有物種，成為高原環(huán)境相關(guān)性研究的理想模型［19］。

本研究利用dd-RAD 測(cè)序技術(shù)深入了解不同地區(qū)喜馬拉雅旱獺間的遺傳變異，并對(duì)其多樣性進(jìn)行分析，青海省3 個(gè)地區(qū)的喜馬拉雅旱獺類(lèi)群中共鑒定出19 279 845 個(gè)SNP 位點(diǎn)，說(shuō)明喜馬拉雅旱獺群體具有豐富的遺傳多樣性，這些多樣性是物種適應(yīng)環(huán)境變化的基礎(chǔ)。

從遺傳變異指標(biāo)中可知，樂(lè)都和玉樹(shù)喜馬拉雅旱獺的觀測(cè)雜合度大于期望雜合度，表明這2 個(gè)地區(qū)的喜馬拉雅旱獺雜合過(guò)度，黃南喜馬拉雅旱獺觀測(cè)雜合度小于期望雜合度，表明種群內(nèi)純合子較多，雜合子缺失［20］。近交系數(shù)越大表示基因純化程度越高，玉樹(shù)喜馬拉雅旱獺近交系數(shù)最大。多態(tài)信息量可判斷等位基因的頻率，等位基因越多，多態(tài)信息量越大，3個(gè)地區(qū)中黃南喜馬拉雅旱獺的多態(tài)信息量最大，相對(duì)于其他2 個(gè)地區(qū)的喜馬拉雅旱獺的多態(tài)性更高。有效等位基因數(shù)可在一定程度上反應(yīng)種群變異程度，黃南喜馬拉雅旱獺有效等位基因數(shù)更接近等位基因數(shù)，表明群體中的等位基因分布較為均勻。等位基因頻率用來(lái)表示基因庫(kù)的豐富程度，黃南的最小等位基因頻率最大，表明黃南樣本基因庫(kù)豐富程度高。3 個(gè)地區(qū)多態(tài)信息含量、有效等位基因數(shù)、最小等位基因頻率和多樣性指數(shù)（Shannon）大小均呈現(xiàn)相同規(guī)律。

從群體結(jié)構(gòu)聚類(lèi)分析可知，44 例喜馬拉雅旱獺分為明顯的兩群。從PCA 分析可知，44 例喜馬拉雅旱獺分成明顯的3組，3個(gè)地區(qū)的喜馬拉雅旱獺各自聚類(lèi)在一起。從系統(tǒng)發(fā)育樹(shù)可見(jiàn)，黃南喜馬拉雅旱獺分為2個(gè)遺傳分支，樂(lè)都和玉樹(shù)樣本匯聚成1個(gè)大的遺傳分支，說(shuō)明這2 個(gè)地區(qū)的喜馬拉雅旱獺親緣關(guān)系更近，且處在相同地區(qū)的喜馬拉雅旱獺親緣關(guān)系更近。從地理距離看，黃南和樂(lè)都種群的距離較近，但沒(méi)有完全聚類(lèi)在一起，沒(méi)有呈現(xiàn)出明顯的系統(tǒng)地理分布格局，同在黃南的喜馬拉雅旱獺分為2 個(gè)遺傳分支的原因可能是青藏高原地勢(shì)復(fù)雜，河流山川眾多，以及近年來(lái)青海省旅游業(yè)發(fā)展迅速，人類(lèi)活動(dòng)等對(duì)喜馬拉雅旱獺種群活動(dòng)范圍起到了限制作用。由以上可知，喜馬拉雅旱獺的種群進(jìn)化較復(fù)雜，群體結(jié)構(gòu)可能與地理環(huán)境有相關(guān)性。

綜上所述，本研究基于簡(jiǎn)化基因組測(cè)序，構(gòu)建系統(tǒng)發(fā)育樹(shù)，從基因水平分析不同地區(qū)喜馬拉雅旱獺個(gè)體之間的遺傳進(jìn)化關(guān)系，有助于對(duì)該物種種群多樣性的保護(hù)，為高原野生動(dòng)物疾病的預(yù)防提供新的方向和思路。