一種快速挖掘雞品種特征性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標記集合篩選方法是實現低成本和快速品種鑒定的基礎。