紀會 ，官久強 ，王會 ，周建旭 ，阿農呷 ，何宗偉 ，樊珍詳 ，邱龍康 ，曹詩曉 ，安添午 ，柏琴 ，鐘金城 *，羅曉林 *

（1. 四川省草原科學研究院，四川 成都 611731；2. 西南民族大學，四川 成都 610041；3. 小金縣科學技術和農業(yè)畜牧水務局，四川 阿壩州 624200；4.甘孜藏族自治州畜牧站，四川 甘孜州 626099；5. 稻城縣農牧農村和科技局畜牧站，四川 甘孜州 627750；6. 新龍縣農牧農村和科技局畜牧站，四川 甘孜州 626800）

牦牛（Bos grunniens）是分布于以青藏高原為中心，及其毗鄰的高山、亞高山地區(qū)的特有珍稀牛種［1］。青藏高原地域環(huán)境復雜，不同區(qū)域的牦牛在生態(tài)適應性方面存在明顯差異，表現為體型外貌和生產性能的差異，形成了各具特色的地域類型和亞型［2］。家牦牛約在7300 年前由野牦牛馴化而來［3］，中國有18 個地方品種（遺傳資源）和2 個培育品種，是世界上擁有牦牛數量和品種類群最多的國家。然而由于長期過度選育和純繁，導致現有家養(yǎng)牦牛出現品種流失的現象，尤其是與其抗病性、肉質等相關性狀的退化已相當嚴重，極大地影響了青藏高原及周邊地區(qū)農牧業(yè)生產的穩(wěn)定性［4］。因此，加強地方牦牛群體遺傳多樣性研究，挖掘新種質資源，對于遺傳資源保護、新品種培育和牦牛生產性能提高等都具有積極的理論和現實意義。

亞丁牦牛和拉日馬牦牛是肉乳兼用型優(yōu)良地方牦牛資源。亞丁牦牛分布于四川省甘孜藏族自治州稻城縣，經產母牛 153 d 平均擠奶量 434 kg，產奶量居中國牦牛之最［5?6］，2005 年亞丁牦牛存欄量約 7.5 萬頭，2015 年降至4.6 萬頭，由于飼草料缺乏，背山封閉，近親交配嚴重，繁殖性能差，其群體有退化趨勢。拉日馬牦牛中心產區(qū)位于四川省甘孜藏族自治州新龍縣拉日馬鎮(zhèn)，2016 年存欄量約8.7 萬頭，具有凈肉率高的特性，公牛凈肉率達45.7%［7］。然而，關于兩個群體分子遺傳背景的研究尚未開展。

限制性內切位點相關的DNA 標記（restriction site-association DNA，RAD）簡化基因組測序，即利用限制性酶切法獲得目的基因組部分片段進行高通量測序，獲得可代表目的基因組信息的高密度SNPs 分子標記，具有簡便、有效及低成本等特點，已廣泛用于群體遺傳多樣性研究，如探究孟加拉虎出現白色條紋的遺傳機理［8］，對鹿苑雞不同保種群保種現狀的評估［9］及中國家兔地方品種的遺傳地位研究［10］。本研究利用RAD 簡化基因組測序對9個地方牦牛群體進行全基因組限制性酶切檢測，以期在DNA 分子水平分析亞丁牦牛和拉日馬牦牛的遺傳多樣性，探明各群體間的進化關系，為今后進一步開展兩種牦牛遺傳資源的保護、遺傳改良和育種提供理論依據。

1 材料與方法

1.1 試驗動物和方法

以亞丁牦牛、拉日馬牦牛、九龍牦牛、麥洼牦牛、金川牦牛、昌臺牦牛、中甸牦牛、玉樹牦牛、類烏齊牦牛為試驗動物。每個群體隨機選擇9～11 頭，采集類烏齊牦牛的組織樣本以及其他8 個牦牛群體的血液樣本（表1）。具體操作：用一次性注射器從牦牛頸外靜脈采集5 mL 血液于EDTA 抗凝管，放入冰盒中帶回，放置?80 ℃冰箱中保存；對類烏齊牦牛進行屠宰后采集組織樣品，經DEPC（diethyl pyrocarbonate，焦碳酸二乙酯）沖洗后，錫箔紙包裝迅速置于液氮保存，帶回實驗室備用。牦牛群體地理分布情況見圖1。

表1 牦牛樣本來源及樣本量Table 1 Source and number of yak sample

圖1 牦牛樣本的地理分布Fig.1 Geographical distribution of yak sample

1.2 基因組電子酶切評估

建庫測序前，為檢測酶切效率和對所有標記進行分析，按照PstI 酶（CTGCA，G）切位點對參考基因組（GCF_000298355.1）進行模擬酶切（電子酶切），統(tǒng)計獲得限制性酶切片段的數量及覆蓋率。共獲得2959206 條限制性酶切片段（restriction fragment，RF），基因組覆蓋率為41.89%，其中有效限制性酶切片段（effective RF：長度大于300 bp）共1558672 條。有效限制性酶切片段在基因組中分布均勻，PstI 酶對牦牛參考基因組的酶切效果較好，可用于建庫測序（圖2）。

圖2 基因組中有效限制性酶切片段分布圖Fig.2 Map of effective restriction fragments（effective-RF）in the genome

1.3 全基因組限制性酶切

采用血液基因組DNA 試劑盒（北京天根）提取血樣DNA，CTAB（cetyltrimethylammonium bromide，十六烷基三甲基溴化銨）法提取組織樣DNA。采用PstI 酶對全基因組進行酶切，88 個樣本共構建88 個簡化基因組測序（restriction site-association DNA，RAD）文庫，在HiSeq3000 測序平臺上用150 bp 雙端序列進行RAD 測序。

1.4 數據處理

測序數據的質量控制和原始數據過濾：通過分析堿基的質量值和組成評估數據的質量。對原始序列（raw reads）中低質量的reads 進行過濾，獲得高質量的凈數據（clean reads），用于后續(xù)的信息分析。過濾條件：去除帶接頭的reads；去除含未知堿基比例大于10%的reads；去除低質量reads（質量值Q≤20 的堿基數占整條reads 的50% 以上）［11］。

單核苷酸多態(tài)性（single nucleotide polymorphisms，SNP）的識別與過濾：使用比對軟件 bwa［12］（0.7.12）采用mem 算法將過濾后的clean reads 比對到參考基因組上，比對參數為-k32-M；比對結果使用picard（1.129）軟件標記 duplicates reads（PCR 過量擴增所形成的 reads）。使用軟件 GATK［13］檢測群體 SNP/indel，參數：-filter“QD＜2.0”，得到可信度較高的SNP/indel 數據集。Indel（insertion-deletion）全稱為插入或缺失突變，指在DNA 或RNA水平上，發(fā)生于個體間的單核苷酸突變位點。為提高分析準確性，對原始SNP 進行過濾，過濾條件為：缺失率（missing rate）不超過20%，第二等位基因頻率（minor allele frequency，MAF）不小于0。

1.5 群體遺傳多樣性分析

雜合度統(tǒng)計：雜合度指同源染色體上的一個位點存在不同等位基因的狀態(tài)，用于衡量群體中的遺傳變異程度［14］。應用 Stacks populations v 2.41 軟件［15］，其參數設置如下：最小等位基因頻率（MAF）＞0.05；標簽最低覆蓋倍數＞3×；在每個等位基因位點處，每個品系至少75%的個體存在序列信息，至少3 個品系存在有效信息［16］。雜合度取值在0～1 之間，雜合度越高，表明群體內遺傳多樣性系數越高。

連鎖不平衡（linkage disequilibrium，LD）分析：連鎖不平衡指群體內不同座位等位基因之間的關聯。采用PopLDdecay 3.31［17］軟件，通過計算一定序列范圍內兩兩位點間的LD 系數（r2），來估算LD 衰減的趨勢。r2代表兩位點的相關性，它在0＜r2＜1 的區(qū)間內變化，如果r2=0，表示兩位點發(fā)生重組的可能性是100%，如果r2=1，則表示位點完全連鎖。通常位點間的關聯性會隨著基因組上位點間距離的增加而不斷下降，稱為LD 衰減速度。使用LD 衰減距離來描述LD 衰減速度的快慢，衰減距離越長，其衰減速度越慢，群內等位基因間連鎖程度越強。

1.6 群體間遺傳關系分析

成對遺傳分化值（Pairwise Fst）［18］：分析不同品種間的遺傳距離。應用Stacks populations v 2.41 軟件，參數設置與雜合度統(tǒng)計一致。根據每個等位基因位點SNP 信息，進行雙向選擇家系內Pairwise Fst 統(tǒng)計，并對各Fst 進行Fisher 精確檢驗，使用P 值對Fst 值進行校正（Corrected Fst），進而根據高斯核函數，計算Smoothed pairwise Fst，利用 40 kb-Slid Windows 對全基因組上的 Fst 進行平均值統(tǒng)計。采用 Wright［19］理論，如群體 Fst 值為 0.00～0.05，說明群體間不存在分化；Fst 值為0.05～0.15，說明群體間處于中度分化；Fst 值為0.15～0.25，則說明群體間分化程度較大；Fst 值大于0.25，說明群體間存在極大分化。

主成分分析（principal component analysis，PCA）：使用軟件GCTA［20］從中獲得各個樣本的主成分值。根據各個樣本在第一主成分（PC1）、第二主成分（PC2）、第三主成分（PC3）3 個綜合指標中的數值大小，使用R 語言繪制PCA 散點圖。散點圖中每一個位點代表一個樣本，圖中兩個樣本距離越遠，則說明兩個樣本遺傳背景差異越大［21］，遺傳背景相似的個體則會在圖中聚為一類。

系統(tǒng)進化樹構建：基于群體SNP 數據運用MEGA 4.0［22］軟件計算樣品間距離，得到個體間的序列差異度矩陣，然后使用IQ-Tree 軟件，最大似然法構建系統(tǒng)進化樹（maximum likelihood tree，ML-tree）。

1.7 群體結構分析

應用Structure［23］軟件進行群體結構分析，用于推算測序樣本分屬幾個群體，群體間基因交流程度以及每個樣本的混血程度?；诙辔稽c基因型將遺傳相似的個體聚類。先預設88 個樣本由若干亞群K（K=1～n）構成，根據每個SNP 位點的基因型采用隱馬科夫?蒙特卡羅鏈模型算法（markov-chain monte carlo，MCMC）［24］，在每個群體均最大程度符合哈代溫伯格平衡檢驗前提下，輸出88 個樣本被分為K個亞群時的最佳分類狀態(tài)，同時每一個K值均會對應產生一個似然值（likelihood），最后對似然值取自然對數后輸出Inlikelihood，Inlikelihood 越大，說明K值越接近于真實情況，一般Inlikelihood 值先隨著K值增大而不斷升高，之后進入平臺期，最佳K值即為ΔK值最大時。

2 結果與分析

2.1 牦?；蚪M的SNP

對RAD 簡化基因組文庫進行Illumina Hiseq 3000 雙端高通量測序，88 個牦牛樣本測序獲得4.2 G reads 原始序列數據，經過濾后有效數據量為3.9 G。各群體原始數據有效率在89.82%～95.54%，測序質量數值大于20、30 的堿基占總體堿基的百分比為Q20≥96.93%、Q30≥91.47%，GC 含量在48.35%～49.07%，說明堿基組成平衡，測序質量較好，數據可靠。將過濾后的高質量reads 與參考基因組比對，比對率均在98.00%以上。平均測序深度6.8×，88 個樣本基因組覆蓋率在30.99%～38.70%，與電子酶切參考基因組預測的覆蓋率相近，預期率達92.48%，說明RAD 簡化測序預測牦?；蚪M效果較好（表2）。

表2 測序數據質量匯總情況Table 2 Summary of sequencing data quality

所有樣本共檢測到20413266 個SNP 位點，亞丁牦牛SNP 數量最多，為7596791，昌臺牦牛SNP 數量最少，為5472579，拉日馬牦牛SNP 數量為6163048，不同群體SNP 信息見圖3。另檢測到插入變異位點（insertion）有 362758 個 ，缺 失 變 異 位 點（deletion）487510 個，插入和缺失位點長度集中在1～10 bp，其中1～3 bp 最多（圖4）。

圖 3 9 個群體 SNP 數量Fig.3 Number of SNP in 9 populations

圖4 Indels 長度統(tǒng)計Fig.4 Indels length statistics

2.2 群體遺傳多樣性

2.2.1 牦牛群體的雜合度 9 個牦牛群體的觀測雜合度差異較小，在0.2172～0.2440 之間，均低于期望雜合度（表3）。九龍牦牛雜合度最低為0.2172，遺傳多樣性最貧乏。亞丁牦牛的觀測雜合度較低為0.2186，遺傳多樣性貧乏；拉日馬牦牛觀測雜合度為0.2233，低于9 個群體平均觀測雜合度0.2276，遺傳多樣性偏低。

表3 9 個牦牛群體的雜合度Table 3 Heterozygosity of 9 yak populations

2.2.2 牦牛群體的連鎖不平衡 9 個群體SNP 距離在0～50 kb 范圍內的LD 衰減速度較快，LD 系數約從 0.5 降至 0.2 以下，SNP 距離在 50～300 kb 范圍內LD 水平較低在 0.1～0.2 之間，顯示在 0～50 kb 范圍內各等位基因連鎖不平衡程度較高。當LD 系數r2=0.2 時，各群體對應LD 衰減距離基本相等，說明各群體LD 衰減速度差異不大（圖5）。

圖5 LD 隨距離增長的衰減大小Fig.5 Attenuation of LD with distance

2.3 牦牛群體間的遺傳關系

通過Pairwise Fst 統(tǒng)計分析不同地方牦牛群體間的遺傳距離。亞丁牦牛、九龍牦牛、中甸牦牛以及類烏齊牦牛平均分化指數高于0.05（表4），遺傳分化程度較大；拉日馬牦牛、麥洼牦牛、金川牦牛、昌臺牦牛以及玉樹牦牛平均分化指數低于0.05，遺傳分化程度較低。亞丁牦牛平均分化指數最高為0.0653，與其他8 個群體存在中等程度（0.05＜Fst＜0.15）的遺傳分化；拉日馬牦牛的平均分化指數最低為0.0443，與玉樹牦牛、麥洼牦牛群體間的遺傳分化程度最小。

表4 群體間Pairwise FstTable 4 Pairwise Fst between population

基于個體基因組SNP 差異程度，按主成分將不同個體聚類成不同亞群（圖6）。亞丁牦牛群體明顯與其他地方牦牛群體分開，而且分布距離較遠，與Fst 統(tǒng)計結果一致，遺傳分化程度較高。在PC2 負值方向，有2個拉日馬牦牛樣本（LRM-4、LRM-7）與九龍牦牛群體聚集，在 PC2 正值方向，金川牦牛 JC-3、JC-4 樣本、麥洼牦牛（MW-4、MW-5、MW-6、MW-8）樣本分別獨立成群，遺傳變異現象明顯，其余個體以及拉日馬牦牛、昌臺牦牛、玉樹牦牛、類烏齊牦牛群體間出現嚴重重疊，親緣關系較近。

圖6 9 個牦牛群體的主成分分析Fig.6 Principal components analysis of 9 yak populations

為了解9 個地方牦牛群體間的親緣關系，通過最大似然法構建系統(tǒng)進化樹（圖7），結果顯示亞丁牦牛與中甸牦牛親緣關系最近，拉日馬牦牛雖有3 個樣本（LRM-4、LRM-7、LRM-6）與九龍牦牛聚為一支，其余7 個樣本聚集明顯，先與麥洼牦牛聚為一支，說明拉日馬牦牛與九龍牦牛存在基因交流，與麥洼牦牛親緣關系近。類烏齊牦牛與玉樹牦牛分群明顯，親緣關系近。金川牦牛和昌臺牦牛聚為一支，無明顯分群，存在嚴重基因交流現象。

圖7 最大似然法構建系統(tǒng)進化樹Fig.7 Phylogenetic tree construction by maximum likelihood method

2.4 牦牛群體的遺傳結構

對牦牛群體進行群體結構分析，確定最佳分群情況，揭示各牦牛群體的血統(tǒng)構成，以及每個個體的血統(tǒng)潛在來自哪些群體以及對應比例（圖8），可追溯其品種的形成歷史及親緣關系。結果顯示從K=2～9，亞丁牦牛始終分為一個獨立群體，分化程度最高，與Fst 統(tǒng)計、PCA 分析結果一致。Inlikelihood 計算結果顯示（圖 9），當K=3 時，ΔK值最高，即 9 個地方牦牛群體88 個個體可分為3 個群體：亞丁牦牛、九龍牦牛、其他7 個地方群體（拉日馬牦牛、昌臺牦牛、金川牦牛、類烏齊牦牛、麥洼牦牛、玉樹牦牛、中甸牦牛）為一個群體，最接近牦牛群體的理論分群情況。亞丁牦牛、九龍牦牛血統(tǒng)純正，僅部分個體出現極少量混雜現象，其他7 個群體血統(tǒng)組成相似。中甸牦牛每個個體均含少量（≤25%）亞丁牦牛和九龍牦牛血統(tǒng)，存在基因交流現象，與亞丁牦牛、九龍牦牛親緣關系近，與系統(tǒng)進化樹分析結果一致。K=3～9，拉日馬牦牛血統(tǒng)中含有一定比例九龍牦牛血統(tǒng)，K=6 時起，麥洼牦牛部分個體（MW-4、MW-5、MW-6、MW-8）被歸為一個獨立群體，與PCA 分析結果一致，拉日馬牦牛含有一定比例麥洼牦牛血統(tǒng)，說明拉日馬牦牛與九龍牦牛、麥洼牦牛均存在基因交流，與麥洼牦牛血統(tǒng)組成較相似，親緣關系較近，與系統(tǒng)進化樹分析結果一致。

圖8 牦牛群體結構堆疊分布Fig.8 Groups structure clustering distribution in yak

圖9 運用Inlikeliood 預估最佳K 值的4 個步驟Fig.9 Four steps to estimate the best K value with inlikelihood

3 討論

關于牦牛群體遺傳多樣性已開展大量的研究，柴志欣等［25］利用RAPD 分子標記分析西藏地區(qū)11 個牦牛類群的遺傳多樣性并進行分類；Zhang 等［26］利用16 種微衛(wèi)星標記分析中國9 個牦牛品種的種群內遺傳變異水平；郭松長等［27］利用線粒體DNAD-loop 區(qū)部分序列對我國10 個家牦牛品種（類群）進行遺傳多樣性及分類研究。亞丁牦牛和拉日馬牦牛是肉乳兼用型優(yōu)良地方牦牛資源，探究其分子遺傳背景，為開展亞丁牦牛和拉日馬牦牛遺傳改良和育種提供理論依據。

本研究利用RAD 測序技術在全基因組水平探究亞丁牦牛和拉日馬牦牛的遺傳多樣性及遺傳結構。群內遺傳多樣性分析結果顯示，9 個牦牛群體的觀測雜合度為0.2172～0.2440，均低于期望雜合度，說明各地方群體內存在一定程度近交率［28］。其中，麥洼牦牛、九龍牦牛、昌臺牦牛、類烏齊牦牛和中甸牦牛的觀測雜合度均低于運用微衛(wèi)星標記所測群體觀測雜合度（0.5650～0.9599）［29?32］，微衛(wèi)星標記所測的雜合度較高與其標記位點多態(tài)性高、分布集中有關，RAD 簡化基因組測序較全面地反映整個基因組多態(tài)位點情況［9］。亞丁牦牛和拉日馬牦牛的觀測雜合度分別為0.2186 和0.2233。亞丁牦牛遺傳多樣性貧乏，可能因為亞丁地區(qū)背山封閉，近親交配嚴重，導致繁殖性能降低，從而引起群體內個體雜合性降低［5］。拉日馬牦牛遺傳多樣性偏低，推測其群內存在一定程度近交率。9 個群體間LD 衰減速度差異較小，可能與各群體間遺傳多樣性差異小有關［33］。

本研究中亞丁牦牛的平均遺傳分化指數Fst 最高為0.0653，與其他8 個地方牦牛群體存在中度分化，群體結構分析顯示從K=2～9 亞丁牦牛始終分為一個獨立群體，血統(tǒng)構成純正，且Inlikelihood 計算結果顯示亞丁牦牛理論上應歸為一個獨立群體。生物的進化以自然環(huán)境為選擇主線，不斷修正和改變基因的頻率，從而引導生物的進化方向［34］。據史料記載，亞丁地區(qū)飼養(yǎng)牦牛歷史悠久，公元619?786 年，藏王將吐蕃國牦牛遷徙至東義片區(qū)即現在的亞丁地區(qū)飼養(yǎng)，該區(qū)域受赤土河、東義河切割，山高谷深，立體氣候明顯，亞丁牦牛在對環(huán)境的適應過程中同時受到人工選擇的作用，亞丁地區(qū)和中甸地區(qū)歷史上曾有相互交換種牦牛和在交界地混牧的習慣［6］，因此中甸牦牛血統(tǒng)中含有少量亞丁牦牛血統(tǒng)，具有親緣關系，經過長期的自群繁育，亞丁牦牛群體血統(tǒng)純正，與其他地方群體間產生了遺傳變異。總之，亞丁牦牛雖體遺傳多樣性貧乏，但血統(tǒng)純正，其遺傳背景與其他8 個地方群體差異明顯，肉奶性能良好，其產奶量居中國牦牛之最，建議亞丁牦牛單獨列為一個牦牛品種進行開發(fā)利用和保護。

拉日馬牦牛的平均遺傳分化指數Fst 最低為0.0443，群體遺傳背景較復雜，與九龍牦牛和麥洼牦牛均存在基因交流，與麥洼牦牛親緣關系最近。拉日馬牦牛與九龍牦牛中心產區(qū)相距較遠，存在基因交流的原因可能是九龍牦牛向拉日馬鎮(zhèn)發(fā)生過遷徙，九龍牦牛屬世界上最好的種牛，存在人工引種的可能。新龍縣志辦記載，在1925 年前后，瓦西洛爾龍作為當時新龍縣拉日馬麥洼部落頭人（稱“麥本”），因部落內部糾紛，帶著牲畜7000 多頭（主要是牦牛）背井離鄉(xiāng)前往現阿壩州藏族羌族自治州紅原縣境內麥洼區(qū)定居。據文獻資料記載及牦?？蒲泄ぷ髡叨啻握{查，麥洼牦牛來自甘孜藏族自治州北部色達、德格、爐霍、新龍等縣［6］。本研究在全基因組水平上進一步證明拉日馬牦牛與麥洼牦牛存在基因交流，屬麥洼牦牛的祖先群之一，拉日馬牦牛較麥洼牦牛在肉奶性能方面存在一定的優(yōu)越性，其凈肉率和產奶量均高于麥洼牦牛［8?9］。

4 結論

本研究以亞丁牦牛和拉日馬牦牛及其他7 個地方牦牛類群為研究對象進行SNP 檢測和群體遺傳關系分析。亞丁牦牛遺傳多樣性貧乏，血統(tǒng)純正，與其他8 個地方牦牛群體存在中度分化，可列為一個獨立遺傳資源進一步研究和保護。拉日馬牦牛遺傳多樣性偏低，屬麥洼牦牛的祖先群之一，其群體遺傳結構較復雜，需結合其他研究方法、增加樣本量進一步開展研究工作。