一种快速挖掘鸡品种特征性SNP标记集合的方法

白 露,王梦杰,马小春,何政肖,谭晓冬,刘 杰,赵桂苹,文 杰,刘冉冉

(中国农业科学院北京畜牧兽医研究所 畜禽营养与饲养全国重点实验室 农业农村部动物遗传育种与繁殖(家禽)重点实验室,北京 100193)

我国是世界上地方鸡种资源最丰富的国家之一,《国家畜禽遗传资源品种名录(2021版)》显示,我国现有地方鸡品种115个。地方鸡种大多具有外貌特征多样、适应性强、肉质风味独特、蛋品质优良等特点[1-3],符合我国传统消费习惯,为培育地方特色肉鸡和蛋鸡新品种提供了丰富的育种素材。然而,大量品种的精准鉴定和保种方法仍然有较大的提升空间。利用特征性分子标记对地方品种和商业化品系进行精准标识,对推动畜禽种质资源保护和商业化利用具有重要意义。

单核苷酸多态性(single nucleotide polymorphism,SNP)位点作为第三代遗传标记与其他分子标记相比具有数量多、分布广泛等优越性[4],已有技术可对其进行快速和规模化筛查,进而实现基因分型[5]。随着二代基因组测序成本降低,全基因组重测序成为种质资源研究[6]、群体进化[7]、基因组育种[8]等研究的常规技术方法,可以挖掘到大量目标性状相关的SNP标记[9-10]。然而,针对家禽具有地方品种和专门化品系众多的特点,需要筛选数量较少的SNP集合进行品种/品系区分,建立简便快速的鉴定方法,辅助种质资源保护和鉴定工作。

目前,中低密度SNP标记集合的检测方法主要包括SNP固相芯片、基于靶向SNP标记集合检测的液相芯片、基于质谱原理的SNP标记集合检测方法等[11-14]。中国农业科学院北京畜牧兽医研究所已研发55K SNP芯片“京芯一号”[15],中国农业大学研发50K SNP芯片 “凤芯一号”[16],江苏省家禽科学研究所研发23K液相芯片“酉芯一号”,山东省农业科学院家禽研究所研发11K液相芯片“鲁芯一号”等,主要服务于经济性状功能基因和分子标记挖掘[17-19]、基因组育种工作[20-21]和种质资源鉴定[22]。

群体分化指数(fixation index,Fst)是检测群体受到自然或人工选择基因组变异的常用方法,是群体间分化程度的衡量指标,可筛选受到选择压力影响的SNP标记[23-26]。连锁不平衡(linkage disequilibrium,LD)分析常应用于独立SNP提取[27-29],快速型白羽肉鸡父系LD衰变距离390 kb,京星黄鸡专门化品系LD衰变距离129 kb,不同品系LD差异较大,可通过各品种/品系的LD分析对SNP标记进行缩减[16,30]。

综上,本研究通过群体分化指数分析和连锁不平衡分析获得少量SNP标记,可以将目标品种与其它代表性品种区分,从而建立目标品种特征性SNP标记集合,为挖掘鸡品种/品系特征性SNP标记集合提供快速有效的方案。

1 材料与方法

1.1 试验数据

试验数据选取来源于中国农业科学院北京畜牧兽医研究所的北京油鸡群体(BJY AC (n=59))、BJY E (n=40)、BJY F (n=40))、京星黄鸡选育系D2系(JXH.D2 (n=49))、京星黄鸡选育系H系(JXH.H (n=59))、茶花鸡(CH (n=30))、大围山微型鸡(DWS (n=24))、武定鸡(WD (n=21)、藏鸡(ZJ (n=10))、大骨鸡(DG (n=7))和瓢鸡(P(n=21))重测序数据;来源于佛山高明区新广农牧有限公司的快速型白羽肉鸡(B (n=60))重测序数据;来源于山东农业科学院家禽研究所的7个山东地方鸡品种A、B、BRG、BRM、D、L和S (SD A (n=10)、SD B (n=10)、SD BRG (n=10)、SD BRM (n=10)、SD D (n=20)、SD L (n=80)、SD S (n=20))重测序数据。试验数据共涉及19个品种/品系共580个个体。

1.2 数据质控

基于10x以上的全基因组重测数据,利用PLINK(V 1.90)[31]软件对SNP进行标准的质量控制,删除缺失率>0.1的个体、删除缺失率>0.1的SNP以及删除次等位基因频率<0.05的个体(--mind 0.1--geno 0.1--maf 0.05),并将19个品种品系测序数据合并为一个37.14 G数据量的vcf文件。保留1～28号染色体上的位点,共保留16 927 197个SNPs标记位点用于后续分析。

1.3 主成分分析和系统进化树

使用GCTA 64(V 1.93.2)[32]软件构建亲缘关系矩阵(--make-grm)后计算PCA,计算每个主成分解释百分比,选择前2个主成分,用RStudio(V1.1.463)绘制主成分分析(principal component analysis,PCA)平面图。

随机提取19个品种/品系各10个个体,大骨鸡7个个体,共187个个体,利用PLINK(V1.90)软件对SNP频率构建遗传距离矩阵(--distance-matrix)。通过MEGA(V 7.0.26)[33]软件,采用领接法(Neighbour Joining,NJ)绘制进化树。

1.4 群体分化指数分析

以1个品种/品系作为目标品种/品系,利用VCFTools(V0.1.13)[34]软件计算采用1对N的方式进行Fst分析,以1 kb为窗口大小、1 kb为步长计算SNP单点Fst值(--fst-window-size 1--fst-window-step 1)。

1.5 连锁不平衡分析

利用PLINK(V1.90)软件提取群体分化指数分析结果中MEAN_FAST≥0.65 SNP位点形成SNP标记集合,对SNP标记位点进行LD分析(--blocks no-pheno-req)。提取全部非LD SNP标记位点以及每个LD中1个SNP标记位点,作为独立SNP标记位点。

2 结 果

2.1 遗传结构分析

对19个品种/品系质控后16 927 197个SNPs进行PCA(图1a)和NJ进化树(图1b)分析,结果表明北京油鸡、快速型白羽肉鸡品系、京星黄鸡H系和京星黄鸡D2系与其它品种遗传距离较远,分层明显。茶花鸡和大围山微型鸡聚成一支,武定鸡和瓢鸡聚成一支,山东地方品种/品系聚在一起。

a.多品种主成分分析;b.NJ进化树结果。B.快速型白羽肉鸡;BJY AC、BJY E、BJY F.北京油鸡;CH.茶花鸡;DG.大骨鸡;DWS.大围山微型鸡;JXH.D2.京星黄鸡选育系D2系;JXH.H.京星黄鸡选育系H系;P.瓢鸡;SD A、SD B、SD BRG、SD BRM、SD D、SD L、SD S.山东地方鸡品种品系A、B、BRG、BRM、D、L、S;WD.武定鸡;ZJ.藏鸡

2.2 目标群体品种/品系特征性SNP标记集合筛选

根据遗传结构分析结果挑选独立于其它群体的北京油鸡、快速型白羽肉鸡、京星黄鸡H系和京星黄鸡D2系,分别通过单位点Fst分析和LD分析筛选特征性SNP标记集合。以快速型白羽肉鸡为目标品种进行单位点Fst分析(图2a),结果表明,与其他品种/品系显著差异的SNP标记主要位于1、5、18和28号染色体上。提取MEAN_FAST≥0.80共346个SNPs标记进行PCA分析,可将快速型白羽肉鸡与其它群体分开(图2b)。对346个SNPs标记进行LD分析,提取所有不连锁的SNPs和每个LD中1个SNP标记,共114个SNPs标记进行群体PCA分析,结果表明114个SNPs标记可将快速型白羽肉鸡与其它品种/品系分开(图2c)。

a.快速型白羽肉鸡vs.其它品种/品系群体分化指数分析结果;b.选择性清除分析筛选MEAN_FAST≥0.80 SNP标记主成分分析结果;c.连锁不平衡分析筛选的SNP标记主成分分析结果

以京星黄鸡选育系H系为目标品系进行单位点Fst分析(图3a),结果表明,与其他品种显著差异的SNP标记主要位于1、2、3和4号染色体上。提取MEAN_FAST≥0.76共356个SNPs标记进行PCA分析,可将京星黄鸡选育系H系与其它群体分开(图3b)。对356个SNPs标记进行LD分析,提取所有不连锁的SNPs和每个LD中1个SNP标记,共220个SNPs标记进行群体PCA分析,结果表明220个SNPs标记可将京星黄鸡选育系H系与其它品种分开(图3c)。

a.京星黄鸡选育系H系vs.其它品种/品系群体分化指数分析结果;b.选择性清除分析筛选MEAN_FAST≥0.76 SNP标记主成分分析结果;c.连锁不平衡分析筛选的SNP标记主成分分析结果

以京星黄鸡选育系D2系为目标品系进行单位点Fst分析(图4a),结果表明,与其他品种/品系显著差异的SNP标记主要位于1、2、4、7、14和23号染色体上。提取MEAN_FAST≥0.76共321个SNPs标记进行PCA分析,可将京星黄鸡选育系D2与其它群体分开(图4b)。对321个SNPs标记进行LD分析,提取所有不连锁的SNPs和每个LD中1个SNP标记,共226个SNPs标记进行群体PCA分析,结果表明226个SNPs标记可将京星黄鸡选育系D2系与其它品种/品系分开(图4c)。

a.京星黄鸡选育系D2系vs.其它品种/品系群体分化指数分析结果;b.选择性清除分析筛选MEAN_FAST≥0.76 SNP标记主成分分析结果;c.连锁不平衡分析筛选的SNP标记主成分分析结果

2.3 遗传距离较近品种/品系特征性SNP标记集合筛选

根据遗传结构分析结果挑选聚集在一个分支的武定鸡和瓢鸡,分别通过单位点Fst分析和LD分析筛选特征性SNP标记集合。以武定鸡为目标品种进行单位点Fst分析(图5a),结果表明,与其他品种显著差异的SNP标记主要位于1、2、4、5和15号染色体上。提取MEAN_FAST≥0.70共368个SNPs标记进行PCA分析,可将武定鸡与其它群体分开(图5b)。对368个SNPs标记进行LD分析,提取所有不连锁的SNPs和每个LD中1个SNP标记,共204个SNPs标记进行群体PCA分析,结果表明204个SNPs标记可将武定鸡与其它品种分开(图5c)。

a.武定鸡vs.其它品种/品系群体分化指数分析结果;b.选择性清除分析筛选MEAN_FAST≥0.70 SNP标记主成分分析结果;c.连锁不平衡分析筛选的SNP标记主成分分析结果

以瓢鸡为目标品种进行单位点Fst分析(图6a),结果表明,与其他品种显著差异的SNP标记主要位于1、2和4号染色体上。提取MEAN_FAST≥0.65共1 178个SNPs标记进行PCA分析,可将瓢鸡与其它群体分开(图6b)。对1 178个SNPs标记进行LD分析,提取所有不连锁的SNPs和每个LD中1个SNP标记,共178个SNPs标记进行群体PCA分析,结果表明178个SNPs标记可将瓢鸡与其它品种分开(图6c)。

a.瓢鸡vs.其它品种/品系群体分化指数分析结果;b.选择性清除分析筛选MEAN_FAST≥0.65 SNP标记主成分分析结果;c.连锁不平衡分析筛选的SNP标记主成分分析结果

3 讨 论

目前SNP标记集可以基于单倍型分析[35]、连锁不平衡分析[36]、全基因组关联分析[37]等方法获得。Judge等[38]基于有系谱和重测序数据的大群体以及中等密度SNP芯片的大群体,利用Delta统计、Fst统计、结合Delta统计和成对Fst值的索引进行计算等方法识别信息量最大的SNPs,通过300个以上SNPs标记精准量化生物样本中安格斯牛和赫里福德牛的比例。Seo等[39]基于鸡600K SNP芯片,GWAS分析后,对病例组和对照组进行LD修剪,得到96个SNPs标记可以将目标鸡与其它鸡群分开。这些基于系谱与GWAS分析的研究可进行少量群体特征性SNP集合的挖掘,但是均需要较大群体的表型信息与测序,而本研究方法仅需要试验群体代表性个体的重测序数据,一般30个左右个体可代表一个品种的遗传多样性,数量较少的个体可作为背景,在试验材料和数据准备上更为简便。

Bertolini等[40]基于奶牛大群体Bovine SNP50 v1 BeadChip 芯片,利用基于Delta、Fst、PCA-chrom 和 PCA-whole等技术,通过品种分配和随机森林筛选出96个SNPs组成的SNP-set可以将品种区分开。Schiavo等[41]基于猪大群体PorcineSNP60 BeadChip 芯片,保留LD分析中所有r2<0.3的SNP和1个LD中任一SNP,然后使用Delta、Fixation指数、主成分分析统计和两种随机森林分类方法筛选到96个SNPs标记位点可进行目标品种鉴定。Cho等[42]基于鸡大群体600K SNP芯片,通过GWAS和LD分析筛选,得到初步的SNP标记集合,然后通过随机森林(RF)和AdaBoost(AB)两种机器学习算法,筛选到Yeonsan Ogye鸡群的38(RF)和43(AB)个共81个最佳SNPs标记集合,在品种区分上显示了100%的准确性。Kumar等[43]基于小等位基因频率连锁不平衡的方法,鉴定到591个品种特异性SNPs组成的集合,适用于鉴别牛的亲缘关系的分配。Gao等[44]基于24个猪品种的62 822个SNPs基因型文件,通过LD、PCA、随机森林及相应的包外误差估计(OOB)和MDA筛选方法获得1 000个SNPs可将目标品种区分。利用随机森林等一系列方法可筛选出较少的品种特征性SNP标记,是未来发展的重点,但需要较强的方法学作为基础。而本方法仅通过一次单点群体分化指数分析和连锁不平衡分析,即可挖掘到114～226个SNPs标记将目标品种与其它代表性品种区分开,更为快捷。

4 结 论

本研究应用19个鸡品种/品系全基因重测序数据进行1对N的单位点群体分化指数分析,以MEAN_FAST≥0.65为筛选标准,对筛选得到的SNP标记进行连锁不平衡分析,在多个品种中确定了114～226个不同染色体上SNPs标记可以将目标品种与其它代表性品种区分开来,从而建立目标品种特征性SNP标记集合。该SNP标记集合筛选方法是实现低成本和快速品种鉴定的基础。