基于简化基因组测序的娄门鸭遗传多样性评价

刘宏祥，沈永杰，张丽华，章双杰，王 靖，朱 杰，陈瑜哲，朱春红，宋卫涛，张 丹，陶志云，徐文娟，刘红林，李慧芳*

(1. 江苏省家禽科学研究所，扬州 225125； 2. 昆山市畜牧兽医站，苏州 215300；3.南京农业大学，南京 210095)

遗传资源是种业科技原始创新、现代种业持续发展的源泉以及国家发展的重要战略资源。我国鸭遗传资源丰富，现有品种数量约占世界的48%左右。其中娄门鸭是江苏省著名的地方品种，原产于昆山市，以肉质优良而闻名。苏州历史悠久的熟食鸭制品大多均选用苏州当地产的膘肥、体大、肉嫩、皮白的娄门鸭为原料。娄门鸭历史上饲养量较大，解放前年产量达到100多万只，远销四周县市。由于受到外来品种的冲击，自上世纪60年代开始饲养量逐年降低，从2008年开始昆山市农业农村局与江苏省家禽科学研究所合作进行娄门鸭抢救性保护工作。从2013年开始新建娄门鸭保种场进行系统保种以来，尚未全面评估娄门鸭的保种效果。因此对娄门鸭的遗传多样性的评价工作迫在眉睫。

遗传多样性的评估可靠性依赖于遗传标记的质量。目前DNA分子遗传标记可以从DNA水平上直接反映个体之间的差异，稳定性高、多态性好，因此是目前最常用、应用范围最广的一类有效标记。其中代表性的DNA分子遗传标记有SSR(简单重复序列的微卫星标记)和SNP(单核苷酸多态性)。微卫星标记由于多态性较高，在遗传多样性的早期研究中得到了大量应用。虽然单个SNP的多态性较微卫星低，通常为二等位标记，但全基因组范围内的SNPs数量多、分布广，检测稳定性和判定准确率都较高，易实现标准化和自动化检测，因此SNP检测方法在遗传多样性评估中越来越受到重视。随着二代高通量技术的发展，基于酶切的简化基因组技术能够得到全基因组范围内足够数量的SNPs位点，且成本较基因组重测序大幅降低，因此成为分析群体遗传多样性和亲缘关系，评估种群近交系数的理想方法，并在不同物种上都得到了大量应用。

本研究通过简化基因组测序技术检测娄门鸭保种群全基因组范围内的SNPs标记，系统分析基于SNPs标记的群体遗传多样性，鉴定保种群的群体结构，评估娄门鸭群体的保种效果，为后续娄门鸭这一宝贵遗传资源保护方案的修正提供科学合理的指导和建议。

1 材料与方法

1.1 试验动物及样品采集

娄门鸭保种群来自江苏省昆山市娄门鸭原种场，所有个体饲养标准和环境条件一致。10周龄时随机选取同一世代163只鸭(39只公鸭，124只母鸭)，翅静脉采血后血样冷冻保存。

1.2 试验方法

1.2.1 基因组DNA的提取 采用酚-氯仿法提取血样基因组DNA。冻存后的每管血样解冻后吸取100 μL，加入900 μL SDS裂解液和20 μL蛋白酶K后摇匀，过夜消化后，加入900 μL酚、氯仿、异戊醇(25∶24∶1)混合液萃取基因组DNA，加入无水乙醇后沉淀基因组DNA，最后加入TE缓冲液溶解基因组DNA，-20 ℃保存备用。

1.2.2 基因组DNA的质量鉴定 通过核酸蛋白浓度测定仪测量基因组DNA的浓度和OD值，采用琼脂糖凝胶电泳检测基因组DNA的完整性。质量合格后进行下一步简化基因组测序。

1.2.3 简化基因组测序 所有基因组DNA样品均委托广州基迪奥生物科技有限公司进行简化基因组测序。具体试验流程为：选择R Ⅰ和Ⅲ限制性内切酶对基因组DNA酶切后加上带barcode的测序接头；混合样本后，构建小片段文库(250～550 bp)，进行PE125双末端测序。

1.3 数据处理

1.3.1 数据质控与基因组比对 为了保证测序数据的质量，首先对下机的clean reads进行过滤：1)去除含adapter的reads；2)去除含N比例大于10%的reads；3)去除低质量reads(质量值Q≤10的碱基数占整条read的50%以上)。使用bwa软件(v0.7.12)中的mem算法将过滤后的reads比对到鸭参考基因组IASCAAS_PekingDuck_PBH1.5(https://www.ncbi.nlm.nih.gov/assembly/GCA_003850225.1)上，比对参数为-k 32 -M。比对完成之后利用plink软件(v1.9)对SNPs位点质控，筛除不合格的样品和SNPs位点。质控标准如下：1)只使用常染色体上的位点；2)最小等位基因频率(MAF)>0.05；3)哈代-温伯格平衡检验的值>10；4)每个样品的SNPs检出率>90%。

1.3.2 群体遗传多样性分析 根据基因组上不同SNPs位点的分型结果，利用R软件计算群体遗传多样性不同指标，包括多态信息含量(polymorphism information content，)、固定指数(-index)、基因多样性指数()、有效等位基因数()以及香农信息指数(Shannon-Weiner index，)。4种遗传多样性指标的计算公式：

=1---2

公式1

公式2

公式3

公式4

=-(ln+ln)

公式5

其中，和分别表示某个SNP位点两个等位基因的频率；表示期望杂合度，表示观测杂合度；表示样品个体数。

对于二等位基因的单个SNP标记而言，当==05时，值最大，其值为0.375。定义单个SNP标记的相对()为以0.375为基准的值，即：

=/0.375

公式6

1.3.3 群体结构分析 利用plink软件(v1.9)和admixture软件(v1.3)分析群体遗传结构。利用筛选后得到的SNPs标记，假设样品的分群数(K值)为1～9，分别聚类分析。根据交叉验证错误率(cross-validation error，CV error)确定最优分群数。使用R包pophelper将每个样本在各个亚群的遗传组成绘制成柱形图。使用gcta软件(v1.92.2)分析群体主成分以及个体间的亲缘关系，获得PCA结果后使用R软件绘制前2个主成分的散点图。

1.3.4 ROH分布及近交系数评估 利用plink软件(v1.9)检测基因组上不同滑动窗口的ROHs，并统计检测到的ROH长度及其不同染色体上的密度。

基于ROH评估近交系数是利用全基因组分子信息评价保种效果的一种方法，即利用ROH计算基因组近交系数，计算方法：

公式 7

其中，当计算个体时，为基因组上ROHs的总长度，为基因组的总长度；当计算染色体时，为单一染色体ROHs的总长度，为该染色体的总长度。

2 结 果

2.1 基因组DNA检测

使用酚-氯仿法对娄门鸭血样提取基因组DNA后，检测所有DNA样品的浓度和OD值。结果显示，每份DNA样品浓度均在50 ng·μL，OD/OD均在1.80以上，表明DNA样品的浓度和纯度均符合要求。对基因组DNA样品进行琼脂糖凝胶电泳，结果显示DNA条带清晰，无拖尾现象，亮度均一，表明基因组DNA样品质量较好，无蛋白质污染以及降解现象(图1)，可用于后续基因组测序。

M. DNA相对分子质量标准；1～8. 基因组DNAM. DL2000 marker; 1-8. Genomic DNA图1 基因组DNA的琼脂糖凝胶电泳图Fig.1 The agarose gel electrophoresis of genomic DNA

2.2 娄门鸭群体基因组数据质量及SNPs分布

对测序下机数据进行clean reads质控后，公鸭群体平均每个个体得到3.59 M高质量reads，母鸭群体平均得到4.38 M高质量reads，与鸭基因组比对后，公、母鸭高质量reads的比对率均达到97%以上(表1)，数据质量达到后续分析的要求。

将高质量reads与鸭基因组比对后，共得到622 205个 SNPs位点。使用plink软件对获得的SNPs进行质控，最终获得374 455个高质量的SNPs位点。这些SNPs位点在不同染色体上的分布呈现明显差异，其中在NC_040046、NC_040047、NC_040048和NC_040049四个染色体上的分布比例达到50.70%(图2)。

表1 测序后的reads在鸭基因组中的比对情况

图2 SNPs在不同染色体上的密度分布图Fig.2 Density distribution of SNPs on different chromosomes

2.3 娄门鸭保种群的遗传多样性分析

基于SNPs位点的测序分型数据，计算不同SNPs位点对应等位基因频率，并统计不同染色体的遗传多样性指标分布情况。多态信息含量在整个基因组上的分布范围在0.001 5～0.375 0之间，其中34.08% SNPs位点的大于0.25，为中度多态性位点。不同染色体上的分布规律相似，均呈现“两端分布集中，中间分布稀少”的状态(图3A)。基因多样性指数(图3C)、有效等位基因数(图3D)、香农信息指数(图3E)的分布规律与相似，从这几个指标之间的相关性(表2)也可以看出，它们的皮尔逊相关系数都在0.97以上，甚至达到0.99。固定指数(-index)与等4个指标的分布规律明显不同，与其他4个 指标的皮尔逊相关系数只有0.20左右。具体到各个染色体，-index除了在NC_040075染色体上呈现两端集中分布外，其他29个染色体均偏向左侧分布(图3B)，这也导致NC_040075染色体上的-index值大于其他染色体。

、、、4个遗传多样性指标在整个基因组上的均值分别为0.154 1、0.192 9、1.336、0.296 9，且不同指标在各染色体上的均值变化规律一致(图4)。NC_040046、NC_040047、NC_040048、NC_040049、NC_040050、NC_040075等6个染色体上、、、4个指标均小于各自的基因组均值，且NC_040046染色体上的4个遗传多样性指标均最小，分别为0.124 0、0.155 1、1.271 0、0.242 0。NC_040062染色体上的4个遗传多样性指标均最大，分别为0.201 6、0.254 6、1.452 1、0.383 7。-index与其他4个指标有较大差异，其基因组上的总体均值为0.320 8，除了-index值最小(0.201 4)的NC_040062染色体以及值最大(0.492 6)的NC_040075染色体外，其他28个染色体的-index值较为接近，均在基因组总体均值上下波动。

图3 不同遗传多样性指标在染色体上的分布规律Fig.3 Distribution of different genetic diversity indexes in chromosomes

图中红色点表示该染色体的遗传多样性指标高于基因组均值，绿色表示低于基因组均值The red dots in the figure indicate that the genetic diversity indexes of the chromosome is higher than the genomic mean, while the green indicates lower图4 不同染色体上的遗传多样性指标比较Fig.4 Comparison of genetic diversity indexes on different chromosomes

表2 不同遗传多样性指标之间的相关性

根据公式(公式1)，SNP标记作为二等位基因，其最大值为0.375。以0.375为基准值，计算SNP标记的，统计娄门鸭保种群的分布特征。结果显示，娄门鸭保种群相对值在0.004 0～1.000 0之间，整个基因组上平均为0.410 9。根据指标划分低度(小于0.25)、中度(0.25～0.50)、高度多态性位点(大于0.50)的标准，娄门鸭保种群9.30%的SNPs位点属于中度多态性位点，38.80%的SNPs位点属于高度多态性位点。

2.4 娄门鸭保种群的群体结构分析

使用admixture软件分析163只鸭群体的遗传结构，根据Cross-validation error最小化原则确定最佳K值为2(图5A)。当K=2时，娄门鸭群体分为两个相对独立的亚群，分别有141个个体和22个个体；当K=3时，娄门鸭群体分为3个相对分明的亚群，分别由16个个体、125个个体和22个个体组成(图5B)。当K=4、5或之后更大数值时，娄门鸭亚群划分更为复杂，且离K=2的最佳K值较远，群体划分与实际结果有较大偏差。

图5 娄门鸭保种群体的admixture遗传结构图Fig.5 Genetic structure of Loumen duck conservation population by admixture

使用GCTA软件分析群体的主成分，结果表明，使用主成分1和主成分2可以将163只娄门鸭个体明显分为3个聚类簇(图6)。其中左上角A亚群包含22个个体(2公,20母)，右上角B亚群包含125个个体(35公,90母)，右下角C亚群包含16个个体(2公,14母)。该3个亚群的个体组成与admixture分析中K=3时的亚群分组结果完全一致。

图6 基于SNPs标记位点的娄门鸭保种群主成分分析图Fig.6 PCA graph of Loumen duck conservation population based on SNPs loci

基于过滤后的标记，使用GCTA软件分析个体间的亲缘关系，获得两两样本间的亲缘关系矩阵。基于亲缘关系矩阵分析发现，娄门鸭保种群大部分个体之间亲缘关系系数较低(图7A)。根据亲缘关系热图(图7B)，可将娄门鸭保种群分为3个亚群，分别为左上角22个个体组成的亚群、右下角16个个体组成的亚群以及剩余125个个体组成的亚群。该亚群分组结果与admixture、PCA两种方法得到的结果一致。

图7 娄门鸭保种群个体间亲缘关系分布(A)及热图(B)Fig.7 Distribution (A) and heat map (B) of genetic relationships among individuals from Loumen duck conservation population

2.5 ROHs片段长度分布及基因组近交系数FROH统计

在163只娄门鸭保种群中共发现2 966条ROHs片段，其中39只公鸭493条(平均每只12.64条)，124只母鸭2 473条(平均每只19.94条)。公鸭和母鸭ROHs片段长度主要集中在0～2 Mb区间，占比分别达到94.93%和92.00%(表3)。

表3 娄门鸭群体中不同ROHs长度的分布

基于ROH计算每个个体不同染色体的基因组近交系数，分析娄门鸭保种群不同染色体近交系数的分布规律。NC_040068、NC_040074染色体的分布较为分散，且值分别达到0.352 4和0.319 3；其他染色体的分布主要集中于左侧，其中NC_040046染色体的值最小，仅为0.006 0(图8A)。对于个体而言，娄门鸭保种群公鸭为0.027 5，母鸭为0.043 3(图8B)，差异没有统计显著性(=0.137 9)。

3 讨 论

经常评估保种群体的遗传多样性及近交系数是保种工作的一项重要内容。在物种资源保护中，尤其在家禽上一般不记录详细的系谱，因此群体内个体间的亲缘关系及近交系数等的评估结果往往与真实值之间存在较大的差距。基因组中SNPs标记不受环境影响，且在基因组中分布广泛和易于检测，是评估种群遗传多样性的理想载体。随着高通量技术的发展，SNP芯片和基因组二代测序方法成为基因组范围内大量SNPs标记批量检测的重要技术。目前，鸭上尚没有一款能够实用的SNP芯片，而简化基因组技术在鸭上的成功开发，促使其成为鸭遗传多样性评估中的常用方法。因此，本课题组利用简化基因组测序技术系统研究了娄门鸭保种群体的遗传多样性，并分析了娄门鸭的群体结构。

A图中每个黑点表示相应染色体的FROH均值，B图中绿点和红点分别表示公鸭、母鸭群体FROH的均值Each black dot in figure A indicates the mean FROH value of the corresponding chromosome, while green and red dots in figure B indicate the mean FROH values of the male and female duck populations, respectively图8 娄门鸭保种群不同个体及染色体基因组近交系数FROH分布Fig.8 The distribution of genomic inbreeding coefficient FROH of different individuals and chromosomes

群体遗传多样性有多种度量方法，各有优缺点和适用范围。本课题组比较了多态信息含量、固定指数-index、有效等位基因数、基因多样性指数和香农信息指数等5个代表性的遗传多样性指标，发现除了-index之外的其他4个 指标之间具有较高的相关性，其在染色体上的分布规律也较为一致，说明这4个遗传多样性指标所蕴含的信息是冗余的，因此选择其中一种即可说明保种群体的遗传多样性情况。

是衡量多样性高低的较好指标，尤其适用于分子标记，因此文献报道中对的研究较为系统，当≥0.5时，该座位为高度多态性位点；当0.25≤<0.5时为中度多态性位点；当<0.25时为低度多态性位点。该判别标准适合于微卫星等复等位基因的分子标记，这些分子标记的特点是单个位点的等位基因数量不固定，当无限多个等位基因数时，最大值为1。SNP标记属于二等位基因，其最大为0.375，因此直接以为0.25和0.50的划分标准来判定娄门鸭保种群结果有失偏颇。为了继续使用该判定标准，课题组以SNP位点最大值0.375为基准，创新性地使用相对()这一指标来判定SNP位点的多态性高低。本研究结果显示，娄门鸭保种群相对值在0.004 0～1.000 0之间，整个基因组上平均相对为0.410 9，中度多态性位点占比9.30%，高度多态性位点达到38.80%，结果表明该群体遗传多样性较为丰富，目前娄门鸭保种方案比较合理。

对娄门鸭保种群每个单SNP位点遗传多样性指标分别分析，并考察这些指标在不同染色体上的分布情况。总体而言，、、、四个遗传多样性指标之间的分布规律相似，而且每个指标在不同染色体上的分布规律也相似；而固定指数-index与其他指标在染色体上的分布规律有较大差异，且-index在不同染色体上的分布规律也不尽相同。该结果进一步表明，、、和四个遗传多样性指标中，可以选择其中一种，比如来说明群体的遗传多样性情况，而-index可作为其他多样性指标信息的有效补充。

固定指数-index是群体中杂合子频率差异的估计值，能够通过杂合子频率偏差反映群体在某个座位的遗传分化程度。当一个群体分化成多个亚群时，整个群体中的杂合子必然会减少，大量位点的-index高于0，表明出现了近亲杂交，群体有分化倾向。娄门鸭保种群整个基因组的-index均值为0.320 8，说明娄门鸭群体可能有多个亚群出现。这从群体结构分析结果中也可以看出。群体结构指遗传变异在物种或群体中的一种非随机分布。处于同一亚群内的不同个体亲缘关系较高或SNPs分布状态接近，而亚群与亚群之间则亲缘关系稍远或SNPs分布状态差异较大。本研究从遗传结构、PCA主成分以及分子亲缘关系三个角度分析娄门鸭保种群的群体结构，均将娄门鸭群体分为3个亚群，其中一个较大的亚群(35公、90母)和两个较小的亚群(分别为2公、20母和2公、14母)。-index和群体结构分析结果都说明娄门鸭保种群出现了一定的分化现象。群体的分化首先从基因组局部区域的遗传多样性变化开始，逐步扩展到基因组大部分区域，分化程度更高时甚至出现外观、性状表型的变化。本研究中,娄门鸭分化成3个亚群，根据-index结果可知其主要表现在染色体局部区域，分化群体间尚没有表现出外观、性状表型的明显差异。群体分化的出现可能由保种群留种时的抽样误差导致的遗传漂变造成，后续保种过程中可以将群体结构分析得到的3个分化的亚群之间的个体进行适当的非随机交配，消除目前的群体分化现象。

利用SNPs位点信息估计分子近交系数的方法主要有基于单个位点纯合性以及长纯合片段(runs of homozygosity，ROH)两种。与传统的基于系谱的近交系数以及SNPs位点杂合度的近交系数相比，基于ROH的基因组近交系数更加接近真实的近交系数。刘家鑫等利用绵羊全基因组ROH比较了10个绵羊群体的近交系数，认为基于ROH评估的近交系数可为绵羊种群的保种提供参考依据。史良玉等利用ROH对不同来源的大白猪产仔数近交衰退进行评估，发现法系群体中只有检测到近交衰退现象，说明基于的方法敏感性更高。曹愉夏等利用简化基因组测序数据计算基于ROH的基因组近交系数，有效校正了狼山鸡由于部分系谱缺失导致的近交系数失真。本研究利用简化基因组测序得到了全基因组水平的SNPs标记信息，分析ROH片段在整个基因组以及不同染色体上的分布情况，并评估基因组近交系数。对于个体而言，娄门鸭公鸭近交系数均值为0.027 5，母鸭为0.043 3，说明该种群近交系数较低，这也与该群体遗传多样性较为丰富有关。对于染色体而言，不同染色体的近交系数差异较大，其中NC_040068、NC_040074的达到0.30以上，远高于其他染色体均值，这可能与娄门鸭在保种过程中这两个染色体上部分区域受到纯化选择有关。后续保种过程中应根据ROH片段鉴定并重点监测这两个染色体上近交系数较高的基因组区域，防止过度近交，避免个别染色体近交系数上升过快。

4 结 论

本研究利用不同遗传多样性指标分析娄门鸭保种群的遗传多样性，显示娄门鸭保种群的遗传多样性较丰富，群体基因组近交系数较低，但个别染色体的近交系数较高，且群体出现了部分分化。后续保种中可以将遗传结构分析得到的3个分化亚群个体之间进行适当的非随机交配，消除目前的群体分化现象，并重点监测染色体近交系数较高的基因组区域，避免个别染色体近交系数上升过快。