王 佟，馬曉明，余道寧，馬超凡，潘和平，梁春年，扎西塔，閻 萍

（1. 中國農(nóng)業(yè)科學(xué)院蘭州畜牧與獸藥研究所，甘肅省牦牛繁育工程重點實驗室，蘭州 730050；2. 農(nóng)業(yè)農(nóng)村部青藏高原畜禽遺傳育種重點實驗室，蘭州 730050；3. 西北民族大學(xué)生命科學(xué)與工程學(xué)院，蘭州 730030；4. 青海省海北州祁連縣畜牧獸醫(yī)工作站，祁連 810400）

牦牛是重要的大型哺乳動物［1］，主要分布在3 000～5 000 m 的高海拔地區(qū)［2］。作為我國高海拔地區(qū)畜牧業(yè)中的優(yōu)勢畜種［3］，牦牛為應(yīng)對高原地區(qū)的極端環(huán)境和天然草場的季節(jié)性變化形成了適應(yīng)能力強和生命力頑強等的各種適應(yīng)機制［4］，牧民的衣食住行都與其緊密相關(guān)［5-6］。由于牧民對牦牛品種繁殖等知識相對缺乏，飼養(yǎng)管理方式落后，造成了牦牛的選育水平普遍不高［7］，進而限制了牦牛的生長發(fā)育，這就使得牦牛的經(jīng)濟價值不能得到應(yīng)有的體現(xiàn)。目前《國家畜禽遺傳資源品種名錄（2021 年版）》收錄的地方牦牛品種共有18 個，培育牦牛品種2 個。當(dāng)前，隨著第三次全國畜禽遺傳資源普查工作的進行，會有更多的地方牦牛品種被鑒定收錄。青海省位于“世界屋脊”青藏高原的東北部，省內(nèi)山脈高聳，地形多樣［8］，這種特殊的地理位置和多變的氣候，孕育了多個牦牛遺傳資源。祁連縣隸屬于青海省海北藏族自治州，位于祁連山中部的群山峻嶺和盆地之中［9］，平均海拔3 169m，屬于高原大陸性氣候，地貌多樣，適合畜牧業(yè)發(fā)展［10］。祁連牦牛作為祁連縣境內(nèi)的牦牛類群之一，具有區(qū)別于其他牦牛類群的獨特外貌特征和生產(chǎn)性能。祁連牦牛體質(zhì)堅實、結(jié)構(gòu)緊湊勻稱、體軀深長、頭粗而長（公牛粗短）、多有角；其肉色深紅、呈明顯的大理石紋。在當(dāng)?shù)鬲毺氐牡乩憝h(huán)境和氣候下，當(dāng)?shù)仫曫B(yǎng)的家牦牛與祁連山型野牦牛進行了較為頻繁的基因交流，從而形成了一個獨特的牦牛遺傳資源——祁連牦牛［11］。

全基因組重測序是在已知基因組序列的前提下，對物種群體內(nèi)不同個體的基因組序列進行測序，然后與參考基因組進行比對，根據(jù)比對信息挖掘群體或個體間的變異信息［12］，進而解析群體結(jié)構(gòu)，挖掘重要基因等。近年來，隨著高通量測序技術(shù)的不斷發(fā)展，全基因組重測序也在各領(lǐng)域得到了廣泛的應(yīng)用，并逐漸代替了其他測序手段［13-14］。目前，在牛［15-16］和羊［17-18］等各類畜禽中已經(jīng)有大量使用基因組重測序的方法與研究。如Kang 等［19］對6個牦牛群體進行了重測序研究，結(jié)果顯示6 個牦牛群體存在明顯差異，并認(rèn)為科才牦牛是中國牦牛的一個獨立群體。Chai 等［20］對我國31 個牦牛群體進行了群體基因組學(xué)研究，發(fā)現(xiàn)牦牛最早的馴化發(fā)生在青藏高原東南部，然后向青藏高原西部和東北部擴散到四川、甘肅和青海。鮑麒等［21］利用全基因組重測序評估了肅南牦牛的遺傳結(jié)構(gòu)，并注釋了多個與組織分化、胚胎發(fā)育和肉品質(zhì)等相關(guān)的選擇區(qū)域。李在文等［22］利用GBS 構(gòu)建了麥洼牦牛3 個保種群的系譜，為麥洼牦牛的后續(xù)保種工作提供了基礎(chǔ)。

當(dāng)前，已有研究者對其他牦牛群體進行了全基因組水平的研究［23-25］。但是對祁連牦牛和其他牦牛資源在基因組水平上還未進行系統(tǒng)的研究。因此，本研究基于全基因組重測序?qū)Σ煌笈Ｈ后w進行比對，同時結(jié)合群體基因組學(xué)分析構(gòu)建牦牛群體基因組特征圖譜，解析祁連牦牛與其他牦牛群體的差異意義重大。

1 材料與方法

1.1 樣本采集

本研究選擇九龍牦牛（Jiulong yak）、天祝白牦牛（Tianzhu White yak）、斯布牦牛（Sibu yak）、帕米爾牦牛（Pamir yak）、木里牦牛（Muli yak）、雪多牦牛（Xueduo yak）、祁連牦牛（Qilian yak）和肅南牦牛（Sunan yak）共8 個牦牛類群的46 個樣本進行測序分析，樣本具體信息詳見表1。本研究采集所有牦牛個體靜脈血10 mL，于含有EDTA的一次性負(fù)壓采血容器內(nèi)，放入-20℃冰箱備用。提取牦牛血液DNA 時使用EasyPure 血液基因組DNA 試劑盒（全式金，北京），使用Qubit 熒光定量儀檢測樣本DNA濃度，1%瓊脂凝膠電泳檢測樣本DNA 完整性，檢測合格的樣本在illuminaPE150 平臺上機測序。

表1 樣本基本信息

1.2 全基因組序列比對與遺傳變異檢測

因為在血液樣本DNA 提取、建庫和測序之后生成的數(shù)據(jù)會存在誤差，并且誤差的存在會對后續(xù)分析產(chǎn)生影響。因此，為了保證后續(xù)分析的準(zhǔn)確性，對上述步驟產(chǎn)生的數(shù)據(jù)誤差進行了質(zhì)量評估和過濾。使用Trimmomatic［26］軟件進行質(zhì)控，GATK 軟件進行SNP 鑒定，并用軟件自帶參數(shù)（Select Variants 和Variant Filtration）對SNP 進行過濾［27］。

1.3 群體基因組學(xué)分析

1.3.1 主成分分析和系統(tǒng)進化樹構(gòu)建

為了將個體按照主成分進行分類，了解研究群體的聚類和差異，采用主成分分析進行聚類，使用PLINK［28］軟件對常染色體SNP 數(shù)據(jù)進行分析。使用bcftools 提取vcf文件子集，將vcf 轉(zhuǎn)換為plink 格式，得到3 個以.map、.nosex、.ped 結(jié)尾的文件，基于.ped 生成bed 二進制文件，然后得到2 個以.bim、.bed 結(jié)尾的文件，之后進行PCA 分析，獲得2 個.eigenval、.eigenvec 結(jié)尾的文件，其中.eigenvec 記錄特征向量，用于坐標(biāo)軸的繪制。最后利用R 語言的ggplot2 包可視化繪制PCA 圖。為了解群體間個體的譜系關(guān)系，使用PLINK 軟件計算遺傳距離矩陣，采用Neighbor-Joining 法構(gòu)建系統(tǒng)進化樹，使用iTOL 在線軟件繪制系統(tǒng)進化樹。

1.3.2 雜合度和近交情況分析

為了評估牦牛的種質(zhì)資源現(xiàn)狀，使用mlRho［29］軟件計算牦牛個體的雜合度以及近交程度。

1.3.3 群體分化指數(shù)分析

為了進一步衡量群體間的分化程度，用VCFtools［30］軟件采用滑動窗口法進行Fst 分析，最后使用R 語言中的corrplot 安裝包對結(jié)果進行可視化。

1.3.4 連續(xù)純合子區(qū)域（ROH）分析

在動物進化過程中，ROH 會在動物基因組上呈現(xiàn)出不同的分布模式，因此動物基因組上ROH 的數(shù)目、長度和分布可以揭示動物種群的進化歷史等。本研究采用PLINK［28］軟件對ROH 進行檢測。

2 結(jié)果與分析

2.1 不同牦牛類群的聚類分化情況

主成分分析作為群體遺傳學(xué)中常用的分析方法，可以將許多相關(guān)性很高的變量轉(zhuǎn)化成彼此相互獨立或不相關(guān)的變量，進而通過結(jié)果了解群體的聚類情況。由圖1A 可知，本研究牦牛群體能夠得到明顯的區(qū)分，說明這幾個牦牛群體之間很少有雜交現(xiàn)象，并且群體內(nèi)樣本相對較聚集，群體間較分散。其中祁連牦牛與肅南牦牛和雪多牦牛2 個群體的遺傳距離較近，與帕米爾牦牛、木里牦牛、九龍牦牛和天祝白牦牛4 個牦牛群體的遺傳距離相對較遠(yuǎn)。為了解不同牦牛的譜系關(guān)系，本研究計算了個體間的遺傳距離矩陣，并采用鄰接法構(gòu)建系統(tǒng)進化樹。圖1B 顯示，所有牦牛群體單獨聚在一起，未出現(xiàn)牦牛個體混雜現(xiàn)象，這與前面的主成分分析結(jié)果相一致。其中斯布牦牛與帕米爾牦牛聚為一類，這也與兩個牦牛群體的地理位置相符合；同處在四川省的木里牦牛和九龍牦牛聚為一類；位于甘肅省境內(nèi)的天祝白牦牛和肅南牦牛聚為一類，然后與祁連牦牛聚為一大類。

圖1 不同牦牛群體的聚類分化情況

2.2 不同牦牛類群的雜合度和近交系數(shù)分析

雜合度和近交系數(shù)分析結(jié)果見表2，幾個牦牛群體的理論雜合度在0.354 4～0.360 0 之間。肅南牦牛的觀測雜合度最大，木里牦牛的觀測雜合度最小，分別為0.346 3和0.311 2。所有牦牛群體的觀測雜合度均小于理論雜合度。通過對牦牛群體的近交系數(shù)計算發(fā)現(xiàn)，8 個牦牛群體的近交系數(shù)在0.022 9～0.136 0 之間，其中木里牦牛最大，肅南牦牛最小。

表2 各牦牛群體的雜合度和近交系數(shù)

2.3 不同牦牛類群的群體分化指數(shù)

在前面的分析中，本研究使用的8 個牦牛群體可以得到明顯的區(qū)分，兩個分析結(jié)果相一致。但為了進一步對群體分化程度數(shù)字化，更加清晰地了解群體間的分化程度，基于SNP 數(shù)據(jù)計算了各牦牛群體間的群體分化指數(shù)。目前認(rèn)為遺傳分化指數(shù)≥0.04 的牦牛群體，可以拆分為新的遺傳品種資源，結(jié)果見圖2。由圖2 可知，九龍牦牛（JL，0.08）、木里牦牛（ML，0.08）、帕米爾牦牛（PMR，0.08）、祁連牦牛（QL，0.05）、斯布牦牛（SB，0.07）、肅南牦牛（SN，0.07）、天祝白牦牛（TZ，0.06）和雪多牦牛（XD，0.07）這8 個牦牛群體的分化指數(shù)均大于0.04，這與上述分析結(jié)果相一致，說明祁連牦牛可以作為一個單獨的遺傳資源，這為祁連牦牛的科學(xué)分類和發(fā)展提供了科學(xué)依據(jù)。

圖2 不同牦牛群體的遺傳分化指數(shù)

2.4 不同牦牛類群的連續(xù)純合子片段分析

在本研究的8 個牦牛群體中共檢測到ROH 15 704個，并發(fā)現(xiàn)ROH 的數(shù)量和長度呈正相關(guān)，介于0.91～1.00之間，結(jié)果見圖3A。其中帕米爾牦牛群體檢測到的ROH數(shù)量最大，肅南牦牛最少，分別為4 566 個和559 個。將ROH 的長度進行分類（＜1 Mb；1～2 Mb；＞2 Mb），ROH數(shù)目＜1 Mb 大小的占絕大部分，其次為1～2 Mb 大小的；8 個牦牛群體（九龍牦牛，天祝白牦牛，斯布牦牛，帕米爾牦牛，木里牦牛，雪多牦牛，祁連牦牛，肅南牦牛）ROH＞2 Mb 類別的數(shù)目分別為0、1、3、3、4、0、0 和0 個，結(jié)果見圖3B。檢測到ROH 的數(shù)量排序為：帕米爾牦牛＞木里牦牛＞天祝白牦牛＞祁連牦牛＞斯布牦牛＞九龍牦牛＞雪多牦牛＞肅南牦牛。

圖3 不同牦牛群體的ROH 分析

3 討論

種群群體遺傳結(jié)構(gòu)的確定對于家養(yǎng)動物與野生動物來說都很重要，這關(guān)系到他們的育種和種群保護［31］。有大量的研究人員使用微衛(wèi)星標(biāo)記研究了牦牛群體的遺傳結(jié)構(gòu)［32-34］。本研究通過群體基因組學(xué)構(gòu)建牦牛群體特征圖譜發(fā)現(xiàn)：所有牦牛群體能夠得到區(qū)分，群體之間較為分散，沒有出現(xiàn)混雜情況。系統(tǒng)進化樹與PCA 分析結(jié)果相一致。通過這兩個分析可以將牦牛群體進行初步的區(qū)分。遺傳多樣性也稱為基因多樣性，是某一生物品種或群體內(nèi)個體之間遺傳變異的總和，是評估種質(zhì)資源現(xiàn)狀的重要指標(biāo)［35］，群體的遺傳多樣性越高，對環(huán)境的適應(yīng)能力就越強［36］。同樣雜合度和近交系數(shù)也可用來評估群體的遺傳多樣性，研究中的理論雜合度是通過計算所得，而觀測雜合度則是兩個樣本等位基因不相同的概率。近交系數(shù)是指基因純化的程度。相關(guān)研究發(fā)現(xiàn)，群體的雜合度與生長速率和繁殖力等呈正相關(guān)［37］。本研究結(jié)果表明，8 個牦牛群體的觀測雜合度略低于期望雜合度，但相差并不大，說明雖然群體內(nèi)受到了一定程度的選擇壓力，但尚未發(fā)生遺傳漂變等。其中肅南牦牛的觀測雜合度最高（0.346 3），但近交系數(shù)最低（0.022 9）；而木里牦牛的觀測雜合度最低（0.311 2），近交系數(shù)最高（0.136 0）；各牦牛群體近交系數(shù)順序為：木里牦牛＞帕米爾牦牛＞祁連牦牛＞斯布牦牛＞九龍牦牛＞天祝白牦牛＞雪多牦牛＞肅南牦牛，說明木里牦牛群體的近交程度最大，肅南牦牛群體的近交程度最小。本研究中8 個牦牛群體的雜合度與王佟等［23］的研究結(jié)果基本相符，但低于李鐸等［33］的研究結(jié)果。推測牦牛群體的遺傳多樣性是波動的，也可能與研究所用的位點數(shù)目、樣本數(shù)量及采樣地區(qū)等多因素相關(guān)，這仍需進一步的深入研究。

群體分化指數(shù)的范圍在0～1 之間，數(shù)值越大，說明群體間的分化程度就越大。本研究中祁連牦牛的群體分化指數(shù)為0.05，與其他牦牛群體間存在明顯的遺傳分歧，可作為潛在的牦牛群體遺傳資源，為今后大力發(fā)展該牦牛群體的資源挖掘和開發(fā)保護工作提供理論基礎(chǔ)。相關(guān)研究表明，ROH 在基因組中的大小和分布受到近親交配和人工選擇等因素的影響［38］，目前ROH 分析常用來計算評估近交系數(shù)，并且這種方法具有較高的準(zhǔn)確度［39］。本研究對8 個牦牛群體的ROH 分析結(jié)果顯示，帕米爾牦牛的ROH 數(shù)目最多，可能是因為帕米爾牦牛經(jīng)過了近交選育，這與其群體的近交系數(shù)較高相一致。畜禽種質(zhì)資源是種業(yè)振興的基礎(chǔ)，對牦牛群體的科學(xué)保護和優(yōu)質(zhì)資源的開發(fā)利用，可以有效地保護種群內(nèi)的基因多樣性，促進品種的改良和創(chuàng)新?？傊?，本研究為牦牛群體的資源保護提供了科學(xué)依據(jù)。

4 結(jié)論

綜上，8 個牦牛群體之間沒有過多的基因交流，其中祁連牦牛是在長期的進化發(fā)育過程中產(chǎn)生的一個獨特類群。本研究從全基因組水平揭示了祁連牦牛的系譜及遺傳結(jié)構(gòu)，為后續(xù)祁連牦牛遺傳資源的保護及開發(fā)利用提供了理論依據(jù)。