基于重測序技術(shù)的塔河馬鹿特異性SNP位點篩選

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

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

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

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

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

1 材料與方法

1.1 材料

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

1.2 方法

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

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

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

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

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

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

2 結(jié) 果

2.1 血液基因組DNA質(zhì)檢結(jié)果

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

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

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

2.4.1 SNP位點的特征分析 將32份樣品測序得到的有效數(shù)據(jù)與試驗前期得到的測序數(shù)據(jù)與參考基因組（MHL_V1）進行比對，過濾后得到20 139 122個高質(zhì)量的SNPs位點。對每條染色體上的SNP位點分布情況進行統(tǒng)計（圖1）。SNP位點基本隨機均勻分布在每條染色體上，其中4號染色體上分布最多，有1 132 360個，26號染色體上分布最少，有317 705個，SNP在染色體上的分布可能與染色體的長度有關(guān)，染色體越長，SNP位點數(shù)量越多。用SNPEff統(tǒng)計所有染色體上變異類型的分布情況表明，1 156 947（6.76%）個SNPs分布在上、下游區(qū)域；基因間的變異最多，有11 118 374個（65.02%）；外顯子區(qū)域有128 027個（0.75%）；內(nèi)含子區(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位點構(gòu)建分子進化樹 基于SNP位點構(gòu)建分子進化樹，結(jié)果顯示， 塔河馬鹿、天山馬鹿和阿爾泰馬鹿區(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在染色體上的分布結(jié)果顯示，剩余位點數(shù)544 717個。其中SNPs位點在4號和5號染色體上分布最多，其余染色體上位點分布較為均勻（圖4）。

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

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

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

4 結(jié) 論

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