鄧伊華，王天嬌，王洪亮，董依萌，劉 欣，邢秀梅

（1.東北林業(yè)大學野生動物與自然保護地學院，哈爾濱 150006；2.中國農(nóng)業(yè)科學院特產(chǎn)研究所,特種經(jīng)濟動物分子生物重點實驗室，長春 130112）

中國養(yǎng)鹿業(yè)歷史悠久，鹿資源豐富，是世界上養(yǎng)鹿最早、鹿產(chǎn)品加工及應用最廣泛的國家之一。中國有8個馬鹿亞種，其中新疆分布有3個馬鹿亞種，分別為塔里木馬鹿（Cervuselaphusyarkandensis）、天山馬鹿（Cervuselaphussongaricus）和阿爾泰馬鹿（Cervuscanadensisasiaticus）[1]。塔里木馬鹿是培育優(yōu)質高產(chǎn)馬鹿的優(yōu)勢資源，在中國近70年的馬鹿養(yǎng)殖業(yè)中做出了巨大貢獻。自20世紀50年代開始，從捕捉野生仔鹿進行馴化和飼養(yǎng)繁殖，到20世紀70年代已經(jīng)完全依靠自繁自育，最終成功選育的第一個馬鹿品種，當?shù)胤Q為塔河馬鹿[2]。現(xiàn)代分子生物學研究表明，塔里木馬鹿為中國馬鹿的祖先，也是世界馬鹿的祖先[3]。塔里木馬鹿獨特的環(huán)境適應性決定它們只能生存在塔里木河流域，其耐粗飼、耐受炎熱和干燥環(huán)境特性一直是學者關注的焦點。由于缺少對塔河馬鹿品種遺傳資源保護的安全意識，養(yǎng)殖者在追求高產(chǎn)的同時，忽視了對塔河馬鹿資源的保護，與天山馬鹿和阿爾泰馬鹿雜交嚴重，受胎率低、幼仔存活率低等問題一直未得到解決。系統(tǒng)開展塔河馬鹿種質資源資源鑒定工作能合理利用塔河馬鹿資源并有效保護塔河馬鹿。

近幾年來，隨著測序技術的不斷發(fā)展，以第二代測序技術（NGS）為主的高通量自動化測序技術在單核苷酸多態(tài)性（SNP）位點的開發(fā)和利用方面逐漸凸顯出巨大優(yōu)勢。在此基礎上，通過全基因組重測序技術對已知基因組序列的不同個體進行基因組重測序，比較參考基因組和全基因組重測序序列，獲得大量的SNP位點、拷貝數(shù)變異（CNV）、插入缺失位點（InDel）和結構變異位點（SV）等遺傳特征[4-6]。全基因組重測序技術在動物變異檢測[7]、性狀定位[8]、基因挖掘[9]、遺傳進化[10]等多個領域均取得了豐碩的成果[11]。Ba等[12]采用ddRAD-seq技術檢測到30個圈養(yǎng)個體（7頭梅花鹿、6頭馬鹿和17頭雜交鹿）的約320 000個全基因組SNPs，通過篩選并觀察雜合度，首次報道了梅花鹿和馬鹿的特異性SNP位點。董世武等[13]通過GBS技術對梅花鹿、馬鹿及其雜交后代（F1和F2） 4個群體共226個樣本進行測序，分別識別出馬鹿特異性SNPs位點474個，梅花鹿特異性SNPs位點558個，為梅花鹿、馬鹿和雜交鹿的鑒別提供了依據(jù)。塔依爾江·麥麥提等[14]利用全基因組重測序技術對塔里木馬鹿干旱環(huán)境適應相關基因進行篩選和研究，成功篩選出塔里木馬鹿過氧化還原酶3（PRDX3）基因。Fan等[15]對500多頭梅花鹿、馬鹿和雜交鹿進行全基因組重測序，篩選得到1 300多萬個SNPs位點，通過生物信息學和芯片設計原則，最終選擇1 000個特異性SNPs位點用于芯片合成，隨機選取266頭鹿進行驗證，結果表明，1K梅花鹿SNP芯片能夠準確地區(qū)分梅花鹿、馬鹿和雜交鹿。從簡化基因組測序篩選特異性位點到利用重測序技術研發(fā)相關基因芯片，測序技術正逐步應用于梅花鹿和馬鹿的鑒別研究、基因挖掘等方面。塔河馬鹿是所有馬鹿種群中規(guī)模最大的群，然而目前種群數(shù)量不足萬只，數(shù)量還在持續(xù)下降，而且人工飼養(yǎng)的馬鹿存在雜交，導致其品質變化，多樣性減少，塔河馬鹿遺傳資源流失嚴重。因此，需要對塔河馬鹿這一寶貴的地方品種資源進行合理的保護與利用，系統(tǒng)地開展塔河馬鹿種質資源鑒定工作。

本研究利用全基因組重測序技術，通過對32份塔河馬鹿進行全基因組測序，結合實驗室前期獲得的塔河馬鹿、天山馬鹿和阿爾泰馬鹿共88份重測序數(shù)據(jù)（未發(fā)表），以組裝到染色體級別的梅花鹿基因組序列作為參考，深入挖掘塔河馬鹿與天山馬鹿、阿爾泰馬鹿之間全基因組水平上的差異，篩選出特異性強，穩(wěn)定性高的塔河馬鹿特異性SNP位點，構建塔河馬鹿分子遺傳標記，為塔河馬鹿種質資源收集、整理、保護和評價工作奠定基礎。

1 材料與方法

1.1 材料

國家特種經(jīng)濟動物資源共享平臺收集并保存的1999年在新疆地區(qū)采集的32份塔河馬鹿血液樣本。

1.2 方法

1.2.1 血液基因組DNA的提取 采用酚-氯仿法提取32份塔河馬鹿血液基因組DNA，提取的DNA用1.0%瓊脂糖凝膠電泳和核酸定量蛋白儀（Agilent 5400，安捷倫）進行質量檢測，經(jīng)檢測合格后-20 ℃保存?zhèn)溆谩?/p>

1.2.2 DNA文庫的構建和測序 檢測合格的DNA樣品用于DNA文庫構建，32份塔河馬鹿DNA分別建庫和測序。建庫和測序由北京諾禾致源科技股份有限公司完成。將DNA樣品通過Covaris超聲波破碎儀隨機打斷，經(jīng)末端修復、加A尾、加測序接頭、純化、PCR擴增等步驟完成整個文庫制備。文庫檢測合格后，把不同文庫按照有效濃度及目標下機數(shù)據(jù)量的需求，利用Illumina HiSeq/MiSeq對32個樣本進行雙末端低深度測序，隨后對原始下機數(shù)據(jù)進行預處理，測序得到原始測序序列（sequenced reads）或者原始數(shù)據(jù)（raw reads：含有帶接頭的、低質量的reads）。將得到的序列進行GC含量分布檢測和測序質量分析。為了保證信息分析質量，對raw reads進行過濾，得到clean reads，基于clean reads進行后續(xù)分析。

1.2.3 測序數(shù)據(jù)比對 使用BWA-MEM軟件[16]將過濾后的clean reads和實驗室前期獲得的重測序數(shù)據(jù)一同與組裝到染色體水平的梅花鹿參考基因組序列（MHL_V1）進行比對，并對本次32份樣品的比對率（mapping rate）、覆蓋度（coverage）和測序深度（depth）進行統(tǒng)計。

1.2.4 基于SNP位點的生物信息學分析

1.2.4.1 SNP位點檢測和注釋 使用SAMTOOLS軟件[17]對BAM文件進行排序，刪除重復讀取的數(shù)據(jù)。利用GATK 4.0.2軟件進行變異檢測，過濾結果輸出到VCF文件保存。使用軟件SNPEff對得到的SNP位點進行注釋，統(tǒng)計SNP位點在33條染色體上的分布。

1.2.4.2 分子進化樹構建 使用SNPhylo軟件對所有樣品進行關系推斷，構建分子進化樹，統(tǒng)計候選特異性SNP位點。

1.2.4.3 候選特異性SNP位點篩選 將1999年收集整理的塔河馬鹿樣品作為一個群體，其余馬鹿作為另一個群體，使用VCFTOOLS[18]計算Fst值，將計算得到的Fst值按照降序排序，設定閾值Fst≥0.25,淘汰不符合條件位點，得到每條染色體剩余SNP數(shù)目用于構建塔河馬鹿遺傳標記。

1.2.4.4 塔河馬鹿特異性遺傳標記構建 統(tǒng)計篩選后每條染色體上的SNP數(shù)目并將所有SNP位點進行編碼，基因型AA編碼為0；基因型AB編碼為1；基因型BB編碼為2，基因型缺失編碼為-9；編碼后的SNP位點數(shù)據(jù)構成矩陣，使用R語言4.0.4版本prcomp函數(shù)[19]進行主成分分析（PCA），計算每條染色體上SNP標記的得分，根據(jù)排名選取每條染色體排名前1 500名的SNP，統(tǒng)計33條染色體篩選出來的SNP位點數(shù)，計算PC1和PC2兩個主成分右奇異向量最大的列。選取前100個SNPs構成特征SNP子集。提取前100個SNPs集合的特征值，計算協(xié)方差矩陣。分別對33條染色體篩選得到的SNP位點和前100個SNPs位點繪制主成分分析聚類圖。 最后統(tǒng)計前100個SNPs在染色體上的分布情況。

2 結 果

2.1 血液基因組DNA質檢結果

1.0%瓊脂糖凝膠電泳檢測結果表明，得到的塔河馬鹿DNA目的條帶單一，無明顯降解和雜質污染。核酸定量蛋白儀檢測其濃度和純度都在要求范圍之內，符合后續(xù)建庫要求。其中DNA樣品總量在0.09182～4.14943 μg之間。

2.2 測序數(shù)據(jù)統(tǒng)計與質量評估

32份塔河馬鹿血液基因組DNA樣品經(jīng)全基因組重測序，共產(chǎn)生raw reads 869 841 762 900 bp，去除低質量序列后得到clean reads 868 791 354 600 bp，平均每個樣品27 149 729 831.25 bp。測序質量Q20≥96.96%、Q30≥92.09%，GC含量在43.29%～46.30%之間，沒有發(fā)生堿基分離的現(xiàn)象。此次建庫測序成功，測序質量較高，獲得的數(shù)據(jù)量、質量和序列長度均符合后續(xù)數(shù)據(jù)分析要求。

2.3 測序數(shù)據(jù)比對情況

以組裝到染色體級別的梅花鹿全基因組為參考序列（MHL_V1），對每個樣品的比對率（mapping rate）、覆蓋度（coverage）和測序深度（depth）進行統(tǒng)計，32份樣品的測序數(shù)據(jù)的平均比對率為98.06%，平均覆蓋度為97.66%，平均深度為6.787。

2.4 基于SNP位點的生物信息學分析

2.4.1 SNP位點的特征分析 將32份樣品測序得到的有效數(shù)據(jù)與試驗前期得到的測序數(shù)據(jù)與參考基因組（MHL_V1）進行比對，過濾后得到20 139 122個高質量的SNPs位點。對每條染色體上的SNP位點分布情況進行統(tǒng)計（圖1）。SNP位點基本隨機均勻分布在每條染色體上，其中4號染色體上分布最多，有1 132 360個，26號染色體上分布最少，有317 705個，SNP在染色體上的分布可能與染色體的長度有關，染色體越長，SNP位點數(shù)量越多。用SNPEff統(tǒng)計所有染色體上變異類型的分布情況表明，1 156 947（6.76%）個SNPs分布在上、下游區(qū)域；基因間的變異最多，有11 118 374個（65.02%）；外顯子區(qū)域有128 027個（0.75%）；內含子區(qū)域有4 701 559個（27.47%）（圖2）。

圖1 每條染色體上的SNP數(shù)目Fig.1 The number of SNP on each chromosome

圖2 變異類型在染色體上的分布Fig.2 Distribution of variant types on chromosomes

2.4.2 基于SNP位點構建分子進化樹 基于SNP位點構建分子進化樹，結果顯示， 塔河馬鹿、天山馬鹿和阿爾泰馬鹿區(qū)分明顯，其中塔河馬鹿單獨匯聚成一支，天山馬鹿和阿爾泰馬鹿聚類位置相近，根據(jù)進化樹聚類位置選擇樣本進行分析，過濾后位點數(shù)為12 050 781（圖3）。

圖3 基于SNPs的分子進化樹Fig.3 Molecular evolutionary tree based on SNPs

2.4.3 候選特異性SNP位點的篩選 將1999年收集整理的塔河馬鹿樣品作為一個群體，其余馬鹿作為另一個群體，使用VCFTOOLS軟件淘汰不符合條件位點，統(tǒng)計每條染色體的候選特異性SNP位點，篩選后的SNP在染色體上的分布結果顯示，剩余位點數(shù)544 717個。其中SNPs位點在4號和5號染色體上分布最多，其余染色體上位點分布較為均勻（圖4）。

圖4 Fst篩選后的SNP在染色體上的分布Fig.4 Distribution of SNP on chromosomes after Fst screening

2.4.4 塔河馬鹿特異性遺傳標記的構建 從每條染色體選取排名前1 500的SNPs位點，33條染色體共篩選得到49 500個SNPs位點。49 500個SNP位點及其前100個SNPs位點主成分分析結果表明，當位點數(shù)下降至100時，1999年的塔河馬鹿依然能準確與天山馬鹿和阿爾泰馬鹿區(qū)分開，區(qū)分效力基本沒有下降。大多數(shù)阿爾泰馬鹿分布PC1的右邊，少量分布與天山馬鹿相近，群體分布較為分散，個體血源可能不純，此外，1999年的塔河馬鹿和近年的塔河馬鹿在聚類圖上區(qū)分明顯，大部分現(xiàn)在的塔河馬鹿與天山馬鹿分布相近，天山馬鹿分布在PC1的左上角，個體聯(lián)系較為緊密（圖5、6）。

圖5 49 500個SNPs的主成分分析Fig.5 PCA of 49 500 SNPs

圖6 前100個SNPs的主成分分析Fig.6 PCA of top 100 SNPs

通過計算篩選后選取前100個SNPs位點作為塔河馬鹿特異性分子遺傳標記，前100個SNPs位點在染色體上的分布可以看到，15和30號染色體的分布最多，其中15號染色體上分別有3個C/C、2個T/T,G/G和A/A基因型各1個，30號染色體有4個T/T，C/C、G/G、A/A基因型各1個。1、4、5和23號染色體上無分布。得到的塔河馬鹿優(yōu)勢基因型在大部分染色體上均有分布，分布情況良好（圖7）。

圖7 前100個SNPs在染色體上的分布Fig.7 Distribution of top 100 SNPs on chromosomes

3 討 論

中國家養(yǎng)馬鹿遺傳背景單一，塔河馬鹿由于高強度的人工選擇和近交導致品種退化現(xiàn)象日益嚴重。近年來,隨著雜交改良的推進，產(chǎn)茸重量上的絕對優(yōu)勢致使雜交鹿在養(yǎng)殖上占據(jù)的市場空間越來越大，養(yǎng)殖者盲目追求鹿茸產(chǎn)量而無序雜交，造成了塔河馬鹿純種資源的匱乏。雜交鹿從肉眼上難以與塔河馬鹿區(qū)分，給塔河馬鹿的保種、育種工作帶來了極大困難。目前二代測序技術憑借其高通量和低成本的優(yōu)勢在動物基因組測序中得到了廣泛的應用，在全基因組水平上篩選特異性SNP位點[20]，結合相關檢測技術和生物信息學分析，推斷動物品種構成，預測雜種或純種[21]。朱濤等[22]選取7個群體的56只鴨進行重測序，篩選出686個存在于特定群體中的SNPs和InDel并隨機選取7個進行了驗證，結果完全符合預期。岳陽等[23]對小梅山豬和其他11個豬品種進行簡化基因組測序，篩選出5個小梅山豬的特異性SNPs位點并進行驗證，建立小梅山豬的SNP分子標記鑒別方法。本研究以大量馬鹿的重測序數(shù)據(jù)為基礎，采用本團隊組裝到染色體的梅花鹿基因組作為參考，更準確地檢測塔河馬鹿SNP的變異情況，利用生物信息學軟件分析SNP位點，提高計算效率，降低研發(fā)成本，為后續(xù)塔河馬鹿鑒別的相關技術開發(fā)及應用提供多樣化的選擇。

本試驗對32份塔河馬鹿DNA樣本進行全基因組重測序，結合試驗前期的56份馬鹿重測序數(shù)據(jù)，與組裝到染色體級別的梅花鹿參考基因組進行比對，通過生物信息學分析并得出結論。結果表明，32個塔河馬鹿樣品測序和比對情況良好，數(shù)據(jù)量符合預期。篩選的SNP位點基本分布在基因間和內含子區(qū)域，與Wang等[24]研究廣西地方鴨中篩選的SNP位點變異分布結果一致，SNP的分布差異可能與DNA在不同基因區(qū)域的功能差異有關[25]。通過構建分子進化樹能進一步分析馬鹿的遺傳多樣性并挖掘個體間的差異[26]，從分子進化樹上可以看到天山馬鹿和阿爾泰馬鹿分支相近，塔河馬鹿單獨聚成一支，與天山馬鹿和阿爾泰馬鹿相距較遠，與蘇瑩等[27]、Shao等[28]研究結果一致。王小鵬[29]利用60K基因芯片結合主成分分析和系統(tǒng)進化樹篩選特異位點構建豬特異性遺傳標簽,劉月麗等[30]利用主成分分析和隨機森林算法結合篩選高質量SNP位點用于羊的品種分類。從聚類圖中可以看到1999年的塔河馬鹿與其他馬鹿區(qū)分明顯，與現(xiàn)在未知血源狀況的塔河馬鹿也能準確區(qū)分，進一步說明現(xiàn)在的塔河馬鹿樣本不純，已經(jīng)進行了雜交，混有天山馬鹿或阿爾泰馬鹿的血源。通過篩選得到的塔河馬鹿特異性SNP位點，不僅將早年的塔河馬鹿與天山馬鹿和阿爾泰馬鹿區(qū)分開，還能與現(xiàn)在雜交的塔河馬鹿區(qū)分。SNP標記作為第3代分子標記被認為是目前應用前景最好的遺傳標記，其在動植物基因組中均分布廣泛，穩(wěn)定性高且富有代表性。SNP標記不再以DNA片段的長度變化作為檢測手段，而是直接以序列變異作為標記，完全摒棄了凝膠電泳而采用測序技術或最新的DNA芯片技術。本研究采用全基因組重測序技術，測序深度和覆蓋度高，得到的塔河馬鹿特異性位點范圍廣、準確性高，另一方面選取的參考群體范圍廣，32份塔河馬鹿重測序數(shù)據(jù)結合之前56份馬鹿的重測序數(shù)據(jù)（11頭塔河馬鹿、11頭天山馬鹿、34頭阿爾泰馬鹿）使特異性位點的篩選結果更加準確。

4 結 論

本研究結果表明，32份塔河馬鹿重測序數(shù)據(jù)過濾后篩選得到了20 139 122個高質量的SNPs位點，建樹過濾后SNP位點數(shù)為12 050 781個，過濾后共得到544 717個SNPs位點用于塔河馬鹿分子遺傳標記構建，最終篩選得到了100個具有高度多態(tài)性、穩(wěn)定性良好的SNPs位點，結果為塔河馬鹿的精準鑒別和核心種質的篩選提供理論依據(jù)，為塔河馬鹿種質資源的收集、整理、保存和評價工作的順利開展奠定堅實基礎。