基于表型和基因组信息评价北京油鸡保种群保种情况

王海龙，王 巧，邢思远，王 杰，李庆贺，郑麦青，崔焕先，刘冉冉，赵桂苹，文 杰

(中国农业科学院北京畜牧兽医研究所 动物营养学国家重点实验室，北京 100193)

畜禽遗传资源是生物多样性的重要组成部分，在农业、经济和科学研究中发挥着重要作用。我国幅员辽阔，不同地理和气候条件孕育了丰富的畜禽遗传资源，约占世界畜禽遗传资源总量的六分之一。据2012年版《中国畜禽遗传资源志》统计，我国有777种动物遗传资源正式命名，包括556个地方品种，其中鸡品种116个，北京油鸡是其中之一。

北京油鸡起源于大约300年前的清代早期，以洼里和大屯两个地区最为集中。典型的北京油鸡具有凤冠、胫羽和五趾等外观特征，兼具地方鸡种肉、蛋品质优良，耐粗饲、抗病、抗逆性强等优良特性[1]，但生产性能相对于国外品种较差。70年代末，为满足市场需求，企业主推高产品种，北京油鸡数量骤降[2]。为保护北京油鸡这一地方鸡种，中国农业科学院北京畜牧兽医研究所等单位开始承担北京油鸡的品种保护工作。2001年北京油鸡被农业农村部列为国家级畜禽品种资源重点保护品种。

亲缘关系较近个体交配的过程称为近交。在畜禽育种中，一般使用近交系数评估个体之间的近交情况。近交系数是指当个体由于近交而杂合基因减少时，纯合基因或纯合子的百分比[3]。近交可以增加种群的纯合性，也会增加纯合致病基因的概率，导致后代繁殖力下降和表型衰退等情况的发生。群体数量是影响保种情况的重要因素，群体数量较小，更容易导致近交发生。

随着二代全基因组重测序技术的发展，测序成本不断降低，重测序技术越来越多的应用到畜禽保种工作中[4-6]。基于SNP标记计算得到基因组近交系数(genomic inbreeding coefficient based on SNPs，FSNP)，成为保种评价的新方法。使用SNP位点信息计算近交系数，可以有效避免系谱信息不完整、误差较大等特点[7]。连续性纯合片段(runs of homozygosity, ROH)[8]是基因组中的一段长纯合片段，长度为数百kb到数Mb不等。自Ferencakovic等[9]首次利用ROH进行近交评估以来，ROH已广泛应用于畜禽近交系数的计算。另有研究发现，FROH的准确性相较于FHOM、FGRM、FUNI计算的结果准确性更高[10]。

本研究以国家级北京油鸡保种场随机交配保种群体为研究对象，基于外貌特征和基因组信息对保种群体的外貌特征、近交系数和有效群体大小进行分析，研究北京油鸡的保种情况。

1 材料与方法

1.1 试验材料与样品采集

本试验随机选取国家级北京油鸡保种场2019年随机交配保种群体40只个体，其中公鸡17只，母鸡23只，采用完全随机交配。同时，对北京油鸡凤冠、五趾、胫羽等外貌特征进行整理，对近年外貌特征的变化趋势进行分析。其中，1979年和1991—1993年的外貌特征记录来自优质黄羽肉鸡品系选育和配套研究论文集(内部资料)，2001年、2005年外貌特征记录来源于现有群体的早期数据[11]。

翅下静脉采血，于EDTA抗凝管中，-20 ℃保存，用于提取DNA。

1.2 鸡基因组DNA提取

采用酚-仿传统法提取DNA，使用两种方法对血液DNA提取效果进行检测：1)琼脂糖凝胶电泳分析DNA降解的程度以及是否存在RNA污染；2)Nanodrop检测DNA的纯度，OD260 nm/OD280 nm合格范围为1.8～2.0。

合格DNA样品送往北京康普森生物技术有限公司进行二代全基因组重测序，基于Illumina NovaSeq技术测序平台，利用双末端测序(Paired-End)的方法，测序策略为Illumina PE150，测序深度为10×。

1.3 基因组变异检测

使用基因组比对软件BWA(v0.7.17)[12]中BWA-MEM算法，将过滤后的Clean Reads比对参考基因组(ftp://ftp.ensembl.org/pub/release-101/fasta/gallus_gallus/dna/)，并去掉未比对上、低质量(MQ<4，MQ为mapping的质量值)、重复的Reads。剩余的高质量Reads用于后续分析。

在比对到参考基因组序列的基础上，通过突变分析软件GATK(v4.1.7.0)[13]检测全基因组中所有SNPs位点，过滤得到高可信度的SNP数据。

1.4 测序数据质控

本研究使用PLINK(v1.90)[14]软件对全基因组测序数据生成的文件进行质量控制。质控标准如下：SNP检出率大于0.9，最小等位基因频率大于0.05，哈迪-温伯格平衡P值大于10-6，同时只保留常染色体，避免性别影响。经过质控后，共剩余40只个体的6 252 214个SNPs用于后续分析。使用BEAGLE.18May20.d20[15]软件对缺失基因型进行填充。

1.5 基于基因组信息估计有效群体大小(Ne)

本研究使用PLINK软件对质控数据进行连锁不平衡修剪。把r2<0.1默认为两个SNPs不连锁。经过连锁不平衡修剪后共剩余43 167个SNPs位点，再使用SNeP(v1.1)[16]软件估计北京油鸡历史世代的有效群体大小，公式[17]：

(1)

使用NeEstimator(v2.1)[18]软件的基于连锁不平衡方法估计当前世代的有效群体大小[19]。

1.6 基于基因组信息估计近交系数(F)

(2)

其中，l为组成ROH的最少SNP数目，a为设定的检测到假阳性ROH的百分率，本研究中设定为0.05，ns为个体的SNP数目，ni为样本数，het为SNP的平均杂合度。本研究中由此公式计算得出组成ROH的最少SNP数目为153个。

根据PLINK软件计算的结果，利用公式计算FROH：

(3)

其中∑iLROHi为常染色体上ROH的总长度，Lauto为常染色体的长度，FROH为基于ROH计算的近交系数。

使用PLINK软件对质控后数据计算基于纯合基因型的近交系数(FHOM)和基于联合配子之间相关性计算的近交系数(FUNI)，公式如下：

FHOM=(O-E)/(L-E)

(4)

其中，O为个体观测纯合基因型的数量，E为期望纯合基因型的数量，L为基因型的SNP数目。

(5)

其中,pi为第i个SNP位点等位基因A的初始频率，xi是基因型中a的个数，当SNP标记的基因型为AA、Aa或aa时，xi分别为0、1或2。

使用GCTA(v1.93.2beta)软件构建基因组关系G矩阵，再使用R(v3.6.3)软件对构建的G矩阵计算基于基因组关系G矩阵的近交系数(FGRM)[21]，公式如下：

(6)

其中，pi同上，xi同上，N为SNP的数目。

1.7 FROH、FHOM、FGRM和FUNI不同算法近交系数间的相关性分析

使用R的PerformanceAnalytics包对2019年北京油鸡保种群体FROH、FHOM、FGRM和FUNI4个近交系数进行相关性分析。

1.8 基于各世代实际留种数量估计近交系数

根据国家级北京油鸡保种场1979—2019年繁育记录，统计本保种群各世代实际留种公母鸡数，估计2019年北京油鸡保种群有效群体大小(Ne)、近交增量(ΔF)、近交系数(Ft)，公式[22]如下：

(7)

(8)

Ft=1-(1-ΔF)t

(9)

其中，Ne为有效群体大小，NS为公鸡，ND为母鸡，ΔF为近交增量，t为世代数，Ft为t世代近交系数。

2 结 果

2.1 北京油鸡保种群体外貌特征的变化情况

对1979年、1991—1993年、2001年、2005年、2017年、2018年和2019年的外貌特征进行统计，结果见图1。胫羽百分比始终在97%以上，表现稳定。凤冠、五趾的百分比在2005—2017年都有不同程度的下降，凤冠百分比从100%下降到93%，五趾百分比从27%下降到15%。2018年凤冠百分比上升到98%，五趾百分比上升到27%，并保持稳定。

图1 北京油鸡外貌特征统计图

2.2 基于各世代实际留种数量估计近交系数

根据国家级北京油鸡保种场1979—2019年各世代实际留种公、母鸡数量(共38个世代)，估计有效群体大小和近交增量，结果如表1所示。北京油鸡保种群38个世代平均近交增量为0.001 9，38个世代繁育后，2019年保种群近交系数为0.070 3。

表1 北京油鸡保种群有效群体大小、近交系数估计

2.3 基于基因组信息的北京油鸡保种效果分析

2.3.1 有效群体大小估计 基于连锁不平衡方法估计2019年北京油鸡保种群体的有效群体大小为193。使用SNeP估计北京油鸡保种群的有效群体大小，结果如图2所示。发现有效群体大小随着世代的减小逐渐呈现平缓下降的趋势，98世代之前的有效群体大小为595，13世代之前的有效群体大小176。北京油鸡保种群有效群体大小从45世代前开始下降加快，北京油鸡世代间隔大约为1年，70年代 刚好是北京油鸡群体数量开始减少的时期，与现实情况相吻合。

图2 北京油鸡历史世代有效群体大小变化图

2.3.2 ROH统计 对40只北京油鸡进行分析，共检测出7 101个ROH。分别统计ROH不同长度所占比例(图3A)、不同染色体ROH数量所占比例(图3B)、不同染色体上ROH的平均长度(图3C)和ROH数量与长度相关统计(图3D)。较短的ROH(0～0.5 Mb)所占比例最多，约占总数的75%，且ROH数量随ROH长度的增加逐渐减少。常染色体ROH数量各不相同，ROH分布不均匀。每条染色体ROH的数量大体随染色体长度的增加而增加。其中1号染色体ROH数量最多，约占ROH总数的19%，16、25、30、31、32号染色体ROH数量最少，而不同染色体上ROH的平均长度没有显著区别。在不同个体中，ROH的数量与ROH覆盖的长度也有很大的不同。随着ROH数量的增加，ROH的总长度也在增加。在这一群体中，ROH数量最多的个体有232个ROH，总长度约为110 Mb，ROH数量最少的个体有34个ROH，总长度约为8.7 Mb。

图3 ROH长度及分布统计图

2.3.3 基因组近交系数比较 比较基于ROH、纯合基因型、基因组关系G矩阵、联合配子之间的相关性4种不同算法下的基因组近交系数(FROH、FHOM、FGRM、FUNI)，结果如表2所示。2019年保种群FROH平均值为(0.079 8±0.017 1)，范围为0.009 1～0.115 5。3种基于SNP位点计算的近交系数FHOM、FGRM、FUNI平均值分别为(0.060 5±0.039 8)、(0.066 2±0.034 7)、(0.063 7±0.035 4)，范围分别为-0.021 8～0.176 1、-0.003 2～0.162 3、0.000 6～0.164 3，分布范围明显大于FROH结果，且3种基于SNP位点计算的近交系数离散程度较大。

表2 基因组信息不同算法计算的近交系数

2.3.4 不同算法近交系数相关性分析 计算2019年保种群体FROH、FHOM、FGRM和FUNI4种不同算法所得近交系数值间的相关性，结果如图4所示。从图4中上三角可以看出，FROH与FGRM显著相关(P<0.01)，相关系数为0.45。其他各组两两之间相关极显著(P<0.001)，相关系数都大于0.5。从图4中下三角可以看出，FHOM与FGRM，FHOM与FUNI以及FGRM与FUNI之间存在较高线性相关。

对角线：近交系数的频率分布直方图，图中曲线是百分位曲线。对角线上方：不同近交系数之间的Spearman相关系数。＊.P<0.05，＊＊.P<0.01，＊＊＊.P<0.001。对角线以下：不同近交系数之间的散点图，散点位置由上侧和右侧F值共同确定，图中曲线是趋势线

3 讨 论

3.1 外貌特征变化对保种工作的指导

北京油鸡保种群体胫羽、凤冠、五趾等品种特有外貌特征得到了良好保留。2005—2017年北京油鸡五趾、凤冠百分比有所下降，自2017年起，繁育过程中适当增加具有凤冠、五趾特征的留种比例，至2018年凤冠、五趾特征百分比上升，恢复到2017年前的水平，并保持稳定。对1979—2019年各世代的外貌特征百分比进行统计，这一数据为国家级北京油鸡保种场利用品种特有外貌特征指导保种工作提供了依据。

3.2 基因组信息对保种工作的指导

传统计算近交系数的方式是基于系谱记录估计的，但经常遇到系谱信息不完整和错误等问题，且忽略初代个体间近交系数也会影响计算准确性[23]。使用系谱信息是基于亲缘关系的血缘同源概率估计期望值，未考虑减数分裂过程中基因重组是随机事件，使用基因组信息估计的是个体间实际近交情况[4]，更能反映保种群体的真实状况。FHOM、FGRM、FUNI是基于状态同源估计的[4, 24]，无法区分血缘同源与状态同源，且某些个体近交系数为负值，因此，分别使用上述几种方法估计近交系数都是不够准确的。FHOM、FGRM、FUNI是通过单点计算得到，结果会受等位基因频率估算的影响，而FROH直接反映了自身长纯合片段占参考基因组长度的比例，受DNA测序样品的质量影响更小[25]，利用ROH估计近交系数可以解决其他近交系数估计存在的缺陷[26-28]。同时，在人类[29-30]、牛[31-32]、猪[33-35]研究中也发现，基于ROH估计近交系数准确性最高。Peripolli等[36]的研究中，系谱计算的近交系数与FROH相关性最强，提出FROH可以作为代替系谱评估保种群体近交系数的有效方式。

本研究中，FROH的平均值大于3种基于SNP位点计算的近交系数FHOM、FGRM、FUNI，可能是由于FHOM、FGRM、FUNI中某些个体值为负数导致的近交系数偏低，造成的结果不准确，这一结果与Alemu等[37]和Shi等[33]的研究一致。相关性分析发现，FHOM与FGRM，FHOM与FUNI以及FGRM与FUNI具有较高的线性相关，且都是基于单点计算的，与曹愉夏等[38]的研究一致。由于FROH是基于血缘同源的估计，而FHOM、FGRM、FUNI是基于状态同源的估计，FROH与FHOM、FROH与FGRM、FROH与FUNI呈中度相关，在近交系数估计时尽量避免使用FHOM、FGRM、FUNI3种计算方法。

4 结 论

北京油鸡1979—2019年40年来外貌特征明显且百分比保持稳定，经过40余年保种工作后FROH为0.079 8，近交系数增长缓慢，说明国家级北京油鸡保种场随机交配保种群体的保种工作是切实可行的。将来对北京油鸡随机交配保种群体每年随机选取一定数量的个体进行二代全基因组重测序检测，有利于以后对保种状况进行动态监控，随时调整保种工作方案。