童雄，羅威，閔力，張志飛，馬新燕，羅成龍，陳衛(wèi)東，徐斌，李大剛

基于全基因組SNPs分析陸豐黃牛和雷瓊牛的群體結構與遺傳多樣性特征

童雄，羅威，閔力，張志飛，馬新燕，羅成龍，陳衛(wèi)東，徐斌，李大剛

廣東省農(nóng)業(yè)科學院動物科學研究所/畜禽育種國家重點實驗室/嶺南現(xiàn)代農(nóng)業(yè)科學與技術廣東省實驗室河源分中心/廣東省畜禽育種與營養(yǎng)研究 重點實驗室，廣州 510640

【目的】研究陸豐黃牛和雷瓊牛與世界不同地域家牛中的系統(tǒng)發(fā)育關系，解析不同家牛群體的遺傳多樣性，為地方家牛資源的鑒定與保護奠定理論基礎?！痉椒ā坎杉?2頭陸豐黃牛和17頭雷瓊牛的組織樣品進行全基因組重測序，整合世界范圍內分布的24個品種92個體的NCBI公共基因組數(shù)據(jù)，共計25個品種121個個體的信息開展群體遺傳學研究。選用Bos taurus ARS-UCD1.2作為參考基因組，經(jīng)基因組序列比對與質量控制獲取高質量Reads，應用GATK軟件檢測基因組SNPs。進一步基于群體SNPs，利用系統(tǒng)進化樹構建、PCA聚類和Admixture評估進行群體結構分析，通過核苷酸多樣性（）、雜合度（）和連鎖不平衡（LD）水平研究群體的遺傳多樣性?！窘Y果】29頭廣東地方牛經(jīng)基因組測序獲取6 905 944 306個Clean Reads，每個樣本平均基因組覆蓋率為97.99%，平均測序深度分別為12.78×。整合NCBI公共基因組數(shù)據(jù)，經(jīng)質控后共檢測出14 664 391個群體SNPs。整合系統(tǒng)進化樹、PCA、Admixture的結果發(fā)現(xiàn)普通牛與瘤牛分化，瘤牛群體存在中國瘤牛與印度瘤牛分化，普通牛群體中東北亞普通牛（韓牛、延邊牛）和西藏牛與歐洲普通牛分化，溫嶺高峰牛和舟山牛從中國瘤牛群體中分化。陸豐黃牛和雷瓊牛均屬于純正的中國瘤牛，但陸豐黃牛與皖南牛之間、雷瓊牛與吉安牛之間呈現(xiàn)最近的親緣關系，說明陸豐黃牛與地域臨近的雷瓊牛屬于兩個獨立的品種。部分陸豐黃牛與雷瓊牛均存在歐洲普通牛和東北亞普通牛的血統(tǒng)混雜，且混雜比例較高，說明這兩個品種牛急需加強群體內的提純復壯。相較于歐洲普通牛和韓牛，中國家牛群體的LD衰減速率更快，核苷酸多樣性和雜合度更高，說明其遺傳多樣性更豐富。相較于其他中國家牛，陸豐黃牛和雷瓊牛的LD水平更低，雜合度更高，且核苷酸多樣性和雜合度的密度分布更為集中，說明兩者受人工選擇強度較低，且維持較高的群體遺傳多樣性。【結論】通過全基因組SNPs標記，系統(tǒng)解析了陸豐黃牛與雷瓊牛的群體遺傳結構和多樣性特征，為這兩個品種獨立分類及其保護利用提供數(shù)據(jù)支撐。

陸豐黃牛；雷瓊牛；全基因組SNPs；群體結構；遺傳多樣性

0 引言

【研究意義】家牛是人類馴化最早的家畜動物之一，距今約10 000年前由世界不同地域的原牛多次馴化而來[1-2]。中國家牛在長期的馴化和培育過程中受到人為選擇和自然選擇的作用，形成了表型豐富多樣的地方品種[3]，目前《國家畜禽遺傳資源品種名錄（2021年版）》收錄的地方品種牛達55個，后續(xù)隨著第三次全國畜禽遺傳資源普查工作的推進，會有更多的地方品種牛被鑒定收錄。陸豐黃牛和雷瓊牛是廣東省特有的地方品種，具有耐熱、耐粗飼、抗病性強等優(yōu)良特性[3-4]，因此，這兩個品種牛具有重要的開發(fā)利用價值和科學研究意義。近年來外來商業(yè)肉牛品種的引入導致這兩個品種牛的種質資源流失，保種形勢日趨嚴峻[4]；同時，由于分子鑒定數(shù)據(jù)的缺乏，陸豐黃牛與雷瓊牛在品種分類上是否歸屬于兩個獨立品種的問題，一直存在爭議。【前人研究進展】《中國畜禽遺傳資源志-牛志》在總結前人成果后，以生態(tài)類型、外貌特征、地理分布區(qū)域、分子生物學、生化遺傳學等因素作為主要的分類依據(jù)，將中國黃牛劃分為4種類型：華北型、中原型、南方型和培育品種[3]，陸豐黃牛與雷瓊牛從地域分布和表型特征上，屬于典型的南方型中國黃牛[4]。前期基于線粒體mt-DNA序列[5]、常染色體基因組[6-7]和Y染色體基因組[7]分子標記開展的分類學研究，世界范圍內家?？煞譃?大類：瘤牛（Indicine, Bos taurues indicus）和普通牛（Taurine, Bos taurus taurus），進一步細分可分為5類型：歐亞普通牛（Eurasian taurine）、東亞普通牛（East Asian taurine）、歐洲普通牛（European taurine）、中國瘤牛（Chinese indicine）和印度瘤牛（Indian indicine）[7]。中國家牛主要來源于3個祖代類群：歐亞普通牛、東亞普通牛和中國瘤牛，中國西北部家牛屬于歐亞普通牛，中國南方家牛屬于中國瘤牛，中國中北部和西南部家牛均屬于普通牛和瘤牛的雜交類群，西藏地區(qū)家牛屬于東亞普通牛[7]。前期常染色體基因組SNPs評估遺傳多樣性研究中均反映出中國家牛的遺傳多樣性高于印度瘤牛，而兩者高于歐洲普通牛；中國家牛中南方瘤牛的遺傳多樣性要高于北方的普通牛[7-8]。與此相反，mt-DNA單倍型多態(tài)性研究發(fā)現(xiàn)中國北方普通牛的遺傳多樣性較南方瘤牛更為豐富[5]。前期血液同工酶多態(tài)性研究報道了雷瓊牛與文山牛和溫嶺高峰牛等南方瘤牛親緣關系相近，且雷瓊牛的群體遺傳多樣性豐富，高于外來商用品種牛[9]?；趍t-DNA序列的母系遺傳學研究則報道雷瓊牛mt-DNA的遺傳多樣性較低[10]，在分類上屬于瘤牛，但其特有的單倍型與印度瘤牛存在顯著差異[11-12]。近年來，通過全基因組SNPs進行群體遺傳學研究中，陸續(xù)報道雷瓊牛具有純正的中國瘤牛血統(tǒng)，且中國家牛中混雜的瘤牛血統(tǒng)均來自雷瓊牛為代表的中國瘤牛而不是印度瘤牛[7-8]。有關陸豐黃牛的群體遺傳學研究較少，目前LIU等在研究8個中國南方家牛品種的群體基因組特征時，發(fā)現(xiàn)陸豐黃牛與其他7個南方牛品種的遺傳距離較遠，且存在一定比例的特異性祖代血統(tǒng)[13]，但該項研究中僅整合南方少數(shù)的幾個牛品種，缺乏其他地域家牛的信息，因此，無法判定采集的陸豐黃牛樣本是否存在遺傳混雜，進而引發(fā)群體結構解析偏差。此外，LIU等先后利用該樣本數(shù)據(jù)，比較中國南方牛群體遺傳多樣性時，發(fā)現(xiàn)陸豐黃牛和海南牛（海南雷瓊牛群體）近交系數(shù)較其他南方牛（文山牛、廣東雷瓊牛群體等）高，且其遺傳多樣性更低[13-14]?！颈狙芯壳腥朦c】近年來隨著高通量測序技術的發(fā)展，世界各地域分布的家牛陸續(xù)開展基因組學研究[6-8, 15-18]，目前通過整合世界各地域代表性牛品種的基因組信息，將是系統(tǒng)解析陸豐黃牛和雷瓊牛群體遺傳結構和遺傳多樣性的契機?！緮M解決的關鍵問題】本研究采集陸豐黃牛和雷瓊牛的群體樣本進行全基因組測序，整合其他地域代表性牛品種的基因組數(shù)據(jù)，解析陸豐黃牛與雷瓊牛之間、兩者與其他地域牛品種之間的親緣關系，指導品種的分類鑒定。同時，不同家牛群體的遺傳多樣性分析，為后續(xù)地方品種保種方案制定提供信息支持。

1 材料與方法

1.1 樣品采集與數(shù)據(jù)庫信息收集

研究于2022年6—7月，分別在廣東省陸豐市潭西鎮(zhèn)陸豐黃牛保種場和廣東省湛江市麻章區(qū)雷瓊牛保種場采集12頭陸豐黃牛（圖1）的耳組織樣和17頭雷瓊牛（圖2）的外周血樣，進行基因組DNA提取、質量檢測、文庫構建與測序。為解析廣東2個地方品種牛與世界其他地域家牛品種之間的系統(tǒng)發(fā)育關系，本研究從NCBI數(shù)據(jù)庫下載24個品種92個個體的全基因組重測序數(shù)據(jù)，所有數(shù)據(jù)庫個體的下載信息見附表1。本研究25個品種121個個體的樣本信息詳見表1。

圖2 雷瓊牛（公）

1.2 試驗方法

1.2.1 基因組DNA提取與質量檢測 采用傳統(tǒng)的苯酚-氯仿法進行組織樣基因組DNA的抽提。根據(jù)Illumina NovaSeq6000系統(tǒng)（Illumina, Inc., San Diego, CA, USA）對基因組DNA樣品的文庫構建質量要求：OD260/280介于1.8—2.0之間，DNA總量≥10μg。

1.2.2 基因組文庫構建與測序 29份組織樣品檢測合格后，每個樣品取用0.5μg DNA模板，根據(jù)Annoroad? Universal DNA Fragmentase kit V2.0（AN200101-L）及Annoroad? Universal DNA Library Prep Kit V2.0（AN200101-L）的方法及流程進行文庫制備。使用NovaSeq 6000 S4 Reagent kit V1.5試劑在NovaSeq 6000 S4平臺進行成簇和測序，運行雙端測序程序（PE），得到150 bp的雙端測序Reads。

1.3 數(shù)據(jù)處理

1.3.1 基因組序列比對與質量控制 將原始數(shù)據(jù)（Raw Data）中帶接頭和低質量的序列進行過濾，從而得到高質量序列（Clean Reads）。過濾步驟如下：（1）去除接頭污染的Reads（Reads中接頭污染的堿基數(shù)大于5 bp）；（2）去除低質量的Reads（超過50%的堿基質量值低于19的Reads）；（3）去除含N比例大于5%的Reads。

選用Bos taurus ARS-UCD1.2[19]作為參考基因組，通過基因組比對軟件BWA v.0.7.17[20]將CleanReads比對到參考基因組上。利用Samtools v.1.9軟件[21]對比對后的序列進行排序，且通過-q 4參數(shù)去除低質量Reads，應用Picard-tools v.2.20.2軟件（http:// broadinstitute.github.io/picard/）去除PCR重復。

1.3.2 基因組SNPs的檢測 應用GATK v.3.8軟件[22]中檢測群體位點變異的方法檢測SNP。通過質量值、深度、重復性等條件進行過濾篩選，過濾參數(shù)設置為：QD<2.0, ReadPosRankSum<-8.0, FS>60.0, QUAL< 30.0, DP<4.0, MQ<40.0, MappingQualityRankSum<-12.5。利用VCFtools v0.1.16軟件[23]提取常染色體SNPs后，以參數(shù)“--min-alleles 2 --max-alleles 2 --maf 0.05 --max- missing 0.9 --minGQ 15.0”進行質控，獲得群體高質量常染色體SNPs。

1.3.3 群體遺傳結構分析 系統(tǒng)進化樹分析：利用Plink v1.9軟件[24]計算群體IBS遺傳距離，通過MEGA v8.0軟件[25]構建Neighbor-joining tree。PCA分析：應用Plink v1.9軟件[24]進行LD-pruning，剔除高連鎖群體SNPs后，利用Eigenstrat v6.1.4軟件[26]進行PCA分析。Admixture分析：基于LD-pruning后的群體SNPs，使用Admixture v1.3.0軟件[27]的 unsupervised mode進行Admixture分析，計算K=2-5的群體遺傳組分。

1.3.4 群體內遺傳多樣性分析 雜合度分析：自制perl腳本運行分析，bin size設為40 000，step window size設為20 000。核苷酸多樣性分析：利用VCFtools v0.1.16軟件[23]進行計算，參數(shù)設置：“--window-pi 40000，--window-pi-step 20000”。連鎖不平衡LD分析：使用PopLDdecay V3.40軟件[28]，對各群體的LD decay進行統(tǒng)計與作圖。

2 結果

2.1 基因組數(shù)據(jù)比對質量與數(shù)據(jù)合并

本研究采集29頭廣東地方牛樣品進行基因組測序，獲得6 905 944 306個Clean Reads，每個樣本平均基因組覆蓋率為97.99%，平均測序深度分別為12.78×。將本研究采集的29個樣本的基因組測序數(shù)據(jù)與NCBI數(shù)據(jù)庫下載的24個品種92個個體的基因組測序數(shù)據(jù)進行合并分析（表1），檢測出14 664 391個群體SNPs?；谫|控條件設置SNP缺失率≤0.15，故LF08個體在后續(xù)研究中被去除。

2.2 不同地域分布家牛的群體遺傳結構

2.2.1 系統(tǒng)進化樹 基于狀態(tài)一致性（Identity By State, IBS）距離構建全群的鄰近進化樹，推斷不同品種牛之間的系統(tǒng)發(fā)育關系（圖3）。本研究采集的雷瓊牛與NCBI數(shù)據(jù)庫下載的3頭雷瓊牛聚為一簇，這說明本研究采集樣本具有品種代表性，且鄰近進化樹構建方法科學合理。進化樹上清晰反映：所有家牛出現(xiàn)普通牛（安格斯牛、海福特牛、利木贊牛、西門塔爾牛、西藏牛、韓牛、延邊牛、魯西牛、渤海黑牛、南陽牛）和瘤牛（錦江牛、皖南牛、廣豐牛、吉安牛、雷瓊牛、陸豐黃牛、文山牛、舟山牛、溫嶺高峰牛、哈里亞娜牛、吉爾牛、沙希華牛、內洛爾牛、塔帕卡牛、婆羅門牛）兩個大類群；在瘤牛群體中，中國瘤牛（錦江牛、皖南牛、廣豐牛、吉安牛、雷瓊牛、陸豐黃牛、文山牛、舟山牛、溫嶺高峰牛）與印度瘤牛（哈里亞娜牛、吉爾牛、沙希華牛、內洛爾牛、塔帕卡牛、婆羅門牛）形成了兩個獨立的類群。此外，在中國瘤牛群體中，浙江省的舟山牛和溫嶺高峰牛形成一個相對獨立的類群。陸豐黃牛和雷瓊牛均屬于中國瘤牛類群，雷瓊牛群體內僅1個個體（Leiqiong12）偏離，其他個體均能聚為一簇；而陸豐黃牛群體內出現(xiàn)2個亞群（亞群1：lufeng01、lufeng04、lufeng07、lufeng10和lufeng12；亞群2：lufeng02、lufeng05、lufeng06和lufeng09），且亞群外的2個個體（Lufeng03和Lufeng11）遺傳距離與中國北方牛（南陽牛和魯西牛）相近。

2.2.2 主成分分析（PCA） 對120個個體進行主成分分析，計算特征值的TW檢驗統(tǒng)計量，并對TW進行顯著性判定。前2個顯著性特征向量PC1和PC2分別解釋16.19%和3.44%的變異。通過PC1和PC2建立的二維PCA圖（圖4）中：PC1分布可解釋普通牛和瘤牛之間存在顯著的遺傳分化；PC2分布解釋了中國瘤牛與印度瘤牛之間的遺傳差異，且陸豐黃牛與雷瓊牛均歸屬于中國瘤牛類群，這與系統(tǒng)進化樹的結果相一致。

表1 本研究中25個品種121個個體的樣本信息

Angus：安格斯牛；Hereford：海福特牛；Limousin：利木贊牛；Simmental：西門塔爾牛；Tibetan：西藏牛；Hanwoo：韓牛；Yanbian：延邊牛；Luxi：魯西牛；Bohai：渤海黑牛；Nanyang：南陽牛；Jinjiang：錦江牛；Wannan：皖南牛；Guangfeng：廣豐牛；Jian：吉安牛；Leiqiong：雷瓊牛；Lufeng：陸豐黃牛；Wenshan：文山牛；Zhoushan：舟山牛；Wenling：溫嶺高峰牛；Hariana：哈里亞娜牛；Gir：吉爾牛；Sahiwal：沙希華牛；Nelore：內洛爾牛；Tharparkar：塔帕卡牛；Brahman：婆羅門牛。品種代號后的數(shù)字代表品種內個體的編號。下同

圖中虛線方框中僅顯示陸豐黃牛和雷瓊牛的PCA。The dotted box in the figure shows only the PCA of Lufeng cattle and Leqiong cattle.

2.2.3 Admixture祖代分析 Admixture評估家牛群體的祖代數(shù)量和遺傳混雜比例，隨著祖代群體數(shù)K的逐步遞增（K=2-5），所研究樣本群體呈現(xiàn)不同的群體結構特征（圖5）。當K=2-3時，普通牛和瘤牛分離，瘤牛群體分化為：中國瘤牛和印度瘤牛，這與PCA和系統(tǒng)進化樹的結果一致。K=4時，西藏牛、韓牛、延邊牛與歐洲普通牛出現(xiàn)明顯的分化。K=5時，溫嶺高峰牛和舟山牛從中國瘤牛中分離出來，這與系統(tǒng)進化樹的結果一致。5個陸豐黃牛和12個雷瓊牛具有單一的祖代成分，而剩余的6個陸豐黃牛和8個雷瓊牛均存在多個祖代成分，且主要受到歐洲普通牛和東亞普通牛的混雜影響。5個單一祖代成分的陸豐黃牛與系統(tǒng)進化樹中亞群1個體完全一致。

2.3 不同地域分布家牛的核苷酸多樣性和雜合度分析

對18個家牛品種（海福特牛、安格斯牛、西門塔爾牛、韓牛、西藏牛、渤海黑牛、魯西牛、南陽牛、舟山牛、婆羅門牛、溫嶺高峰牛、文山牛、雷瓊牛、陸豐黃牛、廣豐牛、吉安牛、皖南牛、錦江牛）進行種群核苷酸多樣性（）和雜合度（）評估（圖6）。除西藏牛外，中國家牛的核苷酸多樣性和雜合度的中位數(shù)與婆羅門牛相近，但均高于歐洲普通牛（海福特牛、安格斯牛、西門塔爾牛）和韓牛。歐洲普通牛和韓牛核苷酸多樣性的四分位距較中國家牛小，說明核苷酸多樣性較中國家牛更為集中。歐洲普通牛和韓牛核苷酸多樣性和雜合度分布不均勻，密度分布集中中位數(shù)與下四分位數(shù)之間，說明歐洲普通牛和韓牛中低的多樣性和雜合度值較集中分布。陸豐黃牛和雷瓊牛和的四分位距較其他中國家牛小，說明和較其他中國家牛更為集中。此外，這兩個品種牛的中位數(shù)高于其他中國家牛，且密度分布的峰值更高寬度更窄，說明其較其他中國家牛更高，且變異程度小。

品種代號后的數(shù)字代表品種內個體的編號 The number after breed code represents the number of the individual within the breed.

2.4 不同地域分布家牛的連鎖不平衡分析

為評估18個家牛群體的連鎖不平衡水平，應用全基因組SNPs繪制不同群體的2衰減圖（圖7）。相比歐洲普通牛（海福特牛、安格斯牛、西門塔爾牛）和韓牛，除西藏牛外，中國家牛（渤海黑牛、魯西牛、南陽牛、舟山牛、溫嶺高峰牛、文山牛、雷瓊牛、陸豐黃牛、廣豐牛、吉安牛、皖南牛、錦江牛）和婆羅門牛群體的LD衰減速率更快。陸豐黃牛和雷瓊牛群體的LD水平（2值）較其他中國家牛群體低。

圖6 不同家牛群體的核苷酸多樣性（a，Pi）和雜合度（b，Hp）

3 討論

3.1 不同地域分布家牛的群體遺傳結構

現(xiàn)代家牛均由共同的祖先-原牛馴化而來，前期mt-DNA遺傳方面的研究發(fā)現(xiàn)普通牛和瘤牛來源于兩個獨立分化的原牛群體，且兩者在母系遺傳分化的時間約為20萬年前[32]。本研究中系統(tǒng)進化樹、PCA和Admixture分析結果一致顯示：普通牛類群和瘤牛類群存在明顯的遺傳分化，這與之前基因芯片[6, 16, 33]和基因組重測序研究[7, 17]中得到的結果一致，從基因組水平上反映出：世界范圍內不同地域分布的家牛主要由普通原牛和瘤原牛分別馴化而來。前期考古學的研究表明，中國南方瘤牛是距今約3 500—2 500年前由印度次大陸的瘤牛遷移擴散而來[34]。CHEN等[7]和WANG等[16]從常染色體和Y染色體水平上均發(fā)現(xiàn)瘤牛群體分化為中國瘤牛和印度瘤牛，且推斷出兩者分歧時間約為3.66萬—4.96萬年，早于瘤牛馴化時間。本研究增補4個中國南方瘤牛品種（陸豐黃牛、雷瓊牛、舟山牛、溫嶺高峰牛）的基因組變異信息，同樣證實了兩個瘤牛類群之間的分化，綜合前期考古和分子遺傳學的研究成果，造成遺傳分化的原因可能是兩個類群具有獨立的祖先。本研究整合了CHEN等[7]研究中部分樣本，因此同樣發(fā)現(xiàn)西藏牛與東北亞的普通牛（韓牛、延邊牛）遺傳距離相近，而與歐洲普通牛存在明顯的遺傳分化。品種形成歷史上，舟山牛是由上海川沙、南匯等地的牛只引入后形成的品種[3, 35]，而溫嶺高峰牛是在當?shù)噩F(xiàn)有牛群的基礎上選育而來[3]。JIANG等在研究舟山牛與其他品種牛之間的親緣關系時，發(fā)現(xiàn)舟山牛與溫嶺高峰牛的親緣關系相近，且兩個品種牛與相鄰地域的中國瘤牛（皖南牛、吉安牛、雷瓊牛）形成明顯的遺傳分化[15]。本研究整合歐洲普通牛、印度瘤牛及更多的中國瘤牛品種（n=9），同樣發(fā)現(xiàn)這兩個品種牛遺傳距離相近，且從中國瘤牛類群中分離，但部分溫嶺高峰牛個體存在祖代血統(tǒng)混雜，這說明舟山牛可能在相對封閉地域環(huán)境下形成了特有的群體，而溫嶺高峰牛在近期可能受到舟山牛和其他中國瘤牛基因交流的共同影響。

圖7 不同家牛群體的連鎖不平衡LD衰減

3.2 陸豐黃牛和雷瓊牛與其他地域家牛品種間的親緣關系

通過采集廣東陸豐黃牛和雷瓊牛的群體樣本，整合NCBI數(shù)據(jù)庫公布的24個品種重測序數(shù)據(jù)，發(fā)現(xiàn)廣東2個地方品種牛均屬于中國瘤牛類群。本研究系統(tǒng)進化樹和Admixture分析發(fā)現(xiàn)5個血統(tǒng)純正的陸豐黃牛個體與地理分布相對較遠的皖南牛呈現(xiàn)最近的親緣關系，這可能是自古徽商與粵商密集商貿往來而引入的基因交流產(chǎn)生的結果，皖南牛核心產(chǎn)區(qū)所在的歙縣等地正是徽商主要的發(fā)源地，而陸豐黃牛主要分布的粵東地區(qū)是粵商商貿的代表性區(qū)域[4, 36]。進一步從Admixture遺傳混雜比例結果來看，未受到遺傳混雜影響的陸豐黃牛與雷瓊牛具有單一純正的祖代，皖南牛則主要來源于廣東來源的瘤牛祖代血統(tǒng)，同時，受到其他中國瘤牛祖代的遺傳混雜影響，因此，可以推測陸豐黃牛在商貿往來介導之下，單向流入皖南地區(qū)進而影響皖南牛的品種形成。CHEN等在研究不同地域家牛祖代特征時發(fā)現(xiàn)皖南牛與雷瓊牛具有相同的單一祖代[7]，而本研究在整合舟山牛、溫嶺高峰牛和陸豐黃牛后，發(fā)現(xiàn)皖南牛存在多種中國瘤牛血統(tǒng)混雜的情況，進而解析皖南牛的品種形成及其與陸豐黃牛的親緣關系。本研究系統(tǒng)進化樹反映出雷瓊牛和吉安牛的親緣關系最近，這與CHEN等[7，15-16]的研究結果一致，結合Admixture遺傳混雜結果可以發(fā)現(xiàn)吉安牛主要血統(tǒng)來自雷瓊牛為代表的祖代，且混雜了少量舟山牛和溫嶺高峰牛來源的祖代及歐亞普通牛祖代，因此，推斷自古廣東與江西的行政管理和貿易流通[37]介導了雷瓊牛與吉安牛相互影響，而吉安牛廣袤的中心產(chǎn)區(qū)（江西省吉安市、贛州市、撫州市、湖南省茶陵、衡山）[3]易受鄰近地區(qū)和外來牛種的基因滲入影響。前期血液同工酶多態(tài)性研究發(fā)現(xiàn)，相比于其他地域家牛，雷瓊牛與中國南方瘤牛（文山牛和溫嶺高峰牛）親緣關系更近[9]；而基于mt-DNA序列開展雷瓊牛母系起源研究中報道了雷瓊牛具有特有的單倍型，與印度瘤牛之間存在遺傳差異[11-12]。CHEN等[7-8]在家牛群體基因組研究中報道中國南方瘤牛可能是獨立馴化而來，雷瓊牛具有純正的中國瘤牛血統(tǒng)，且中國家牛中瘤牛血統(tǒng)來自雷瓊牛為代表的中國瘤牛而不是印度瘤牛，這與本研究系統(tǒng)進化樹、PCA和Admixture分析結果一致。本研究中部分陸豐黃牛與雷瓊牛均存在歐洲普通牛和東亞普通牛的血統(tǒng)混雜，且混雜比例較高，這說明廣東地方品種牛在近期可能受到外來品種牛的基因滲入影響，急需加強品種內的提純復壯。

3.3 不同地域家牛的遺傳多樣性比較

通過核苷酸多樣性和雜合度的中位數(shù)及數(shù)據(jù)分布結果均反映出中國家牛的遺傳多樣性與印度瘤牛相近，而均高于歐洲普通牛（海福特牛、安格斯牛、西門塔爾牛）和韓牛，這可能與歐洲普通牛和韓牛作為商業(yè)品種牛經(jīng)過高度選育有關。同時，本研究發(fā)現(xiàn)中國家牛和婆羅門牛群體的LD衰減速率快于歐洲普通牛和韓牛，這與XIA等[38-39]報道的LD衰減趨勢一致，同時也驗證了中國家牛和婆羅門牛具有更高的遺傳多樣性。本研究Admixture的結果反映出中國北方家牛和部分南方家牛存在廣泛的多祖代遺傳混雜，這可能是維持中國家牛群體較高的遺傳多樣性的重要來源。前期中國黃牛基因組研究中也陸續(xù)報道中國家牛中廣泛的基因交流，進而影響了秦川牛[8]、滇中牛[7]、潿洲黃牛[16]、大別山牛[40]等群體的遺傳多樣性。北方家牛中歐洲普通牛的血統(tǒng)比例較中國南方家牛更高，這與CHEN等[7,16]報道的一致。廣東2個地方品種牛較其他中國家牛具有更為集中的核苷酸多樣性和雜合度，筆者結合Admixture的結果推斷這可能是廣東2個品種受人工選擇作用的影響較低。此外，廣東2個地方品種牛具有較高的雜合度，高于文山牛等其他中國家牛品種，這可能是本研究選取的參考基因組與其遺傳距離較遠所產(chǎn)生或是群體高遺傳多樣性的體現(xiàn)。從中國家牛群體的LD水平來看，廣東2個地方品種牛的LD水平明顯低于其他中國家牛，這進一步說明這2個品種牛具有高的遺傳多樣性，且受人工選擇的強度較低。與本研究中雜合度和LD的趨勢正好相反，LIU等先后利用芯片技術研究中國南方家牛的群體遺傳多樣性時，發(fā)現(xiàn)相比文山牛、雷瓊牛（廣東雷瓊牛群體）、昭通牛等5個南方牛品種，陸豐黃牛與海南牛（海南雷瓊牛群體）具有更低的雜合度（和）和更高的近交系數(shù)（）[13-14]，這可能與不同研究中采集樣本的群體結構與多樣性有關。楊關福等在研究雷瓊牛血液蛋白的遺傳多樣性時，發(fā)現(xiàn)雷瓊牛的多態(tài)座位百分比和平均雜合度較高，且均高于荷斯坦牛[41]，這與本研究中雷瓊牛的遺傳多樣性高于歐洲普通牛的趨勢是一致的。與此相反，前期mt-DNA的遺傳多樣性研究中，雷瓊牛D-loop區(qū)的核苷酸多樣性僅為0.15%，低于秦川牛、蒙古牛、中國荷斯坦牛等品種[10]。

4 結論

本研究采集29頭廣東地方牛（12頭陸豐黃牛和17頭雷瓊牛），整合24個品種92個體的NCBI公共基因組數(shù)據(jù)，開展陸豐黃牛和雷瓊牛的遺傳多樣性與群體結構特征研究。陸豐黃牛和雷瓊牛均屬于純正的中國瘤牛，陸豐黃牛與皖南牛之間、雷瓊牛與吉安牛之間呈現(xiàn)最近的親緣關系。中國家牛群體相較于歐洲普通牛和韓牛，具有更快的LD衰減速率和更高的核苷酸多樣性和雜合度相比于其他中國家牛，陸豐黃牛和雷瓊牛LD水平更低，雜合度更高，核苷酸多樣性和雜合度的密度分布更為集中。部分陸豐黃牛與雷瓊牛均存在較高比例的歐洲普通牛和東北亞普通牛的血統(tǒng)混雜，在后續(xù)的品種保護中急需進行群體內的提純復壯。

[1] LOFTUS R T, MACHUGH D E, BRADLEY D G, SHARP P M, CUNNINGHAM P. Evidence for two independent domestications of cattle. Proceedings of the National Academy of Sciences of the United States of America, 1994, 91(7): 2757-2761.

[2] HESSE B. The first steps of animal domestication. Journal of Ethnobiology, 200626(1): 171-174.

[3] 國家畜禽遺傳資源委員會組. 中國畜禽遺傳資源志-牛志. 北京: 中國農(nóng)業(yè)出版社, 2011.

China National Commission of Animal Genetic Resources. Animal genetic resources in China. Beijing: China Agriculture Press, 2011. (in Chinese)

[4] 廣東省農(nóng)業(yè)廳, 廣東省畜牧獸醫(yī)局, 廣東省畜禽遺傳資源委員會. 廣東省地方畜禽遺傳資源志. 廣州: 廣東科技出版社, 2018.

Department of Agriculture of Guangdong Province. Records of local livestock and poultry genetic resources in Guangdong Province. Guangzhou: Guangdong Science & Technology Press, 2018. (in Chinese)

[5] XIA X, QU K, ZHANG G, JIA Y, MA Z, ZHAO X, HUANG Y, CHEN H, HUANG B, LEI C. Comprehensive analysis of the mitochondrial DNA diversity in Chinese cattle. Animal Genetics, 2019, 50(1): 70-73.

[6] GAO Y H, GAUTIER M, DING X D, ZHANG H, WANG Y C, WANG X, FARUQUE M O, LI J Y, YE S H, GOU X, HAN J L, LENSTRA J A, ZHANG Y. Species composition and environmental adaptation of indigenous Chinese cattle. Scientific Reports, 2017, 7(1): 1-14.

[7] CHEN N B, CAI Y D, CHEN Q M, LI R, WANG K, HUANG Y Z, HU S M, HUANG S S, ZHANG H C, ZHENG Z Q, et al. Whole-genome resequencing reveals world-wide ancestry and adaptive introgression events of domesticated cattle in East Asia. Nature Communications, 2018, 9(1): 1-13.

[8] MEI C G, WANG H C, LIAO Q J, WANG L Z, CHENG G, WANG H B, ZHAO C P, ZHAO S C, SONG J Z, GUANG X M, et al. Genetic architecture and selection of Chinese cattle revealed by whole genome resequencing. Molecular Biology and Evolution, 2018, 35(3): 688-699.

[9] 俞英, 文際坤, 朱芳賢, 趙開典, 宿兵, 王文, 張亞平, 劉愛華, 季維智, 施立明. 云南文山黃牛和迪慶黃牛遺傳多樣性的蛋白電泳研究. 動物學研究, 1997, 18(3): 333-339.

YU Y, WEN J K, ZHU F X, ZHAO K D, SU B, WANG W, ZHANG Y P, LIU A H, JI W Z, SHI L M. Genetic diversity of Wenshan and Diqing yellow cattle in Yunan province assayed by protein electrophoresis. Zoological Research, 1997, 18 (3): 333-339. (in Chinese)

[10] 王朝鋒, 雷初朝, 陳宏, 蔡欣, 黨瑞華, 歐江濤. 雷瓊牛mtDNA D—loop遺傳多態(tài)性研究. 黃牛雜志, 2005(5): 14-15.

WANG C F, LEI C C, CHEN H, CAI X, DANG R H, OU J T. Study on mtDNA D-loop genetic diversity in Leiqiong cattle. Journal Yellow Cattle Science, 2005(5): 14-15. (in Chinese)

[11] YU Y, NIE L, HE Z Q, WEN J K, JIAN C S, ZHANG Y P. Mitochondrial DNA variation in cattle of South China: origin and introgression. Animal Genetics, 1999, 30(4): 245-250.

[12] 耿榮慶, 常洪, 冀德君, 王蘭萍, 常國斌, 李世平, 馬國龍, 陳宏宇, 常春芳, 李永紅. 雷瓊牛母系起源的遺傳學證據(jù). 畜牧獸醫(yī)學報, 2008, 39(7)849-852.

GENG R Q, CHANG H, JI D J, WANG L P, CHANG G B, LI S P, MA G L, CHEN H Y, CHANG C F, LI Y H. Genetic evidence of maternal origin of Leiqiong cattle. Acta Veterinaria et Zootechnica Sinica, 2008, 39(7)849-852. (in Chinese)

[13] LIU Y Q, XU L Y, YANG L, ZHAO G Y, LI J Y, LIU D W, LI Y K. Discovery of genomic characteristics and selection signatures in southern Chinese local cattle. Frontiers in Genetics, 2020, 11: 533052.

[14] LIU Y Q, ZHAO G Y, LIN X J, ZHANG J H, HOU G Y, ZHANG L P, LIU D W, LI Y K, LI J Y, XU L Y. Genomic inbreeding and runs of homozygosity analysis of indigenous cattle populations in Southern China. PLoS ONE, 2022, 17(8): e0271718.

[15] JIANG L H, KON T, CHEN C Y, ICHIKAWA R, ZHENG Q Y, PEI L Y, TAKEMURA I, NSOBI L H, TABATA H, PAN H, OMORI Y, OGURA A. Whole-genome sequencing of endangered Zhoushan cattle suggests its origin and the association of MC1R with black coat colour. Scientific Reports, 2021, 11(1): 1-11.

[16] WANG X G, JU Z H, JIANG Q, ZHONG J F, LIU C K, WANG J P, HOFF J L, SCHNABEL R D, ZHAO H, GAO Y P, et al. Introgression, admixture, and selection facilitate genetic adaptation to high-altitude environments in cattle. Genomics, 2021, 113(3): 1491-1503.

[17] KIM K, KWON T, DESSIE T, YOO D, MWAI O A, JANG J, SUNG S, LEE S, SALIM B, JUNG J, et al. The mosaic genome of indigenous African cattle as a unique genetic resource for African pastoralism. Nature Genetics, 2020, 52(10): 1099-1110.

[18] 宋娜娜, 鐘金城, 柴志欣, 汪琦, 何世明, 吳錦波, 蹇尚林, 冉強, 蒙欣. 三江黃牛全基因組數(shù)據(jù)分析. 中國農(nóng)業(yè)科學, 2017, 50(1): 183-194.

SONG N N, ZHONG J C, CHAI Z X, WANG Q, HE S M, WU J B, JIAN S L, RAN Q, MENG X. The whole genome data analysis of Sanjiang cattle. Scientia Agricultura Sinica, 2017, 50(1): 183-194. (in Chinese)

[19] ROSEN B D, BICKHART D M, SCHNABEL R D, KOREN S, ELSIK C G, TSENG E, ROWAN T N, LOW W Y, ZIMIN A, COULDREY C, et al. De novo assembly of the cattle reference genome with single-molecule sequencing. GigaScience, 2020, 9(3): giaa021.

[20] LI H, DURBIN R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 2009, 25(14): 1754-1760.

[21] LI H, HANDSAKER B, WYSOKER A, FENNELL T, RUAN J, HOMER N, MARTH G, ABECASIS G, DURBIN R, SUBGROUP 1 G P D P. The sequence alignment/map format and SAMtools. Bioinformatics, 2009, 25(16): 2078-2079.

[22] MCKENNA A, HANNA M, BANKS E, SIVACHENKO A, CIBULSKIS K, KERNYTSKY A, GARIMELLA K, ALTSHULER D, GABRIEL S, DALY M, DEPRISTO M A. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Research, 2010, 20(9): 1297-1303.

[23] DANECEK P, AUTON A, ABECASIS G, ALBERS C A, BANKS E, DEPRISTO M A, HANDSAKER R E, LUNTER G, MARTH G T, SHERRY S T, MCVEAN G, DURBIN R, GROUP 1 G P A. The variant call format and VCFtools. Bioinformatics, 2011, 27(15): 2156-2158.

[24] PURCELL S, NEALE B, TODD-BROWN K, THOMAS L, FERREIRA M A R, BENDER D, MALLER J, SKLAR P, DE BAKKER P I W, DALY M J, SHAM P C. PLINK: a tool set for whole-genome association and population-based linkage analyses. The American Journal of Human Genetics, 2007, 81(3): 559-575.

[25] KUMAR S, STECHER G, LI M, KNYAZ C, TAMURA K. MEGA X: molecular evolutionary genetics analysis across computing platforms. Molecular Biology and Evolution, 2018, 35(6): 1547-1549.

[26] PATTERSON N, PRICE A L, REICH D. Population structure and eigenanalysis. PLoS Genetics, 2006, 2(12): e190.

[27] ALEXANDER D H, NOVEMBRE J, LANGE K. Fast model-based estimation of ancestry in unrelated individuals. Genome Research, 2009, 19(9): 1655-1664.

[28] ZHANG C, DONG S S, XU J Y, HE W M, YANG T L. PopLDdecay: a fast and effective tool for linkage disequilibrium decay analysis based on variant call format files. Bioinformatics, 2019, 35(10): 1786-1788.

[29] STOTHARD P, LIAO X P, ARANTES A S, DE PAUW M, COROS C, PLASTOW G S, SARGOLZAEI M, CROWLEY J J, BASARAB J A, SCHENKEL F, MOORE S, MILLER S P. A large and diverse collection of bovine genome sequences from the Canadian Cattle Genome Project. GigaScience, 2015, 4(1): s13742-15.

[30] SHIN D H, LEE H J, CHO S, KIM H J, HWANG J Y, LEE C K, JEONG J, YOON D, KIM H. Deleted copy number variation of Hanwoo and Holstein using next generation sequencing at the population level. BMC Genomics, 2014, 15: 240.

[31] BICKHART D M, XU L Y, HUTCHISON J L, COLE J B, NULL D J, SCHROEDER S G, SONG J Z, GARCIA J F, SONSTEGARD T S, VAN TASSELL C P, et al. Diversity and population-genetic properties of copy number variations and multicopy genes in cattle. DNA Research, 2016, 23(3): 253-262.

[32] ACHILLI A, BONFIGLIO S, OLIVIERI A, MALUSà A, PALA M, KASHANI B H, PEREGO U A, AJMONE-MARSAN P, LIOTTA L, SEMINO O, et al. The multifaceted origin of taurine cattle reflected by the mitochondrial genome. PLoS One, 2009, 4(6): e5753.

[33] DECKER J E, MCKAY S D, ROLF M M, KIM J, MOLINA ALCALá A, SONSTEGARD T S, HANOTTE O, G?THERSTR?M A, SEABURY C M, PRAHARANI L, et al. Worldwide patterns of ancestry, divergence, and admixture in domesticated cattle. PLoS Genetics, 2014, 10(3): e1004254.

[34] CAI D W, SUN Y, TANG Z W, HU S M, LI W Y, ZHAO X B, XIANG H, ZHOU H. The origins of Chinese domestic cattle as revealed by ancient DNA analysis. Journal of Archaeological Science, 2014, 41: 423-434.

[35] CHEN Y C. Four interesting endangered breeds of animals in China. Animal Genetic Resources Information, 1995, 16: 29-35.

[36] 魏霞, 劉正剛. 明清安徽與廣東的貿易往來. 安徽史學, 2001(4): 18-21.

WEI X, LIU Z G. On the trade of Anhui Province with Guangdong Province during the Ming and Qing dynasties. Historiography Anhui, 2001(4): 18-21. (in Chinese)

[37] 王海. 明清粵贛通道與兩省毗鄰山地互動發(fā)展研究[D]. 廣州: 暨南大學, 2008.

WANG H. Study on the interactive development of Guangdong- jiangxi passage and the adjacent mountainous areas of the two provinces in Ming and Qing dynasties[D]. Guangzhou: Jinan University, 2008. (in Chinese)

[38] XIA X T, ZHANG S J, ZHANG H J, ZHANG Z J, CHEN N B, LI Z G, SUN H X, LIU X, LYU S J, WANG X W, LI Z M, YANG P, XU J W, DING X T, SHI Q T, WANG E Y, RU B R, XU Z J, LEI C Z, CHEN H, HUANG Y Z. Assessing genomic diversity and signatures of selection in Jiaxian Red cattle using whole-genome sequencing data. BMC Genomics, 2021, 22(1): 43.

[39] ZHANG S J, YAO Z, LI X M, ZHANG Z J, LIU X, YANG P, CHEN N B, XIA X T, LYU S J, SHI Q T, et al. Assessing genomic diversity and signatures of selection in Pinan cattle using whole-genome sequencing data. BMC Genomics, 2022, 23(1): 460.

[40] GUAN X W, ZHAO S P, XIANG W X, JIN H, CHEN N B, LEI C Z, JIA Y T, XU L. Genetic diversity and selective signature in Dabieshan cattle revealed by whole-genome resequencing. Biology, 2022, 11(9): 1327.

[41] 楊關福, 張細權, 李加琪, 聶龍, 王文. 徐聞黃牛和海南黃牛血液蛋白的遺傳多樣性. 華南農(nóng)業(yè)大學學報, 1996, 17(2): 23-27.

YANG G F, ZHANG X Q, LI J Q, NIE L, WANG W. Genetic diversity in Xuwen and Hainan yellow cattle based on blood protein electrophoresis. Journal of South China Agricultural University, 1996, 17(2): 23-27. (in Chinese)

Population Structure and Genetic Diversity of Lufeng Cattle and Leiqiong Cattle Based on Genome-Wide SNPs

TONG Xiong, LUO Wei, MIN Li, ZHANG ZhiFei, MA XinYan, LUO ChengLong, CHEN WeiDong, XU Bin, LI DaGang

Institute of Animal Science, Guangdong Academy of Agricultural Sciences/State Key Laboratory of Livestock and Poultry Breeding/ Heyuan Branch of Guangdong Laboratory for Lingnan Modern Agriculture/Guangdong Key Laboratory of Animal Breeding and Nutrition, Guangzhou 510640

【Objective】Phylogenetic relationship among Lufeng cattle, Leiqiong cattle and domestic cattle in different regions worldwide, and genetic diversity of different domestic cattle populations were studied, so as to lay a theoretical foundation for the identification and protection of domestic cattle resources. 【Method】Tissue samples of 12 individuals of Lufeng cattle and 17 ones of Leiqiong cattle were collected for whole genome resequencing. Combined with another 92 cattle genomes from 24 breeds worldwide available in the NCBI database, a panel of cattle genomes comprising 121 individuals were generated from 25 breeds to carry out population genetics study. Bos taurus ARS-UCD1.2 assembly was selected as the reference genome. High-quality reads were obtained by genome alignment and quality control. Genomic SNPs were detected by GATK software. Population structure was analyzed by phylogenetic tree construction, PCA clustering, and Admixture evaluation. Genetic diversity of the populations was studied by estimating nucleotide diversity (), heterozygosity (), and linkage disequilibrium (LD). 【Result】A total of 6 905 944 306 clean reads were obtained by genome sequencing from 29 individuals of the two cattle breeds in Guangdong. Average genome coverage and average sequencing depth of each sample were 97.99% and 12.78×, respectively. After integrating the NCBI cattle genome data, 14 664 391 population SNPs were identified. The results of phylogenetic tree, PCA and Admixture showed that a primary division was found betweencattle from taurine and indicine. Moreover, indicine cattle split on Chinese and Indian cattle, and Northeast Asian cattle (Hanwoo and Yanbian) and Tibetan cattle separated from European taurine group, while Wenling cattle and Zhoushan cattle differentiated from Chinese indicine group. Lufeng cattle and Leiqiong cattle belonged to pure Chinese indicine cattle. Lufeng cattle and Wannan cattle, Leiqiong cattle and Ji'an cattle showed the closest relationship, respectively. The relationships indicated that Lufeng cattle and Leiqiong cattle in the adjacent areas belonged to two independent breeds. Some Lufeng and Leiqiong individuals were genetically admixed with European taurine and Northeast Asian taurine cattle, and the admixed proportion was high, indicating that these two breeds needed to strengthen the purification and rejuvenation within the population. Compared with European taurine cattle and Korean cattle, for Chinese domestic cattle, LD decay rate was faster, and nucleotide diversity () and heterozygosity () were higher, indicating that genetic diversity of Chinese domestic cattle was richer. Compared with other Chinese domestic cattle, LD levels of Lufeng cattle and Leiqiong cattle were lower, heterozygosity () was higher, and the density distribution of nucleotide diversity () and heterozygosity () was more concentrated, indicating that the two cattle breeds were less subject to artificial selection and maintain higher genetic diversity. 【Conclusion】Population structure and genetic diversity of Lufeng cattle and Leiqiong cattle were analyzed by genome-wide SNPs, which provided data support for independent classification and conservation of these two cattle breeds.

Lufeng cattle; Leiqiong cattle; whole-genome SNPs; population structure; genetic diversity

10.3864/j.issn.0578-1752.2023.14.014

2023-01-04；

2023-04-23

廣東省現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術體系創(chuàng)新團隊項目（2022KJ114）、嶺南現(xiàn)代農(nóng)業(yè)科學與技術廣東省實驗室河源分中心自主科研項目（DT20220023）、2022年省級鄉(xiāng)村振興戰(zhàn)略專項資金種業(yè)振興項目（2022-XBH-00-008）

童雄，E-mail：tongxiong@gdaas.cn。羅威，E-mail：luowei@gdaas.cn。童雄和羅威為同等貢獻作者。通信作者徐斌，E-mail：xubin@gdaas.cn。通信作者李大剛，E-mail：lidagang@gdaas.cn

（責任編輯 林鑒非）