广益黑猪繁殖性状的全基因组关联分析

朱吉，任慧波，吴发平，崔清明，刘莹莹，邓缘，李华丽，胡雄贵，左剑波，邹亿文，陈晨*，彭英林

广益黑猪繁殖性状的全基因组关联分析

朱吉1，任慧波1，吴发平2，崔清明1，刘莹莹1，邓缘1，李华丽1，胡雄贵1，左剑波1，邹亿文3，陈晨1*，彭英林1

(1.湖南省畜牧兽医研究所，湖南 长沙 410131；2.长沙县农业农村局，湖南 长沙 410100；3.湖南广益农业开发集团股份公司，湖南 新化 417600)

选取441头广益黑猪经产母猪为研究对象，利用猪50K SNP芯片对猪耳组织DNA进行基因型分型，PLINK 1.9质控后，采用GMAT中的重复力模型进行猪繁殖性状相关的全基因组关联分析，确定显著位点。结果表明：441头经产广益黑猪母猪的耳组织DNA基因分型共获得50 898个SNPs，经质控剩余46 165个SNPs位点用于关联分析；平均亲缘关系系数为–0.002 2，平均亲缘关系较远，不存在明显的群体分层；总产仔数性状有1个SNP在全基因组范围内达到显著相关，4个SNPs达到潜在显著关联，候选基因包括、ENSSSCG00000003753、、；产活仔数性状有6个SNPs潜在显著关联，候选基因包括、ENSSSCG00000022411、、ENSSSCG00000029007、、、ENSSSCG00000023343；弱仔性状有1个SNP达潜在显著关联，候选基因包括和；经基因功能分析推测，可能是影响猪总产仔数的候选基因。

广益黑猪；繁殖性状；全基因组关联分析；SNP；候选基因

提高猪的繁殖性能是其遗传育种研究的重点和难点之一[1]。大量的研究[2–4]表明，利用常规育种手段提高猪产仔性能的遗传进展较为缓慢，全基因组关联分析(GWAS)在动植物遗传育种上则发挥了重要作用，取得了快速发展。

已有GWAS用于猪总产仔数、产活仔数、排卵率、产死仔数、木乃伊等相关繁殖性状分析[5–8]。汤香[9]利用50K芯片对宁乡猪9个繁殖性状进行GWAS分析，发现32个显著的SNP位点与总产仔数、死胎数、妊娠周期、成活率和断乳死亡率等5个繁殖性状显著相关，经基因注释挖掘到15个可能与宁乡猪繁殖性状有关的候选基因；进一步的关联分析结果发现，基因第6外显子的55 bp处存在T/C突变，TT基因型的总产仔数和产活仔数均显著高于CC型的。蔡薇等[10]基于50K SNP芯片技术，搜寻到、、是最可能影响金华猪和嵊县花猪总产仔数和产活仔数的候选基因。宋志芳[11]利用GWAS分析深县猪繁殖性能，初步推断、、等12个基因可能是潜在的候选基因。莫家远等[12]利用GWAS分析，在巴马香猪上筛选到4个基因(、、、)影响产活仔数。陈俊峰等[13]采用GGP porcine 50K芯片对淮南母猪开展全基因组关联分析，在总产仔数、有效产仔数、窝重3个繁殖性状共发现了49个SNPs和64个候选基因，筛出了、、和等作为总产仔数性状的候选基因。SCHNEIDER等[14]利用全基因组关联分析猪的死胎性状，定位到5个QTL，筛选得到候选基因。目前，对培育猪种繁殖性状的GWAS研究较为鲜见。

广益黑猪是以湘西黑猪为母本、杜洛克猪为父本，历经6个世代培育而成的优良新猪种。广益黑猪综合体现了中国猪种的早熟多产、耐粗饲、肉质香嫩和外来猪种体型大、生长快、瘦肉率高的特点。为进一步挖掘广益黑猪的繁殖性能的遗传潜力，笔者选取了441头广益黑猪母猪，利用中芯一号50K SNP芯片开展GWAS分析，筛选主效基因，以期为分子育种手段提升繁殖性能提供依据。

1　材料与方法

1.1　供试广益黑猪

441头经产广益黑猪母猪均来源于湖南广益农业开发集团股份公司，涵盖12个家系。所有个体的营养标准一致、饲养环境相同。

1.2　耳组织样采集及DNA提取

于2021年12月采集414头母猪的耳组织样，保存于75%乙醇冻存管，置泡沫箱，交北京康普森公司提取DNA，经质检合格后用猪中芯一号50K SNP芯片分型并输出结果成plink格式，以备开展后续分析。

1.3　芯片数据质控和基因型填充

运用PLINK 1.9进行芯片数据质控，删除掉基因型缺失率>10%、最小等位基因频率<1%、哈迪–温伯格平衡检验<0.000 001的位点，质控后运用Beagle5.0基因型填充[15]。

1.4　全基因组关联分析模型

运用PLINK 1.9计算广益黑猪群体的主成分矩阵(PCA)、亲缘系数矩阵，并运用GMAT进行重复力模型GWAS分析。模型如下：

=++++

其中，为表型向量；为固定效应的设计矩阵；为固定效应向量，包括胎次和年份效应；为检验标记效应的设计矩阵；为检验标记效应向量；为多基因效应向量，服从正态分布(0,2)，为标记推断的亲缘关系矩阵，2为多基因效应的方差；为永久环境效应向量，服从正态分布(0,σ2)，为单位矩阵，σ2为永久环境效应的方差；为残差效应向量，服从正态分布(0,2)，2为残差方差。采用1/作为全基因组显著阈值，5/作为潜在显著性阈值，其中，为有效SNP位点数。

1.5　候选基因注释

针对显著SNP位点的名称和所在基因组位置，运用ANNOVAR，参照猪基因组序列版本(Sus scrofa 10.2)，在其上下游进行基因注释。

1.6　富集分析

运用R包clusterProfiler，对注释候选基因进行GO和KEGG功能富集分析。采用R语言的enrich plot包进行可视化绘图。结合文献查阅、综合分析，确定生猪繁殖性状候选基因。

2　结果与分析

2.1　供试广益黑猪的表型数据分析结果

供试广益黑猪的总产仔数、产活仔数、弱仔数、死胎数等繁殖性状列于表1。从表1可看出，每窝总产仔数、产活仔数、死胎数、弱仔数的平均值分别为9.91、9.36、0.22、0.24头，变异系数分别为19.86%、23.50%、282.96%、286.06%。该结果表明广益黑猪繁殖性状的变异系数较高，具有较大的遗传改良潜力。

总产仔数数据为722窝样本统计所得；产活仔数、死胎数、弱仔数数据为436窝样本统计所得。

2.2　供试广益黑猪的基因分型及数据质控

基因分型后共获得50 898个SNPs。质控时，剔除SNP缺失率>10%的位点531个、最小等位基因频率<1%的位点3891个、哈迪–温伯格平衡检验<0.000 001的位点311个，剩余46 165个SNPs位点用于后续关联分析。

2.3　供试广益黑猪的群体结构分析结果

图1体现的是猪只之间的亲缘关系远近，414头广益黑猪母猪的平均亲缘关系系数为–0.002 2，最小值为–0.211 3，最大值为0.829 9，平均亲缘关系较远。从图2可看出，供试广益黑猪群体不存在明显的群体分层。

图1　供试广益黑猪的亲缘关系矩阵热图

图2　供试广益黑猪PCA的群体结构

2.4　供试广益黑猪的全基因组关联分析结果

从图3可看出，4个图中值的观测与期望值均显示出强烈偏差，提示与性状显著关联的SNP可能会比预期多。

a　总产仔数；b　产活仔数；c　死胎数；d　弱仔数。

从图4可看出，总产仔数、产活仔数、弱仔等3个性状均检测到SNPs达到全基因组显著或潜在显著水平，死胎性状未检测到达全基因组潜在显著水平的SNP；总产仔数的显著或潜在显著水平的SNPs位于6和12号染色体上，产活仔数的潜在显著水平的SNPs位于6和13号染色体上，弱仔的潜在显著水平的SNPs位于13号染色体上。

从表2可知，3个性状共计有12个SNPs达到全基因组潜在显著水平，关联基因13个。与总产仔数显著关联的SNP有1个，位于6号染色体上，与总产仔数潜在显著关联的SNPs有4个，分别位于6号和12号染色体上。6号染色体上的1个显著SNP(CNC10062413)、1个潜在显著SNP(CNC10062410)处于1个未知功能片段和基因之间的位置，1个潜在显著SNP(CNC10062432)处于和ENSSSCG00000003753基因之间的位置；12号染色体的2个潜在显著SNPs(CNC10120322、CNC10120325)分别处于、的内含子上。与产活仔数潜在显著关联的SNPs有6个，分别位于6号和13号染色体上。6号染色体上的2个潜在显著SNPs(CNC10061224、CNC10061225)处于和ENSSSCG00000022411基因之间的位置，1个潜在显著SNP(CNC10061169)处于和ENSSSCG00000029007基因之间的位置，2个潜在显著SNPs(CNC10061151、CNC10061152)处于和基因之间的位置；13号染色体的1个潜在显著SNP(CNC10130001)处于1个未知功能片段和ENSSSCG00000023343基因之间的位置。与弱仔潜在显著关联的SNP有1个，位于13号染色体上，处于和基因之间的位置。

a　总产仔数；b　产活仔数；c　死胎数；d　弱仔数。虚线水平线为全基因组显著关联阈值线(1/n)；实线水平线为全基因组水平潜在显著阈值线(5/n)。

表2　供试广益黑猪繁殖性状全基因组关联显著的SNPs和相关基因

“*”示SNP与性状显著关联。

2.5　GO和KEEG富集分析结果

GO和KEGG富集分析结果(表3和表4)表明，总产仔数相关基因主要在磷酸转移酶活性、磷脂酰肌醇介导的信号、胞质分裂的调控、胞质运输等方面起作用，并参与磷酸肌醇代谢、磷脂酰肌醇信号系统、自体吞噬等；主要在p38–MAPK级联、应激激活MAPK级联反应的正向调节、血管新生的细胞迁移调控、细胞核因子酉乙蛋白信号的正调控、细胞衰老的负调控起作用，并在下丘脑分泌促性腺激素释放激素信号通路、神经营养素信号通路、MAPK信号通路等富集。产活仔数相关基因主要在造血祖细胞分化、血细胞生成、造血或淋巴器官发育等生物学进程富集。弱仔相关基因主要在细胞紧密连接、趋化因子信号通路、Rap1信号通路、肌动蛋白骨架调节、环磷腺苷信号通路、Ras信号通路等富集。

表3　供试广益黑猪注释候选基因的GO富集分析结果

表4　供试广益黑猪注释候选基因的KEGG富集分析结果

3　结论与讨论

本研究中，与总产仔数显著关联的SNP关联的候选基因，也名为34，是唯一一个第3类磷脂酰肌醇3–激酶(PI3K)，在自噬、细胞内吞、Mtor/S6K1信号通路、内吞体–高尔基体膜运输等生物过程中起着重要调节作用[16]。以往研究[17–18]发现，在小鼠神经、肝脏、心脏等多个组织器

官的正常生理功能起重要作用，小鼠胚胎急性删除成纤维细胞后，会导致内吞作用受阻，自噬退化。孙妍[19]研究发现，可能通过自噬、内吞等生物过程的介导，对小鼠的卵母细胞发育起重要作用。通过有效激活包含JNK、p38、ERK5在内的MAPKs、NF–κB信号通路，调控胚胎心血管的早期发育、肌细胞形成、炎症和免疫应答等[20]。YANG等[21]在小鼠上发现MEKK3(有丝分裂原活化蛋白激酶)，MAP3K家族的丝氨酸/苏氨酸蛋白激酶，在小鼠的早期胚胎发育中有着至关重要的作用，敲除可造成小鼠早期胚胎因血管发育不全以及卵黄囊结构破损而死亡。另有研究表明，ERK信号通路参与正常妊娠的维持，在羊的子宫动脉收缩力的调节中发挥重要作用[22]；p38–MAPK活化会导致人胎儿组织发生衰老和早衰[23]；对小鼠胚胎血管生成至关重要，并且在胚胎血管发育中发挥内在作用[24]。HWANG等[25]研究显示，的mRNA表达量与产仔数呈负相关；同一SNP，也被报道与德国大白猪的产活仔数有关[26]。综合分析可知，基因可能是总产仔数的重要候选基因。

与产活仔数潜在显著关联的SNP邻近的，属于锌指蛋白家族，对其研究较少，仅局限于基因结构、转录调控等信息的研究[27]。锌指蛋白在生物体内分布广泛，作为转录调控因子，在基因表达调控、细胞分化、胚胎发育等方面发挥着重要的作用[28]。基因有可能是猪产活仔数相关的候选基因，仍需进一步验证。

与弱仔数潜在显著关联的SNP关联的候选基因，调节细胞骨架重组、基因转录、细胞迁移和黏附、细胞增殖凋亡等过程[29]。现有的研究基本围绕其“促癌基因”特性开展的[30]。本研究结果的基因富集结果显示，在Rap1信号通路和Ras信号通路富集。Rap1是一种小GTP酶，以高度依赖细胞类型的方式调节MAP激酶(MAPK)活性。Ras蛋白是GTP酶，激活的RAS(RAS–GTP)通过效应因子调节PI3K活性。这与相关研究[31]的与MAPK/P38信号通路相关、磷酸化后可间接激活MAPK/Erk信号通路的结果相印证，也提示了该基因作为繁殖性状相关基因的可能，后续应进一步验证。

[1] 赵谦，浦亚斌，关伟军，等．猪重要性状全基因组关联分析的研究进展[J]．畜牧兽医学报，2015，46(6)：873–881．

[2] 林晋军，姚威，彭沙莎，等．水稻抽穗期性状QTL的全基因组关联定位分析[J]．湖南农业大学学报(自然科学版)，2021，47(3)：283–290．

[3] 张雁明，邢国芳，刘美桃，等．全基因组关联分析：基因组学研究的机遇与挑战[J]．生物技术通报，2013(6)：1–6．

[4] 肖瑜，马海明．猪的全基因组测序研究进展[J]．中国畜牧杂志，2019，55(5)：15–20．

[5] CHEN Z T，YE S P，TENG J Y，et al．Genome-wide association studies for the number of animals born alive and dead in duroc pigs[J]．Theriogenology，2019，139：36–42．

[6] VERARDO L L，SILVA F F，LOPES M S，et al. Revealing new candidate genes for reproductive traits in pigs：combining Bayesian GWAS and functional path- ways[J]．Genetics Selection Evolution，2016，48：9．

[7] HE L C，LI P H，MA X，et al．Identification of new single nucleotide polymorphisms affecting total number born and candidate genes related to ovulation rate in Chinese Erhualian pigs[J]．Animal Genetics，2017，48(1)：48–54．

[8] WU P X，WANG K，ZHOU J，et al．A genome wide association study for the number of animals born dead in domestic pigs[J]．BMC Genetics，2019，20(1)：4．

[9] 汤香．宁乡猪繁殖性状候选基因GWAS筛选及其关联分析[D]．长沙：湖南农业大学，2020．

[10] 蔡薇，罗才玉，项云，等．基于50K SNP 芯片技术对金华猪和嵊县花猪繁殖性状的全基因组关联分析[J]．中国畜牧杂志，2020，56(7)：52–58．

[11] 宋志芳．利用GWAS对深县猪繁殖性状的研究[D]．保定：河北农业大学，2018．

[12] 莫家远，高九昱，奉玲丽，等．巴马香猪产活仔数性状全基因组关联分析[J]．中国畜牧兽医，2020，47(12)：3965–3975．

[13] 陈俊峰，王璟，张家庆，等．全基因组关联分析筛选淮南猪繁殖性状候选基因[J]．中国畜牧杂志，2023，59(7)：111–114．

[14] SCHNEIDER J F，MILES J R，BROWN-BRANDL T M，et al．Genomewide association analysis for average birth interval and stillbirth in swine[J]．Journal of Animal Science，2015，93(2)：529–540．

[15] MAREES A T，DE KLUIVER H，STRINGER S，et al．A tutorial on conducting genome-wide association studies：quality control and statistical analysis[J]．International Journal of Methods in Psychiatric Research，2018，27(2)：e1608．

[16] BACKER J M．The regulation and function of Class Ⅲ PI3Ks：novel roles for Vps34[J]．The Biochemical Journal，2008，410(1)：1–17．

[17] SIMONSEN A，LIPPÉR，CHRISTOFORIDIS S，et al. EEA1 links PI3K function to Rab5 regulation of endosome fusion[J]．Nature，1998，394：494–498．

[18] CALLAGHAN J，SIMONSEN A，GAULLIER J M，et al. The endosome fusion regulator early endosomal autoantigen 1(EEA1) is a dimer[J]．The Biochemical Journal，1999，338(Pt 2)：539–543．

[19] 孙妍．PIK3C3在卵母细胞中的作用研究[D]．南京：南京医科大学，2017．

[20] 赵爽，吴丽莎，陈翔．MEKK3蛋白激酶的研究进展[J]．复旦学报(医学版)，2010，37(6)：742–746．

[21] YANG J H，BOERM M，MCCARTY M，et al．Mekk3 is essential for early embryonic cardiovascular develop- ment[J]．Nature Genetics，2000，24(3)：309–313．

[22] XIAO D， ZHANG L．ERK MAP kinases regulate smooth muscle contraction in ovine uterine artery：effect of pregnancy[J]．American Journal of Physiology．Heart and Circulatory Physiology，2002，282(1)：H292–H300．

[23] MENON R，PAPACONSTANTINOU J．p38 Mitogen activated protein kinase(MAPK)：a new therapeutic target for reducing the risk of adverse pregnancy outcomes[J]. Expert Opinion on Therapeutic Targets，2016，20(12)：1397–1412．

[24] FISHER O S，DENG H Q，LIU D，et al．Structure and vascular function of MEKK3-cerebral cavernous malformations 2 complex[J]．Nature Communications，2015，6：7937．

[25] HWANG J H，AN S M，YU G E，et al．Association of single-nucleotide polymorphisms in NAT9 andgenes with litter size traits in Berkshire pigs[J]．Archives Animal Breeding，2018，61(4)：379–386．

[26] SPÖTTER A，HAMANN H，MÜLLER S，et al．Effect of polymorphisms in four candidate genes for fertility on litter size in a German pig line[J]．Reproduction in Domestic Animals 2010，45(4)：579–584．

[27] HUNTLEY S，BAGGOTT D M，HAMILTON A T，et al. A comprehensive catalog of human KRAB-associated zinc finger genes：insights into the evolutionary history of a large family of transcriptional repressors[J]．Genome Research，2006，16(5)：669–677．

[28] CLARY L M，OKKEMA P G．The EGR family gene egrh-1 functions non-autonomously in the control of oocyte meiotic maturation and ovulation in[J]. Development，2010，137(18)：3129–3137．

[29] BOLLAG G，CROMPTON A M，PEVERLY-MITCHELL D，et al．Activation of Rac1 by human Tiam 1[J]. Methods in Enzymology，2000，325：51–61．

[30] LEE S H，KUNZ J，LIN S H，et al．16-kDa prolactin inhibits endothelial cell migration by down-regulating the Ras-Tiam1-Rac1-Pak1 signaling pathway[J]．Cancer Research，2007，67(22)：11045–11053．

[31] WOODCOCK S A，ROONEY C，LIONTOS M，et al. Src-induced disassembly of adherens junctions requires localized phosphorylation and degradation of the Rac activator Tiam1[J]．Molecular Cell，2009，33(5)：639–653．

Genome-wide association study of reproductive traits in Guangyi black pigs

ZHU Ji1，REN Huibo1，WU Faping2，CUI Qingming1，LIU Yingying1，DENG Yuan1，LI Huali1，HU Xionggui1，ZUO Jianbo1，ZOU Yiwen3，CHEN Chen1*，PENG Yinglin1*

(1.Hunan Institute of Animal and Veterinary Science, Changsha, Hunan 410131, China; 2.Changsha County Agricultural and Rural Bureau, Changsha, Hunan 410100，China; 3.Hunan Guangyi Agricultural Development Corporation, Xinhua, Hunan 417600, China)

This experiment involved the collection of ear tissue samples from 441 multiparous sows of Guangyi black pigs. DNA was extracted and genotyped using a porcine 50K SNP beadchip. Quality control(QC) was performed using PLINK 1.9 software. Genome-wide association analysis of porcine reproductive traits was conducted based on QC results using the mixed linear model of GMAT. Results revealed that among 414 Guangyi black sows, a total of 50 898 SNPs were obtained, but after quality control, 46 165 SNPs remained suitable for GWAS. The population’s kinship coefficient was –0.002 2, indicating a relatively low average genetic relationship and negligible population stratification. One SNP exhibited significant correlation, and four SNPs showed potential significant correlation with total number born at a genome-wide level. The associated candidate genes included, ENSSSCG00000003753,and. Additionally, six SNPs potentially displayed significant correlation with the number born alive, involving associated candidate genes such as, ENSSSCG00000022411,, ENSSSCG00000029007,,and ENSSSCG00000023343. Furthermore, one SNP showed potential significant correlation with the number born weak, and the associated candidate genes wereand. Gene functional analysis suggested thatandmight be candidate genes influencing the total number born in Guangyi black pigs.

Guangyi black pigs; reproductive traits; genome-wide association analysis; SNP; candidate gene

S828.2

A

1007–1032(2023)04–0453–08

10.13331/j.cnki.jhau.2023.04.012

2023–01–06

2023–08–06

湖南省自然科学基金项目(2021JJ30386)；湖南省重点研发计划项目(2021NK1009、2020NK2024、2019NK2193)；湖南省重点实验室开放研究基金项目(2017TP1030)

朱吉(1981—)，女，湖南望城人，硕士，副研究员，主要从事猪的遗传育种研究，24625546@qq.com；*通信作者，陈晨，博士，副研究员，主要从事猪的遗传育种研究，404071654@qq.com

朱吉，任慧波，吴发平，崔清明，刘莹莹，邓缘，李华丽，胡雄贵，左剑波，邹亿文，陈晨，彭英林．广益黑猪繁殖性状的全基因组关联分析[J]．湖南农业大学学报(自然科学版)，2023，49(4)：453–460．

ZHU J，REN H B，WU F P，CUI Q M，LIU Y Y，DENG Y，LI H L，HU X G，ZUO J B，ZOU Y W，CHEN C，PENG Y L．Genome-wide association study of reproductive traits in Guangyi black pigs[J]．Journal of Hunan Agricultural University(Natural Sciences)，2023，49(4)：453–460．

http://xb.hunau.edu.cn

责任编辑：邹慧玲

英文编辑：柳正