鮑 麒，包鵬甲，馬曉明，吳曉云，孟光耀，褚 敏，曹紅梅，郎建英，安玉峰，梁春年，閻 萍

(1.中國農(nóng)業(yè)科學(xué)院蘭州畜牧與獸藥研究所，蘭州 730050；2.農(nóng)業(yè)農(nóng)村部青藏高原畜禽遺傳育種重點(diǎn)實(shí)驗(yàn)室，蘭州 730050；3.甘肅省牦牛繁育工程重點(diǎn)實(shí)驗(yàn)室，蘭州 730050；4.肅南裕固族自治縣畜牧獸醫(yī)服務(wù)中心，肅南 734400)

牦牛(Bosgrunniens)在生物學(xué)分類上屬于偶蹄目(Artiodactyla)?？?Bovidae)牛亞科(Bovinae)牛屬(Bos)，主要分布于青藏高原及其毗鄰高山、亞高山高寒地區(qū)[1]。目前，全世界現(xiàn)有牦牛約2 200萬頭，其中95%以上分布在中國[1]。據(jù)國家畜禽遺傳資源品種名錄(2021年版)，中國牦?，F(xiàn)有18個(gè)地方品種以及2個(gè)培育品種，主要分布在中國西藏、青海、甘肅、四川等地[2]。肅南牦牛是青藏高原型牦牛的地方類群，主要分布于張掖市肅南裕固族自治縣，群體數(shù)量約6萬頭。肅南牦牛毛色以黑褐色為主，體型外貌上帶有野牦牛的特征，是肅南地區(qū)唯一在培的牦牛品種。從全基因組層面解析肅南牦牛的遺傳變異，有助于進(jìn)一步了解這一品種牦牛的特性以及挖掘經(jīng)濟(jì)性狀相關(guān)的候選基因。

群體結(jié)構(gòu)的鑒定是群體基因組學(xué)的研究基礎(chǔ)，探索不同群體之間的系統(tǒng)發(fā)育關(guān)系、聚類分析結(jié)果和祖先成分組成，可以了解不同物種的分類地位和分化程度，在群體親緣關(guān)系鑒定、物種起源分析、種群歷史等研究中發(fā)揮作用[3-4]，為下游分析提供背景信息。選擇信號(hào)(selection signal)包括選擇性清除和搭車效應(yīng)，是挖掘與表型相關(guān)基因的重要手段[5]。隨著高通量測(cè)序技術(shù)成本的不斷降低，利用選擇信號(hào)解析物種之間性狀差異的研究日益普遍。劉真[6]利用遺傳分化指數(shù)(Fst)法對(duì)短脂尾的蒙古羊和瘦尾的藏羊群體進(jìn)行選擇信號(hào)分析，篩選出50個(gè)與脂類代謝相關(guān)的基因，如過氧化物酶體增生激活受體γ(PPARG)、視黃酸受體γ(RXRG)、溶質(zhì)載體家族27成員2(SLC27A2)、?；o酶A合成酶長鏈家族成員6(ACSL6)等。Shen等[7]利用跨群體復(fù)合似然比(XP-CLR)、Fst和θπ比值(θπratio)方法對(duì)延邊牛和非洲熱帶牛的重測(cè)序數(shù)據(jù)進(jìn)行選擇信號(hào)分析，篩選出3個(gè)與冷適應(yīng)性相關(guān)的候選基因(皮質(zhì)醇穩(wěn)定蛋白(CORT)、成纖維細(xì)胞生長因子5(FGF5)、脂肪酸合酶(CD36))，從基因組層面解析了延邊牛的耐寒機(jī)制。吳林慧等[8]在對(duì)恩施黑豬的選擇信號(hào)研究中發(fā)現(xiàn)一系列基因分別與胴體長度、精子形成、骨骼肌分化和總產(chǎn)仔數(shù)相關(guān)，包括溶血磷脂酸受體2(LPAR2)、NADH 脫氫酶1α復(fù)合體亞基(NDUFA13)、肌細(xì)胞增強(qiáng)因子2B(MEF2B)和芳烴受體基因(AHR)。連續(xù)純合片段(ROH)是指二倍體生物基因組中出現(xiàn)的連續(xù)純合區(qū)域[9]，能夠用來評(píng)估家畜近交系數(shù)、推斷種群歷史、輔助解析某些復(fù)雜性狀的形成機(jī)制[10]。Ferenakovi等[11]利用4個(gè)家牛(瑞士褐牛、德國弗萊維赫牛、挪威紅牛、蒂羅爾灰色牛)SNP數(shù)據(jù)，基于ROH長度和系譜估計(jì)分別計(jì)算近交系數(shù)FROH與FPED，結(jié)果發(fā)現(xiàn)ROH長度>8和16 MB時(shí)估計(jì)的FROH和FPED接近，且利用FROH評(píng)估各牛種之間近交水平的結(jié)果與FPED評(píng)估的結(jié)果一致。Wu等[12]在對(duì)20頭皖南黑豬ROH分析中鑒定到1 913個(gè)ROHs，利用ROH長度計(jì)算的平均基因組近交系數(shù)為0.5234，呈現(xiàn)高度近交的狀態(tài)，揭示了皖南黑豬目前出現(xiàn)繁殖性能衰退的原因。莫家遠(yuǎn)等[13]對(duì)包括東山豬在內(nèi)的4個(gè)廣西地方豬種進(jìn)行ROH檢測(cè)，注釋的候選基因主要富集在與免疫、繁殖和肉品質(zhì)相關(guān)的通路。Xu等[14]對(duì)金華豬群體篩選出的高頻ROH區(qū)域進(jìn)行注釋，鑒定到與生殖相關(guān)的基因(同源框A3(HOXA3)、HOXA7、HOXA10和HOXA11)，肉質(zhì)相關(guān)基因(成肌細(xì)胞決定基因1(MYOD1)、脂素家族蛋白3(LPIN3)和β-連環(huán)蛋白樣蛋白1(CTNNBL1))、食欲相關(guān)基因(核連蛋白2(NUCB2))和抗病力相關(guān)基因(黏蛋白抗原4(MUC4)、MUC13、MUC20、Leishmanolysin-like肽酶(LMLN)、整合素蛋白β5(ITGB5)、心臟發(fā)育蛋白1(HEG1)、溶質(zhì)載體家族12成員8(SLC12A8)和肌球蛋白輕鏈激酶(MYLK))。

與基因芯片技術(shù)相比，全基因組測(cè)序檢測(cè)精度高，對(duì)基因組完整度依賴性小，能發(fā)現(xiàn)更多單核苷酸多態(tài)性(SNP)或稀有位點(diǎn)，因此能找到更短的ROH片段，有利于檢測(cè)親緣關(guān)系較遠(yuǎn)的物種之間的雜合度，避免出現(xiàn)ROH假陽性情況[15]。目前，關(guān)于肅南牦牛在全基因組范圍內(nèi)的群體遺傳結(jié)構(gòu)解析以及功能基因挖掘的研究鮮見報(bào)道?；诖耍狙芯繉?duì)肅南牦牛以及巴州牦牛、九龍牦牛、斯布牦牛、天祝白牦牛4個(gè)著名牦牛品種進(jìn)行10×全基因組重測(cè)序，對(duì)其進(jìn)行遺傳結(jié)構(gòu)、ROH以及選擇信號(hào)分析，比較不同牦牛之間的異同，以期為肅南牦牛的育種改良工作提供基礎(chǔ)信息。

1 材料與方法

1.1 樣本采集、DNA提取與全基因組重測(cè)序

本研究選擇肅南牦牛(Sunan)、巴州牦牛(Bazhou)、九龍牦牛(Jiulong)、斯布牦牛(Sibu)和天祝白牦牛(Tianzhu)5個(gè)牦牛品種進(jìn)行全基因組重測(cè)序分析。樣本分別采集于甘肅省肅南裕固族自治縣、新疆維吾爾自治區(qū)巴音郭楞蒙古自治州、四川省九龍縣、西藏自治區(qū)墨竹工卡縣和甘肅省天祝藏族自治縣，樣本數(shù)分別為8、8、12、10和10頭。所有個(gè)體采集頸靜脈血液10 mL于EDTA抗凝管中，置于-20 ℃待測(cè)。使用EasyPure血液基因組DNA提取試劑盒(北京全式金生物技術(shù)有限公司)提取血樣DNA，通過測(cè)定A260 nm/A280 nm比值和1.0%瓊脂糖凝膠電泳篩選，檢測(cè)提取DNA的質(zhì)量和完整性。用Covaris超聲波DNA破碎儀將合格的基因組DNA樣本隨機(jī)分成片段。對(duì)片段DNA進(jìn)行末端修復(fù)，在3′-端加入堿基A連接測(cè)序接頭，吸附富集后進(jìn)行PCR擴(kuò)增形成測(cè)序文庫。構(gòu)建好的文庫在BGISEQ-500平臺(tái)進(jìn)行測(cè)序。測(cè)序得到的原始圖像數(shù)據(jù)通過BGISEQ-500 Base Calling軟件轉(zhuǎn)換為原始數(shù)據(jù)，以FASTQ文件格式存儲(chǔ)。

1.2 全基因組序列比對(duì)及變異檢測(cè)

使用Trimmomatic[16]對(duì)原始讀長(raw reads)進(jìn)行質(zhì)控，去除適配器序列、污染和低質(zhì)量reads，將高質(zhì)量的reads與最新的牦牛參考基因組(GCA_005887515.1)比對(duì)得到SAM文件(Sequence Alignment Map)。使用SAMtools(v1.9)將SAM文件排序并轉(zhuǎn)換為二進(jìn)制格式(BAM)[17]。對(duì)上游得到的每個(gè)樣本的BAM文件采用GATK v4.1.8.0進(jìn)行SNP鑒定，采用MarkDuplicates、HaplotypeCaller和SelectVariants模塊獲取SNP的原始文件[18]。通過VariantFiltration模塊，設(shè)置“QUAL<30.0、QD<2.0、FS>60.0、MQ<40.0、HaplotypeScore>13.0”的過濾條件，對(duì)上游得到的初始VCF文件執(zhí)行硬過濾，最終得到含有高質(zhì)量SNP的VCF文件。采用samtools flagstat[17]命令計(jì)算每個(gè)樣本的統(tǒng)計(jì)信息，包括平均覆蓋率、原始讀長計(jì)數(shù)、映射讀長和正確配對(duì)的讀長比率?？紤]到基因型數(shù)據(jù)的低等位基因頻率以及性別因素可能會(huì)影響后續(xù)的分析，使用PLINK v1.9[19]剔除最小等位基因頻率(MAFs)<0.05、個(gè)體檢出率<0.95、SNP檢出率<0.99、無染色體分配和性染色體上的SNP，參數(shù)設(shè)置為“--maf 0.05 --mind 0.05 --geno 0.01 --chr 1-29”。

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

基于SNP信息，采用PLINK進(jìn)行主成分分析(PCA)以確定群體間的遺傳結(jié)構(gòu)。PCA的可視化基于R包ggplot2完成。使用ADMIXTURE軟件[20]中最大似然估計(jì)方法進(jìn)行祖先成分分析，最佳K值的確定參考Evanno等[21]的方法。采用PLINK軟件過濾強(qiáng)連鎖SNP位點(diǎn)(--indep-pairwise 50 10 0.3)，隨后構(gòu)建遺傳距離矩陣，使用腳本轉(zhuǎn)換為MEGA格式后采用鄰接法構(gòu)建系統(tǒng)進(jìn)化樹。

1.4 親緣關(guān)系系數(shù)計(jì)算

為了消除親緣關(guān)系較近的樣品對(duì)下游分析的影響，采用KING軟件[22]計(jì)算兩兩個(gè)體之間的親緣關(guān)系系數(shù)，剔除存在三代以內(nèi)親緣關(guān)系的個(gè)體，即關(guān)系判斷閾值>0.0442的樣本。親緣關(guān)系矩陣可視化使用R包c(diǎn)orrplot完成。

1.5 連鎖不平衡分析

采用PopLDdecay軟件[23]計(jì)算每個(gè)品種的連鎖不平衡(LD)程度，并使用軟件自帶的Plot_MultiPop.pl腳本繪制LD衰減曲線圖。

1.6 ROH分析及近交系數(shù)計(jì)算

每個(gè)個(gè)體的ROH使用PLINK軟件計(jì)算，參數(shù)設(shè)置為“--homozyg --homozyg-density 50 --homozyg-gap 100 --homozyg-kb 300 --homozyg-snp 50 --homozyg-window-het 3 --homozyg-window-snp 50 --homozyg-window-threshold 0.05”。對(duì)得到的ROH數(shù)量和長度進(jìn)行統(tǒng)計(jì)和分類，使用GraphPad Prism 8.0軟件進(jìn)行展示。參考Mcquillan等[24]方法計(jì)算近交系數(shù)(FROH)，計(jì)算公式為：FROH=常染色體上ROH片段長度之和/常染色體總長度之和。將各自群體內(nèi)檢測(cè)到ROH的SNP進(jìn)行頻次統(tǒng)計(jì)，以頻率Top 5%的SNP定義為高頻SNP，將其構(gòu)成的ROH區(qū)域作為高頻ROH區(qū)域。使用snpEff和BEDtools對(duì)各品種高頻ROH區(qū)域進(jìn)行基因注釋，并對(duì)定位到的基因及重疊基因分別進(jìn)行GO功能和KEGG信號(hào)通路富集分析。

1.7 選擇信號(hào)分析和候選基因定位

ROH的形成受到自然選擇和人工選擇的影響，檢測(cè)ROH片段中位點(diǎn)的選擇強(qiáng)度可能會(huì)鑒定到與表型相關(guān)的關(guān)鍵基因[12]。采用rehh包[25]分別計(jì)算各高頻ROH區(qū)域內(nèi)SNP的綜合單倍型評(píng)分(iHS)，將iHS評(píng)分Top 1%的位點(diǎn)定義為強(qiáng)受選擇位點(diǎn)。針對(duì)iHS法得到的受選擇位點(diǎn)，使用1.6的方法對(duì)其進(jìn)行注釋篩選出候選基因。為了縮短計(jì)算時(shí)間，使用R包parallel進(jìn)行并行計(jì)算。

2 結(jié) 果

2.1 5個(gè)牦牛群體變異檢測(cè)結(jié)果

本研究對(duì)巴州牦牛、九龍牦牛、斯布牦牛、肅南牦牛以及天祝白牦牛共48個(gè)樣本進(jìn)行重測(cè)序，共獲得816 G高質(zhì)量數(shù)據(jù)，平均每個(gè)樣本獲得了17 G的數(shù)據(jù)，平均測(cè)序深度為7.14×。將高質(zhì)量測(cè)序數(shù)據(jù)比對(duì)到參考基因組之后，平均比對(duì)率為98.864%。為了進(jìn)行后續(xù)群體遺傳分析，對(duì)5個(gè)群體一同進(jìn)行SNP變異檢測(cè)，經(jīng)過濾后得到15 092 883個(gè)SNP位點(diǎn)。此外，5個(gè)群體SNP的注釋結(jié)果顯示，各牦牛群體中的SNP位點(diǎn)分布相似，主要都位于內(nèi)含子區(qū)域，其次是基因間區(qū)、基因下游區(qū)域、基因上游區(qū)域、外顯子區(qū)域，平均占比分別為51.47%、37.96%、4.57%、4.49%、1.09%(表1)。

2.2 牦牛群體結(jié)構(gòu)分析

對(duì)5個(gè)地方牦牛群體進(jìn)行主成分分析、親緣關(guān)系分析、群體祖先成分分析和系統(tǒng)進(jìn)化樹分析的結(jié)果見圖1。由圖1可知，5個(gè)牦牛群體按照第一主成分(2.5397%)和第二主成分可以區(qū)分開(1.8990%)，說明5個(gè)群體間沒有雜交現(xiàn)象，與品種所處的地理位置相符合。除肅南牦牛群體內(nèi)較為聚集外，其他4個(gè)牦牛品種群體內(nèi)較分散(圖1A)。親緣關(guān)系結(jié)果顯示，5個(gè)牦牛群體內(nèi)親緣關(guān)系較低，關(guān)系判斷閾值均<0，不存在三代以內(nèi)親緣關(guān)系的個(gè)體，滿足后續(xù)分析的要求(<0.0442)(圖1B)。祖先成分分析顯示，最佳分群數(shù)為4群；K=2時(shí)，巴州牦牛首先被鑒定出來，肅南牦牛與天祝白牦牛祖先成分相似；K=3時(shí)，巴州牦牛、九龍牦牛、天祝白牦?？梢悦黠@區(qū)分；K=4時(shí)，肅南牦牛擁有其他4個(gè)牦牛祖先成分，且天祝白牦牛祖先成分所占比例最高(圖1C)。系統(tǒng)進(jìn)化樹顯示，巴州牦牛、九龍牦牛、斯布牦牛、肅南牦牛和天祝白牦牛群體各自聚集在一個(gè)分支上，肅南牦牛聚集在進(jìn)化樹一端，其他4個(gè)牦牛品種集中在進(jìn)化樹另一端(圖1D)。

A，5個(gè)牦牛群體的主成分分析；B，5個(gè)牦牛群體的親緣關(guān)系鑒定；C，祖先成分分析；D，系統(tǒng)進(jìn)化樹A,Principal component analysis of five yak populations;B,Genetic relationship identification of five yak populations;C,Ancestral component analysis;D,Phylogenetic tree圖1 群體基因組學(xué)分析Fig.1 Population genomics analysis

2.3 LD分析

為了評(píng)估不同牦牛群體的LD情況，采用PopLDdecay分別計(jì)算各個(gè)群體的成對(duì)r2值，用于衡量LD水平(圖2)。LD分析顯示，5個(gè)牦牛LD衰減的順序?yàn)椋壕琵堦笈?斯布牦牛>天祝白牦牛>巴州牦牛>肅南牦牛。

圖2 LD分析Fig.2 Linkage disequilibrium analysis

2.4 ROH基本統(tǒng)計(jì)及基因組近交系數(shù)結(jié)果

在5個(gè)牦牛品種中共檢測(cè)到8 426個(gè)ROHs，不同群體ROH數(shù)量和長度高度相關(guān)，介于0.98～0.99之間(圖3A)。按照ROH長度進(jìn)行分類，可將ROH分為3類(0～0.5 Mb；0.5～1 Mb；>1 Mb)，ROH數(shù)目在0～0.5 Mb大小的占絕大部分，其次是0.5～1 Mb的ROH。肅南牦牛、天祝白牦牛、巴州牦牛、九龍牦牛和斯布牦牛分別存在26、4、3、1和1個(gè)>1 Mb的ROH(圖3B)。每條染色體都有ROH分布，高頻ROH主要位于22、16、4、2和9號(hào)染色體(圖3C)。所有品種中ROH平均長度最長為肅南牦牛(133.33 MB)，最短為天祝白牦牛(52.02 Mb)。ROH數(shù)量和高頻ROH數(shù)量排序?yàn)椋好C南牦牛>九龍牦牛>天祝白牦牛>巴州牦牛>斯布牦牛；平均近交系數(shù)排序?yàn)椋好C南牦牛>巴州/九龍牦牛>天祝白牦牛>斯布牦牛(表2)。

A,群體內(nèi)每頭個(gè)體ROH總數(shù)目與ROH總長度的散點(diǎn)圖；B,短、中和長片段ROH數(shù)量直方圖；C,不同染色體ROH分布直方圖A,Scatter plot of the total number of ROH and the total length of each individual in the population;B,Histogram of the number of short,medium and long ROH;C,Histogram of the number of ROH on different chromosomes圖3 ROH分析Fig.3 ROH analysis

表2 ROH和近交系數(shù)統(tǒng)計(jì)

2.5 高頻ROH區(qū)域基因注釋與富集分析

對(duì)高頻區(qū)域內(nèi)SNP進(jìn)行注釋，其中巴州牦牛、九龍牦牛、斯布牦牛、肅南牦牛和天祝白牦牛分別定位到58、73、52、465和35個(gè)基因，其中包含核內(nèi)不均一性核糖核蛋白K(HNRNPK)、驅(qū)動(dòng)蛋白27(KIF27)、G激酶錨定蛋白1(GKAP1)等3個(gè)重疊基因，均位于7號(hào)染色體。對(duì)肅南牦牛篩選出的候選基因進(jìn)行富集分析，以Q-value<0.05作為篩選條件，結(jié)果顯示肅南牦牛富集到17條GO條目(RNA聚合酶對(duì)轉(zhuǎn)錄的負(fù)調(diào)控Ⅱ、ATP結(jié)合、胚胎骨骼系統(tǒng)形態(tài)發(fā)生、脂肪氧合酶通路等)和2條KEGG信號(hào)通路(Notch信號(hào)通路和Hippo信號(hào)通路)(表3)。此外，對(duì)3個(gè)重疊基因的富集分析結(jié)果顯示，這些基因參與了16個(gè)GO條目和3條KEGG信號(hào)通路，包括基因表達(dá)調(diào)控、細(xì)胞連接、胰島素受體信號(hào)通路的正調(diào)控、體元投射通路等(圖4)。

表3 肅南牦牛候選基因的富集分析結(jié)果

續(xù)表

圖4 重疊基因富集分析結(jié)果Fig.4 Enrichment analysis results of overlap gene

2.6 選擇信號(hào)分析

為了解高頻ROH區(qū)域受到的選擇作用，采用基于單體型(haplotype)的iHS方法對(duì)5個(gè)牦牛群體共有的高頻ROH區(qū)域進(jìn)行選擇信號(hào)檢測(cè)，結(jié)果顯示，候選基因主要定位在16號(hào)染色體3個(gè)高頻ROH區(qū)域上(Chr16：64 780 362-65 091 013；Chr16：75 041 712-75 388 110；Chr16：68 520 840-68 824 137)。N-α乙酰轉(zhuǎn)移酶25(NAA25)在每個(gè)品種基因組上都篩選出來(圖5)。巴州牦?；蚪M上篩選到內(nèi)質(zhì)網(wǎng)分子伴侶29(ERP29)、跨膜蛋白16(TMEM116)、TRAF型鋅指結(jié)構(gòu)域蛋白1(TRAFD1)(圖5A)，斯布和肅南牦?；蚪M上篩選到ERP29、TMEM116(圖5C、5D)。天祝白牦牛基因組上除了NAA25和TRAFD1外，還篩選到HECT結(jié)構(gòu)域E3泛素連接酶(HECTD4)(圖5E)。

A～E，分別代表巴州牦牛、九龍牦牛、斯布牦牛、肅南牦牛和天祝白牦牛共享高頻ROH區(qū)域iHS分析A-E，The iHS analysis in shared ROH islands of Bazhou,Jiulong,Sibu,Sunan and Tianzhu White yak,respectively圖5 5個(gè)牦牛群體共有高頻ROH區(qū)域iHS分析Fig.5 iHS analysis in shared ROH islands of five yak populations

3 討 論

3.1 5個(gè)牦牛品種群體結(jié)構(gòu)分析

肅南牦牛是肅南地區(qū)在培的牦牛品種，其與國內(nèi)其他牦牛品種之間的分化程度、分類地位以及進(jìn)化歷史尚不清楚，群體基因組學(xué)研究有助于在基因組水平解決這個(gè)問題。相比于其他牦牛群體，肅南牦牛群體內(nèi)部存在最小的變異幅度，可能是因?yàn)槊C南牦牛由于種質(zhì)資源開發(fā)時(shí)間較晚與地理環(huán)境原因，因此其變異幅度較小。肅南牦牛與其他牦牛群體在基因組水平上差異較大，其在系統(tǒng)發(fā)育樹分支長度遠(yuǎn)長于其他4個(gè)牦牛品種。祖先成分分析能夠推斷每個(gè)材料的基因組變異來源于假設(shè)亞群的可能性[20]。肅南牦牛包含了其他4種牦牛群體祖先血統(tǒng)，這可能與肅南牦牛培育及改良過程中引入其他牦牛品種優(yōu)良個(gè)體有關(guān)，這也與前期文獻(xiàn)報(bào)道的一致[26]。在各個(gè)成分中，其包含天祝白牦牛祖先成分的比例最大，表明肅南牦牛群體與天祝白牦牛群體存在廣泛交流，這也與二者所處的地理位置較近相符合[26-27]。根據(jù)變異程度、進(jìn)化距離和祖先成分推斷，肅南牦牛可能是在長期外部交流下產(chǎn)生的一個(gè)獨(dú)特類群，與地理距離最近的天祝白牦牛交流廣泛。

3.2 ROH分析

家畜經(jīng)過較強(qiáng)人工選育后，會(huì)導(dǎo)致特定區(qū)間純和度的增加，使得群體內(nèi)ROH富集的頻率增加，基因組多樣性也隨之降低[28]。對(duì)5個(gè)品種的ROH分類顯示，牦牛ROH片段普遍較短(<2 Mb)，統(tǒng)計(jì)發(fā)現(xiàn)0～0.5 Mb的ROH數(shù)目最多，其次是0.5～1 Mb片段，說明小片段ROH對(duì)整個(gè)基因組的近交貢獻(xiàn)最大。相比于其他4個(gè)牦牛品種，肅南牦牛的ROH片段總長度和數(shù)目最多，其原因可能是肅南牦牛經(jīng)過高度的近交選育，核苷酸多樣性最低。近交會(huì)增加純和等位基因頻率，從而增加ROH的長度和數(shù)量，因此ROH也常被用來評(píng)估家畜群體的近交配程度[10]。本研究計(jì)算結(jié)果顯示肅南牦牛平均近交系數(shù)要高于其余4個(gè)牦牛群體，其原因可能與其育種模式有關(guān)，即使用有限數(shù)量的高品質(zhì)公牛作為父本進(jìn)行人工授精，造成群體的近交水平的顯著增加。由于LD狀態(tài)受到家畜近交水平的影響，馴化程度越高，選擇強(qiáng)度越大的群體，LD衰減越慢[10]。肅南牦牛和巴州牦牛相對(duì)其他3種牦牛具有較高的連鎖程度，這可能與其受到較強(qiáng)的選擇強(qiáng)度有關(guān)的，這同樣解釋了肅南牦牛群體具有較高近交水平的結(jié)果。

肅南牦牛高頻ROH區(qū)域注釋到465個(gè)候選基因，富集到22條顯著通路，其中胚胎骨骼系統(tǒng)形態(tài)發(fā)生、前后模式規(guī)范、甲狀腺發(fā)育通路與機(jī)體發(fā)育和大腦形態(tài)形成有關(guān)，且參與3條通路的主要為同源異形盒(Homeobox，HOX)基因家族成員(HOXA3、HOXA5、HOXD3)。HOX蛋白是一種轉(zhuǎn)錄因子，在胚胎發(fā)育過程參與組織分化過程[29]。此外研究還發(fā)現(xiàn)HOX家族成員與毛囊形成和毛發(fā)生長密切相關(guān)，參與調(diào)控角蛋白和角蛋白關(guān)聯(lián)蛋白的表達(dá)，維持毛囊正常形態(tài)等[30]。甲狀腺發(fā)育通路包含F(xiàn)GF8基因，該基因與中腦發(fā)育有關(guān)，在肢體發(fā)育的誘導(dǎo)、啟動(dòng)和維持中起作用[31-32]，且有研究顯示FGF8在小鼠表皮中過表達(dá)會(huì)抑制毛囊的發(fā)育[33]。HOXA3、HOXA5、HOXD3和FGF8這4個(gè)基因可能參與調(diào)控肅南牦牛毛發(fā)生長過程，值得進(jìn)一步研究。脂肪氧合酶通路包含4個(gè)基因，其中3個(gè)為脂氧合酶(LOX)基因家族成員編碼基因(ALOX12B、ALOX15B、ALOXE3)，該家族主要參與氧化不飽和脂肪酸如花生四烯酸等特定的過氧化物[34]。研究顯示添加花生四烯酸會(huì)影響動(dòng)物肉類的適口性[35]，該類基因可能與肅南牦牛肉品質(zhì)有關(guān)。參與顯著通路的基因主要以基因家族的形式出現(xiàn)，這一特點(diǎn)可能與肅南牦牛品種培育過程中特定功能的強(qiáng)化有關(guān)。另外還富集到一條Hippo通路，該通路最初在果蠅中被發(fā)現(xiàn)，被認(rèn)為是一種保守的信號(hào)通路，在控制器官大小、限制細(xì)胞增殖和誘導(dǎo)細(xì)胞凋亡方面發(fā)揮著重要作用[36]，可能與肅南牦牛機(jī)體發(fā)育有關(guān)。

對(duì)各高頻ROH區(qū)域篩選的候選基因取交集，發(fā)現(xiàn)HNRNPK、KIF27、GKAP1 3個(gè)基因在5個(gè)牦牛品種中是共享的，均位于7號(hào)染色體1個(gè)高頻ROH區(qū)域中(Chr7：12 661 870-13 045 935)。對(duì)共享基因的富集分析結(jié)果顯示，3個(gè)基因參與了胰島素受體信號(hào)通路的正調(diào)控通路。Bao等[37]研究表明胰島素相關(guān)通路是牦牛能夠適應(yīng)高原環(huán)境的主要調(diào)節(jié)機(jī)制之一，胰島素的動(dòng)態(tài)調(diào)節(jié)有助于高原動(dòng)物更好地利用能量，抵御惡劣的自然環(huán)境。參與腦室系統(tǒng)開發(fā)通路的KIF27是哺乳動(dòng)物中COS2的同源基因之一，其功能是作為Hedgehog信號(hào)通路的調(diào)節(jié)因子，參與胚胎期發(fā)育、組織修復(fù)、器官形成等過程[38]，也有小鼠上的研究證明該通路在低氧條件下肺血管重建過程中發(fā)揮作用[39]。HNRNPK是一種對(duì)細(xì)胞周期有廣泛影響的RNA結(jié)合蛋白。小鼠研究中發(fā)現(xiàn)，血紅素氧合酶的破壞可能增強(qiáng)hnRNPK介導(dǎo)的蛋白翻譯抑制，進(jìn)而削弱β-catenin/hnRNPK調(diào)控的基因表達(dá)，這是協(xié)同肺修復(fù)和再生所必需的[40]。相比平原地區(qū)的牛，牦牛生活在高原低氧極端環(huán)境中，其肺部結(jié)構(gòu)發(fā)生了一系列適應(yīng)性變化，如牦牛的囊狀期較短，肺泡期較長，肺臟發(fā)育時(shí)間較早等[41]，推測(cè)hnRNPK基因與牦牛低氧氣濃度影響下的肺功能修復(fù)過程有關(guān)。本研究鑒定的3個(gè)共享基因可能通過胰島素代謝、肺血管重建和肺功能修復(fù)途徑參與高原適應(yīng)性的調(diào)節(jié)，這一結(jié)果有助于解析牦牛高原適應(yīng)性進(jìn)化的分子機(jī)制和遺傳原理。

3.3 受選擇基因分析

iHS統(tǒng)計(jì)量可以靈敏地檢測(cè)正選擇信號(hào)及近期選擇，可用來揭示選擇信號(hào)背后的功能基因[42]。本研究用iHS方法對(duì)共有高頻ROH區(qū)域內(nèi)的位點(diǎn)進(jìn)行選擇信號(hào)分析，發(fā)現(xiàn)5個(gè)牦牛16號(hào)染色體上3個(gè)高頻ROH區(qū)域(Chr16：64 780 362-65 091 013；Chr16：75 041 712-75 388 110；Chr16：68 520 840-68 824 137)內(nèi)部分SNP受到強(qiáng)烈選擇，這3個(gè)相鄰區(qū)段可能是5個(gè)牦牛品種在染色體序列發(fā)生不斷重組的進(jìn)化過程中保留下來的受選擇熱點(diǎn)區(qū)域。注釋結(jié)果顯示，該區(qū)段得到5個(gè)不同的基因(NAA25、ERP29、TMEM116、TRAFD1、HECTD4)，這些基因部分或全部牦牛品種中檢測(cè)到，可能與滑動(dòng)窗口和步長的大小設(shè)置有關(guān)。NAA25在每個(gè)品種中都鑒定出，該基因編碼N-端乙酰轉(zhuǎn)移酶B復(fù)合體(NatB)的輔助亞基，能夠乙酰化蛋氨酸殘基，可能在正常細(xì)胞周期進(jìn)程中起作用[43]。除九龍和天祝白牦牛外，ERP29在其他牦牛都鑒定出，該基因編碼蛋白是內(nèi)質(zhì)網(wǎng)腔內(nèi)的一種駐留蛋白，參與內(nèi)質(zhì)網(wǎng)分泌蛋白的加工過程[44]。TMEM116在巴州牦牛、斯布牦牛和肅南牦牛上鑒定出，該基因編碼蛋白屬于跨膜蛋白(TMEM)家族，與細(xì)胞間、細(xì)胞內(nèi)的信號(hào)傳導(dǎo)、免疫相關(guān)疾病相關(guān)[45]。TRAFD1在巴州牦牛和天祝白牦牛中鑒定出來，該基因負(fù)向調(diào)控人類和哺乳動(dòng)物的Toll樣受體(TLR)和RIG-Ⅰ樣受體(RLR)信號(hào)轉(zhuǎn)導(dǎo)[46]。TLR介導(dǎo)的信號(hào)通路與免疫反應(yīng)中病原體的消除、獲得性免疫體系的建立有關(guān)[47]。RLR是識(shí)別胞漿中病毒RNA的主要受體，在抵抗病毒感染的過程起關(guān)鍵作用[48]。TMEM116、TRAFD1受到強(qiáng)烈選擇可能與牦牛品種良好的抗病性有關(guān)。在天祝白牦牛中鑒定到HECTD4，該基因在人上可能編碼E3泛素蛋白連接酶[49]。其在牦牛中的作用尚不清楚，需要進(jìn)一步的驗(yàn)證。本研究受選擇分析中顯著位點(diǎn)集中在16號(hào)染色體3個(gè)相鄰區(qū)段上，推測(cè)是由于牦牛之間分化程度較低，牦牛共有高頻ROH區(qū)段內(nèi)位點(diǎn)受到較為一致的選擇壓力。篩選到受選擇基因主要與抗病性、內(nèi)質(zhì)網(wǎng)分泌蛋白加工、以及細(xì)胞周期調(diào)節(jié)相關(guān)，這一結(jié)果有助于進(jìn)一步了解不同牦牛類群整體的進(jìn)化方向。

4 結(jié) 論

研究結(jié)果表明，肅南牦牛是在長期外部交流下產(chǎn)生的一個(gè)獨(dú)特類群，受天祝白牦牛血統(tǒng)影響較大。相比于其他4個(gè)成熟牦牛品種，具有較高的近交水平，可能與品種培育過程受到的強(qiáng)化選擇有關(guān)。在肅南牦牛高頻ROH區(qū)域中鑒定到與胚胎發(fā)育與組織分化過程(HOXA3、HOXA5、HOXD3)，肉品質(zhì)(ALOX12B、ALOX15B、ALOXE3)和中腦發(fā)育(FGF8)相關(guān)的基因。此外，與其他4個(gè)牦牛品種共享3個(gè)與高原適應(yīng)性相關(guān)的基因hnRNPK、KIF27、GKAP1。在共有高頻ROH區(qū)域鑒定到的受選擇基因主要與抗病性、內(nèi)質(zhì)網(wǎng)分泌蛋白加工及細(xì)胞周期調(diào)節(jié)相關(guān)，包括NAA25、ERP29、TMEM116、TRAFD1、HECTD4基因。本研究結(jié)果有助于了解肅南牦牛的品種特性，為正在進(jìn)行的基因組選擇育種和下一步育種計(jì)劃的制定提供參考。