基于高密度SNP 芯片評估“科爾沁肉牛”遺傳背景

林雨濃，馬萬欣，包玲玲，何曙光，石順利，張秋生，趙澈勒格日，高會江，李俊雅*，王澤昭*

（1.中國農業科學院北京畜牧獸醫研究所，北京 100193；2.內蒙古自治區通遼市畜牧業發展中心，內蒙古通遼 028000）

肉牛產業是內蒙古自治區經濟發展的支柱產業，是農牧民增收的重要途徑。“科爾沁牛”是以西門塔爾牛為父本，蒙古牛及三河牛和蒙古牛的雜種牛為母本而培育成的乳肉兼用品種。作為通遼地區地方畜牧業的重要資源，具有適應性強，宜牧，耐粗飼，耐寒，抗病力強等特點。然而，隨著肉牛市場需求的增加和養殖業的發展，“科爾沁牛”在生產性能和經濟效益方面遇到了瓶頸。為了提升群體肉用性能，當地農業部門利用北美肉用型西門塔爾牛為父本，本地“科爾沁牛”為母本進行雜交改良，歷經30 年，經過“雜交改良、橫交固定和選育提高”三個階段，培育出肉用型“科爾沁肉牛”群體。

近年來，群體遺傳背景分析成為評估和研究動物群體遺傳結構的重要工具。隨著SNP 芯片和測序技術的發展，有助于人們更全面地了解不同群體之間的遺傳多樣性和遺傳聯系。Eva M.Strucken 等[1]采用777K 芯片對印度本地的15 個牛種進行遺傳多樣性與群體結構研究，發現印度瘤牛比印度本地普通牛有更高的遺傳多樣性，他們推斷本地普通牛較小的有效群體導致了群體基因頻率受到遺傳漂變的影響更大，而低有效群體數量增加了近交衰退發生的概率，使其在選擇的過程中遺傳多樣性降低速度快于印度瘤牛。在動物育種過程中，保持一定的遺傳多樣性有利于提高群體的環境適應性[2]。Alejandro Amaya 等[3]對約2.7 萬頭哥倫比亞西門塔爾牛（1975—2017年）的數據進行了分析，發現該種群的近交系數由5.06%下降到了2.25%，認為是得益于歐洲和美系西門塔爾牛之間的基因交流，才使該種群的遺傳多樣性得以提升。

我國疆域遼闊、資源豐富，但目前由于經濟效益等問題致使國外種源占有率居高不下。物種作為戰略性資源、是農業的基石，但是目前關于國內地方品種的研究存在著許多空白和滯后。馬均等[4]研究發現，我國自主培育的專門化肉牛新品種—華西牛與蒙古牛、三河牛和其他家系西門塔爾牛出現不同的比例的群體結構組成，形成了自己獨特的遺傳結構特征。劉正喜等[5]在2022 年使用全基因組重測序數據對渤海黑牛的群體結構、遺傳多樣性、親緣關系分析后發現，渤海黑牛與郟縣紅牛和魯西牛的關系密切，該研究為逐年減少的渤海黑牛群體的種質資源保護和品種發展提供了關系的理論基礎。

本研究目的是分析“科爾沁肉牛”、華西牛和美系西門塔爾牛這三個群體之間的遺傳關聯，評估“科爾沁肉牛”改良育種在遺傳層面的影響，旨在為“科爾沁肉牛”新品種培育提供遺傳背景支持。

1 材料與方法

1.1 試驗群體及SNP 分型

“科爾沁肉牛”（Kerqin beefcattle，“KBC”），來自內蒙古科爾沁牛業股份有限公司的437 頭核心群牛只；美系西門塔爾牛（American Simmental cattle，ASC）來自于通遼京緣種牛繁育有限責任公司的25 頭種牛；華西牛（HuaXi cattle，HXC），來自內蒙古奧科斯牧業有限公司的55 頭核心群牛只。

每個樣本使用含有EDTA 抗凝劑的真空采血管采血5mL 送樣檢測，并取2mL 血樣進行DNA提取。使用Illumina BovineHD BeadChip 芯片進行變異檢測，該芯片含有覆蓋牛基因組的777962個標記位點。

1.2 數據質控

使用Genome Studio 軟件進行基因分型。通過PLINK v1.9[6]軟件提取出常染色體后進行數據質量控制（Quality control，QC）。QC 順序如下：①刪除等位基因檢出率（Call rate）小于90%的位點；②刪除最小等位基因頻率（Minor allele frequency，MAF）小于0.01 的位點；③刪除哈代溫伯格平衡檢驗P 值小于1×10 箒6 的SNP 位點；④刪除SNP 缺失率大于10%的樣本個體。通過QC 后數據用于后續分析。單獨使用PLINK 中的--indep-pairwise 參數對所有個體的連鎖位點進行過濾，剩余位點進行群體結構分析。過濾后的個體數見表1。

表1 試驗樣本數量

1.3 遺傳多樣性分析

對QC 后數據統計群體期望雜合度（Expected heterozygosity，He）、群體觀察雜合度（Observed heterozygosity，Ho）、MAF 分布情況。 MAF 分布情況是按照不同MAF 大小將SNP 標記分為5 組，分別是（MAFSNP<0.1）、（0.1≤MAFSNP＜0.2）、（0.2≤MAFSNP＜0.3）、（0.3≤MAFSNP＜0.4）和（0.4≤MAFSNP），統計3 個群體SNP 標記的分布情況。

1.4 群體分析

使用GCAT 軟件進行了群體主成分（Principal component analysis，PCA）計算，輸出數據導入到R 中，使用ggplot2 包[7]進行可視化。隨后使用鄰接法（Neighbor-Joining，NJ）構建進化樹[8]：從成對FST 值的距離矩陣生成鄰居網絡樹，并在Splitstree 4.14.5[9]中實現建樹。使用Admixture v.1.3 軟件的默認設置進行群體結構分析。設置了10 倍交叉驗證選擇最合適的K 值[10]。使用PopLDdecay 軟件[11]計算了所有數據的LD 衰減距離，并繪制了LD 衰減圖。

2 結果

2.1 群體遺傳多樣性

由表2 可知，經過質控后的“KBC”群體具有最高的平均MAF 值（0.28），ASC 群體的平均MAF 值最低（0.23）。平均He 的范圍取值從0.34（HXC）到0.37（“KBC”）。平 均Ho 最高為“KBC”（0.37），最低的是ASC（0.36），HXC和ASC 的平均Ho 略大于He，“KBC”的Ho 和He 幾乎相等。

表2 群體遺傳多樣性參數

由圖1 可知，“KBC”群體被過濾掉了最多的遺傳標記。“KBC”群體過濾掉了81778 個MAF 小于0.01 的標記。QC 后，ASC 群體擁有最多的低頻標記（MAF＜0.1）。隨著MAF 的增大，HXC 和“KBC”群體內中高MAF 的SNPs 多于ASC 群體。在三個群體中，“KBC”具有最多的MAF≥0.3 的SNP 個數，具體數量為（413400），其次是HXC（363989），最少的是ASC（344174）。

圖1 不同群體的MAF 分布

2.2 PCA 和進化樹

由圖2（a）可知，三個群體主要聚集成為兩個簇，將“KBC”群體和HXC 群體分開。此外，“KBC”群體與ASC 群體中部分個體出現明顯的聚集和重疊。根據群體在圖中的空間位置，可以判斷出“KBC”和HXC 的遺傳距離較遠，和ASC 之間遺傳距離較近。

圖2 三個牛群體的主成分分析結果和鄰居連接系統發育樹

為了觀察三個群體在進化距離上的相對位置，本研究構建了環形進化樹（圖2b）。“KBC”群體和HXC 群體分別位于進化樹的起點位置和末端位置，表明這兩個群體的進化距離相對較遠。ASC 群體和HXC 群體處于兩個比較靠近的分支，暗示這兩個群體經歷了不同的進化路徑，但是其分化路徑相對較近。ASC 群體處于“KBC”群體的一個分支下，表明這兩個群體可能具有共同的祖先，使得在基因層面上具備一定的相似性。

2.3 種群結構分析

通過假設群體結構數的先驗值K，推斷群體內亞群個數。通過對比不同祖先群體遺傳成分比例觀察“KBC”、ASC 和HXC 群體間遺傳相似性。由圖3 可知，在K=2 時，“KBC”的祖先比例和HXC 相似，少數個體祖先比例與ASC 相似；當K=3 時，“KBC”群體和ASC 群體部分個體祖先比例相似，此時HXC 群體出現與“KBC”群體和ASC 群體截然不同的祖先成分比例；隨著K 值不斷增大，“KBC”群體和HXC群體差異的趨勢逐漸增大；當K=5 時，“KBC”出現了更為復雜的遺傳比例；K=2～4 時，ASC群體的祖先組成比例變化不明顯，“KBC”的祖先比例變化復雜。

圖3 當K=2～5 時，基于LD 過濾的三種牛種群的SNPs 的結構分析

2.4 連鎖不平衡

如圖4 所示，三個群體的LD-decay 曲線中“KBC”群體衰減速度最快，ASC 的LD 衰減速度較慢，HXC 群體的LD 衰減速度位于兩者之間。三個群體的遠端SNP 之間的連鎖水平出現差異，在≥200Kb 之后并未衰減到相近的r2范圍。隨著距離逐漸增大，三個群體的LD 衰減曲線逐漸平緩，r2結果始終保持是r2ASC≥r2HXC≥r2KBC。

圖4 三個群體的連鎖不平衡衰減圖

3 討論

3.1 遺傳多樣性參數

遺傳多樣性是生物多樣性的核心維度之一，反映了群體內部基因的多樣性水平。它在維持種群的健康、適應性和進化潛力方面起著關鍵作用[12,13]。目前，通過基因芯片和測序技術輔助分析種群遺傳進展成為熱點研究內容[14]，通過評估種群遺傳多樣性對于制定未來的育種目標至關重要[2]。本研究對“KBC”群體、HXC 群體和ASC群體的基因型數據進行分析，通過對群體遺傳多樣性的統計，以期發現“KBC”在育種過程中的遺傳變化。研究發現，“KBC”的MAF、He 和Ho 分別為0.28、0.37 和0.37，在三個群體中均為最高，并且“KBC”群體He 和Ho 的一致性在三個群體中為最高。這一結果反映出“KBC”群體內的遺傳多樣性比ASC 和HXC 群體較高，同時暗示著近親交配的水平較低[15]。種群遺傳多樣性越高，意味著進化潛力越大，對環境變化的反應能力越強[16]。此外，人工選擇也會極大的影響群體的多樣性[17]。“KBC”是以“科爾沁牛”（三河牛×西門塔爾牛）為母本，以引進的美系西門塔爾牛為父本雜交培育而成的群體，與ASC 相比，其培育期間的選擇強度相對較低，使其還保持著一定的遺傳改良潛力。Maiorano AM 等[18]研究發現，親本自交系雜交產生的雜種優勢賦予了后代更高的遺傳多樣性，并加強了雜種個體優勢性狀表現[19,20]。“KBC”群體在遺傳改良的過程中，保留了“科爾沁牛”群體的適應性強的同時，借助ASC 群體的基因交流提升了肉用性能。在哈迪-溫伯格平衡下，“KBC”的Ho 大于He，說明雜交系統仍然是擴大品種遺傳變異的重要方法。

3.2 PCA 和進化樹

PCA 是利用解釋群體最大遺傳方差的特征值向量對群體進行分離[21]。本結果表明，“KBC”和HXC 群體各自聚集成了兩個簇并保持了一定的距離。同時，“KBC”群體中部分個體與ASC群體出現明顯的聚集，這些“重疊”，標志著個體之間的基因型具有一定的相似性。這種相似性是由于群體之間存在相似的遺傳背景或近親關系導致了其共享某些遺傳特征或具有相似的基因型組合。

通過構建進化樹，發現“KBC”群體和ASC群體進化距離較近，與HXC 群體距離較遠。需要注意的是，環形進化樹是無根樹，僅提供了群體之間遺傳和進化距離的相對關系。群體在進化樹上的位置不代表起源和進化方向，而是反映了在進化距離上的相對位置。最初“KBC”群體在培育中使用了ASC 改良其肉用性能，在基因組上留下了選擇足跡（Selective Footprint）[22]。經過長期的人工選擇，增加了“KBC”群體和ASC 群體的遺傳相似性，因此其在進化樹上的位置處于同一分支。后續研究如果要全面地理解群體間的關系和演化過程，需要分析表型特征信息或分子鐘估計[23]等進行更深入的研究。

3.3 群體結構分析

群體遺傳結構分析基于遺傳標記信息將種群分為若干亞群，對于研究群體和個體間的遺傳關聯具有重要意義[24]。揭示了群體間的遺傳交流程度、遷移歷史以及潛在的遺傳壁壘[25]。研究發現，在K=2～3 時“KBC”群體內部分個體與ASC 祖先組成比例十分相似，在K=2～4 時，ASC 群體的組成比例比較穩定，而“KBC”群體祖先比例變化明顯。Ocampo 等[26]發現不同品種之間的分化模式和遺傳關系與每個品種的進化歷史高度一致。“KBC”群體和HXC 群體在雜交改良階段都引入過西門塔爾的血統，因此，當K值比較小的時候，更容易地捕捉到群體之間的主要遺傳差距。在“KBC”群體內，不同個體之間的祖先組成比例也在不斷變化，一方面是其原始群體本身的多態性所引起；另一方面，經過不斷的雜交改良，致使“KBC”群體存在多祖先的基因滲入。

本研究中，在假設K 值范圍內并未出現交叉驗證（Cross-Validation，CV）錯誤率的拐點，可能是研究中ASC 群體和HXC 群體使用的樣本數量較小，捕捉三個群體間的遺傳分化效果較差。當K 大于4 時，三個群體祖先比例出現明顯不同，且個體祖先比例變得相對混亂。這是因為，隨著K 值增大，模型捕捉群體之間細微的遺傳差異的靈敏度也隨之提高。

3.4 LD 衰減

通過對比群體之間的連鎖不平衡水平，對于合理利用種質資源至關重要。三個群體中，“KBC”群體衰減的最快，ASC 群體的LD-decay衰減的最慢。而人工選擇的效果可以反映在每個群體的基因組連鎖不平衡水平上[27]。在物種遺傳改良過程中，在高強度的選擇下，家畜群體中有利的染色體片段頻率迅速增加[28]，并促進群體中LD 水平的增加[29,30]。當群體中LD 水平增加，提高后續選擇的準確性，可以更快地改善目標性狀。HXC 群體已經于2021 年12月通過國家畜禽遺傳資源委員會審定，被認定為肉牛新品種。而“KBC”群體目前還處于新品種培育的過程中。因此，“KBC”群體內遠端標記的連鎖效果比其他兩個群體較差。連鎖不平衡分析可以反映相應的獨特育種選擇歷史和種群結構[27,31,32]，同時，也對后續基因組信息挖掘提供指示。值得注意的是，在研究中隨著標記之間的距離逐漸增大，三個群體的LD 衰減逐漸平緩，但是一直存在差距。后續如要進行GWAS 分析時，為了保證準確性，“KBC”群體可能會需要比HXC 和ASC 更多的基因組信息。

4 結論

本研究使用遺傳信息數據評估了“KBC”、HXC 和ASC 群體的遺傳多樣性和種群結構。“KBC”群體內中高頻（MAF≥0.3）的標記數量最多，且具有最高的遺傳多樣性和最低的連鎖不平衡衰減速度。在進化方向上，“KBC”群體與HXC 群體的距離更遠，但部分個體的種群組成與ASC 群體分化較差。通過進行群體遺傳背景分析，發現“科爾沁肉牛”經過雜交改良、橫交固定、選育提高三個階段的持續選育，在遺傳方面（分子水平）與華西牛和美系西門塔爾牛品種存在明顯差異，已具備形成新培育品種的先決條件。