猪产仔数的全基因组关联分析

付言峰，赵为民，李 辉，李碧侠，程金花，徐小波，任守文

(江苏省农业科学院 畜牧研究所；
江苏现代农业(生猪)产业技术体系阜宁推广示范基地， 江苏省农业种质资源保护与利用平台，农业部种养结合重点实验室，江苏 南京 210014)

产仔数是养猪生产中一个重要的育种性状、繁殖性状和经济性状[1-2]，而总产仔数和产活仔数一直是种猪选育中最常记录的2个产仔数指标。由于产仔数遗传力较低且又是伴性性状[3]，利用常规的育种方法选育难度较大，且进展缓慢[4]。随着分子生物学的发展，全基因组SNPs分子标记逐渐应用于猪产仔数的选择育种中。

利用全基因组SNPs分子标记，可以进行全基因组关联分析(Genome-wide association study，GWAS)[5]和全基因组选择育种(Genomic selection)[6]。GWAS是利用全基因组范围内筛选出的高密度分子标记(现大多为SNP)进行试验群体的全基因组扫描，再通过统计学方法将扫描所得的分子标记基因型与表型性状关联分析，最后鉴定出影响表型性状的相关分子标记和候选基因[7]，近年来已在奶牛和生猪育种生产上得到广泛的应用。其中在猪上，研究人员利用全基因组关联分析技术在嵊县花猪、金华猪群体分别筛选了160，124 个达到显著性水平的SNP[8]，还利用全基因组SNPs芯片在杜洛克猪中筛选了乳腺癌抗雌激素耐药基因3(BCAR3)、N-乙基马来酰亚胺敏感因子(NSF)等初情期日龄潜在的候选基因[9]。

本研究将利用Illumina全基因组猪50K SNPs芯片检测了186头加系大白猪的基因组，并与这些猪的总产仔数、产活仔数性状进行了关联分析，进一步对获得的显著关联SNPs染色体物理位置100 kb范围的基因进行筛选挖掘，最后对这些候选基因进行GO功能和KEGG通路注释，拟解决试验猪高产仔数关键位点和基因的筛选，研究结果将对进一步揭示猪高产仔性状的遗传基础，提高猪育种效率和经济效益具有较大的借鉴意义。

1.1 试验动物和样品采集

试验动物来自江苏天兆实业有限公司。首先，在1 800头加系大白猪基础群母猪中选取了300头组成核心群，其次，又在核心群中挑选了186头能繁母猪为研究对象，按猪场常规进行饲喂和管理，同时在自由采食及饮水条件下，注意进行疾病的防治和行为体况的观察。分2次进行耳组织样的采集，然后置于干冰中冷冻保存，回到实验室后置于-20 ℃冰箱保存，用于提取基因组DNA。

1.2 DNA提取和全基因组SNPs芯片检测

DNA提取方法：将约2 mg 猪组织样剪碎，加入600 μL 组织DNA 提取液(组织裂解液)，上下颠倒10次混匀；
加入蛋白酶K(20 mg/mL)6 μL，使其至终浓度200 μg/mL，50～55 ℃水浴锅内消化过夜；
加入等体积(606 μL)苯酚抽提，12 000 r/min 离心10 min，取上清(≤700 μL)；
加入等体积苯酚∶氯仿(1∶1)抽提，12 000 r/min 离心10 min，取上清(≤700 μL)；
加入等体积氯仿抽提，12 000 r/min 离心10 min，取上清(≤500 μL)；
加入2倍体积的无水乙醇沉淀DNA，12 000 r/min 离心10 min，弃上清。用至少500 μL 70%乙醇快速洗涤沉淀，重复2次(为防止DNA溶解，此步操作应快速进行)；
常温晾干，用30，100～200 μL，1 000～2 000 μL灭菌双蒸水(或TE溶液)溶解DNA；
0.7%的琼脂糖电泳初步检测浓度和纯度，紫外分光光度仪测定DNA 含量。

DNA经过质控检测合格之后，进行Illumina全基因组猪50K SNPs芯片检测，质控后共得到1.75 Gb的有效数据，用于后期的变异检测。

1.3 全基因组关联分析(Genome-wide association studies，GWAS)

本研究使用GEMMA软件[10]中的混合线性模型，进行猪产仔表型性状与全基因组SNPs之间的GWAS分析，GWAS分析模型为：Y=μ+Xb+u+e，其中，Y表示加系大白猪产仔性状表型值；
μ表示产羔数平均值；
X表示固定效应矩阵，b为固定效应向量；
u表示剩余多基因效应；
e表示产仔表型值的随机残差，u和e均服从正态分布。

本研究对湖羊GWAS 结果中的差异显著性水平(P值)采用Bonferroni 校正，并使用R软件中的”qqman”[11]进行曼哈顿图和Q-Q图(Quantile-Quantile Plot)绘制。

1.4 候选基因筛选

利用Ensembl网站(http：//www.ensembl.org/index.html)的大白猪基因组序列Sscrofa 11.1(GCA_001700135.1)进行每个SNP前后40 kb范围内定位的候选基因筛选，再利用其基因注释工具和NCBI网站(http：//www.ncbi.nlm.nih.gov/)Pubmed文献数据库进行基因的功能注释。

1.5 基因型质控

凡是符合以下3个条件中的任何一个，基因型质控过程将删除该个体或位点：分型缺失数>10%的位点、最小等位基因频率低于0.01的位点、哈代-温伯格检验P值<0.000 001 的位点。

2.1 表型数据统计

选用的186头母猪均有详细的产仔数记录，这些记录包括：耳号、样品号、猪只类型、父亲耳号、出生日期、出生体质量、乳头数、达100 kg日龄、备注、配种日期、产仔日期、总产仔数、产活仔数、健仔数、断奶仔猪数、出生窝质量、出生个体均质量。这些母猪均含有第1，2胎产仔数记录，最高胎次为第8胎(表1)。

表1 试验猪产仔数记录Tab.1 Litter size records of experimental pigs

2.2 全基因组关联分析

在经SNP质控之后，在18对(36条)常染色体上一共得到36 867个SNP标记，用于试验猪产仔数(总产仔数、产活仔数)的全基因组关联分析研究。这些加系大白猪所有胎次总产仔数共显著关联到2个SNP，分别为seq-rs323899658和seq-rs329781338，其中seq-rs323899658 SNP没有注释到基因，而seq-rs329781338 SNP又注释到3个基因，分别为SSBP1(线粒体 DNA 单链结合蛋白基因)、WEE2和KIAA1147(表2)。

表2 猪总产仔数显著关联SNP的候选基因注释结果Tab.2 Genes associated with total number born(TNB)of pigs

这些试验猪产活仔数共显著关联到7个SNP，分别为：seq-rs81238474、seq-rs321377412、seq-rs340736313、seq-rs80782154、seq-rs81244816、seq-rs81467772和seq-rs329781338。这7个SNP中有5个SNP有基因注释，2个没有基因注释；
这5个有注释的SNP共注释到22个基因，分别为TAS2R60、EPHA1、FAM131B、CLCN1、CASP2、TMEM139、TRBV21OR9-2、PRSS2、TRBV19、TRBV24-1、U6、SSBP1、WEE2、KIAA1147、NKX2-8、ALKBH3、HSD17B12、U6、HOMER2、WHAMM、FSD2和SCARNA15(表3)。

表3 猪产活仔数显著关联SNP的候选基因注释结果Tab.3 Genes associated with number born alive(NBA)of pigs

加系大白猪所有胎次总产仔数和产活仔数的GWAS结果如图1，2所示，从曼哈顿图结果(图1-A、图2-A)来看，总产仔数有2个SNP达到Bonferroni校正的阈值线，产活仔数有7个SNP达到Bonferroni校正的阈值线，表明加系大白猪总仔数和产活仔数分别有2，7个显著的SNP位点。从总产仔数和产活仔数的QQ图结果来看(图1-B、图2-B)，总产仔数有2个P值的观测值明显超过期望值，剩下P值的观测值都低于期望值或与其相当。产活仔数有较多的P值观测值都高于期望值，剩下的2个值相当，基本没有低于期望值的观测值。由此可见，QQ图与曼哈顿图结果一致。

图1 猪总产仔数曼哈顿图(A)和QQ-Plot图(B)Fig.1 Manhattan chart(A)and QQ plot chart(B)of total number born(TNB)in pigs

图2 猪产活仔数曼哈顿图(A)和QQ-Plot图(B)Fig.2 Manhattan chart(A)and QQ plot chart(B)of number born alive(TNA)in pigs

2.3 候选基因功能富集分析

在ENSEMBL网站上下载家猪的参考基因组信息，将显著SNP位点注释到距离其100 kb以内的基因上，总产仔数的2个显著SNP注释到3个基因，产活仔数的7个SNP注释到21个基因。然后使用clusterProfiler包将注释到的候选基因根据GO(Gene Ontology)和KEGG数据库进行基因功能富集分析。

猪总产仔数候选基因GO富集到的生物学过程一共42个(P<0.05)按显著性由高到低(P值由低到高)分别为：多生物过程的调控、减数分裂细胞周期的负调控、减数分裂核分裂的调控、参与卵母细胞成熟的减数分裂、细胞周期过程的调控、参与卵母细胞成熟的减数分裂细胞周期过程的调控、卵母细胞成熟、卵母细胞发育、卵母细胞分化、生殖过程的负调控和减数分裂细胞周期的调控。这些生物学过程均来自WEE2基因(图3-A)。

猪总产仔数候选基因GO富集到的分子功能一共5个(P<0.01)，按显著性由高到低(P值由低到高)分别为：非跨膜蛋白酪氨酸激酶活性、跨膜受体蛋白酪氨酸激酶活性、蛋白酪氨酸激酶活性、跨膜受体蛋白激酶活性和镁离子结合，这些分子功能均来自WEE2基因(图3-B)。猪总产仔数候选基因生物通路一共5条(P<0.05)，按显著性由高到低(P值由低到高)分别为：错配修复、DNA复制、同源重组、细胞周期和人类免疫缺陷病毒1型感染。这些生物通路均来自WEE2和SSBP1基因(图3-C)。

猪产活仔数候选基因GO富集到的分子功能一共3个(P<0.01)，按显著性由高到低(P值由低到高)分别为：跨膜受体蛋白酪氨酸激酶活性、蛋白酪氨酸激酶活性和跨膜受体蛋白激酶活性。这3种富集到的分子功能均同时来自EphA1和WEE2基因(图4)。

3.1 群体分层分析

GWAS研究是基于全基因组范围内不同遗传变异位点(多为SNPs)之间会出现连锁不平衡的现象来确定影响某些表型性状或数量性状的基因[12]，而影响连锁不平衡的因素很多，如遗传漂变、群体分层等，群体分层会导致GWAS结果假阳性[12-13]。根据本研究得到的QQ图结果可以看出，P值观察值和预测值基本相同，说明试验所用的加系大白猪样本未出现群体分层现象，本研究所得到的GWAS结果是可靠的，所使用的混合线性模型是合理的。

A.生物学过程；
B.分子功能；
C.KEGG通路。

A.Biological process；
B.Molecular function；
C.KEGG pathway.

图4 猪产活仔数候选基因GO功能富集分析-分子功能Fig.4 Go function enrichment analysis of candidate genes associated with total number born(TNB) in pigs，molecular function

3.2 全基因组关联分析

利用SNP 生物芯片技术和全基因组关联分析(GWAS)技术能够使用较低的成本和较高的准确率筛选出与繁殖性状关联性强的遗传标记，从而进一步筛选出相应的候选基因，为后期的功能验证以及选种选育、培养更高繁殖力的商品猪种打下基础[8]。GWAS曼哈顿图显示，本研究所筛选到的所有SNP都均匀分布，QQ图也显示，本研究选择的混合线性模型是适合的，前面的研究结果也佐证了这一结论[12]。

3.2.1 总产仔数 本研究所用的加系大白猪的总产仔数性状鉴定到2个显著性的SNP位点，其中1个SNP又注释到3个基因SSBP1、WEE2和KIAA1147，而这3个基因中只有WEE2基因GO功能富集到生物学过程和分子功能，只有WEE2和SSBP1基因KEGG通路能富集到，且前3位分别为错配修复、DNA复制和同源重组。有研究表明，WEE2基因(WEE1同源体，也称为WEE1B)属于母源基因，仅在卵巢的卵母细胞中高表达。在小鼠GV 期卵子中降调WEE2基因表达后，小鼠卵母细胞能够完成减数第一次分裂后排出极体，由MI 卵子发育至PB1 卵子，但通过IVF受精后，卵子不能形成原核，表现为受精障碍[14]。由此可见，WEE2基因在卵母细胞受精过程中发挥至关重要的作用。WEE2基因纯合突变c.G585C(p.K195N)纯合错义突变将导致受精失败和女性不孕[15]。本研究结果表明，WEE2基因与加系大白猪总产仔数相关。

另一个基因SSBP1(又称mtSSB)，中文名称为线粒体 DNA 单链结合蛋白，具有促进线粒体 DNA复制和损伤修复进而调控线粒体生成的重要作用[16]，近年来，SSBP1在肿瘤中作用受到越来越多的关注，有研究表明，SSBP1在卵巢癌组织中的表达显著低于正常卵巢组织[17]。此基因的生物学过程、分子功能和信号通路与总产仔数的紧密相关。

3.2.2 产活仔数 本研究所用的加系大白猪的产活仔数性状鉴定到7个显著关联SNP，其中5个SNP有基因注释，共注释到22个基因，这些基因GO功能富集到3种分子功能，跨膜受体蛋白酪氨酸激酶活性、蛋白酪氨酸激酶活性和跨膜受体蛋白激酶活性，对应的基因为EphA1和WEE2。GO功能中的生物学过程、细胞学分区和KEGG通路均未富集到。

目前已有许多研究表明，EphA1(促红细胞生成素产生肝细胞受体A1)和ephrin A1(促红细胞生成素产生肝细胞配体A1)在人[18]、小鼠和猪[19]的胚胎附植过程中发挥着重要的作用。EphA1基因第15 外显子上有一个G→C的突变，且该SNP突变的GG基因型可作为大白猪和长白猪中的优势基因型[20]。猪胚胎附植前、中、后期ephrin A1和Eph A1在其子宫内膜附植点mRNA表达量均呈先升高后降低的趋势，且妊娠猪表达量极显著高于空怀猪(P<0.01)[21]，EphA1的CpG岛甲基化很可能会影响此基因的表达[22]。梅山猪胚胎附植前、中、后期ephrin A1和Eph A1在其子宫内膜附植点mRNA和蛋白表达量均呈先升高后降低的趋势[23]。另外，胚胎附植期在人的子宫内膜上皮细胞和子宫内膜癌细胞中均能检测到Eph A1的mRNA和蛋白表达[18]。

综上所述，加系大白猪所有胎次母猪总产仔数和产活仔数的全基因组范围内的显著关联SNP分别为seq-rs329781338和seq-rs81244816，这2个SNP注释到的功能候选基因分别为WEE2和EphA1。本研究结果可为下一步提高猪产仔数和产仔数的基因编辑精准育种提供一定参考。