全基因组关联分析鉴定法系大白猪繁殖性状基因组区域和候选基因

李秀领，陈健梅，吴姿仪，韩雪蕾，乔瑞敏，王克君，杨 峰，王献伟，李新建*

（1.河南农业大学动物科技学院，河南郑州 450046；2.河南省畜牧总站，河南郑州 450008）

母猪的繁殖能力在养猪业的经济效益构成中起着重要作用。为了获取更高的收益，提高母猪的繁殖性能尤为重要。总产仔数（Total Number Born，TNB）、产活仔数（Number Born Alive，NBA）和出生窝重（Litter Weight Born Alive，LWB）是评估母猪繁殖性能的关键指标。然而，繁殖性状的遗传结构较为复杂，由于其低遗传力特性（约0.1），传统的育种计划获得的遗传改良非常缓慢。因此，标记辅助选择（Marker-Assisted Selection，MAS）和基因组选择（Genomic Selection，GS）被广泛应用于繁殖性状的遗传改良。

近几十年来，随着分子育种技术的发展，全基因组关联分析（Genome-Wide Association Study，GWAS）已被广泛用于识别复杂性状中的数量性状基因座（Quan titative Trait Loci，QTL）和数量性状核苷酸（Quantita tive Trait Nucleotides，QTN）。在猪相关的研究中，GWAS 已经成为定位与复杂性状相关的SNPs 和候选基因的有力方法之一。在GWAS 的统计方法中，单步全基因组关联分析（Single-Step GWAS，ssGWAS）已被用于研究猪的经济性状。与传统GWAS 相比，ssGWAS 可以将个体的基因型、表型和系谱信息都整合进混合线性模型中。因此，ssGWAS 可以提高低遗传力性状候选基因定位的精度。

此外，繁殖性状受许多微效基因控制，但到目前为止，只有少数候选基因被确定，而且只是解释了一小部分遗传变异。因此，繁殖性状的具体遗传机制仍然需要进一步研究。本研究的目的是通过在法系大白猪群体中进行ssGWAS 检测影响繁殖性状（TNB、NBA 和LWB）的遗传变异并确定候选基因。

1 材料与方法

1.1 动物和表型数据 在本研究中，收集了河南某核心养猪场392 头法系大白猪的系谱信息以及2020 年1—10 月的表型信息。其中，表型包括TNB、NBA 和LWB，所有表型记录均为母猪第1 胎次产仔记录。所有个体均被分组在标准的商业猪舍中并且每个个体都是在相同的饲料和饲养管理条件下饲养的。TNB、NBA和LWB 的表型记录的平均值、标准差、最小值和最大值见表1。

表1 繁殖性状的描述性统计

1.2 基因型数据和质控 使用组织基因组DNA 提取试剂盒从猪耳组织样中提取DNA，并使用纽勤公司GeneSeek Genomic Profiler（GGP）猪50K SNP 芯 片进行基因分型。数据质量控制由PLINK v1.90进行。剔除基因型检出率<0.90、次等位基因频率<0.05 和Hardy-Weinberg 平衡检验<10的SNP 标记，删 除SNP 标记检出率<0.90 的个体。所有位于未映射区域的SNPs 也被排除。最终保留了43 098 个SNPs 和392 个个体用于后续分析。

1.3 遗传参数估计 使用GCTA v1.93.2进行限制性最大似然算法来估计所有性状的遗传参数。所有的常染色体SNPs 被用来计算基因组关系矩阵（GRM），然后计算主成分。用单变量动物模型来估计遗传力，用双变量GREML 分析来估计性状间的遗传相关。统计模型为：

其中，是表型观测值的向量；代表固定效应，包括前3 个主成分、产仔季节和妊娠天数，是的相关矩阵；是随机多基因效应的向量；是残差效应的向量。

1.4 ssGWAS GWAS 是使用BLUPF90 软件家族进行的。RENUMF90 模块用于准备表型和系谱数据，BLUPF90 模块用于计算基因组估计育种值（Genomic Estimated Breeding Values，GEBV），POSTGSF90 模块用于将GEBVs 转换为SNPs 效应，并计算相邻窗口所解释的方差。

在本研究中，单性状动物模型被用于GWAS 分析。统计学模型可以写成：

其中，是基于所有动物系谱的分子关系矩阵，是基因分型个体的分子关系矩阵的子集，是基因组关系矩阵。上标-1 表示它是相关矩阵的逆数。VanRaden（2008）构建G 矩阵的方法如下：

其中，Z 是根据等位基因频率编码的标记矩阵，D 是SNP 效应方差的权重对角矩阵（最初D=I）。ssGWAS的计算方法如下：

①构建G 矩阵；

②计算所有动物的GEBVs；

⑤归一化SNP 权重以保持总方差不变；

⑥根据SNP 权重重新计算矩阵；='；

⑦循环步骤②。

从第②步到第⑦步重复一次迭代过程，按照Wang等的描述，计算第i 个连续SNP（SNP 窗口）解释的遗传变异百分比：

Beissinger 等报告说，与其他较大的窗口尺寸相比，5 个或10 个SNP 组成的滑动窗口的检测率与假阳性率都是最有利的。因此，本研究计算了由10 个相邻SNP组成的滑动窗口所解释的加性遗传变异的比例。

1.5 候选基因组区域注释 解释的遗传变异大于1%的基因组窗口被认为是与性状相关的。在11.1基因组（http://www.ensembl.org/）上搜索了候选窗口内的基因或部分重叠的基因。使用大鼠基因组数据库（RGD，https://rgd.mcw.edu/rgdweb/homepage/）进行人类表型本体论（Human Phenotype Ontology，HPO）分析，选择Bonferroni 校正后<0.05 的条目进行后续分析。随后，使用KOBAS 在线网站对与生殖系统相关条目中涉及的候选基因进行分析，以确定基因参与的KEGG 信号通路。<0.05 的KEGG 通路被选中。使用PigQTLdb 数据库（http://www.animalgenome.org/cgi-bin/QTLdb/SS/index）对猪基因组中先前确定的QTL 进行评估。

2 结果

2.1 表型统计 表1 列出了3 个繁殖性状的描述性统计。其中，TNB 和NBA 的范围是0～24，平均值分别为14.91 和13.11。LWB 平均值为15.14。此外，所有性状的变异系数值都较高，这表明数据在某种程度上是离散的。另外，3 个性状之间高度正相关，表型相关系数分别为0.86、0.67 和0.81（表2）。

2.2 遗传力和遗传相关 本研究对3 个繁殖性状的遗传力和遗传相关都进行了估计（表2）。3 个性状的估计遗传力在0.06～0.12 之间，其中，NBA 遗传力最高，为0.12，TNB 遗传力最低，为0.06。3 个性状两两间的遗传相关从0.62 到0.83，都呈现高度正相关。

表2 繁殖性状的遗传力（对角线，粗体）、遗传相关（上三角）和表型相关系数（下三角）

2.3 GWAS 本研究将连续10 个SNPs 划分为一个窗口后，共有43 098 个窗口用于后续的GWAS 分析，其中最后9 个窗口包含的SNPs 少于10 个。

本研究中3 个性状分别鉴定到95、72、109 个相关的窗口，它们分别解释了TNB、NBA 和LWB 总遗传变异的1% 或以上（图1）。共有40 个基因组区域与繁殖性状相关，其中TNB、NBA 和LWB 分别包括15、14、11 个区域。对于TNB，共有15 个基因组区域，分布在染色体（SSC）1、3、5、7、13、15、16、17 和X 上，其中最高的区域解释了4.794%的遗传变异。有14 个基因组区域与NBA 显著相关，它们分布在SSC 1、2、3、6、7、16、17 和X 上，最高的区域解释了2.214%的遗传变异。对于LWB，109 个窗口串联成11 个基因组区域，这些基因组区域分布在SSC 2、3、6、10、13、16 和X 上，它们所解释的最大遗传变异为2.364%（表3）。

表3 繁殖性状相关基因组区域

图1 繁殖性状ssGWAS 曼哈顿图

2.4 候选基因鉴定 经过注释，共得到532 个基因，包括381 个编码蛋白的基因。随后，381 个基因被用来进行HPO 分析，结果富集在以下3 个条目：泌尿生殖系统表型、心血管系统表型和神经系统表型（图2A）。在泌尿生殖系统表型中，共得到81 个候选基因（图2B，2C）。有趣的是，大多数候选基因都位于2 号染色体上。随后，对候选基因进行了进一步的KEGG 分析，得到了82 条显著的通路。在82 条显著的通路中共发现7 条与生殖相关的KEGG 通路，包括5 条与生殖激素相关的信号通路和2 条卵细胞减数分裂及成熟相关的通路（表4）。对于所有的候选基因，本研究查询了以前的研究是否报道过与生殖性状有关的文献，结果列于表5。

表4 候选基因参与的显著的繁殖相关KEGG 通路

表5 候选基因中先前报道过的与繁殖相关的基因

图2 候选基因鉴定

3 讨 论

产仔性状包含仔猪的TNB、NBA 和LWB，与母猪的繁殖性能密切相关，影响猪的生产效率和经济利益。因此，探索猪产仔性状的遗传结构是极为重要的。此外，探索产仔性状的遗传机制非常复杂，整个繁殖过程受多种激素和生理途径的调节。GWAS 是检测与目标性状相关的候选基因最有效方法之一。本研究利用ssGWAS方法探索了392 头大白猪的TNB、NBA 和LWB 的遗传变异，通过HPO 和KEGG 分析进一步挖掘了与生殖相关的基因以及这些基因涉及的通路。

繁殖性状的遗传力通常在0.1 左右。在本研究中，3个繁殖性状的估计遗传力为0.06～0.12，与前人研究结果一致。此外，3 个性状之间的遗传和表型相关都呈现出高度正相关，这些结果与Wolf 等的研究一致。

本研究发现了40 个解释1%以上遗传变异的基因组区域。在这些区域中，还发现了一些以前报道的与繁殖有关的QTL。在3 个性状中SSC13 上的48.65～52.85 Mb 区域包含了最多解释超过1% 的遗传变异的窗口，并且在这个基因组区域中有9 个先前报道的相关QTLs。这些QTLs 包括3 个黄体数量QTLs（493、8 826、24 285），2 个卵巢重量QTLs（18 354、18 356），1 个青春期年龄QTLs（8 824），1 个子宫角重量QTLs（18 409），1 个总产仔数QTLs（15 1333）和1 个死胎数QTLs（57 525）。此外，SSC14 的104.44～106.90 Mb 区域解释了2.87%的TNB 性状的遗传变异，并且在此区域注释到了基因。据报道，基因在胚胎发育中起着不可或缺的作用。

基于关联分析鉴定到的显著基因组区域内的蛋白编码基因，本研究进行了HPO 分析，发现这些基因聚成了3 类：泌尿生殖系统、心血管系统和神经系统。共有7 种表型涉及泌尿生殖系统，包括原发性闭经、闭经、月经周期异常、女性生殖系统生理功能异常、泌尿生殖系统异常、生殖系统异常以及生殖系统生理学异常。这7 个表型共包含81 个候选基因，且大多聚集在SSC 2 上。经过搜索文献，在SSC2 上获得了10 个先前报道的与生殖有关的基因。其中，畸形表皮自动调节因子1（DEAF1）编码一个多域转录因子。相关结合位点的变化可能与山羊的产仔数有关。基因参与磷酸戊糖途径（PPP）的第2 阶段，并在糖酵解和PPP 之间提供一个链接。缺陷的患者在新生儿期出现低出生体重，甚至在婴儿期早期死亡。胰岛素样生长因子2（）基因是一个父系印记基因，已知它能增加猪的肌肉生长并减少脂肪沉积。然而，也有报道称对母猪的繁殖力有影响。此外，基因已被确认为人类、小鼠和牛的一个印记基因。印记基因与早期胚胎和胎盘的发育有关系。这些基因的表达水平的变化可以导致异常的生长和胚胎死亡。据报道，在猪的心脏、卵巢等部位有母体表达，在肺和肾有双侧表达。因此，基因在猪身上的印记还有待验证。基因是sirtuin家族的一个重要成员，在调节细胞生理和衰老方面发挥了广泛的作用。Gorczyca 等发现基因在猪卵巢卵泡细胞中表达。

为进一步了解候选基因的功能和这些基因之间的联系，本研究还进行了KEGG 分析，在所有显著的KEGG 信号通路中得到了7 条与繁殖性状相关的通路。这7 条通路分别是催产素信号通路、催乳素信号通路、GnRH 信号通路、松弛素信号通路、雌激素信号通路、母细胞减数分裂和孕激素介导卵细胞成熟通路。其中催产素、催乳素、促性腺激素释放激素（GnRH）、松弛素和雌激素都属于生殖激素，与生殖活动如排卵、怀孕、胚胎发育、分娩和泌乳直接相关。参与生殖激素相关途径的共有7 个基因，即和CREB3L3。其中，MAP2K2编码的蛋白质属于MAPK 激酶家族，基因对胎盘形成期间合体滋养层的正常发育有贡献。基因属于Ras 肿瘤基因家族，其编码的产物在信号转导途径中发挥作用。有报道称miR-210 通过和调节水牛排卵前的颗粒细胞功能。基因参与调控血管生长和基础血压。Zou 等研究表明基因在胎盘血管的形成中有重要作用。基因编码AMP 依赖性转录因子家族，它可以调节脂质代谢并介导对感染的反应。相关研究表明通过控制GH 信号传导，对身体的生长有负向调节作用。此外，有2 个基因参与卵母细胞减数分裂和成熟相关通路，分别是和。Aurora-A（AURKA）是一种丝氨酸/苏氨酸蛋白激酶。Yao 等研究结果显示基因在猪卵母细胞减数分裂成熟、受精和早期胚胎有丝分裂过程中对微管的组装起着关键的调节作用。基因编码胰岛素，这是一种多肽激素，可调节生理学的各个方面，是决定后生动物雌性生殖的关键因素。据报道，胰岛素参与了奶牛的卵泡成熟。

4 结 论

本研究使用50K SNP 芯片对法系大白猪群体中的TNB、NBA 和LWB 表型进行了ssGWAS 分析，鉴定到了40 个基因组区域。区域内注释到的所有编码蛋白的基因被用于进行HPO 分析，发现了81 个与繁殖有关的候选基因，这些基因主要涉及GnRH 信号通路、雌性激素信号通路和卵细胞成熟通路等。其中，本研究还发现一些与产仔数有关的基因，包括、、、、等。尽管筛选出了一些与繁殖相关的候选基因，但部分基因在畜禽中还没有被研究过，因此这些候选基因的功能还需要进一步探讨。下一步，将利用高密度测序数据对这些区域进行精细定位，以进一步确定与产仔数相关的基因并筛选出致因突变，为猪的遗传机制分析和育种提供参考。