基于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。总体突变率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等。解析了策勒黑羊常年发情、多胎的分子遗传机理，对抗逆性相关候选基因分析，发现策勒黑羊群体携带抗肺炎相关的分子标记，这些发现为策勒黑羊保种和改良提供了分子水平的参考。