基于簡化基因組測序技術的喜馬拉雅旱獺遺傳多樣性分析

李 優，南新營，趙坤鈺，李耀東

（青海大學，西寧，810016）

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

1 研究方法

1.1 樣本采集

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

1.2 基因組DNA提取與質量控制

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

1.3 dd-RAD文庫構建

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

1.4 全基因組SNP篩查、分型與質控

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

1.5 數據統計與分析

利用BWA 軟件［13］、GATK 軟件、ADMIXTURE 軟件［14］、GCTA 軟件［15］和Fast Tree 軟件［16］分別對SNP注釋、遺傳分化系數統計，并對SNP位點觀測雜合度（Ho）、期望雜合度（He）、近交系數（Fis）、多樣性指數（Shannon）、多態信息含量（PIC）、有效等位基因數（Ne）和最小等位基因頻率（MAF）7 個遺傳多樣性參數進行分析。基于得到的SNP，使用ADMIXTURE軟件計算群體結構，以分群數（K值）為2～10 進行聚類，做出群體結構聚類分析。使用GCTA 軟件，通過主成分分析（PCA）獲得樣本遺傳距離遠近關系，輔助進化分析。

2 結果

2.1 測序數據統計

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

表1 喜馬拉雅旱獺dd-RAD測序數據Tab.1 dd-RAD sequencing data of Marmota himalayana

2.2 SNP檢測與位點開發

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

表2 喜馬拉雅旱獺SNP統計結果Tab.2 Statistical results of SNP of Marmota himalayana

2.3 遺傳進化分析

2.3.1 遺傳統計量分析

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

表3 遺傳統計量信息Tab.3 The results of genetic statistics

2.3.2 群體結構聚類

群體結構分析結果表明，44 例喜馬拉雅旱獺遺傳背景一致的個體都聚類在一起，說明聚類結果準確。擁有最低交叉驗證錯誤率（CV error）的分群數被認為是最優分群，由圖1 可知，交叉驗證錯誤率曲線呈現上升趨勢，當分群數（K）為2 時，最低交叉驗證錯誤率值最小，即K=2 為最佳分群數，3 個地區的喜馬拉雅旱獺被劃分為2 個不同的類群，說明樣本之間共享2個基因庫。

圖1 群體結構（A）及不同K值對應的交叉驗證錯誤率（B）Fig.1 Group structure（A）and cross validation error rate corresponding to different K values（B）

2.3.3 PCA分析

主成分分析結果顯示，選取的喜馬拉雅旱獺樣本存在3 個明顯的分組，樂都、玉樹和黃南的喜馬拉雅旱獺各自聚在一起（圖2），部分旱獺親緣關系較為緊密，而部分旱獺親緣關系較遠。

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

2.3.4 系統進化分析

系統進化樹表明，44 例喜馬拉雅旱獺樣本匯聚成3 個較大的遺傳分支。第1 類、第2 類是黃南的喜馬拉雅旱獺，第3類包括樂都和玉樹2個地區的喜馬拉雅旱獺（圖3）。

圖3 基于鄰接法的44例喜馬拉雅旱獺的系統發育樹Fig.3 The phylogenetic tree of 44 Marmota himalayana based on adjacency method

3 討論

青藏高原特有的地理環境使之成為高原自然界的原始“本底”，保存了許多珍稀瀕危生物，為開展有關青藏高原生物學、高原醫學、動植物遺傳多樣性及保護等學科的研究，提供了理想的基地和天然實驗室［18］。喜馬拉雅旱獺作為青藏高原特有物種，成為高原環境相關性研究的理想模型［19］。

本研究利用dd-RAD 測序技術深入了解不同地區喜馬拉雅旱獺間的遺傳變異，并對其多樣性進行分析，青海省3 個地區的喜馬拉雅旱獺類群中共鑒定出19 279 845 個SNP 位點，說明喜馬拉雅旱獺群體具有豐富的遺傳多樣性，這些多樣性是物種適應環境變化的基礎。

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

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

綜上所述，本研究基于簡化基因組測序，構建系統發育樹，從基因水平分析不同地區喜馬拉雅旱獺個體之間的遺傳進化關系，有助于對該物種種群多樣性的保護，為高原野生動物疾病的預防提供新的方向和思路。