馬鈞，樊安平，王武生，張金川，江曉軍，馬瑞軍，賈社強，劉飛，雷初朝，黃永震

研究報告

全基因組重測序解析秦川牛保種群遺傳多樣性和遺傳結(jié)構(gòu)

馬鈞1，樊安平2，王武生2，張金川2，江曉軍2，馬瑞軍2，賈社強2，劉飛2，雷初朝1，黃永震1

1. 西北農(nóng)林科技大學動物科技學院，楊凌 712100 2. 陜西省農(nóng)牧良種場，寶雞 722203

在畜禽資源保護中，群體遺傳多樣性和遺傳結(jié)構(gòu)是決定保種效果的重要因素。本研究采用全基因組重測序技術(shù)檢測100頭秦川牛(30頭公牛、70頭母牛)的基因組變異，通過分析群體遺傳多樣性、連續(xù)純合片段(runs of homozygosity，ROH)分布特征、親緣關(guān)系和家系結(jié)構(gòu)，對秦川牛的保種效果進行了綜合評估。結(jié)果顯示，100頭秦川牛共檢測到20,968,017個高質(zhì)量SNPs位點，平均最小等位基因頻率為0.191±0.124，平均多態(tài)信息含量為0.279±0.131，平均觀察雜合度為0.275±0.131，平均期望雜合度為0.279±0.131，表明秦川保種群遺傳多樣性較為豐富。秦川牛保種群體平均狀態(tài)同源(identity by state，IBS)遺傳距離為0.243±0.020，其中公牛為0.242±0.021，親緣關(guān)系G矩陣結(jié)果與IBS距離矩陣結(jié)果一致，均顯示秦川牛保種群部分個體間親緣關(guān)系較近。100頭秦川牛個體共檢測到8258個基因組ROH，ROH總長度為9.64 GB，平均ROH長度為1.167±1.203 Mb，69.35%的ROH是長度為0.5～1 Mb的短ROH。個體平均ROH總長度為96.40 Mb?；赗OH的平均近交系數(shù)為0.039±0.039，其中30頭秦川公牛的平均近交系數(shù)為0.044±0.035，表明部分公牛個體存在一定程度的近交積累。進化樹結(jié)果顯示，秦川牛保種群所測個體可分為8個家系，包括7個含公牛家系和1個不含公牛家系。本研究表明，秦川牛保種群的遺傳多樣性較為豐富，未出現(xiàn)較大程度近交積累，但部分個體間存在近交風險，應(yīng)強化選配以確保秦川牛資源的可持續(xù)發(fā)展。

秦川牛；全基因組重測序；遺傳多樣性；遺傳結(jié)構(gòu)；親緣關(guān)系

根據(jù)《國家畜禽遺傳資源品種名錄(2021年版)》，我國目前擁有55個地方普通牛品種，是世界上牛種資源最多的國家之一。然而，近年來引進國外優(yōu)良牛種雜交改良本土黃牛，雖然在一定程度上提高了國內(nèi)牛肉生產(chǎn)水平，但也加劇了本土純種黃牛數(shù)量的減少。畜禽種質(zhì)資源是國家的戰(zhàn)略性資源，做好地方品種的保種工作對于落實種業(yè)振興戰(zhàn)略至關(guān)重要。秦川牛是我國優(yōu)良地方黃牛品種，具有體軀高大、適應(yīng)性好、抗病力強、肉質(zhì)細致、瘦肉率高及大理石紋品質(zhì)優(yōu)良等特點，是一個種質(zhì)特色鮮明的優(yōu)良地方品種[1]。然而，目前秦川牛保種群體的遺傳多樣性及其家系結(jié)構(gòu)尚不清楚，嚴重阻礙了秦川牛種質(zhì)資源的保護和利用，亟需相關(guān)研究為該群體保種及選配工作的科學開展提供理論指導。

遺傳多樣性是生物適應(yīng)環(huán)境與進化的基石，也是評估種質(zhì)資源現(xiàn)狀的重要參考指標[2]。群體的遺傳多樣性越高，適應(yīng)環(huán)境變化的能力就越強[3]。衡量群體遺傳多樣性有助于了解群體的種質(zhì)資源現(xiàn)狀和未來趨勢，是畜禽種質(zhì)資源保護和開發(fā)利用的關(guān)鍵[4]。評估基因組親緣關(guān)系及近交系數(shù)有助于了解群體內(nèi)親緣關(guān)系和近交水平，為制定合理的選配計劃提供參考。連續(xù)純合片段(runs of homozygosity，ROH)是生物基因組中純合基因型的連續(xù)片段，是子代繼承了親本的相同單倍型而形成的[5]。研究表明，基因重組、近親交配、自然和人工選擇會影響ROH在基因組中的大小和分布[6]。ROH分析可用于評估群體遺傳多樣性[5,7]、解析近交水平[8]、推測近交歷史[9]、人工選擇的追蹤[10]及優(yōu)化育種規(guī)劃等[11]。

單核苷酸多態(tài)性(single nucleotide polymer-phism，SNP)是一種廣泛存在于基因組中的遺傳變異，被廣泛應(yīng)用于種質(zhì)資源評估、遺傳多樣性分析和系統(tǒng)發(fā)育進化等研究[12～14]。目前，常見的SNP分型方法包括全基因組重測序(whole genome sequ-en-cing，WGS)、簡化基因組測序、SNP芯片等。在馬身豬保種群中使用50K SNP芯片進行遺傳多樣性及親緣關(guān)系分析，發(fā)現(xiàn)該群體的遺傳多樣性豐富，但近交程度高[15]。在青峪豬保種群體的遺傳多樣性和家系結(jié)構(gòu)評估中，發(fā)現(xiàn)該群體在閉鎖繁育中出現(xiàn)了遺傳多樣性損失[16]。在合川黑豬保種群中，發(fā)現(xiàn)該群體遺傳多樣性較豐富，但部分個體間存在較大的近交風險[17]。利用GBS(genotyping by sequencing)簡化基因組測序技術(shù)評估麥洼牦牛3個保種群的保種效果，表明現(xiàn)有保種策略是可行的，研究結(jié)果為麥洼牦牛的保護和利用提供了一定遺傳學參考[18]。師睿等[19]基于SNP芯片數(shù)據(jù)分析新疆近交牛的基因組純合程度，為該種質(zhì)的育種規(guī)劃及后續(xù)開發(fā)利用提供了一定參考依據(jù)。李隱俠等[20]基于綿羊SNP 50K v3芯片對柯爾克孜羊進行遺傳多樣性和遺傳結(jié)構(gòu)分析，結(jié)果顯示柯爾克孜羊遺傳多樣性較豐富，群體近交程度低。馬克巖等[21]基于簡化基因組測序技術(shù)，發(fā)現(xiàn)永登七山羊的遺傳多樣性較低，與其他3個甘肅地方綿羊品種(蘭州大尾羊、灘羊和岷縣黑裘皮羊)有明顯分化。

本研究以陜西省農(nóng)牧良種場—國家秦川牛保種場的100頭秦川牛保種個體為研究對象，基于全基因組重測序數(shù)據(jù)評估了該群體的遺傳多樣性、近交系數(shù)、親緣關(guān)系及家系結(jié)構(gòu)，以解析秦川牛保種群的種質(zhì)資源現(xiàn)狀，旨在為該群體優(yōu)質(zhì)保種、科學選配及開發(fā)利用提供遺傳學依據(jù)。

1 材料與方法

1.1 實驗材料

本研究選取來自陜西省農(nóng)牧良種場(國家秦川牛保種場)100頭年齡在2～4歲之間的種牛為研究對象，其中公牛30頭，母牛70頭(表1)。為了方便后續(xù)數(shù)據(jù)分析和對比研究，采用“真實牛號＋性別”的方式標識牛只個體號。所有牛只使用一次性含EDTA抗凝劑的真空采血管頸靜脈采血，每頭牛采集血液6 mL，與管內(nèi)EDTA抗凝劑充分混勻后，低溫運輸至實驗室，于-20℃冰箱冷凍保存。血液基因組DNA使用天根生化科技(北京)有限公司的血液基因組DNA提取試劑盒(DP348-03)進行提取，將提取的DNA樣品用瓊脂糖凝膠電泳分析DNA的純度和完整性，Nanodrop檢測DNA的純度(260/280比值)，并使用Qubit 2.0(美國Invitrogen公司)對DNA濃度進行精確定量。根據(jù)美國Illumina公司推薦的基因組重測序文庫構(gòu)建方法進行文庫制備，使用Illumina NovaSeq平臺生成150 bp雙末端reads。

1.2 基因組SNPs檢測

下機的原始數(shù)據(jù)利用Trimmomatic軟件去除接頭和低質(zhì)量reads[22]，然后用BWA-MEM軟件比對至牛參考基因組(ARS-UCD 1.2)[23]。隨后，利用 SAMtools軟件對比對后的bam文件進行排序[24]，排序后使用Picard (https://broadinstitute.github.io/picard/)的MarkDuplicates模塊去除重復的reads。利用GATK v4.3.0.0的BaseRecalibrator 模塊進行堿基質(zhì)量重校正[25]。對于校正后的bam文件，利用 Qualimap 2軟件統(tǒng)計樣本的測序深度[26]。使用GATK的HaplotypeCaller模塊分別對每條常染色體進行變異檢測，并用GenotypeGVCFs模塊合并得到總VCF文件，最后使用SelectVariants模塊處理得到SNP數(shù)據(jù)集。

表1 秦川牛樣品采集信息

個體號中的M和F分別表示公牛和母牛。

1.3 基因型數(shù)據(jù)質(zhì)控

為得到高質(zhì)量基因組變異數(shù)據(jù)，首先利用GATK v4.3.0.0的VariantFiltration模塊對SNP進行過濾，過濾參數(shù)采用GATK官方推薦的標準：①Variant Confidence/Quality by Depth(QD)<2.0；②RMS Mapping Quality(MQ)< 40.0；③Phred-scaled-value using Fisher’s exact test to detect strand bias(FS)> 60；④Z-score according to the Wilcoxon rank sum test of Alt. Ref read mapping qualities(MQRankSum)< –12.5；⑤Z-score according to the Wilcoxon rank sum test of Alt. Ref read position bias(ReadPos-Rank-Sum)<–8；⑥Symmetric Odds Ratio of 2×2 contin-gency table to detect strand bias(SOR)>3.0[25]。然后，使用BCFtools v1.8過濾得到雙等位SNP集合[27]。最后，使用PLINK v1.9軟件去除不符合哈代溫伯格平衡、MAF值小于0.05和性染色體上的SNP位點，以及基因型缺失大于10%的個體[28]，獲得20,968,017個常染色體SNPs位點用于后續(xù)分析。最后使用PLINK基于位點間的連鎖不平衡進行過濾，參數(shù)為“--indep-pairwise 50 5 0.2”，通過過濾后的1,052,677個SNPs位點用于親緣關(guān)系和家系結(jié)構(gòu)分析。

1.4 遺傳多樣性分析

為了解現(xiàn)有秦川牛保種群體的基因組遺傳多樣性現(xiàn)狀。本研究利用PLINK v1.9軟件通過計算期望雜合度(expected heterozygosity，He)、觀察雜合度(observed heterozygosity，Ho)、最小等位基因頻率(minor allele frequency，MAF)和多態(tài)信息含量(polymorphism information content，PIC)，對秦川牛保種群體的遺傳多樣性進行評估[28]。

1.5 全基因組ROH統(tǒng)計分析

本研究利用PLINK v1.9軟件統(tǒng)計每個樣本的基因組純合性片段，參數(shù)設(shè)置為：滑窗大小設(shè)為50個SNP(--homozyg-window-snp 50)，一段ROH中每50 kb必須有1個SNP(--homozyg-density 50)，滑動窗口允許3個SNP位點雜合(--homozyg-window- het 3)，窗口允許缺少數(shù)設(shè)為5(--homozyg-window- missing 5)[28]。對ROH在秦川牛群體中的分布、長度和數(shù)目情況等進行統(tǒng)計分析。為更加了解群體歷史，進一步將其根據(jù)物理長度分為短片段(0.5～1 Mb)、中片段(1～5 Mb)和長片段(>5 Mb)三類進行分類統(tǒng)計?；赗OH的基因組近交系數(shù)(F)由每個個體檢測到的常染色體ROH片段總長度占基因組常染色體總長度(2.5 Gb)的比值進行計算[29]，公式為：

1.6 秦川牛保種群體親緣關(guān)系分析

為了評估個體間的親緣關(guān)系，采用IBS遺傳距離和G矩陣兩種方法進行分析[30]。使用PLINK v1.9軟件計算個體間的遺傳距離，構(gòu)建IBS距離矩陣。使用GCTA v1.94.1軟件構(gòu)建個體間的親緣關(guān)系G矩陣。最后，使用R語言繪制熱圖可視化個體之間的IBS遺傳距離和G矩陣結(jié)果，并分析保種群個體間的親緣關(guān)系。

1.7 秦川牛群體家系構(gòu)建

為了解析試驗群體的家系結(jié)構(gòu)，采用PHYLIP v3.69軟件基于個體之間IBS距離矩陣使用鄰接法(neighbor-joining，NJ)構(gòu)建系統(tǒng)發(fā)育樹[31]，并使用FigTree v1.3.1軟件進行可視化[32]。然后，結(jié)合親緣關(guān)系分析結(jié)果，以個體間親緣關(guān)系系數(shù)大于等于0.1為標準劃分家系，分析秦川牛保種群體的家系結(jié)構(gòu)[17,20,30]。

2 結(jié)果與分析

2.1 秦川?；蚪M遺傳多樣性分析結(jié)果

本研究對100頭秦川牛進行基因組重測序，產(chǎn)生22,024,263,265個高質(zhì)量reads，其中99.84%的reads可比對到參考基因組上，所有樣本的平均測序深度達到11.68×。通過GATK v4.3.0.0、BCFtools v1.8和PLINK v1.9等一系列軟件對實驗樣本的基因組變異進行質(zhì)控，最終獲得一個包含20,968,017個SNPs位點的秦川牛高質(zhì)量基因組SNP變異集合。在常染色體中，1號染色體上鑒定到的SNPs數(shù)目最多，為1,316,538個；25號染色體的SNPs數(shù)目最少，為373,259個(圖1)。功能注釋結(jié)果如圖2所示，這些SNPs主要分布在基因間區(qū)(12,424,579，59.25%)或內(nèi)含子區(qū)(7,924,936，37.80%)；外顯子僅占總SNPs數(shù)量的0.74%(156,185)，包括51,581個非同義SNPs和100,071個同義SNPs。

所有個體的基因型平均檢出率為0.997?；讷@得的高質(zhì)量SNPs評估秦川牛保種群體的遺傳多樣性水平，結(jié)果顯示，秦川牛群體基因組SNP變異的MAF介于0.05～0.1的比例最高，為31.54%，介于0.4～0.5的占比最低，為9.2%，平均MAF為0.191±0.124(圖3A，表2)?；蚪MSNPs的多態(tài)信息含量范圍介于0.098～0.500之間，平均多態(tài)信息含量為0.279±0.131，具體分布如圖3B所示。其中，多態(tài)信息含量在0～0.25之間的SNPs有10,145,311個，占48.38%；多態(tài)信息含量在0.25～0.5之間的SNPs有10,822,706個，占比51.62%。秦川牛保種群平均觀察雜合度為0.275±0.131，平均期望雜合度為0.279±0.131(表2)，平均觀察雜合度略小于平均期望雜合度，但二者相差不大，推測其可能受到了選擇或近交。

圖1 秦川?；蚪MSNPs在每條染色體上的分布情況

圖2 秦川?；蚪MSNPs功能注釋

圖3 秦川牛保種群體遺傳多樣性分析結(jié)果

A：最小等位基因頻率占比圖；B：多態(tài)信息含量占比圖。

表2 秦川牛群體遺傳多樣性

2.2 秦川牛保種群體IBS距離矩陣分析結(jié)果

利用Plink v1.9軟件計算秦川牛個體間IBS遺傳距離，結(jié)果顯示100頭秦川牛個體間IBS距離值介于0.128～0.278之間，平均IBS值為0.243±0.020，表明個體間的平均遺傳距離較遠，且差異較大；30頭種公牛之間IBS遺傳距離在0.128～0.269之間，平均遺傳距離為0.242±0.021，表明種公牛之間的遺傳距離均較遠。秦川牛保種群體IBS距離矩陣可視化結(jié)果如圖4所示。大部分秦川牛個體間IBS遺傳距離較遠，呈中等程度的親緣關(guān)系(圖4顏色較深的方格)；部分個體間的IBS遺傳距離較近，存在較高的親緣關(guān)系(圖4顏色較淺的方格)，表明該保種群存在較大的潛在近交風險。

2.3 秦川牛保種群體G矩陣分析結(jié)果

基于全基因組SNPs位點構(gòu)建的G矩陣能真實反映個體間的親緣關(guān)系。利用GCTA軟件構(gòu)建基因組關(guān)系G矩陣，進一步分析秦川牛保種群體的親緣關(guān)系。結(jié)果如圖5所示，顏色越接近紅色表示親緣關(guān)系越近。同IBS距離矩陣結(jié)果相似，大部分個體間的親緣關(guān)系呈中等程度(圖5中顏色較淺的方格)，但部分牛之間的親緣關(guān)系較近(圖5中顏色較深的方格)，均表明秦川牛保種群體存在較大的近交風險，必須強化選配措施。

橫、縱坐標為本研究所使用的100頭秦川牛個體，每一個小方格代表兩個個體間的遺傳距離值，顏色越接近紅色表明兩個個體間的遺傳距離越大，反之亦然。

2.4 秦川牛保種群體基因組ROH基本統(tǒng)計及近交系數(shù)結(jié)果

本研究在100頭秦川牛保種群體中共檢測到8258個ROH片段，ROH總長度為9.64 GB；平均每頭秦川牛有82.58個ROH片段，個體平均ROH總長度為96.40 Mb，單個ROH片段平均長度為1.167±1.203 Mb。圖6A展示了秦川牛群體單個ROH長度分布情況，發(fā)現(xiàn)大部分ROH長度都在4 Mb以下，長度為0.5～1 Mb的ROH最多，占比69.35%，其次是長度為1～2 Mb的ROH，比例為17.89%，大于4 Mb的占比很小。在整個研究群體中，不同牛只檢測出的ROH數(shù)目變化較大，從最少的25個ROH片段到最多的328個ROH片段。個體基因組ROH總長度范圍在21.50～880.90 Mb之間，在100頭秦川牛中，79頭牛的ROH長度在50～100 MB之間(圖6B)。ROH在染色體上的數(shù)量分布如圖6C所示，在秦川牛的29條常染色體上，1號染色上分布數(shù)目最多，為637個，28和29號染色體上分布數(shù)目最少，均為85個。從染色體上ROH分布長度看，1號染色體ROH總長度最長，為775.71 Mb，而29號染色體分布的ROH長度最短，為110.07 Mb(圖6D)。為深入了解群體歷史，進一步將ROH依據(jù)物理長度將其分成長、中、短三類，結(jié)果如圖6E所示。0.5～1 Mb之間的短ROH片段數(shù)目和長度分別是5727個和3764.76 Mb，1～5 Mb之間的中ROH片段數(shù)目和長度分別為2344個和4594.06 Mb，大于5 Mb的長ROH片段數(shù)目和長度分別為187個和1281.38 Mb。短ROH數(shù)目在總ROH數(shù)目中比例最大，但其長度在總長度中占比并非最大，在ROH總長度中占比最多的是1～5 Mb的ROH?；赗OH的近交系數(shù)F分析結(jié)果顯示(圖6F)，現(xiàn)有保種群的個體近交系數(shù)范圍在0.009～0.352之間，平均近交系數(shù)為0.039± 0.039，其中母牛個體號為1294F的個體近交系數(shù)最大為0.352。30頭秦川公牛F為0.044±0.035，略高于保種群整體平均水平，表明部分公牛個體已出現(xiàn)一定程度的近交積累。

圖5 秦川牛保種群G矩陣分析可視化結(jié)果

橫、縱坐標為本研究所使用的100頭秦川牛個體，每一個小方格代表兩個個體間的親緣關(guān)系，顏色越接近紅色表明兩個個體間的親緣關(guān)系越近，反之亦然。

圖6 秦川牛基因組ROH統(tǒng)計分析

A：秦川?；蚪M單個ROH占比統(tǒng)計；B：秦川牛個體ROH長度的樣本數(shù)分布；C：ROH在不同染色體上的數(shù)量分布情況；D：ROH在不同染色體上的長度分布情況；E：短、中和長三類ROH的數(shù)量和長度統(tǒng)計圖；F：基于ROH的近交系數(shù)可視化。

2.5 秦川牛保種群體家系結(jié)構(gòu)分析

由于種公牛在保種群中的重要性，本研究使用NJ法對30頭公牛進行了聚類分析。結(jié)合基因組親緣關(guān)系分析結(jié)果，以公牛間基因組親緣系數(shù)大于等于0.1為標準進行聚類分析[17,20,30]。結(jié)果如圖7所示，現(xiàn)有種公牛樣本可分為7個家系，分別命名為家系A(chǔ)、家系B、家系C、家系D、家系E、家系F和家系G。進化樹中用同一顏色標注的樣本被評估為同一家系(表3)。其中，存在部分交叉?zhèn)€體。例如，樣本1329M既與家系E的個體親緣關(guān)系較近，又與家系D的部分個體(1278M)親緣關(guān)系較近，所以同時被歸并到這兩個家系中。最終，這7個家系分別含有10、7、2、3、3、2和4個公牛。根據(jù)所有個體間的G矩陣親緣關(guān)系分析結(jié)果，將與公牛的親緣關(guān)系值大于0.1的母牛與公牛劃分為同一家系，最終將保種群所測個體分為7個含有公牛的家系。在這7個家系外，還有13頭母牛與所檢測的公牛親緣關(guān)系都比較遠，將其劃分至1個不包含公牛血緣的家系H(表3)。此外，部分家系中的母牛存在交叉現(xiàn)象，同時被歸并到不同的家系中。其中母牛1109F同時被劃分到B和E家系中，母牛1409F同時被劃分到家系D和E，母牛1158F同時被劃分到家系E和F，母牛1185F同時被劃分到家系A(chǔ)和G，母牛1144F同時被劃分到家系A(chǔ)、B、C和F中。

圖7 秦川牛保種群體公牛進化樹結(jié)果

3 討論

畜禽種質(zhì)資源是保障國家畜產(chǎn)品供給的戰(zhàn)略性資源，是種業(yè)振興的物質(zhì)基礎(chǔ)。保種工作旨在維護種群內(nèi)的基因多樣性，減少遺傳漂變和基因缺失等的不良影響，保障群體的健康發(fā)展和適應(yīng)環(huán)境能力，促進畜禽品種改良和創(chuàng)新[33]?？茖W評估保種效果可以確定現(xiàn)有群體是否存在遺傳風險以及采取何種措施進行改善。秦川牛是我國優(yōu)良的地方黃牛品種，在我國畜牧業(yè)發(fā)展中發(fā)揮著極其重要的作用，但現(xiàn)有保種群的全基因組遺傳多樣性和群體結(jié)構(gòu)尚未得到系統(tǒng)性研究。本研究采用全基因組重測序技術(shù)檢測秦川牛保種場的100頭秦川牛，鑒定到20,968,017個SNPs位點，并基于這些高質(zhì)量SNPs從遺傳多樣性、親緣關(guān)系、近交水平和家系結(jié)構(gòu)等方面對保種群的保種效果進行評估，為秦川牛保種群制定保種計劃提供科學依據(jù)。

遺傳多樣性反映了物種或種群在進化過程中的適應(yīng)性和生存能力，種群內(nèi)的基因多樣性越高，意味著群體更具適應(yīng)性和生存能力[21,34]。雜合度水平對生物體的機能表現(xiàn)有深遠的影響，有研究發(fā)現(xiàn)生物體的雜合度與其適應(yīng)性(例如生長速率、存活率或繁殖力)之間存在正相關(guān)[35]，而低雜合度與近親繁殖抑郁癥和適應(yīng)性降低有關(guān)[36]。Gao等[37]基于BovineSNP50芯片發(fā)現(xiàn)包括32頭秦川牛在內(nèi)45個世界范圍牛品種的Ho為0.142～0.327，其中32頭秦川牛的Ho為0.297。本研究中，秦川?；蚪M的Ho和He分別為0.275±0.131和0.279±0.131，屬于中等雜合度水平，但雜合度有所降低。同時，與Zhang等[38]基于7003個SNPs在30頭秦川牛的研究結(jié)果(Ho：0.370；He：0.364)也相差較大，這可能與研究所用位點數(shù)目、位點選擇、采樣地區(qū)以及樣本數(shù)目有關(guān)。本研究采用全基因組重測序進行分析，相對于商業(yè)芯片，獲得的基因組SNP變異更加全面。然而，試驗所用保種群個體經(jīng)長期閉群繁育，可能存在遺傳多樣性有所降低的問題。此外，秦川?；蚪M的平均觀察雜合度略低于平均期望雜合度，但二者相差不大，表明該保種群具有較低的雜合子喪失，推測其可能受到了一定程度選擇壓力，但未出現(xiàn)顯著的遺傳漂移或近交等。最小等位基因頻率是對每個群體中不常見等位基因頻率的估計[39]。秦川牛群體MAF分布結(jié)果顯示，該群體低等頻率的等位基因數(shù)目較多，而較高MAF值的SNPs相對較少，表明該群體的遺傳多樣性程度相對較高?；蚪M多態(tài)信息含量反映某一位點在群體中的多態(tài)程度，常用來評估群體遺傳多樣性。根據(jù)Botstein等[40]報道，對于微衛(wèi)星位點，PIC值大于0.5時，可被視為高多態(tài)性位點；反之，當PIC值小于0.25時，歸類為低多態(tài)性位點；而當PIC值處于0.25～0.5之間時，則屬于中度多態(tài)性位點。本研究發(fā)現(xiàn)秦川牛保種群基因組SNPs位點的PIC值分布相對均勻，平均PIC值為0.279±0.131，大部分位點PIC值都分布在0.1～0.5之間，其中51.62%的SNPs位點多態(tài)信息含量在0.25～0.5之間，說明該群體的遺傳多樣性較為豐度，具有中度遺傳變異，尤其在較高的PIC值區(qū)間內(nèi)。以上結(jié)果表明，近年來秦川牛保種群體所采用的保種措施有效，群體遺傳多樣性并未出現(xiàn)明顯丟失。

表3 秦川牛保種群體家系劃分結(jié)果

個體號中的M和F分別表示公牛和母牛。

ROH在基因組上的分布情況反映了群體遺傳多樣性、近交歷史及近交水平等[5,7～9]，并且研究發(fā)現(xiàn)基于ROH計算的近交系數(shù)(F)與基于系譜的近交系數(shù)(F)具有較高相關(guān)性，且比F、F等近交系數(shù)計算方法更能反映群體真實近交水平[41]。短ROH表示發(fā)生在較遠時期的親緣關(guān)系，長ROH反映近期發(fā)生的近交[42]。長ROH數(shù)目越多，家系內(nèi)存在近交的可能性越高。本研究在100頭秦川牛保種個體中共檢測到8258個ROH片段，其中長度為0.5～1 Mb的短ROH數(shù)目最多，長度為5 Mb以上的ROH數(shù)目最少，揭示秦川牛保種群在較近時期并未出現(xiàn)較大程度近交?；赗OH的近交系數(shù)分析顯示，秦川牛保種群的近交系數(shù)為0.039±0.039，處于相對較低的近交水平，但高于先前研究報道的麥洼牦牛保種群的平均近交系數(shù)(0.0168)[18]。30頭公牛的平均近交水平略高于整體平均水平，表明部分公牛個體已出現(xiàn)一定程度的近交積累，在后續(xù)種公牛選留中建議保留近交水平較低的個體，以避免保種群遺傳多樣性喪失。此外，研究還發(fā)現(xiàn)部分個體具有較高的近交水平，其中母牛1294F個體的近交系數(shù)達到0.352，出現(xiàn)較高程度的近交。在后續(xù)保種工作中，對于這些近交系數(shù)相對較高的個體要適當引入外血，或使用與其關(guān)系較遠家系的公牛配對以降低其后代的近交系數(shù)。

目前，保種場普遍采用閉群繁育方式來保護畜禽遺傳資源，種群遺傳多樣性和群體結(jié)構(gòu)極易受到人工選擇和雜交等因素的影響[43]。因此，理清保種群體個體間的親緣關(guān)系和群體遺傳結(jié)構(gòu)對于保存群體遺傳多樣性及保障保種群體的可持續(xù)發(fā)展至關(guān)重要?；贗BS遺傳距離和G矩陣的親緣關(guān)系分析結(jié)果顯示，部分秦川牛個體間的親緣關(guān)系相對較近，這些個體間具有較大的近交風險，在后續(xù)配種工作中，要避免這些親緣關(guān)系較近的公母牛進行配對。此外，本研究中秦川牛保種群的平均IBS遺傳距離為0.243±0.020，公牛平均IBS遺傳距離為0.242±0.021，公牛的平均IBS遺傳距離略低于整個秦川牛保種群，這與先前蔡春波等[15]和劉彬等[16]分別在馬身豬和青峪豬保種群的研究結(jié)果有所不同，這表明部分秦川公牛之間親緣關(guān)系相對較近，因此在后期選留種公牛時需控制好親緣關(guān)系系數(shù)。由于家系構(gòu)成在保存畜禽遺傳多樣性中的重要作用，基于30頭種公牛的IBS遺傳距離構(gòu)建了系統(tǒng)進化樹，并以個體間的親緣關(guān)系遠近程度，將本研究所測秦川牛劃分了7個擁有種公牛的家系和1個不含公共公牛的家系，并且部分家系牛只存在交叉現(xiàn)象。

為了驗證基因組分析在保種實踐中的可靠性，本研究查閱了保種場的系譜記錄。結(jié)果發(fā)現(xiàn)，母牛1294F的父親和爺爺是同一頭公牛，證實了該牛只基因組近交系數(shù)過大是配種失誤產(chǎn)生的近交造成的。本研究進一步將家系劃分結(jié)果與保種場系譜記錄進行比較，發(fā)現(xiàn)在同一家系中，親緣關(guān)系熱圖上位置相近且顏色相似的個體大部分是同一頭種公?；蛟摷蚁灯渌５暮蟠?。此外，一些個體與兩個或多個家系的部分牛只都具有較高的親緣關(guān)系值，是由于其擁有相同的母親或祖代。因此，基于基因組的家系研究結(jié)果可有效填補個體系譜記錄不全的缺陷，優(yōu)化保種群的家系構(gòu)成，并發(fā)現(xiàn)畜禽保種中的潛在問題。在后續(xù)保種選配過程中，可以結(jié)合家系結(jié)構(gòu)和保種場的系譜記錄，制定合理的配種方案，減少同一家系公母個體之間的配種，從而科學地開展保種工作，避免種群出現(xiàn)近交和品種內(nèi)分化。

本研究基于高通量全基因組重測序數(shù)據(jù)從遺傳多樣性、ROH分布特征、親緣關(guān)系和家系結(jié)構(gòu)等方面對秦川牛的保種效果進行綜合評估，發(fā)現(xiàn)秦川牛保種群的遺傳多樣性較為豐富；保種群體未出現(xiàn)大面積近交，但部分個體間存在較大的近交風險。本研究結(jié)果可為秦川牛的保種選配工作提供一定的科學依據(jù)。

[1] Wang B, Ma GJ, Zhang Y, Zhou FL, Cheng ZG, Hao P. External factors and improvement model of Qinchuan cattle hybrid improvement., 2020, (21): 39–40.王勃, 馬國際, 張巖, 周福來, 程助國, 郝鵬. 秦川牛雜交改良外部因素及改良模式. 畜牧獸醫(yī)科學(電子版), 2020, (21): 39–40.

[2] Yan BY. Population genomics research on Chaidamu yellow cattle[Dissertation]. Lanzhou University, 2019.閆碧瑤. 柴達木黃牛的群體基因組學研究[學位論文]. 蘭州大學, 2019.

[3] Bravo S, Larama G, Qui?ones J, Paz E, Rodero E, Sepúlveda N. Genetic diversity and phylogenetic relationship among araucana creole sheep and Spanish sheep breeds., 2019, 172: 23–30.

[4] Jia XB, Ding P, Chen SY, Zhao SK, Wang J, Lai SJ. Analysis of MC1R, MITF, TYR, TYRP1, and MLPH genes polymorphism in four rabbit breeds with different coat colors., 2021, 11(1): 81.

[5] Liu JX, Wei X, Deng TY, Xie R, Han JL, Du LX, Zhao FP, Wang LX. Genome-wide scan for run of homozygosity and identification of corresponding candidate genes in sheep populations., 2019, 50(8): 1554–1566.劉家鑫, 魏霞, 鄧天宇, 謝銳, 韓建林, 杜立新, 趙福平, 王立賢. 綿羊全基因組ROH檢測及候選基因鑒定. 畜牧獸醫(yī)學報, 2019, 50(8): 1554–1566.

[6] Broman KW, Weber JL. Long homozygous chromosomal segments in reference families from the centre d'Etude du polymorphisme humain., 1999, 65(6): 1493–1500.

[7] Al-Mamun HA, Clark SA, Kwan P, Gondro C. Genome- wide linkage disequilibrium and genetic diversity in five populations of Australian domestic sheep., 2015, 47: 90.

[8] Purfield DC, Berry DP, McParland S, Bradley DG. Runs of homozygosity and population history in cattle., 2012, 13: 70.

[9] Browning SR, Browning BL. Identity by descent between distant relatives: detection and applications., 2012, 46: 617–633.

[10] Liu G, Sun FZ, Zhu FX, Feng HY, Han X. Runs of homozygosity and its application on livestock genome study., 2019, 41(4): 304–317.劉剛, 孫飛舟, 朱芳賢, 馮海永, 韓旭. 連續(xù)性純合片段在畜禽基因組研究中的應(yīng)用. 遺傳, 2019, 41(4): 304–317.

[11] Marras G, Gaspa G, Sorbolini S, Dimauro C, Ajmone- Marsan P, Valentini A, Williams JL, Macciotta NPP. Analysis of runs of homozygosity and their relationship with inbreeding in five cattle breeds farmed in Italy., 2015, 46(2): 110–121.

[12] Ma J, Gao X, Li JY, Gao HJ, Wang ZZ, Zhang LP, Xu LY, Gao H, Li HW, Wang YH, Zhu B, Cai WT, Wang CY, Chen Y. Assessing the genetic background and selection signatures of Huaxi cattle using high-density SNP array., 2021, 11(12): 3469.

[13] Van der Nest MA, Hlongwane N, Hadebe K, Chan WY, van der Merwe NA, De Vos L, Greyling B, Kooverjee BB, Soma P, Dzomba EF, Bradfield M, Muchadeyi FC. Breed ancestry, divergence, admixture, and selection patterns of the Simbra crossbreed., 2020, 11: 608650.

[14] Zhang SJ, Yao Z, Li XM, Zhang ZJ, Liu X, Yang P, Chen NB, Xia XT, Lyu SJ, Shi QT, Wang EY, Ru BR, Jiang Y, Lei CZ, Chen H, Huang YZ. Assessing genomic diversity and signatures of selection in Pinan cattle using whole-genome sequencing data., 2022, 23(1): 460.

[15] Cai CB, Zhang XL, Zhang WF, Yang Y, Gao PF, Guo XH, Li BG, Cao GQ. Evaluation of genetic structure in Mashen pigs conserved population based on SNP chip., 2021, 52(4): 920–931.蔡春波, 張雪蓮, 張萬峰, 楊陽, 高鵬飛, 郭曉紅, 李步高, 曹果清. 運用SNP芯片評估馬身豬保種群體的遺傳結(jié)構(gòu). 畜牧獸醫(yī)學報, 2021, 52(4): 920–931.

[16] Liu B, Shen LY, Chen Y, Li Q, Liao K, Guo ZY, Zhang SH, Zhu L. Analysis of genetic structure of conservation population in Qingyu pig based on SNP chip., 2020, 51(2): 260–269.劉彬, 沈林園, 陳映, 李強, 廖坤, 郭芝忺, 張順華, 朱礪. 基于SNP芯片分析青峪豬保種群體的遺傳結(jié)構(gòu). 畜牧獸醫(yī)學報, 2020, 51(2): 260–269.

[17] Long X, Chen L, Wu PX, Zhang TH, Pan HM, Zhang L, Wang JY, Guo ZY, Chai J. Evaluation of the genetic structure and selection signatures in Hechuan Black pigs conserved population., 2023, 54(5): 1854–1867.龍熙, 陳力, 吳平先, 張廷煥, 潘紅梅, 張亮, 王金勇, 郭宗義, 柴捷. 合川黑豬保種群遺傳結(jié)構(gòu)及選擇信號分析. 畜牧獸醫(yī)學報, 2023, 54(5): 1854–1867.

[18] Ma SL, Li XW, Li X, Xie SQ, Liu YL, Tang J, Jiang MF. Assessment of genetic structure of 3 Maiwa yak preserved populations based on genotyping-by-sequencing techno-logy., 2022, 31(9): 183–194.馬士龍, 李小偉, 李響, 謝書瓊, 劉益麗, 唐嬌, 江明鋒. 基于GBS簡化基因組測序評估3個麥洼牦牛保種群的遺傳結(jié)構(gòu)研究. 草業(yè)學報, 2022, 31(9): 183–194.

[19] Shi R, Zhang Y, Wang YC, Huang T, Lu GC, Yue T, Lu ZX, Huang XX, Wei XP, Feng ST, Chen J, Kagedeer Wulan, Abulizi Ruxianguli, Abulizi Ruxianguli. The evaluation of genomic homozygosity for Xinjiang inbred population by SNP panels., 2020, 42(5): 493–506.師睿, 張毅, 王雅春, 黃濤, 盧國昌, 岳濤, 盧振西, 黃錫霞, 衛(wèi)新璞, 馮書堂, 陳軍, 烏蘭·卡格德爾, 茹先古麗·阿不力孜, 努爾胡馬爾·木合塔爾. 利用SNP芯片信息評估新疆近交牛基因組純合度. 遺傳, 2020, 42(5): 493–506.

[20] Li YX, Yashengjiang Nasier, Sai Like Duman, Qian Y, Cao SX, Wang WL, Meng CH, Zhang J, Zhang JL. Evaluation of genetic diversity and genetic structure in Kirgiz sheep population cased on SNPs chip., 2023, 54(2): 572–583.李隱俠, 牙生江·那斯爾, 賽里克·都曼, 錢勇, 曹少先, 王偉列, 孟春花, 張俊, 張建麗. SNP芯片評估柯爾克孜羊群體遺傳多樣性和遺傳結(jié)構(gòu). 畜牧獸醫(yī)學報, 2023, 54(2): 572–583.

[21] Ma KY, Han JT, Bai YQ, Li TT, Ma YJ. Genetic diversity analysis of Yongdeng seven goats based on specific-locus amplified fragment sequencing., 2023, 54(5): 1939–1950.馬克巖, 韓金濤, 白雅琴, 李討討, 馬友記. 基于簡化基因組測序的永登七山羊遺傳多樣性分析. 畜牧獸醫(yī)學報, 2023, 54(5): 1939–1950.

[22] Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data., 2014, 30(15): 2114–2120.

[23] Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM., 2013.

[24] Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, 1000 Genome Project Data Processing Subgroup. The Sequence Alignment/Map format and SAMtools., 2009, 25(16): 2078–2079.

[25] McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, DePristo MA. The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data., 2010, 20(9): 1297– 1303.

[26] Okonechnikov K, Conesa A, García-Alcalde F. Qualimap 2: advanced multi-sample quality control for high- throughput sequencing data., 2016, 32(2): 292–294.

[27] Li H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data., 2011, 27(21): 2987–2993.

[28] Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, Maller J, Sklar P, de Bakker PIW, Daly MJ, Sham PC. PLINK: a tool set for whole-genome association and population-based linkage analyses., 2007, 81(3): 559–575.

[29] Zhao GY. Association and prediction of traits based on genomic homozygous segments in beef cattle[Disserta-tion]. Chinese Academy of Agricultural Sciences, 2021.趙國耀. 基于肉?；蚪M純合片段的性狀關(guān)聯(lián)與預測 [學位論文]. 中國農(nóng)業(yè)科學院, 2021.

[30] Liu CL, Lu D, Zhou QY, Wan MC, Hao XD, Zhang XH, Zheng RQ, Ji HY. Analysis of population genetic structure of Hang pigs by high density SNP chips., 2022, 53(8): 2502–2513.劉晨龍, 盧丹, 周泉勇, 萬明春, 郝曉東, 張獻賀, 鄭瑞強, 季華員. 利用高密度SNP芯片分析杭豬的群體遺傳結(jié)構(gòu). 畜牧獸醫(yī)學報, 2022, 53(8): 2502–2513.

[31] Baum BR. PHYLIP: phylogeny inference package. Version 3.2. Joel felsenstein., 1989, 64(4): 539–541.

[32] Xu LY, Bickhart DM, Cole JB, Schroeder SG, Song JZ, Tassell CPV, Sonstegard TS, Liu GE. Genomic signatures reveal new evidences for selection of important traits in domestic cattle., 2015, 32(3): 711–725.

[33] Edwards CE, Tessier BC, Swift JF, Bassüner B, Linan AG, Albrecht MA, Yatskievych GA. Conservation genetics of the threatened plant species(Missouri bladderpod) reveals strong genetic structure and a possible cryptic species., 2021, 16(3): e0247586.

[34] Liu F, Qu YK, Geng C, Wang AM, Zhang JH, Li JF, Chen KJ, Liu B, Tian HY, Yang WP, Yu YB. Analysis of the population structure and genetic diversity of the red swamp crayfish () in China using SSR markers., 2020, 47: 59–71.

[35] Reed DH, Frankham R. Correlation between fitness and genetic diversity., 2003, 17(1): 230–237.

[36] Eldridge MDB, King JM, Loupis AK, Spencer PBS, Taylor AC, Pope LC, Hall GP. Unprecedented low levels of genetic variation and inbreeding depression in an island population of the Black-Footed Rock-Wallaby., 1999, 13(3): 531–541.

[37] Gao YH, Gautier M, Ding XD, Zhang H, Wang YC, Wang X, Faruque MO, Li JY, Ye SH, Gou X, Han JL, Lenstra JA, Zhang Y. Species composition and environmental adap-tation of indigenous Chinese cattle., 2017, 7(1): 16196.

[38] Zhang WG, Gao X, Zhang Y, Zhao YM, Zhang JB, Jia YT, Zhu B, Xu LY, Zhang LP, Gao HJ, Li JY, Chen Y. Genome-wide assessment of genetic diversity and popu-lation structure insights into admixture and introgression in Chinese indigenous cattle., 2018, 19(1): 114.

[39] Yuan J, Xu GQ, Zhou X, Xu SP, Li S, Li WM, Liu B. SNP chip-based monitoring of population conservation effect of Tongcheng pigs., 2022, 53(08): 2514–2523.袁嬌, 徐國強, 周翔, 徐三平, 李勝, 黎望明, 劉榜. 基于SNP芯片監(jiān)測通城豬的保種效果. 畜牧獸醫(yī)學報, 2022, 53(08): 2514–2523.

[40] Botstein D, White RL, Skolnick M, Davis RW. Constru-ction of a genetic linkage map in man using restriction fragment length polymorphisms., 1980, 32(3): 314–331.

[41] Forutan M, Ansari Mahyari S, Baes C, Melzer N, Schenkel FS, Sargolzaei M. Inbreeding and runs of homozygo-sity before and after genomic selection in North American Holstein cattle., 2018, 19(1): 98.

[42] Keller MC, Visscher PM, Goddard M. Quantification of inbreeding due to distant ancestors and its detection using dense single nucleotide polymorphism data., 2011, 189: 237–249.

[43] Wang Y, Dong RL, Li X, Cui C, Yu GH. Analysis of the genetic diversity and family structure of the Licha Black pig population on Jiaodong Peninsula, Shandong Province, China., 2022, 12(8): 1045.

Analysis of genetic diversity and genetic structure of Qinchuan cattle conservation population using whole-genome resequencing

Jun Ma1, Anping Fan2, Wusheng Wang2, Jinchuan Zhang2, Xiaojun Jiang2, Ruijun Ma2, Sheqiang Jia2, Fei Liu2, Chuchao Lei1, Yongzhen Huang1

In the conservation of livestock and poultry resources, population genetic diversity and genetic structure of the conservation population are important factors affecting the effectiveness of conservation. In this study, whole-genome resequencing technology was used to detect genomic variation in 100 Qinchuan cattle (30 bulls and 70 cows). By analyzing population genetic diversity, runs of homozygosity (ROH) distribution features, kinship relationships, and family structure, the conservation effectiveness of Qinchuan cattle was comprehensively evaluated. The results showed that a total of 20,968,017 high-quality SNPs were detected in 100 Qinchuan cattle, the average minimum allele frequency was 0.191±0.124, the average polymorphic information content was 0.279±0.131, and the average observed heterozygosity was 0.275±0.131, the average expected heterozygosity is 0.279±0.131, indicating that the genetic diversity of the Qinchuan cattle conservation population is relatively rich. The average identity by state (IBS) distance of the Qinchuan conservation population was 0.243±0.020, with a value of 0.242±0.021 for the bulls. The results of the kinship G-matrix were consistent with the results of the IBS distance matrix, both showing that some individuals in the conservation population had close kinship. A total of 8258 genomic ROH were detected in 100 Qinchuan cattle, with a total length of 9.64 GB. The average length of ROH fragments was 1.167±1.203 Mb, 69.35% of the ROH were short ROH with a length of 0.5～1 Mb, and the average total length of ROH per individual was 96.40 Mb. The average inbreeding coefficient based on ROH was 0.039±0.039, with a value of 0.044±0.035 for the bulls, indicating that some bulls had a certain degree of inbreeding accumulation. The results of the phylogenetic tree combined with kinship analysis showed that the individuals in the Qinchuan cattle conservation population could be divided into eight families, including seven families with bulls and one family without bulls. This study demonstrated that the genetic diversity of the Qinchuan conservation population is relatively rich, with no significant inbreeding accumulation, but there is a risk of inbreeding among some individuals. Therefore, it is necessary to strengthen selection and mating to ensure the sustainable development of Qinchuan cattle resources.

Qinchuan cattle; whole-genome resequencing; genetic diversity; genetic structure; kinship

2023-04-25;

2023-06-11;

2023-06-30

陜西省重點研發(fā)計劃項目(編號：2023-YBNY-141)和財政部與農(nóng)業(yè)農(nóng)村部-國家現(xiàn)代農(nóng)業(yè)產(chǎn)業(yè)技術(shù)體系(編號：CARS-37)資助[Sup-portedby the Key R&D project of Shaanxi Province (2023-YBNY-141) and Ministry of Finance and Ministry of Agriculture and Rural Affairs-National Modern Agricultural Industry Technology System (CARS-37)]

馬鈞，在讀博士研究生，專業(yè)方向：動物遺傳育種與繁殖。E-mail: Junma96@163.com

黃永震，博士，副教授，研究方向：動物遺傳育種與繁殖。E-mail: hyzsci@nwafu.edu.cn

10.16288/j.yczz.23-115

(責任編委: 趙要風)