基于RAD-seq静原鸡保种群体的遗传变异分析

马丽霞，曹国伟，朱红芳，邓占钊，蔡正云，周成浩，韩 威，顾亚玲，张 娟*

(1. 宁夏大学农学院，银川 750021； 2. 彭阳县畜牧技术推广服务中心，固原 756599；3. 彭阳县动物卫生监督所，固原 756599； 4. 江苏省家禽科学研究所 国家级地方鸡种基因库，扬州 225125)

静原鸡原名固原鸡，也称静宁鸡，原产地和中心产区位于宁夏回族自治区固原市和甘肃省静宁县。2006年被列入《国家级畜禽遗传资源保护名录》，是宁夏回族自治区5个区级畜禽遗传资源保护品种之一。静原鸡根据羽色可以分为白羽、麻羽和黑羽，具有典型的地方禽类品种属性，如：耐粗饲、抗逆性强、环境适应性好以及良好的山区放养性等。同时，静原鸡兼具屠宰率高、肉质细嫩、味道鲜美、脂肪酸含量高、胆固醇含量低等特点，已逐渐成为当地绿色食品品牌之一。但由于引入品种鸡的冲击，使得这一优良地方品种杂化现象严重，并出现生活力降低、生长势减弱、生产性能下降、抗病能力差等一系列不良现象。因此，加快静原鸡的保种对保护国家及地方畜禽种质资源具有重要意义。现阶段，静原鸡的保护主要以保种场为主体，以家系等量留种的原则开展活体保护，虽然在极大程度上保护了静原鸡群体的遗传多样性，但也存在因初始群体系谱不完善造成的祖源不清问题。

在物种资源保护中，系谱记录依赖于基础群，当系谱数据记录存在问题时，会导致个体间的亲缘关系及近交系数出现错误，因此会出现近交衰退和适应能力下降等现象，而种群内的遗传多样性对于适应能力和避免近交衰退是必要的。目前，单核苷酸多态性(single nucleotide polymorphism, SNP)已成为最理想的群体遗传多样性和基因组学的分子标记，而限制性酶切位点相关的DNA(re-striction-site associated DNA，RAD)测序技术能鉴定出广泛覆盖基因组范围的SNP位点，被广泛应用于资源群体的遗传多样性研究。在家禽中，韩威等基于RAD-seq技术在19个地方鸡种中鉴定出400 562个SNPs标记，揭示了19个品种的遗传多样性和群体结构并发掘了与免疫系统调节、生殖机能调控、应激响应等生物学过程相关的种质特性基因。在草食动物中，Liu等利用RAD-seq技术研究了6个中国地方兔品种和2个引进兔品种的遗传多样性和群体结构，鉴定出了与黑素生成相关的基因。Liu等采用RAD-seq技术在内蒙古、甘肃、青海和新疆4个地区的47只家养双峰驼中检测到1 568 087个SNPs，且新疆骆驼的遗传多样性最高，连锁不平衡(LD)系数的衰减率最快。

鉴于此，本研究利用RAD-seq技术从分子水平研究静原鸡群体的遗传多样性和分子遗传进化，比较白羽、麻羽和黑羽静原鸡的遗传多样性差异以及群体结构间的关系，筛选与静原鸡羽色性状相关的候选基因，以期为静原鸡不同羽色的品系化培育提供理论依据。

1 材料与方法

1.1 样品收集和DNA的提取

本试验所用静原鸡样品来自宁夏回族自治区固原市彭阳县静原鸡繁育中心。在相同饲养条件下选取180日龄品种特征明显、发育正常且健康的白羽(white feather, WF)、麻羽(hemp feather, HF)、黑羽(black feather, BF)静原鸡各60只(40只母鸡，20只公鸡)，翅下静脉采血1.5 mL，EDTA(乙二胺四乙酸)抗凝，通过酚-氯仿法抽提血液基因组DNA。

1.2 RAD文库的构建和测序

利用ddRAD建库方式构建长度范围在300～500 bp的pair-end文库，进行R I(G^AATTC)和III(Hin1IICATG^)双酶切简化基因组测序。

1.3 原始数据的质控和过滤

利用GATK和samtools程序对原始测序数据进行质量控制，截除Reads中的测序接头以及引物序列，过滤平均质量值小于Q5的Reads、过滤N个数 > 5的Reads和长度过短序列，获取SNP质量值Q≥30，最小等位基因频率(minor allele frequency, MAF)≥ 0.05，SNP信息完整度≥0.70的高质量数据。其它过滤参数：QD < 2.0，MQ < 40.0，FS > 60.0，SOR > 6.0，MQRankSum < -12.5，ReadPosRankSum < -8.0，符合参数即过滤。对SNP进行过滤后，使用得到的SNP标记进行PCA、系统进化树、群体结构、选择性清除等分析。

1.4 SNP的检测

利用BWA软件将获得的Clean Reads比对到参考基因组上(http://ftp.ensembl.org/pub/release-104/fasta/gallus_gallus)，根据Clean Reads在参考基因组的定位结果，使用GATK软件进行SNP的检测。

1.5 群体遗传多样性分析

通过PopGen软件计算观察杂合度()、群体内核酸多态性()和群体的平均近交系数()。利用软件PopLDdecay进行连锁不平衡(linkage disequilibrium，LD)分析，参数设置为：-OutPairLD 5。

1.6 群体结构分析

通过软件GCTA进行主成分分析(principal components analysis，PCA)，然后利用R包绘制PCA散点图。对过滤后筛选得到的SNPs标记，采用软件FastTree中的极大似然法(maximum likelihood，ML)构建系统进化树(phylogenetic tree)，并使用Shimodaira-Hasegawa test方法判断每个节点的可信度。使用admixture软件进行群体遗传结构分析，预先设定K=2～15的遗传簇数。

1.7 选择性清除分析

对经过质量控制后的SNP进行选择性清除分析。采用遗传分化系数()和核苷酸多样度(π)相结合的方法，以100K为一个window，10K为步长取一个区域(群体中任意两条不同序列的碱基差异数取平均值)，值按照大小排列，取大小排在前5%的区域信息，π的比值取位于前5%与后5%的区域，根据与π比值筛选出的区域取交集作为关联结果，该过程使用软件vcftools进行分析。

1.8 全基因组关联分析

使用TASSEL(https://tassel.bitbucket.io)软件的MLM模型对不同群体的羽色性状进行关联分析，计算模型为：=α+β+μ+，其中，为表型性状，为基因型(固定效应)的指示矩阵，为固定效应的估计参数；为群体遗传结构的指示矩阵，为SNP的效应；为个体亲缘关系指示矩阵，为预测的随机个体；是随机残差，服从e～(0，δe2)。通过关联的显着度(-value)筛选出潜在的候选SNPs，利用ANNOVAR对位点进行注释。

1.9 基因注释和富集分析

将选择性清除分析和全基因组关联分析得到的显着SNPs位点映射至ENSEMBL数据库(http://ftp.ensembl.org/pub/release-104/fasta/gallus_gallus)公布的鸡参考基因组上，进行显着位点的注释。利用KOBAS(http://bioinfo.org/kobas/annotate/)对受选择的基因进行KEGG分析。

2 结 果

2.1 RAD测序和数据质量

利用RAD-seq技术对静原鸡3个类群180只鸡进行序列测定。通过Illumina平台简化基因组测序获得198.83 Gb 的Clean Data，白羽、麻羽、黑羽静原鸡的Reads数分别为3 731 394、4 100 199和3 674 557。 Q30达到93%以上，Q20达到97%以上，GC含量达到39%以上(表1)。总体而言，测序数据显示出很高的PHRED质量。

表1 样品测序数据评估统计表Table 1 Statistical table for evaluation of sample sequencing data

2.2 群体遗传多样性分析

静原鸡白羽、麻羽群体内鉴定到的SNPs依次为238 533和233 562个，黑羽的SNPs数最多为240 820个。 3种羽色的静原鸡、和分别在0.273 2～0.278 2、0.304 9～0.309 6和0.096 1～0.109 8之间(表2)。LD衰减分析表明，不同类群LD系数的衰减率差异较小。3个类群SNPs距离在0～50 kb范围内衰减较快，衰减系数从0.5降至0.1以下。SNP距离在50～300 kb范围LD水平较低，当LD系数r=0.1时，各类群对应的LD衰减距离基本相等(图1)。

图中横坐标表示连锁不平衡发生的距离，纵坐标是连锁不平衡相关系数The abscissa in the figure represents the distance of linkage disequilibrium, and the ordinate is the correlation coefficient of linkage disequilibrium图1 LD随距离增长的衰减图Fig.1 LD decay with distance

表2 静原鸡3种羽色遗传参数的比较Table 2 Comparison of genetic parameters of 3 feather colors of Jingyuan chicken

2.3 群体结构分析

主成分分析的结果中PC1和PC2这两个特征向量分别可以解释8.98%和4.41%的遗传变异，结果表明除了麻羽静原鸡144号和黑羽静原鸡180号两个个体以外，静原鸡其他个体形成了白羽、麻羽和黑羽3个互不重叠的类群，且类群之间分型明显(图2)。本研究利用极大似然法构建系统发育树，结果显示相同羽色的静原鸡亲缘关系较近，不同羽色的静原鸡之间亲缘关系相对较远，该结果与PCA结果一致(图3)。本研究利用admixture软件分析不同遗传物质在每个个体中的比例，在对每个K值进行群体结构分析时，进行了交叉验证，K=3的交叉验证最适合于静原鸡的真实分化史(图4A)。结果显示，当K=2时，白羽与其他羽色分离，表明该类群与麻羽和黑羽系统发育上存在相对较远的距离。当K=3时，180个个体根据羽色分成了3个类群，与系统进化树分析结果一致。当K=4时，白羽的祖先背景较为纯正，具有主要的遗传祖先。虽然黑羽和麻羽有多个遗传祖先，但一个主要的遗传祖先是显而易见的(图4B)。

图中PC1为主成分1，PC2为主成分2。图中每个点代表一个样品，每组样品用不同颜色表示PC1 is the principal component 1, PC2 is the principal component 2. Each point in the figure represents a sample, and each group of samples is represented by each color图2 PCA分析Fig.2 PCA analysis

树形拓扑结构直观展示了不同种类之间的进化关系，亲缘关系较近物种的进化分枝往往聚成一簇The tree-shaped topological structure intuitively shows the evolutionary relationship between different species, and the evolutionary branches of species with close kinship are often clustered together图3 系统进化树Fig.3 Phylogenetic tree

A.CV-error图: 横坐标表示不同的K值，纵坐标表示不同的K值所对应的CV-error值。B.群体遗传结构图：横坐标是每个样品的结构，其中每一列表示一个个体，其中不同颜色片段的长度表示该个体基因组中某个祖先所占的比例。图片上侧K=2～4表示假定的祖先群体个数从2一直到4，K值是样本所包含的亚群或者祖先数A. CV-error plot: The horizontal coordinates indicate different K values, and the vertical coordinates indicate the CV-error values corresponding to different K values. B. Genetic structure of the population: The abscissa is the structure of each sample, where each column represents an individual, and the length of the different color fragments represents the proportion of an ancestor in the individual’s genome. On the upper side of the picture, K=2 to 4 means that the number of assumed ancestor groups ranges from 2 to 4, and the K value is the number of subgroups or ancestors contained in the sample图4 群体结构分析图Fig.4 population structure analysis diagram

2.4 选择性清除分析

通过结合π和来选择前5%的区域。以静原鸡白羽为试验组，黑羽为正向选择组，共获得819个正向选择区域(图5A)。以白羽为试验组，麻羽为正向选择组，共获得正向选择区域730个(图5B)。以黑羽为对照组，麻羽为正向选择组，检测到719个正向选择区域(图5C)。对筛选区域进行基因注释之后利用KOBAS对检测的候选基因进行KEGG富集分析，KEGG富集结果表明候选基因显着富集到37条通路(<0.05)，其中与羽色相关的通路中富集最显着的是黑色素生成途径，此外还有酪氨酸代谢和Wnt信号传导等(图6)。在这些通路上最终筛选出11个与静原鸡羽色相关的候选基因(表3)。

A.白羽vs.黑羽的选择性清除分析；B.白羽vs.麻羽的选择性清除分析；C.黑羽vs.麻羽的选择性清除分析。横坐标为π的比值和纵坐标为Fst值分别对应上面的频率分布图和右侧的频率分布图，中部的点图则代表不同窗口内的相应的Fst和π比值。其中最上方蓝色和红色区域为π选择出来的top 5%区域，绿色区域为Fst所选择top 5%区域，中间蓝色和红色区域为Fst和π的交集，即为候选的位点A. Selective sweep analysis of white feathers vs. black feathers; B. Selective sweep analysis of white feathers vs. hemp feathers; C. Selective sweep analysis of black feathers vs. hemp feathers. The frequency distribution diagram above and the frequency distribution diagram on the right correspond to the π ratio on the abscissa and Fst value on the ordinate, respectively, the dot diagram in the middle represent the corresponding Fst and π ratios in different windows. The top blue and red areas are the top 5% area selected by π, the green areas are the top 5% area selected by Fst, and the middle blue and red areas are the intersection of Fst and π, which is the candidate sites图5 选择性清除分析Fig.5 Selective sweep analysis

表3 选择性清除分析的选择区域及候选基因Table 3 Selection regions and candidate genes for selective sweep analysis

图6 KEGG通路富集分析Fig.6 KEGG pathway enrichment analysis

2.5 全基因组关联分析

通过对静原鸡羽色进行全基因组关联分析后发现，白羽性状关联到66个显着相关的SNPs位点(<1.48×10)，分布在1、4、8、15号4条染色体上(图7A)，临近或坐落于、、2等30个基因，单个关联标记位点解释的表型变异(R)范围为3.16%～9.18%，前10个基因见表4；麻羽性状鉴定到259个显着相关SNPs位点(<1.48×10)，这些位点主要分布在1～9号染色体上(图7C)，临近或坐落于2、135、3等116个基因，单个关联标记位点解释的表型变异(R)为6.73%～19.53%，前10个基因见表5，黑羽性状关联到200个显着相关的SNPs位点(<1.48×10)，与黑羽性状相关的SNPs主要位于1～9号染色体上(图7E)，临近或坐落于4、4、135等114个基因，单个关联标记位点解释的表型变异(R)为5.85%～14.85%，前10个基因见表6。QQ图是关联分析结果重要的质控图，λ能够判断群体是否分层，λ在1.05以上说明群体存在分层，本研究的λ值均小于1.05说明个体均匀分布，不存在群体分层现象(图7B、7D、7F)。通过文献查阅和KOBAS对关联的基因进行KEGG富集分析(图8)，在白羽和麻羽静原鸡中分别筛选出2个(、)与羽色相关的基因，在黑羽静原鸡中筛选出两个基因(6、2)与羽色相关，这些基因主要在黑色素合成和酪氨酸代谢中发挥重要的调控作用(表7)。

A.白羽的曼哈顿图； C.麻羽的曼哈顿图； E.黑羽的曼哈顿图，横坐标为染色体。纵坐标-lgP，P值越小关联性越强，表现为纵坐标越大，红线表示0.05/全部SNPs的阈值，蓝线表示1/全部SNPs的阈值。B.白羽的QQ图；D.麻羽的QQ图；F.黑羽的QQ图。QQ图的纵坐标是SNP位点的P-value值(这是实际得到的结果)，与曼哈顿图一样也是表示为-lg(P-value)；横轴是则是均匀分布的概率值(这是期望的结果)，同样也是换算为-lg(P-value)A. Manhattan figure of white feathers; C. Manhattan figure of hemp feathers; E. Manhattan figure of black feathers. The horizontal coordinate is chromosomes, the vertical coordinate is -lgP, the smaller the P value, the stronger the association, which is expressed as a larger vertical coordinate, the red line indicates the threshold of 0.05/all SNPs and the blue line indicates the threshold of 1/all SNPs. B. White feather’s QQ picture; D. Hemp feather’s QQ picture; F. Black feather’s QQ picture. The vertical coordinate of the QQ plot is the P-value of the SNP loci (this is the actual result obtained), which is also expressed as -lg(P-value) as in the Manhattan plot; the horizontal axis is the probability value of the uniform distribution (this is the expected result), which is also converted to -lg (P-value)图7 静原鸡羽色性状的曼哈顿图和QQ图Fig.7 Manhattan and QQ maps of feather color traits of Jingyuan chicken

A. 白羽相关基因的KEGG通路富集分析；B. 麻羽相关基因的KEGG通路富集分析；C.黑羽相关基因的KEGG通路富集分析。x轴表示富集度，y轴表示通路名称，气泡大小表示富集到通路中的基因数量，颜色表示富集到该通路富集显着性A. KEGG pathway enrichment analysis of genes related to white feathers; B. KEGG pathway enrichment analysis of genes related to hemp feathers; C. KEGG pathway enrichment analysis of genes related to black feathers. x-axis indicates the degree of enrichment, the y-axis indicates the pathway names, bubble size indicates the number of genes enriched to the pathway, and color indicates the enrichment significance enriched to pathway图8 KEGG富集分析Fig.8 KEGG enrichment analysis

表4 白羽关联显着SNPs统计信息表Table 4 Statistical information table of significant SNPs associated with white feathers

表5 麻羽关联显着SNPs统计信息表Table 5 Statistical information table of significant SNPs associated with hemp feathers

表6 黑羽关联显着SNPs统计信息表Table 6 Statistical information table of significant SNPs associated with black feathers

表7 GWAS筛选的候选基因Table 7 Candidate genes screened by GWAS

3 讨 论

3.1 静原鸡遗传多样性和群体结构多样性

随着生产技术及生产力的发展，畜禽经过长期的自然选择、遗传漂变和人工选择，其遗传多样性随之发生改变。静原鸡作为国家级地方保护资源，其遗传多样性同样受到了严重威胁。2013年，固原市彭阳县建立了静原鸡原种保种场，对静原鸡进行系统保护和开发利用，但是在这期间由于保种知识的欠缺，导致系谱数据的大量缺失和错误，因此本研究通过分子保种技术保护静原鸡群体的遗传多样性并持续监测保种群的状态，评估和制定有效的保种技术手段，对静原鸡以及我国地方鸡种的保护与开发利用具有重要意义。是反映等位基因频率的一个指标，其数值越接近，表明群体的分布越均匀。可以明确近交程度，对于选种选配防止近交衰退和品系繁育具有重要的指导意义。越高，说明物种的多样性越高。本研究利用RAD-seq在静原鸡共检测出712 915个SNPs，3个类群的、和分别在0.273 2～0.278 2、0.304 9～0.309 6和0.096 1～0.109 8之间，与韩威等所研究的19个地方鸡种相比，静原鸡的和均高于19个地方鸡种，低于瓢鸡、河南斗鸡、金湖乌凤鸡等8个地方鸡种，由此表明静原鸡遗传多样性丰富，品系选育良好，资源活力好，开发利用潜力较大。LD是评价遗传多样性的一个重要指标。本研究中，3个类群LD系数的衰减率差异较小，说明该群体没有经过较强的选择，受选择的强度低。主成分分析中除麻羽静原鸡144号和黑羽静原鸡180号两个个体之外，其他个体根据羽色聚在一起，这与进化树和群体结构分析的结果一致，说明静原鸡根据羽色分成了不同的类群，这为静原鸡的分群保种提供了理论依据。

3.2 选择性清除分析

物种的进化离不开选择，在选择的过程中，群体内的有利突变发生后，这个突变基因的适合度越高，越容易被固定，与此基因座连锁的染色体区域，由于搭车效应也被固定下来，大片紧密连锁的染色体区域因此失去多态性，这种由于搭车效应引起多态性下降的现象，遗传上称为选择清除。通过选择清除分析，可以鉴定出差异基因组区域，结合功能注释进一步预测表型相关的候选基因。本研究通过选择性清除分别获得819、730、719个选择区域，这些区域内筛选到与静原鸡羽色相关的候选基因共11个(4、16、、、、1、、1R、2A、、)。TYR是一种酪氨酸酶，它的表达缺乏将导致白羽球茎黑色素生物合成不足，是白羽形成的直接原因。李洪林鉴定出控制兴义矮脚鸡隐性白羽的主效基因为。在色素细胞的增殖和分化中起重要作用，包括黑素母细胞和黑素细胞，研究表明2在番鸭黑羽中的表达量高于白羽。FZD4是一种Wnt蛋白的受体，通过Wnt信号途径调控黑色素生成。16是WNT信号通路中的重要基因，在神经嵴来源的黑素细胞发育和迁移等过程中发挥重要作用。1或3表达缺陷的小鼠表现出神经嵴黑色素细胞形成的丧失。具有内在的GTPase活性，在调节细胞增殖中起重要作用，研究表明基因突变会造成黑色素细胞增生，从而导致黑色素瘤。是色素沉着的主要调节剂，是Wnt信号通路的靶标，降低1的表达可能会下调的表达，最终减少黑色素的合成。在黑色素合成途径中，前体多巴(3, 4-二羟基苯丙氨酸)和多巴胺由和合成。1R编码一个七跨膜结构域G蛋白偶联受体，主要在发育中的羽毛和头发的黑色素细胞中表达。Nam等研究发现，位置和表达水平的差异对韩国本土鸡品种的羽毛色素沉着的影响比1R表达水平的变化更大。CAMK2A是钙调蛋白依赖性蛋白激酶，研究表明，2是调控中国五指山黑猪毛色的候选基因。PRKCA和PRKCB是两种蛋白激酶，PRKCA和PRKCB的激活是抑制黑色素合成的原因，研究已证实了的表达与小鼠皮肤色素沉着大致相关。综上所述，调控静原鸡羽色的基因主要与黑色素合成相关。

3.3 全基因组关联分析

全基因组关联分析是利用自然群体探测物种的遗传变异，进而分析复杂性状遗传结构的有力工具，近年来已在鸡、牛、猪、羊等物种中均有广泛研究。本研究针对静原鸡的3种羽色性状开展全基因组关联分析，初步鉴定到影响静原鸡群体白、麻、黑3种羽色性状相关的致因突变位点和候选基因。其中鉴定到与白羽性状显着关联的位点7个，对这些基因进行功能注释，发现富集在酪氨酸酶通路，该通路在合成黑色素和其他多酚化合物等色素时起重要作用，关联分析结果表明，有2个 显着关联位点189168614、189176762位于基因的第3内含子上。该基因被选择性清除和全基因组关联分析同时筛选到，由此推测是控制静原鸡白羽的关键候选基因，具体调控过程还需进一步验证。另外，对黑羽显着关联水平SNPs分析发现919883和129265651位点与黑色素形成有关，这两个位点分别位于2基因的第2内含子和6的第2内含子上，2是一种磷酸肌醇磷脂酶，6是一种Wnt蛋白的受体。Zhang等通过逆转录聚合酶链式反应(RT-PCR)分析了2基因在各种癌症系中的mRNA表达情况，同时结合其他生物学手段发现2是影响黑色素瘤增殖和凋亡的关键因素。Korstanje等基于生物信息学分析确定了两个导致色素沉着缺陷的寡聚体复合物6和5R，研究发现6敲除将导致色素沉着障碍，由此得出6能够间接影响色素沉着。但有关2和6对羽色的调控研究尚未见报道。除此之外与麻羽性状显着关联的位点有178个，其中与羽色相关的只有13923834位点，该位点位于的9号内含子上，TH是酪氨酸3-羟化酶，Fernstrom等研究表明，在TH酶的催化下，能生成DOPA，而DOPA 是黑色素形成的关键物质。由此说明可能是调控羽色的基因。综上所述，静原鸡3种羽色形成的分子机制不同，遗传基础复杂，以现有的资源群体并不能准确关联到所有致因突变和潜在候选基因，在后续研究中将扩大资源群体规模、优化模型和方法，进一步对相关突变位点进行挖掘和验证。

4 结 论

本研究利用RAD-seq技术对白羽、麻羽和黑羽静原鸡群体进行了全基因组范围内的遗传变异检测，基于鉴定到的SNPs计算群体遗传学统计指标，系统评价了静原鸡群体的遗传多样性和遗传结构，通过选择性清除和全基因组关联分析鉴定了与羽色相关的4、16、、、、1、、1R、2A、、基因。综上所述，本研究结果将对静原鸡这一地方种质资源群体的保护奠定理论基础，同时为静原鸡不同羽色的品系化培育提供新的遗传标记和基因靶点。