基于重测序数据的湖羊产羔数性状全基因组关联分析

谢锐，朱墨，童世锋，赵福平，刘杨*

(1. 南京农业大学动物科技学院，江苏 南京 210095； 2. 农业农村部都市农业重点实验室，上海 200240；3. 中国农业科学院北京畜牧兽医研究所，北京 100193)

繁殖性状是羊的重要经济性状，繁殖性能的高低直接决定了整个养殖经济效益的高低。如何提高羊产羔数一直是科学家们致力于研究解决的重要课题。提高母羊产羔数的前提是提高母羊的排卵数，而哺乳动物的排卵过程受到复杂的机制调控，而且不同品种羊之间还存在着显著差异，如特塞尔绵羊和萨福克羊每个排卵周期只排1个卵，而高产的布鲁拉美利奴羊每个排卵周期可以排10个卵[1]。除了国外的羊具有优秀繁殖性能外，我国本地绵羊品种湖羊也因具有“一年多胎，一胎多羔”的优秀繁殖性能而闻名世界。但影响湖羊产羔数的关键功能基因以及高繁殖力性状的遗传机理还有待进一步研究。

随着二代基因组测序技术的快速发展，越来越多的研究利用基因组测序方法来寻找与羊繁殖性状有关的候选基因，这为从基因组水平上研究目标经济性状提供了便利。全基因组关联分析(GWAS)是在全基因组范围内寻找与动植物重要经济性状相关联的遗传变异，是一种复杂性状功能基因鉴定的分析策略[2]。GWAS概念是由Risch在1996年提出[3]，最早是被应用在医学上研究复杂疾病的遗传模式[4]。GWAS方法因具有速度快、通量高和高精度等显著优点,其在畜牧学遗传育种中的作用日益凸显，目前已经成为挖掘与家畜重要经济性状相关候选基因最重要的方法之一[5]。狄江等[6]采用单标记回归方法对中国美利奴(新疆型)周岁细毛羊进行全基因组关联分析，检测到3个与羊毛产量和2个与羊毛长度显著相关的候选SNP位点，并根据这些SNP位点信息鉴定到3个与羊毛生长、发育相关的生物学过程相关基因，分别为IBIN、HSD17B11和PIAS1基因。王珊等[7]利用全基因组重测序数据对多浪羊产羔性状进行了全基因组关联分析，经过基因注释和富集分析筛选出了2个与多浪羊高繁殖力相关的基因。张军霞等[8]通过多态性位点与性状的关联分析发现小尾寒羊中高繁殖力群体产羔数与KITLG基因g.47468C>T位点显著相关，可作为产羔性状的候选基因。兰蓉等[9]利用全基因组关联分析方法对云南黑山羊产羔性状进行研究，发现2个与山羊产羔数显著相关的SNP位点。刘世佳等[10]利用关联分析的方法发现BMP15基因B2突变(C718T)与产羔数显著相关，可作为提高绵羊产羔数性状的分子标记应用于育种。本研究旨在通过重测序数据，结合GWAS分析技术对影响湖羊产羔数性状的候选基因进行挖掘，这有利于进一步揭示湖羊高繁殖性状的遗传基础，缩短湖羊的育种年限，降低育种成本，提高经济效益。

1 材料与方法

1.1 样本采集及测序

从甘肃某羊场随机挑选206只湖羊，统一采集耳组织样，置于-20℃的干冰中冷冻保存，用于提取基因组DNA。样品采集的DNA经过质控检测合格之后，统一送到北京贝瑞和康生物技术有限公司进行双端150 bp全基因组重测序，平均测序深度为6×。质控后共得到3 301 Gb的有效数据(clean reads)，用于后期的变异检测。

1.2 变异检测及SNP质控

使用绵羊4.0(Ovis_aries 4.0)参考基因组，并使用BWA[11](version 0.7.17)的mem模块测序得到的clean reads 与参考基因组进行比对；使用GATK4[12](version 4.0.2.1)中的Picard 模块标记PCR重复序列，变异检测流程使用GATK4软件中的Haplotype Caller模块和Variant Filtration模块进行SNP的鉴定和初步过滤，将结果输出到VCF文件保存。最后使用VCFtools[13](version 0.1.17)对SNP进行质控，并保留符合以下标准的SNP：①所有SNP位点的平均测序深度大于2×；②只保留常染色体上的SNP；③SNP最大缺失率5%；④只保留二等位基因；⑤最小等位基因频率大于0.05；⑥哈迪温伯格平衡P值大于1×10-6。

1.3 GWAS分析

使用EMMAX[14](version beta-07)中的混合线性模型，利用湖羊产羔数表型性状与全基因组上SNP来进行GWAS分析，混合线性模型如下：

y=μ+Xb+u+e，

其中，y表示湖羊产羔数性状表型值；μ表示产羔数平均值；X表示固定效应矩阵，b为固定效应向量；u表示剩余多基因效应；e表示产羔数表型值的随机残差，u和e均服从正态分布。

本研究中对湖羊GWAS结果中的P值采用Bonferroni校正，并使用R软件中的CMplot包进行曼哈顿图和Q_Q图(Quantile-Quantile Plot)绘制。

2 结果与分析

2.1 表型数据统计

选用的206只母羊均有第一胎产羔数数据，其中181只母羊有第二胎产羔数数据，其余25只母羊第二胎产羔数数据缺失。在湖羊第一胎产羔数数据中，产羔数最多的为4只，产羔数最少的为1只，平均产羔数为2.1只；在湖羊第二胎产羔数数据中，产羔数最多的为5只，产羔数最少的为1只，平均产羔数为2.36只；此外，还统计了湖羊2胎平均产羔数数据，2胎平均产羔数最多的为4只，产羔数最少的为1只。第二胎产羔数和2胎平均产羔数表型数据统计时已剔除数据缺失的个体。统计学结果显示湖羊第一胎平均产羔数和第二胎平均产羔数差异显著。

2.2 全基因组关联分析

经SNP质控之后，在26条常染色体上一共得到1 604 526个SNP标记用于全基因组关联分析研究。湖羊第一胎产羔数、第二胎产羔数和2胎平均产羔数GWAS结果分别如图1、图2和图3所示。从曼哈顿图结果来看，没有SNP位点达到Bonferroni校正的阈值线，表明湖羊产羔数的全基因组关联分析未发现显著SNP位点。从第一胎产羔数、第二胎产羔数以及2胎平均产羔数的Q_Q图结果来看，所有的P值的观测值都没有明显超过期望值，也未找到与湖羊产羔数显著相关的SNP位点，这与曼哈顿图结果一致。

图1 湖羊第一胎产羔数全基因组关联分析曼哈顿图(左)和Q_Q图(右)

图2 湖羊第二胎产羔数全基因组关联分析曼哈顿图(左)和Q_Q图(右)

图3 湖羊2胎平均产羔数全基因组关联分析曼哈顿图(左)和Q_Q图(右)

3 讨论

3.1 群体分层分析

GWAS研究是基于全基因组范围内不同位点之间会出现连锁不平衡的现象来确定影响某些表型性状或数量性状的基因，而影响连锁不平衡的因素很多，如遗传漂变、群体分层等，若群体出现分层现象会导致关联分析的结果出现假阳性而不可靠[15-16]。根据本研究得到的Q_Q图结果可以看出，P值观察值和预测值基本相同。同时本研究还根据湖羊不同胎次产羔数GWAS结果的P值计算了实际的膨胀系数，其中，第一胎产羔数GWAS结果的实际膨胀系数为0.999；第二胎产羔数GWAS结果的实际膨胀系数为1.007；2胎平均产羔数GWAS结果的实际膨胀系数为0.999，实际膨胀系数与预期膨胀系数(1)很接近，说明湖羊样本未出现群体分层的现象，本研究所得到的GWAS结果是可靠的，所使用的混合线性模型是合理的。

3.2 全基因组关联分析

通过查阅现有的文献得知，目前还未有与湖羊产羔数性状的GWAS研究报道。GWAS已经是应用非常普遍的功能基因筛选方法，但仍存在一些缺陷，比如检测结果中得不到显著性的SNP位点，或者得到的假阳性的显著性结果[17]。出现这样的结果可能是由于以下因素导致：①GWAS研究群体的规模小，导致检验功效较低；②SNP标记的密度低或者标记分布不均匀；③选择的统计模型不适合；④表型值检测不准确；⑤目标性状遗传力低且受微效多基因控制。

本研究使用的是湖羊二代重测序数据，SNP密度很高，从曼哈顿图中的SNP密度分布图可以看出，鉴定到的所有SNP都均匀分布，未出现分布不均匀的情况。Q_Q图结果也显示选择的混合线性模型是适合的。然而，湖羊产羔数的全基因组关联分析的结果未鉴定到显著性的SNP位点，推测可能由以下原因所导致：①湖羊的繁殖性状的遗传力低，只有0.125，并且繁殖性状是一种数量性状，受微效多基因控制[18]，导致单个基因的效应对湖羊整个产羔数贡献不大，达不到GWAS研究检测的显著性阈值，从而不能通过GWAS的分析方法挖掘出来；②GWAS分析中使用的群体规模应该越大越好，群体大小是制约GWAS分析检验功效的第一要素，群体越大有助于阳性结果的检出。