基于重測(cè)序數(shù)據(jù)的昌都黑山羊遺傳多樣性及群體結(jié)構(gòu)分析

阿旺措吉，仁青措姆，黃舒泓，周先輝，德 吉，索朗達(dá)，王小龍，吳玉江，巴 貴*

(1.西藏自治區(qū)農(nóng)牧科學(xué)院畜牧獸醫(yī)研究所，西藏 拉薩 850000;2.西北農(nóng)林科技大學(xué) 動(dòng)物科技學(xué) 院，陜西 楊凌 712100)

藏山羊(Caprahircus)是中國(guó)獨(dú)特的種質(zhì)資源，在青藏高原的農(nóng)業(yè)、文化、歷史、經(jīng)濟(jì)乃至宗教等方面都起到了重要的作用，其主要分布在西藏自治區(qū)、四川省阿壩和甘孜自治州、青海省玉樹(shù)自治州、甘肅省甘南、果洛自治州及新疆部分地區(qū)等[1]，是中國(guó)最古老的羊品種之一，最早記錄于公元前3 300-2 000年[2]，外國(guó)商人稱(chēng)謂的克什米爾山羊(Cashmere goat)即為藏山羊。長(zhǎng)期以來(lái)，由于復(fù)雜多樣的環(huán)境、高海拔和極端氣候，藏山羊具備了獨(dú)特的遺傳特性[3]，加之其群體大、山羊絨產(chǎn)量高、肉品質(zhì)及皮革品質(zhì)優(yōu)良，現(xiàn)已成為中國(guó)寶貴的基因庫(kù)[4]。

全基因組重測(cè)序基于高通量基因測(cè)序技術(shù)，并依托計(jì)算機(jī)對(duì)研究群體基因組與已發(fā)表參考基因組進(jìn)行比對(duì)，可以在基因水平對(duì)不同目標(biāo)間的基因差異進(jìn)行對(duì)比。目前，已有大量使用基因組重測(cè)序的方法與研究。如李旭靜等[5]利用全基因組重測(cè)序?qū)d羊經(jīng)濟(jì)性狀基因進(jìn)行篩選，注釋了多個(gè)與毛、乳、肉產(chǎn)量相關(guān)的選擇區(qū)域。王統(tǒng)苗等[6]使用全基因組重測(cè)序技術(shù)對(duì)鴨群的遺傳資源進(jìn)行結(jié)構(gòu)分析，將多個(gè)不同地區(qū)的鴨種資源進(jìn)行分類(lèi)與合并。

前人通過(guò)應(yīng)用分子遺傳學(xué)標(biāo)記、基因組學(xué)、生物信息學(xué)等[7]方法，積累了大量關(guān)于藏山羊遺傳資源多樣性的資料。如王杰等[8]采用SSR分子標(biāo)記對(duì)高原型藏山羊、山谷型藏山羊進(jìn)行遺傳多態(tài)性研究，發(fā)現(xiàn)藏山羊群體遺傳多態(tài)性豐富且對(duì)生存環(huán)境的適應(yīng)能力較強(qiáng)。王永等[9]應(yīng)用ISSR標(biāo)記分析了西藏日土藏山羊的遺傳多樣性，結(jié)果顯示該品種藏山羊的遺傳多態(tài)性較豐富，個(gè)體間雖有差異但同質(zhì)性較好。鄧娟等[10]通過(guò)對(duì)細(xì)胞色素 b基因(Cytb)全序列進(jìn)行擴(kuò)增和測(cè)序，研究了西藏地區(qū)藏山羊的母系起源及遺傳多樣性。

本研究收集西藏昌都3個(gè)不同地區(qū)的藏山羊群體血液DNA并測(cè)序。用SNP位點(diǎn)分析、遺傳多樣性分析、群體遺傳結(jié)構(gòu)分析、遺傳關(guān)系分析及主成分分析等多種分析方式，從基因組水平解釋了藏山羊的種質(zhì)特性，也為后續(xù)藏山羊重要生產(chǎn)性狀的基因定位、分子育種及遺傳資源的保護(hù)和利用提供重要的理論依據(jù)。

1 材料與方法

1.1 試驗(yàn)動(dòng)物

試驗(yàn)樣本采自西藏昌都市邊壩縣邊壩鎮(zhèn)洛亞碼羊場(chǎng)(LYM，n=9)、普玉一村(PY，n=9)及芒康縣朱八龍鄉(xiāng)(MK，n=10)，共28個(gè)樣本。采集健康個(gè)體的靜脈抗凝全血樣本，-80 ℃冰箱中保存?zhèn)溆谩?/p>

1.2 試驗(yàn)方法

1.2.1 DNA提取及建庫(kù) 血液基因組DNA樣品制備：使用購(gòu)買(mǎi)自華大基因(BGI)的試劑盒，按照說(shuō)明書(shū)將-80 ℃藏山羊靜脈血解凍后提取基因組DNA，并冷藏保存。

DNA文庫(kù)的構(gòu)建：使用Nanodrop超微量核酸分析儀檢測(cè)血液基因組DNA樣品的濃度，并使用凝膠電泳檢測(cè)DNA樣品的純度。使用超聲波儀打斷0.001 mg基因組DNA，選擇0.2～0.4 mg大小基因片段。在DNA片段3'端鏈接DNA接頭。使用PCR擴(kuò)增加入DNA接頭的片段并提純回收，將雙鏈PCR產(chǎn)物解鏈，破壞未被環(huán)化的DNA分子，獲得DNA文庫(kù)。

1.2.2 測(cè)序及參考基因組比對(duì) 使用華大基因的DNBSEQ-T7基因測(cè)序儀器平臺(tái)進(jìn)行測(cè)序，測(cè)序工作由華大基因公司完成。下機(jī)后的Raw reads使用trimmomatic軟件進(jìn)行質(zhì)控分析，測(cè)序數(shù)據(jù)過(guò)濾得到Clean reads。使用BWA軟件將測(cè)序數(shù)據(jù)比對(duì)到參考基因組，結(jié)合Samtools軟件轉(zhuǎn)換數(shù)據(jù)格式，Picard軟件對(duì)比對(duì)文件去重，得到最終分析所用文件。

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

1.3.1 SNP位點(diǎn)分析及質(zhì)控 處理后比對(duì)文件，通過(guò)GATK軟件的HaplotypeCaller模塊進(jìn)行SNP位點(diǎn)鑒定，之后經(jīng)過(guò)Bcftools及VariantFiltration模塊進(jìn)行質(zhì)控和過(guò)濾，再通過(guò)SelectVariants模塊以及Bcftools合并各個(gè)體的GVCF文件[11]。

1.3.2 遺傳多樣性分析 利用PLINK軟件，設(shè)定參數(shù)滑動(dòng)窗口大小為50 kb，步長(zhǎng)為20 kb，計(jì)算3個(gè)群體的核苷酸多樣性PI值[12]。通過(guò)PopLDdecay軟件計(jì)算3個(gè)群體的LD衰減，得到群體內(nèi)的LD連鎖不平衡信息[13]，最后利用R語(yǔ)言腳本以及perl語(yǔ)言腳本處理結(jié)果生成可視化圖形。

1.3.3 群體遺傳結(jié)構(gòu)分析 將包含SNP信息的VCF文件通過(guò)PLINK軟件轉(zhuǎn)換為BED及PED文件[12]，輸入到Admixture軟件中推斷群體遺傳結(jié)構(gòu)，設(shè)置群體亞群數(shù)目K=2到12計(jì)算交叉驗(yàn)證誤差(Cross-validation error)，通過(guò)繪制Cross-validation error圖確定最佳分群數(shù)[14]，利用R語(yǔ)言腳本及ggplot2軟件包對(duì)結(jié)果生成可視化圖形[15]。

1.3.4 遺傳關(guān)系分析及主成分分析 利用PLINK軟件生成PED、MAP及BED文件，根據(jù)LD進(jìn)行過(guò)濾，通過(guò)編寫(xiě)的perl腳本計(jì)算遺傳距離，生成Neighborhood-Join法構(gòu)建的進(jìn)化樹(shù)文件，得到的meg文件輸入到MEGA軟件生成可視化的NJ樹(shù)[16]，利用在線網(wǎng)站Interactive Tree Of Life(https://itol.embl.de/)進(jìn)行進(jìn)化樹(shù)可視化的美化[17]。

使用EIGENSOFT的smartpca程序進(jìn)行主成分分析[18]，得到的主成分信息利用R語(yǔ)言腳本及ggplot2軟件包生成可視化圖形[15]。

2 結(jié)果與分析

2.1 藏山羊重測(cè)序數(shù)據(jù)質(zhì)量分析及參考基因組比對(duì)

本次研究利用華大基因DNBSEQ-T7平臺(tái)對(duì)28個(gè)藏山羊個(gè)體血液樣本進(jìn)行全基因組重測(cè)序(表1)。所得下機(jī)數(shù)據(jù)初始Raw reads數(shù)量范圍為155 118 313～272 533 434，經(jīng)過(guò)軟件質(zhì)控過(guò)濾，Clean reads數(shù)量范圍在151 935 132～267 865 582，過(guò)濾比率在97.83%～98.35%之間，過(guò)濾比率是指數(shù)據(jù)清洗后reads數(shù)量與原始數(shù)據(jù)(Raw reads)之間的比值，比值越高則測(cè)序質(zhì)量越高。Q30堿基比率均高于95%，測(cè)序質(zhì)量Q值指的是測(cè)序過(guò)程堿基識(shí)別(Base Calling)過(guò)程中，對(duì)所識(shí)別的堿基給出的錯(cuò)誤概率，Q30即為堿基錯(cuò)誤率0.1%，正確率99.9%，因此Q30是評(píng)估測(cè)序質(zhì)量標(biāo)準(zhǔn)的重要一環(huán)，且GC含量也位于40%～50%之間，說(shuō)明樣本的建庫(kù)良好，測(cè)序質(zhì)量達(dá)到重測(cè)序分析標(biāo)準(zhǔn)。

表1 28個(gè)藏山羊個(gè)體重測(cè)序數(shù)據(jù)質(zhì)量控制信息

使用bwa程序mem模塊將28個(gè)藏山羊重測(cè)序數(shù)據(jù)比對(duì)到山羊參考基因組ARS1(NCBI accession: GCF_001704415.1)，比對(duì)情況如表2所示。使用測(cè)序數(shù)據(jù)質(zhì)控后reads數(shù)與參考基因組reads數(shù)的匹配比對(duì)率(Mapping rate)以及平均覆蓋度(Mean coverage)作為評(píng)估數(shù)據(jù)可信度、以及分析正確度的指標(biāo)，通常比對(duì)率與覆蓋度越高，則分析數(shù)據(jù)的可信度越高。比對(duì)上reads數(shù)在303 831 381～535 839 901之間，比對(duì)率在99.61%～99.93%之間，平均覆蓋度在15.2262×～26.5768×之間，說(shuō)明比對(duì)情況良好，可用于后續(xù)SNP位點(diǎn)分析及群體結(jié)構(gòu)分析。

表2 28個(gè)藏山羊個(gè)體重測(cè)序數(shù)據(jù)與參考基因組比對(duì)情況

2.2 藏山羊全基因組SNP位點(diǎn)鑒定

經(jīng)前期質(zhì)控比較及SNP篩選，在3個(gè)采樣群體28個(gè)個(gè)體中共檢測(cè)到13 726 331個(gè)SNP位點(diǎn)，不同的染色體上SNP位點(diǎn)數(shù)目見(jiàn)表3。由于SNP位點(diǎn)是DNA序列上發(fā)生堿基改變的位點(diǎn)，這些位點(diǎn)由于突變的不確定性而在不同群體與個(gè)體中具有鮮明特點(diǎn)?？傮w來(lái)說(shuō)，鑒定出的SNP位點(diǎn)數(shù)量與染色體的長(zhǎng)度呈正相關(guān)。

表3 28個(gè)藏山羊相對(duì)于參考基因組的SNP位點(diǎn)信息

2.3 藏山羊遺傳多樣性分析

對(duì)3個(gè)群體進(jìn)行核苷酸多樣性分析，最終得到LYM、MK、PY3個(gè)采樣群體核苷酸多樣性(PI)分別為0.00182、0.00178、0.00183，除個(gè)別位點(diǎn)存在差異外，整體分布情況接近，差異極小，表明3個(gè)采樣群體SNP情況接近，如圖1所示。同核苷酸多樣性分析結(jié)果類(lèi)似，3個(gè)采樣群體的連鎖不平橫衰減速度也極為相近，整體趨勢(shì)及進(jìn)化速度較為一致，證明采樣群體間關(guān)系緊密。

圖1 昌都黑山羊采樣群體核苷酸多樣性

2.4 藏山羊群體遺傳結(jié)構(gòu)

使用admixture分析對(duì)3個(gè)采樣群體，以及阿里、班戈、措勤、林芝、尼瑪、那曲、日土、山南共8個(gè)西藏其它地區(qū)藏山羊品種或群體進(jìn)行分析。由于群體遺傳結(jié)構(gòu)是由群體的表型種類(lèi)、頻率、基因種類(lèi)、基因型種類(lèi)及頻率等組成的。而地理環(huán)境因素之間的差異，導(dǎo)致群體間產(chǎn)生不同的變異與分化，根據(jù)變異的數(shù)量與方式不同，可以將群體按照親緣關(guān)系劃分為不同亞群。對(duì)最優(yōu)亞群數(shù)K值的評(píng)估，結(jié)果如圖2所示，當(dāng)交叉驗(yàn)證誤差到達(dá)最小值時(shí)，K=5為假定的最佳分群數(shù)目。K=5時(shí)群體遺傳結(jié)構(gòu)如圖3所示，3個(gè)采樣群體的群體遺傳結(jié)構(gòu)類(lèi)似，但LYM采樣群體及PY采樣群體存在一定的雜交特征，這可能是由于羊群組群時(shí)混合其他群體導(dǎo)致的。

圖2 不同K值下的交叉驗(yàn)證誤差值分布

圖3 K=5時(shí)藏山羊種群遺傳結(jié)構(gòu)圖

2.5 藏山羊遺傳關(guān)系及主成分分析

通過(guò)加入尼瑪(NM)、措勤(CQ)、林芝(LZ)、日土(RT)、山南(SN)、阿里(AL)、班戈(BG)、那曲(NQ)共8個(gè)西藏其它地區(qū)藏山羊品種或群體，同本研究的3個(gè)采樣群體進(jìn)行Neighborhood-Join法構(gòu)建進(jìn)化樹(shù)以及主成分分析。NJ樹(shù)的拓?fù)浣Y(jié)構(gòu)如圖4所示，本研究中3個(gè)采樣群體所在分支離西藏其他地區(qū)藏山羊群體遺傳分支距離較遠(yuǎn)，且3個(gè)采樣群體共同聚為一支，群體間親緣關(guān)系更近，可能為同一群體。

圖4 藏山羊NJ進(jìn)化樹(shù)分析結(jié)果

使用主成分分析，3個(gè)采樣群體間相互距離較近，并與其他藏山羊群體明顯分離，如圖5所示。進(jìn)化樹(shù)分析以及主成分分析結(jié)果與遺傳結(jié)構(gòu)分析結(jié)果一致，說(shuō)明本研究3個(gè)采樣群體親緣關(guān)系以及遺傳距離較近，為共同品種藏山羊。

圖5 藏山羊主成分分析結(jié)果

3 討 論

在公元前3 000-公元前2 000年便有了關(guān)于藏山羊的記載[2]，由于西藏地區(qū)海拔高、溫度低的環(huán)境特點(diǎn)，西藏地區(qū)山羊的適應(yīng)性與低海拔地區(qū)山羊品種差異明顯，因此藏山羊?qū)τ谥袊?guó)其他山羊品種的生產(chǎn)性能以及相關(guān)環(huán)境適應(yīng)性研究具有重要作用。同時(shí)，西藏高原地區(qū)的地理差異導(dǎo)致不同地區(qū)之間存在長(zhǎng)期的地理隔離[19]，可能在藏山羊種群內(nèi)產(chǎn)生不同亞種。王杰等[8]、王永等[9]和秦國(guó)慶等[1]采用分子標(biāo)記法對(duì)西藏地區(qū)山羊的種群多樣性進(jìn)行分析，發(fā)現(xiàn)不同地區(qū)藏山羊群體差異較大，具有較高的種群內(nèi)多樣性。

遺傳多樣性是生物多樣性的一個(gè)層次，是某地區(qū)內(nèi)所有生物遺傳信息的集合，是物種在長(zhǎng)期進(jìn)化過(guò)程中不斷發(fā)生變異所產(chǎn)生的特征。遺傳多樣性越豐富，自然環(huán)境的抗逆性也越強(qiáng)。本研究利用了全基因組SNP篩選獲得了高質(zhì)量的SNP集合，并利用此集合進(jìn)行了核苷酸多樣性分析、連鎖不平橫衰減計(jì)算。并加入了西藏地區(qū)其他品種山羊SNP數(shù)據(jù)集進(jìn)行群體遺傳結(jié)構(gòu)分析以及進(jìn)化樹(shù)的構(gòu)建，研究結(jié)果對(duì)于中國(guó)的遺傳資源保護(hù)具有重要意義。

昌都黑山羊是西藏地區(qū)東部地區(qū)特有的高原型山羊群體，具有抗逆性強(qiáng)、板皮厚滿(mǎn)、皮毛質(zhì)軟細(xì)小、瘦肉率較高、肉質(zhì)鮮嫩、膻味小等優(yōu)良性狀，在當(dāng)?shù)匮虍a(chǎn)業(yè)發(fā)展過(guò)程中發(fā)揮著重要作用。本研究使用昌都黑山羊的3個(gè)群體(LYM、PY和MK)共28個(gè)藏山羊個(gè)體的重測(cè)序數(shù)據(jù)，經(jīng)遺傳多樣性分析，3個(gè)群體的SNP分布與數(shù)量均相似，且連鎖不平衡的整體趨勢(shì)與衰減速度也較為一致，可以認(rèn)定昌都地區(qū)3個(gè)群體為同一品種。后續(xù)通過(guò)加入8個(gè)西藏地區(qū)其他品種的SNP數(shù)據(jù)進(jìn)行群體遺傳結(jié)構(gòu)分析、NJ法進(jìn)化樹(shù)構(gòu)建以及主成分分析，發(fā)現(xiàn)昌都黑山羊的3個(gè)群體與西藏地區(qū)其他品種遺傳距離較遠(yuǎn)且分化明顯。

綜上所述，3地采樣的昌都黑山羊群體為同一品種，且品種內(nèi)存在一定的雜交。遺傳進(jìn)化關(guān)系表明，該品種與西藏其他地區(qū)的山羊品種遺傳分化明顯，可被認(rèn)為是獨(dú)立的西藏山羊群體，為后續(xù)新遺傳資源的鑒定提供了重要的理論支撐。本研究通過(guò)28個(gè)昌都黑山羊重測(cè)序數(shù)據(jù)獲得了高精度的SNP數(shù)據(jù)集，可用于后續(xù)昌都黑山羊分子育種及生物信息學(xué)分析，如BLUP育種值計(jì)算、種間基因滲入研究、環(huán)境適應(yīng)性研究等，為發(fā)展昌都地區(qū)優(yōu)勢(shì)特色畜牧業(yè)、提高昌都黑山羊的良種化程度提供理論依據(jù)。

4 結(jié) 論

本研究發(fā)現(xiàn)昌都黑山羊具有豐富的遺傳多樣性；昌都地區(qū)3個(gè)黑山羊群體內(nèi)親緣關(guān)系與遺傳結(jié)構(gòu)較近，且相較西藏其他8個(gè)地區(qū)藏山羊群體遺傳分化較大。研究結(jié)果對(duì)于中國(guó)家畜遺傳資源保護(hù)具有重要意義。