基于ovine SNP 50 K芯片對策勒黑羊繁殖和抗逆性狀的遺傳性分析

鄭浪漫 張璐璐 王若桐 周 穩 韓志鵬 劉春潔 王彩軍 張成龍 宮昌海 劉書東*

(1.塔里木大學 動物科學與技術學院，新疆維吾爾自治區 阿拉爾 843300;2.新疆生產建設兵團塔里木畜牧科技重點實驗室，新疆維吾爾自治區 阿拉爾 843300;3.新疆巴州畜牧工作站，新疆維吾爾自治區 巴州 841205)

策勒黑羊主產于塔克拉瑪干沙漠南緣、昆侖山北坡，長期受環境因素如干旱少雨、紫外輻射強等的影響很大。因策勒黑羊羔皮制作的衣物深受當地人們喜歡，經當地農牧民長期精心選育，具有繁殖率高、性成熟早、常年發情、抗逆性強等特點[1]。這些特點受基因和環境共同調控，其遺傳規律由分子標記的顯隱性和互作等共同作用，使得基因或分子標記間存在的變異或固定在環境的影響下得以呈現。經過人工選擇和自然選擇后的策勒黑羊群體所存在的基因組中留下具有繁殖和抗逆性狀的印跡。這些基因印跡是揭示策勒黑羊遺傳規律的基礎，也是解析該綿羊品種特異分子標記的必要條件。

在群體中優異基因可以通過基因組選擇的方法進行篩選，還可采用候選分子標記法篩選，如PCR、核酸雜交和免疫篩選等[2]。伴隨著基因組測序技術的快速發展，生物復雜性狀相關分子標記和候選基因的篩選更加便利[3]。Ovine SNP 50 K bead chip由Illumina公司根據綿羊基因組大小設計，可以均勻覆蓋羊全基因組[4]，涵蓋了如羊角、體重、肉質、繁殖和疾病等重要性狀位點，適用于基因的篩選和定位、基因組育種、品種鑒定等。影響綿羊繁殖、抗逆性狀的基因有很多，但是策勒黑羊關于繁殖和抗逆性狀的相關研究尚還不夠深入，需要大量挖掘與繁殖、抗逆性狀相關分子標記?；诓呃蘸谘蚍敝?、抗逆性狀研究現狀和Illumina ovine SNP 50 K bead chip的優點，本研究選用Illumina ovine SNP 50 K bead chip挖掘策勒黑羊優異基因。

哈代-溫伯格(Hardy-Weinberg)定律是指在理想狀態下，各等位基因的頻率在遺傳中是穩定不變的，即保持著基因平衡。如果生物種群的基因型和等位基因頻率在世代中保持不變，則生物種群處于哈代-溫伯格平衡狀態[5]。而經選擇后的群體，基因或分子標記受到影響[6]。哈代-溫伯格平衡(HWE)遺傳變異的統計分析是遺傳數據集分析的重要組成部分，精確的測試程序是測試雙等位基因變異最先進技術，使用快速遞歸程序，在全基因組范圍內對HWE的雙等位基因變異進行精確測試[7]。在綿羊產仔數相關研究中，通過3種基因組方法(GWAS，XP-NSL和ROH)發現GTF2A1與綿羊的繁殖性狀相關[8]。GNAQ調節GnRH的表達和分泌，影響哈薩克羊季節性發情[9]。針對綿羊抗逆性相關研究中，TMEM154與綿羊慢病毒抗性相關[10-11]。

為探析策勒黑羊繁殖和抗逆性狀相關基因，本研究選擇100只無親緣關系的策勒黑羊，基于Illumina ovine SNP 50 K bead chip，通過哈代-溫伯格值分析，篩選策勒黑羊繁殖和抗逆性相關的分子標記，解析其分子遺傳機理。該研究為策勒黑羊優異分子標記挖掘提供參考，同時也為該綿羊種質資源保種和改良提供數據基礎。

1 材料與方法

1.1 試驗樣本

本研究于和田地區策勒縣昆侖綠源羊管家種畜場隨機選取無親緣關系的策勒黑羊100只，同一營養水平舍飼，健康成年母羊。頸靜脈采血，采樣時間為2020年。

1.2 DNA的提取、SNP分型及質量控制

天根試劑盒提取DNA，將合格的DNA樣品送新疆康普森生物技術有限公司，公司基于Infinium HD全基因基因型分型系統實現基因組DNA分型。Infinium HD全基因組基因型分型系統包括商業化綿羊芯片(Illumina ovine SNP 50 K bead chip)制備、iScan芯片掃描儀掃描和GenomeStudio軟件對原始信號文件標準均一化轉換成可讀的基因型數據并導出保存[12]。從公司獲得SNP分型數據。使用Plink 1.90軟件[13]對SNP分型數據進行質量控制：剔除檢出率小于90%、性染色體、最小等位基因頻率(MAF)<5%、哈代溫伯格平衡P<1×10-6的SNPs。

1.3 哈代-溫伯格值計算

本研究中，使用Plink軟件進行哈代溫伯格值計算，對哈代溫伯格值進行卡方適合性檢驗[14]，卡方適合性檢驗公式為：

式中：O(Observed)為每個SNP的觀測值；E(Expected)為每個SNP的期望值。根據卡方值和自由度得到P值，如P<0.05，則該SNP位點不符合哈迪-溫伯格定律；如P>0.05，則該SNP位點符合哈迪-溫伯格定律，基于R語言qqman軟件包繪制曼哈頓圖。

1.4 候選基因的富集分析

對計算出的哈代-溫伯格值進行統計分析，繪制曼哈頓圖。將篩選出從小到大的前1% SNPs作為受選擇位點，參照NCBI數據庫(https:∥www.ncbi.nlm.nih.gov/genome/?term=Oar_v4.0)中的Oar_v4.0基因組數據，對篩選出的SNPs基因注釋。利用David(https:∥david.ncifcrf.gov/)在線軟件對候選基因進行GO和KEGG富集分析[15]。

1.5 候選基因的蛋白互作網絡分析

利用String 11.0(https:∥string-db.org/)在線軟件對所獲得的候選基因進行蛋白互作網絡分析，其中包括蛋白質之間直接物理作用和蛋白質之間間接功能作用2個方面。

2 結果與分析

2.1 策勒黑羊哈代-溫伯格值檢測

使用Illumina ovine SNP 50 K bead chip對SNP分型后得到51 646個SNP位點，質控后剩余47 909個SNP位點。對47 909個SNPs進行哈代-溫伯格檢測，及對染色體上SNPs哈代溫-伯格檢測情況統計分析。在數學上,一般將發生概率小于0.05的事件稱為小概率事件[16]。突變率P<0.05為小概率，可以認為2、5、6、9、10、11、15、20、21、23和26號染色體符合哈代-溫伯格定律(突變率P<0.05)，1、3、4、7、8、12、13、14、16、17、18、19、22、24和25號染色體不符合哈代-溫伯格定律(突變率P>0.05)，突變率情況如表1?？傮w突變率P=0.052>0.05，此群體不符合哈代-溫伯格定律。

表1 策勒黑羊哈代-溫伯格檢測情況Table 1 Qira black sheep Hardy-Weinberg test

取哈代-溫伯格檢測從小到大的P值前1%[17]的SNPs，以P=0.006 3為閾值，得到SNPs 480個。將數據結果使用R語言中qqman包進行繪制曼哈頓圖(圖1)，圖中呈現候選SNPs和所對應的染色體分布情況。

黑線為閾值(-lg(P)=2.2)線，黑線以上點為篩選出的SNP。在-lg(P)中，P=0.006 3(哈代-溫伯格檢測從小到大排列前1%)。The black line is the threshold (-lg(P)=2.2) line, and the points above the black line are the selected SNPs. In -lg(P), P=0.006 3 (Hardy-Weinberg test ranks the top 1% from smallest to largest).

2.2 基因注釋及生物學功能分析

使用綿羊參考基因組(Oar_v4.0)，對篩選出來的480個SNPs上下游各延伸50 kb進行基因注釋，得到417個基因。在轉錄激活因子活性，RNA聚合酶II核心啟動子近端區域序列特異性結合條目(Transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding)基因MEF2A、SOX9、TFAP2B、EGR2、FUBP3、NEUROD6、CREB3L2、NFATC2、EBF2、MITF和POU2F3達到極顯著水平(P>0.01)。在軟骨內骨生長條目(Endochondral bone growth)內的基因BNC2、MSX2和OSTN達到極顯著水平(P>0.01)。在mRNA監測通路(mRNA surveillance pathway)內的基因PPP1CB、UPF1、PPP1CC、PABPC4L、PAPOLA、PPP2R2A、MSI2和PELO到極顯著水平(P>0.01)。將分析結果使用R語言繪制條形圖與KEGG氣泡圖，得到如下圖(圖2和3)所示結果。通過蛋白網絡互作圖(圖4)可知，除了基因LTBP1、ID3、ZNF70、PAPSS1、PCSK5、ARHGAP10、OSTN、GTF2I、HTR1B、SALL1、PDE4B、ZNF135、RCE1和ZMYM4外，其他33個基因都通過相互作用連接在一起。參與生物信號傳遞、基因表達調節、能量和物質代謝及細胞周期調控等過程。

圖2 策勒黑羊優異基因GO富集分析Fig.2 GO enrichment analysis of superior genes in Qira black sheep

圖3 策勒黑羊優異候選基因KEGG富集分析Fig.3 KEGG enrichment analysis of excellent candidate genes in Qira black sheep.

圖4 候選蛋白互作網絡Fig.4 Candidate protein interaction network

將顯著SNP位點與綿羊參考基因組(Oar_v4.0)比對，通過GO富集分析、KEGG通路分析、綿羊QTL數據庫(https:∥www.animalgenome.org/cgi-bin/QTLdb/OA/index)、String數據庫(https:∥string-db.org/)等手段分析候選基因生物學功能。與生長相關的基因有OSTN、TFAP2B和CTNNA2等。與抗逆相關的基因有ZNF70、CREB3L2和TMEM154等。與繁殖相關的基因有SOX9、GTF2A1和GNAQ等。經生物學功能分析后，得到54個候選基因，生物學功能及蛋白質氨基酸數、等電點等預測如下表(表2)。

3 討 論

3.1 哈代-溫伯格檢測

本研究中以0.05為界限，1、3、4、7、8、12、13、14、16、17、18、19、22、24、25號染色體不符合哈代-溫伯格定律，這些染色體受選擇程度高，環境對于策勒黑羊的選擇主要發生在這些染色體上。從整體看，此群體不符合哈代-溫伯格定律，說明策勒黑羊群體有發生突變、遷移、選擇、無隨機交配、群體縮小的情況。然而2、5、6、9、10、11、15、20、21、23和26號染色體符合哈代-溫伯格定律，但在這些染色體上存在哈代-溫伯格不平衡的SNPs。證明了這些染色體和其染色體上的分子標記還可以進一步被選育，來提高策勒黑羊生產性能。這一現象為策勒黑羊種質資源進一步開發利用提供了新的研究思路。例如，在OAR(Ovis aries chromosome，OAR)2存在與產肉性狀相關的分子標記；在OAR11上存在繁殖性狀相關分子標記，針對這些重要經濟性狀染色體和分子標記可以進一步進行分子選育，提高該品種的生產性能。并未獲得BMPR1B多胎性狀候選基因，可能因人工選育后到達平衡位點或者策勒黑羊存在更多的多胎基因。

3.2 策勒黑羊繁殖性狀的遺傳解析

繁殖性狀主要包括多胎、初生活仔數、初生重等，受微效多基因控制。本研究共發現13個基因和繁殖性狀相關，主要有GTF2A1、SOX9、PPARG、PPP1CC、MEF2A和GNAQ等。

一般轉錄因子IIA亞基1(General transcription factor IIA subunit 1，GTF2A1)基因位于OAR7，物理距離為89 482 053～89 516 754 bp。GTF2A1與羅德島紅母雞的產蛋性能相關[18]。狄冉等[19]研究發現：GTF2A1基因g.89505005G>A位點存在3種基因型：AA、AG和GG，且基因型頻率和等位基因頻率在單、多羔品種間差異均達到顯著水平，A等位基因是提高綿羊產羔數的一個潛在有效的DNA標記。ATP合成酶F1亞基α(ATP synthase F1 subunit alpha，ATP5A1)位于OAR18，物理距離為8 551 683～8 552 189 bp。編碼催化線粒體H(+)-ATP合酶的組分[20]，是公認的6種卵母細胞質量標志物之一[21]。

SRY盒轉錄因子9(SRY-box transcription factor 9，SOX9)基因位于OAR11，物理距離為57 682 208～57 799 964 bp。SOX9抑制對于雌性性腺發育至關重要。SOX9是一種轉錄因子，在發育過程中具有多種作用，包括介導睪丸分化的關鍵作用等，SOX9在人類和小鼠中的表達喪失會導致XY雌性發育，而在XX胚胎性腺中的不當激活可以使雄性發育[22]。在目前發現的基因中，DAX-1是SOX9的抑制因子，在雌性通路中，它們直接或間接抑制SOX9的表達；而在雄性通路中，SRY能解除DAX1對SOX9的抑制。

過氧化物酶體增殖物激活受體γ(Peroxisome proliferator activated receptor gamma，PPARG)位于OAR19，物理距離為56 552 358～56 616 707 bp。Ferst等[23]的研究中發現，PPARG參與細胞凋亡、細胞周期、雌二醇和孕酮合成以及代謝的調節。PPARG的相對mRNA豐度在優勢卵泡和次級卵泡中相似，在LH激增后降低，在排卵前升高。此外，發現顆粒細胞中PPARGmRNA水平與卵泡液中孕酮濃度之間存在二次相關性。

蛋白磷酸酶1催化亞基γ(Protein phosphatase 1 catalytic subunit gamma，PPP1CC)，位于OAR17，物理距離為49 020 928～49 022 373 bp。PPP1CC基因轉錄后，經過組織特異性剪接，產生普遍表達的同種型PPP1CC1和睪丸富集和精子特異性同種型PPP1CC2，這對于完成精子發生至關重要[24]。PPP1CC基因的靶向破壞可消除剪接變體PPP1CC1和PPP1CC2，導致男性因精子發生受損而不育[25]。

肌細胞增強因子-2(Myocyte enhancer factor 2D，MEF2)，位于OAR1，物理距離為5 979 082～6 074 239 bp。MEF2家族轉錄因子的成員是各種細胞類型和組織中細胞增殖，分化和侵襲的關鍵調節因子，并且調節人類胎盤發育過程中的滋養層增殖，侵襲和分化[26]。

G蛋白亞基alpha q(G protein subunit alpha q，GNAQ)，位于OAR2，物理距離為58 446 152～58 656 152 bp。GNAQ表達可調節 Kisspeptin表達，并通過Kisspeptin-GPR54信號通路間接調節GNRH基因，這些間接調控加強了GnRH效應并最終通過HPGA調節綿羊發情，GNAQ是控制綿羊季節性發情的主要基因[27]。

骨形態發生蛋白2(Bone morphogenetic protein 2，BMP2)，位于OAR13，物理距離為48 877 314～48 889 861 bp。骨形態發生蛋白(BMP)家族是TGF-β超家族中最大的配體亞群，與參與生殖細胞增殖、分化、凋亡和組織重塑的生物過程和細胞事件密切相關[28-29]。BMP2/4是20多種BMP亞型中的一組。越來越多的研究表明，BMP2影響子宮內膜細胞的增殖、分化和凋亡[30]。免疫組化結果顯示,BMP2蛋白在整個發情周期差異表達,并主要定位于子宮腔上皮細胞和腺上皮細胞[31]。蛋白磷酸酶2A調節亞基B55α(Protein phosphatase 2 regulatory subunit Balpha，PPP2R2A)，位于OAR2，物理距離為39 087 769～39 171 020 bp。Li等[32]發現PPP2R2A作為一種新的調節因子，通過體外調節湖羊子宮內膜基質細胞的增殖和凋亡來影響胚胎植入。

蛋白磷酸酶1催化亞基β(Protein phosphatase 1 catalytic subunit beta，PPP1CB)，位于OAR3，物理距離為35 555 357～35 598 569 bp。PPP1CB基因在牛早期性別分化過程中起重要作用[33]。

ATP合成酶F1亞基α(ATP synthase F1 subunit alpha，ATP5A1)，位于OAR23，物理距離為8 551 683～8 552 189 bp。通過比較蛋白質組學發現ATP5A1與水牛精子活力相關[34]。白定平等研究發現，ATP5A1與莫斯科鴨睪丸的分化和發育有關。

Erb-b2受體酪氨酸激酶4(Erb-b2 receptor tyrosine kinase 4，ERBB4)，位于OAR2，物理距離為212 081 308～212 592 819 bp。ERBB4在卵巢中具有重要作用，小鼠顆粒細胞中的ERBB4缺失造成粒細胞和卵母細胞的細胞間連接缺陷，導致調節卵泡局部微環境的基因表達發生變化[35]。

蛋白轉化酶枯草桿菌蛋白酶/kexin 5型(Proprotein convertase subtilisin/kexin type 5，PCSK5)，位于OAR2，物理距離為59 899 679～60 402 674 bp。PCSK5表達在小鼠卵巢卵泡發育中的作用[36]。采用實時熒光定量PCR對妊娠大鼠研究發現，PCSK5在妊娠早期和產后顯示出相對較高的水平，但在妊娠中期則顯示出較低的水平[37]。

ISL LIM 同源盒1(ISL LIM homeobox 1，ISL1)，位于OAR16，物理距離為27 720 710～27 730 148 bp。ISL1和JMJD3合作改變心臟表觀基因組，指導驅動心臟分化的基因表達變化[38]。ISL1在外生殖器發育過程中也起著關鍵作用[39]，是人類女性生殖道發育過程中上皮和間充質分化標志物之一。

本研究發現的與繁殖性狀相關基因與相關研究報道初步解析了策勒黑羊繁殖性狀，相關候選基因在策勒黑羊繁殖性狀中的驗證和作用機制有待進一步深入研究。確定的分子標記和候選基因的獲得，有助于我們下一步對策勒黑羊進行候選分子標記驗證。

3.3 策勒黑羊抗逆性狀的遺傳解析

策勒黑羊生存在塔克拉瑪干沙漠南緣，經自然和人工選擇具有抗逆性強。本研究共發現30個抗逆性相關的基因，針對常見的抗逆性相關候選基因進行闡述，其中包括ZNF70、TMEM154和CREB3L2等。

鋅指蛋白70(Zinc finger protein 70，ZNF70)，位于OAR17，物理距離為70 357 561～70 482 842 bp。經拷貝數變異(CNV)分析被確定為Nguni牛適應性增強的基因之一[40]。

綿羊跨膜蛋白154(Transmembrane protein 154，TMEM154)，位于OAR17，物理距離為4 832 841～4 906 239 bp。其多態性與進行性肺炎的易感性或抗性相關[41]。在北美綿羊中已經鑒定出與進行性肺炎感染有關的主要基因為TMEM154[42]。全基因組關聯研究已將TMEM154E35K多態性描述為選擇某些美國和歐洲品種中抗性動物的良好遺傳標記[43]。陳磊等[44]研究發現綿羊TMEM154基因在第2外顯子103 bp處發生了G-A突變，3種基因型AA、AG和GG中G103A位點的GG(E35)基因型與綿羊進行性肺炎的高敏感性相關；在策勒黑羊群體中屬于優勢基因型。

cAMP反應性元素結合蛋白3樣2(cAMP responsive element binding protein 3 like 2，CREB3L2)基因，位于OAR4，物理距離為 101 365 537～101 490 569 bp。該基因可在內質網應激傳感器BBF2H7/CREB3L2分泌管腔區域調控軟骨細胞增殖,可調控軟骨細胞增殖[45]。CREB3L2在垂體部位的表達是小鼠促進治療性蛋白生成的有效手段[46]。

免疫球蛋白λ樣多肽5(Immunoglobulin lambda like polypeptide 5，IGLL5)，位于OAR17，物理距離為70 387 677～70 429 652 bp。IGLL5與透明細胞腎細胞癌中腫瘤浸潤的免疫細胞相關[47]。洪東等通過全基因組測序發現，基因IGLL5為亞洲野驢免疫性狀候選基因[48]。

LIM結構域4(LIM domain only 4，LMO4)，位于OAR1，物理距離為63 829 917～64 013 649 bp。Kenny等[49]在小鼠胚胎中發現LMO4在胸腺原始細胞中被激活，且在成熟的T細胞中表達。

激活素A受體I型(Activin a receptor type I，ACVR1)編碼TGFβ受體超家族的骨形態發生蛋白I型受體。它參與各種生物過程，包括骨骼，心臟，軟骨，神經和生殖系統的發育和調節。此外，ACVR1還與纖維發育性骨化癥、心臟畸形和生殖系統改變等疾病相關[50]。

除了基因ZNF70、TMEM154、CREB3L2、IGLL5、LMO4和ACVR1有較多的相關研究和報道外，剩余24個基因SALL1、LTBP1和ID3等富集在與抗逆性相關的上皮細胞受細菌侵襲(Bacterial invasion of epithelial cells)、抗生素的生物合成(Biosynthesis of antibiotics)等功能通路。

以上與抗逆性狀相關基因的發現初步解析了策勒黑羊抗逆性狀。在此基礎上，將對策勒黑羊進行候選分子標記篩選的方法驗證，深入解析候選基因在策勒黑羊抗逆性狀中的作用機制。

3.4 策勒黑羊其他性狀的遺傳解析

轉錄因子 AP-2 β(Transcription factor AP-2 beta，TFAP2B)，位于OAR20，物理距離為23 153 665～23 178 184 bp。Ibrahim等[51]進行全基因組掃描發現基因TFAP2B與Barki綿羊的斷奶重相關。前人研究發現，在能量短缺時TFAP2B可以調節膳食脂肪攝入來減緩消瘦[52]。

鋅指MYM型蛋白4(Zinc finger MYM-type containing 4，ZMYM4)，位于OAR1，物理距離為10 059 208～10 173 385 bp。Ga?lle Marenne等[53]發現基因ZMYM4與肥胖相關。

骨細胞素(Osteocrin，OSTN)，位于OAR1，物理距離為194 361 786～194 378 063 bp參與負荷誘導的長骨生長，骨膜成骨細胞產生的OSTN通過增強C型利鈉肽(CNP)依賴性增殖和軟骨細胞的成熟來調節生長板的生長，從而導致長骨的伸長，此外，OSTN與CNP合作調節骨骼形成[54]。

Rho GTP酶激活蛋白10(Rho GTPase activating protein 10，ARHGAP10)，位于OAR17，物理距離為9 826 962～10 204 650 bp。ARHGAP10基因參與日本黑牛肌內脂肪的形成[55]。

鈣粘蛋白alpha 2(Catenin alpha 2，CTNNA2)，位于綿羊OAR3，物理距離為51 755 958～51 866 901 bp。CTNNA2與血漿維生素K水平相關[56]，CTNNA2也被確定為與首次產犢年齡和犢牛生長率顯著相關[57]。

生長分化因子6(Growth differentiation factor 6，GDF6)，位于綿羊OAR9，物理距離為80 657 044～80 740 501 bp。GDF6在最初因其序列同源性被認為是BMPs家族成員之一，不過其活性在軟骨中更突出，而且GDF6已被證明在體外和體內均能刺激軟骨細胞相關蛋白的表達[58]。

4 結 論

本研究共獲得480個SNPs，417個候選基因，其中與繁殖性狀相關的基因有SOX9、GTF2A1和GNAQ等；與抗逆性狀相關的基因有TMEM154、ZNF70和CREB3L2等。解析了策勒黑羊常年發情、多胎的分子遺傳機理，對抗逆性相關候選基因分析，發現策勒黑羊群體攜帶抗肺炎相關的分子標記，這些發現為策勒黑羊保種和改良提供了分子水平的參考。