肅南牦牛群體遺傳結(jié)構(gòu)、選擇信號分析和ROH檢測

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

(1.中國農(nóng)業(yè)科學(xué)院蘭州畜牧與獸藥研究所，蘭州 730050；2.農(nóng)業(yè)農(nóng)村部青藏高原畜禽遺傳育種重點實驗室，蘭州 730050；3.甘肅省牦牛繁育工程重點實驗室，蘭州 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個地方品種以及2個培育品種，主要分布在中國西藏、青海、甘肅、四川等地[2]。肅南牦牛是青藏高原型牦牛的地方類群，主要分布于張掖市肅南裕固族自治縣，群體數(shù)量約6萬頭。肅南牦牛毛色以黑褐色為主，體型外貌上帶有野牦牛的特征，是肅南地區(qū)唯一在培的牦牛品種。從全基因組層面解析肅南牦牛的遺傳變異，有助于進一步了解這一品種牦牛的特性以及挖掘經(jīng)濟性狀相關(guān)的候選基因。

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

1 材料與方法

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

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

1.2 全基因組序列比對及變異檢測

使用Trimmomatic[16]對原始讀長(raw reads)進行質(zhì)控，去除適配器序列、污染和低質(zhì)量reads，將高質(zhì)量的reads與最新的牦牛參考基因組(GCA_005887515.1)比對得到SAM文件(Sequence Alignment Map)。使用SAMtools(v1.9)將SAM文件排序并轉(zhuǎn)換為二進制格式(BAM)[17]。對上游得到的每個樣本的BAM文件采用GATK v4.1.8.0進行SNP鑒定，采用MarkDuplicates、HaplotypeCaller和SelectVariants模塊獲取SNP的原始文件[18]。通過VariantFiltration模塊，設(shè)置“QUAL<30.0、QD<2.0、FS>60.0、MQ<40.0、HaplotypeScore>13.0”的過濾條件，對上游得到的初始VCF文件執(zhí)行硬過濾，最終得到含有高質(zhì)量SNP的VCF文件。采用samtools flagstat[17]命令計算每個樣本的統(tǒng)計信息，包括平均覆蓋率、原始讀長計數(shù)、映射讀長和正確配對的讀長比率?？紤]到基因型數(shù)據(jù)的低等位基因頻率以及性別因素可能會影響后續(xù)的分析，使用PLINK v1.9[19]剔除最小等位基因頻率(MAFs)<0.05、個體檢出率<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進行主成分分析(PCA)以確定群體間的遺傳結(jié)構(gòu)。PCA的可視化基于R包ggplot2完成。使用ADMIXTURE軟件[20]中最大似然估計方法進行祖先成分分析，最佳K值的確定參考Evanno等[21]的方法。采用PLINK軟件過濾強連鎖SNP位點(--indep-pairwise 50 10 0.3)，隨后構(gòu)建遺傳距離矩陣，使用腳本轉(zhuǎn)換為MEGA格式后采用鄰接法構(gòu)建系統(tǒng)進化樹。

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

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

1.5 連鎖不平衡分析

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

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

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

1.7 選擇信號分析和候選基因定位

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

2 結(jié) 果

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

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

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

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

A，5個牦牛群體的主成分分析；B，5個牦牛群體的親緣關(guān)系鑒定；C，祖先成分分析；D，系統(tǒng)進化樹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分析

為了評估不同牦牛群體的LD情況，采用PopLDdecay分別計算各個群體的成對r2值，用于衡量LD水平(圖2)。LD分析顯示，5個牦牛LD衰減的順序為：九龍牦牛>斯布牦牛>天祝白牦牛>巴州牦牛>肅南牦牛。

圖2 LD分析Fig.2 Linkage disequilibrium analysis

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

在5個牦牛品種中共檢測到8 426個ROHs，不同群體ROH數(shù)量和長度高度相關(guān)，介于0.98～0.99之間(圖3A)。按照ROH長度進行分類，可將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個>1 Mb的ROH(圖3B)。每條染色體都有ROH分布，高頻ROH主要位于22、16、4、2和9號染色體(圖3C)。所有品種中ROH平均長度最長為肅南牦牛(133.33 MB)，最短為天祝白牦牛(52.02 Mb)。ROH數(shù)量和高頻ROH數(shù)量排序為：肅南牦牛>九龍牦牛>天祝白牦牛>巴州牦牛>斯布牦牛；平均近交系數(shù)排序為：肅南牦牛>巴州/九龍牦牛>天祝白牦牛>斯布牦牛(表2)。

A,群體內(nèi)每頭個體ROH總數(shù)目與ROH總長度的散點圖；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)計

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

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

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

續(xù)表

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

2.6 選擇信號分析

為了解高頻ROH區(qū)域受到的選擇作用，采用基于單體型(haplotype)的iHS方法對5個牦牛群體共有的高頻ROH區(qū)域進行選擇信號檢測，結(jié)果顯示，候選基因主要定位在16號染色體3個高頻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)在每個品種基因組上都篩選出來(圖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個牦牛群體共有高頻ROH區(qū)域iHS分析Fig.5 iHS analysis in shared ROH islands of five yak populations

3 討 論

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

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

3.2 ROH分析

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

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

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

3.3 受選擇基因分析

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

4 結(jié) 論

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