表 1(Table 1)
图 1(Fig. 1)
基于RAD-seq静原鸡保种群体的遗传变异分析
引用本文
马丽霞, 曹国伟, 朱红芳, 邓占钊, 蔡正云, 周成浩, 韩威, 顾亚玲, 张娟. 基于RAD-seq静原鸡保种群体的遗传变异分析[J]. 畜牧兽医学报, 2022, 53(7): 2104-2117.
MA Lixia, CAO Guowei, ZHU Hongfang, DENG Zhanzhao, CAI Zhengyun, ZHOU Chenghao, HAN Wei, GU Yaling, ZHANG Juan. Analysis of Genetic Variation in a Conserved Population of Jingyuan Chickens Based on RAD-seq[J]. Acta Veterinaria et Zootechnica Sinica, 2022, 53(7): 2104-2117.
基于RAD-seq静原鸡保种群体的遗传变异分析
1. 宁夏大学农学院, 银川 750021;
2. 彭阳县畜牧技术推广服务中心, 固原 756599;
3. 彭阳县动物卫生监督所, 固原 756599;
4. 江苏省家禽科学研究所 国家级地方鸡种基因库, 扬州 225125
2. 彭阳县畜牧技术推广服务中心, 固原 756599;
3. 彭阳县动物卫生监督所, 固原 756599;
4. 江苏省家禽科学研究所 国家级地方鸡种基因库, 扬州 225125
摘要：旨在分析静原鸡白羽、麻羽、黑羽保种群之间的遗传变异，挖掘调控特征性状的关键候选基因。本研究利用RAD-seq技术对180日龄特征明显、发育正常且健康的白羽、麻羽、黑羽静原鸡(每个羽色选取60个个体，40只母鸡，20只公鸡)进行测序，基于SNP标记计算观察杂合度(Ho)、群体内核酸多态性(Pi)、群体的平均近交系数(Fis)等指标，分析静原鸡3个类群的遗传多样性和群体结构，并通过选择性清除分析和全基因组关联分析(GWAS)筛选出候选基因，利用KOBAS对候选基因进行KEGG (Kyoto Encyclopedia of Genes and Genomes)通路富集, 最终筛选出调控静原鸡羽色的候选基因。测序结果表明，静原鸡180个样品共产生198.83 Gb Clean Data，Q30达到93%以上。遗传多样性分析表明，白羽、麻羽、黑羽静原鸡分别鉴定出238 533、233 562和240 820个SNPs标记，Ho、Pi和Fis分别在0.273 2~0.278 2、0.304 9~0.309 6和0.096 1~0.109 8之间。群体结构分析表明，静原鸡根据不同羽色分为不同类群。通过选择性清除分析和全基因组关联分析共筛选出11个(FZD4、WNT16、EDNRB、TYR、KRAS、CTNNB1、DDC、MC1R、CAMK2A、PRKCB、PRKCA)与静原鸡羽色相关的候选基因，富集结果显示这些基因主要与黑色素生成、酪氨酸代谢和Wnt信号传导等通路相关。综上所述，本研究利用SNPs标记信息可以全面地评价静原鸡的保种现状，为静原鸡的遗传资源保护奠定理论基础。同时，筛选出了与静原鸡羽色性状相关的候选基因，为静原鸡不同羽色的品系化培育提供新的遗传标记和基因靶点。
关键词：静原鸡 RAD-seq 遗传多样性分析 选择性清除分析 全基因组关联分析
Analysis of Genetic Variation in a Conserved Population of Jingyuan Chickens Based on RAD-seq
MA Lixia1
, CAO Guowei2, ZHU Hongfang3, DENG Zhanzhao2, CAI Zhengyun1, ZHOU Chenghao4, HAN Wei4, GU Yaling1, ZHANG Juan1
, CAO Guowei2, ZHU Hongfang3, DENG Zhanzhao2, CAI Zhengyun1, ZHOU Chenghao4, HAN Wei4, GU Yaling1, ZHANG Juan1
1. College of Agriculture, Ningxia University, Yinchuan 750021, China;
2. Pengyang County Animal Husbandry Technology Popularization Service Center, Guyuan 756599, China;
3. Pengyang County Animal Health Supervision Institute, Guyuan 756599, China;
4. National Chicken Genetic Resources, Jiangsu Institute of Poultry Science, Yangzhou 225125, China
2. Pengyang County Animal Husbandry Technology Popularization Service Center, Guyuan 756599, China;
3. Pengyang County Animal Health Supervision Institute, Guyuan 756599, China;
4. National Chicken Genetic Resources, Jiangsu Institute of Poultry Science, Yangzhou 225125, China
Abstract: The aim of this study was to analyze the genetic variation among the white-, hemp- and black-feathered breeding populations of Jingyuan chickens, and to explore the key candidate genes regulating the characteristic traits. In this study, the RAD-seq technology was used to sequence 180-day-old white-, hemp- and black-feathered Jingyuan chickens (60 individuals of each plumage color, 40 hens and 20 roosters) with distinctive features, normal development and health, and the observed heterozygosity (Ho), nucleic acid polymorphism (Pi) within the population and the average inbreeding coefficient (Fis) of the population based on SNPs markers were calculated to analyze the genetic diversity and population structure of Jingyuan chickens. The candidate genes were screened by selective sweep analysis and genome-wide association analysis (GWAS), and enriched by KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways using KOBAS. The candidate genes regulating plumage color of Jingyuan chickens were finally screened. The sequencing results showed that a total of 198.83 Gb of clean data was obtained from 180 samples of Jingyuan chickens, with Q30 reaching more than 93%. Genetic diversity analysis showed that 238 533, 233 562 and 240 820 SNPs markers were identified for the white-, hemp-, and black- feathers of Jingyuan chickens, with Ho, Pi and Fis range of 0.273 2-0.278 2, 0.304 9-0.309 6 and 0.096 1-0.109 8, respectively. The analysis of population structure showed that Jingyuan chickens were divided into different groups according to different feather colors. A total of 11 candidate genes (FZD4, WNT16, EDNRB, TYR, KRAS, CTNNB1, DDC, MC1R, CAMK2A, PRKCB, PRKCA) associated with plumage color in Jingyuan chickens were screened by selective sweep analysis and genome-wide association analysis, and the enrichment results showed that these genes were related to the pathways of melanogenesis, tyrosine metabolism, and Wnt signaling. In summary, the breeding status of the Jingyuan chicken can be comprehensively evaluated using SNPs marker information and lay a theoretical foundation for the conservation of the genetic resources of the Jingyuan chicken in this study. At the same time, the candidate genes related to the plumage color trait of Jingyuan chicken were screened to provide new genetic markers and gene targets for the breeding of different plumage colors of Jingyuan chicken.
Key words:
Jingyuan chicken RAD-seq genetic diversity analysis selective sweep analysis genome-wide association analysis
静原鸡原名固原鸡，也称静宁鸡，原产地和中心产区位于宁夏回族自治区固原市和甘肃省静宁县。2006年被列入《国家级畜禽遗传资源保护名录》，是宁夏回族自治区5个区级畜禽遗传资源保护品种之一[1]。静原鸡根据羽色可以分为白羽、麻羽和黑羽，具有典型的地方禽类品种属性，如：耐粗饲、抗逆性强、环境适应性好以及良好的山区放养性等。同时，静原鸡兼具屠宰率高、肉质细嫩、味道鲜美、脂肪酸含量高、胆固醇含量低等特点，已逐渐成为当地绿色食品品牌之一[2]。但由于引入品种鸡的冲击，使得这一优良地方品种杂化现象严重，并出现生活力降低、生长势减弱、生产性能下降、抗病能力差等一系列不良现象。因此，加快静原鸡的保种对保护国家及地方畜禽种质资源具有重要意义。现阶段，静原鸡的保护主要以保种场为主体，以家系等量留种的原则开展活体保护，虽然在极大程度上保护了静原鸡群体的遗传多样性，但也存在因初始群体系谱不完善造成的祖源不清问题。
在物种资源保护中，系谱记录依赖于基础群，当系谱数据记录存在问题时，会导致个体间的亲缘关系及近交系数出现错误，因此会出现近交衰退和适应能力下降等现象，而种群内的遗传多样性对于适应能力和避免近交衰退是必要的[3]。目前，单核苷酸多态性(single nucleotide polymorphism, SNP)已成为最理想的群体遗传多样性和基因组学的分子标记，而限制性酶切位点相关的DNA(re-striction-site associated DNA，RAD)测序技术能鉴定出广泛覆盖基因组范围的SNP位点，被广泛应用于资源群体的遗传多样性研究[4]。在家禽中，韩威等[5]基于RAD-seq技术在19个地方鸡种中鉴定出400 562个SNPs标记，揭示了19个品种的遗传多样性和群体结构并发掘了与免疫系统调节、生殖机能调控、应激响应等生物学过程相关的种质特性基因。在草食动物中，Liu等[6]利用RAD-seq技术研究了6个中国地方兔品种和2个引进兔品种的遗传多样性和群体结构，鉴定出了与黑素生成相关的基因。Liu等[7]采用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文库，进行EcoR I(G^AATTC)和Nla 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软件计算观察杂合度(Ho)、群体内核酸多态性(Pi)和群体的平均近交系数(Fis)。利用软件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进行选择性清除分析。采用遗传分化系数(Fst)和核苷酸多样度(π)相结合的方法，以100K为一个window，10K为步长取一个区域(群体中任意两条不同序列的碱基差异数取平均值)，Fst值按照大小排列，取大小排在前5%的区域信息，π的比值取位于前5%与后5%的区域，根据Fst与π比值筛选出的区域取交集作为关联结果，该过程使用软件vcftools进行分析。
1.8 全基因组关联分析使用TASSEL(https://tassel.bitbucket.io)软件的MLM模型对不同群体的羽色性状进行关联分析，计算模型为：y=Xα+Zβ+Wμ+e，其中，y为表型性状，X 为基因型(固定效应)的指示矩阵，α为固定效应的估计参数；Z为群体遗传结构的指示矩阵，β为SNP的效应；W为个体亲缘关系指示矩阵，μ为预测的随机个体；e是随机残差，服从e~(0，δe2)。通过关联的显著度(P-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 |
静原鸡白羽、麻羽群体内鉴定到的SNPs依次为238 533和233 562个，黑羽的SNPs数最多为240 820个。3种羽色的静原鸡Ho、Pi和Fis分别在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系数r2=0.1时，各类群对应的LD衰减距离基本相等(图 1)。
|
|
表 2 静原鸡3种羽色遗传参数的比较 Table 2 Comparison of genetic parameters of 3 feather colors of Jingyuan chicken |
|
图中横坐标表示连锁不平衡发生的距离，纵坐标是连锁不平衡相关系数 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 |
主成分分析的结果中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 |
通过结合π和Fst来选择前5%的区域。以静原鸡白羽为试验组，黑羽为正向选择组，共获得819个正向选择区域(图 5A)。以白羽为试验组，麻羽为正向选择组，共获得正向选择区域730个(图 5B)。以黑羽为对照组，麻羽为正向选择组，检测到719个正向选择区域(图 5C)。对筛选区域进行基因注释之后利用KOBAS对检测的候选基因进行KEGG富集分析，KEGG富集结果表明候选基因显著富集到37条通路(P < 集最显著的是黑色素生成途径，此外还有酪氨酸代谢和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 |
|
图 6 KEGG通路富集分析 Fig. 6 KEGG pathway enrichment analysis |
|
|
表 3 选择性清除分析的选择区域及候选基因 Table 3 Selection regions and candidate genes for selective sweep analysis |
通过对静原鸡羽色进行全基因组关联分析后发现，白羽性状关联到66个显著相关的SNPs位点(P＜1.48×10-7)，分布在1、4、8、15号4条染色体上(图 7A)，临近或坐落于TYR、GABRE、SRPX2等30个基因，单个关联标记位点解释的表型变异(R2)范围为3.16%~9.18%，前10个基因见表 4；麻羽性状鉴定到259个显著相关SNPs位点(P＜1.48×10-7)，这些位点主要分布在1~9号染色体上(图 7C)，临近或坐落于RP2、TMEM135、IQSEC3等116个基因，单个关联标记位点解释的表型变异(R2)为6.73%~19.53%，前10个基因见表 5，黑羽性状关联到200个显著相关的SNPs位点(P＜1.48×10-7)，与黑羽性状相关的SNPs主要位于1~9号染色体上(图 7E)，临近或坐落于FGD4、DOCK4、TMEM135等114个基因，单个关联标记位点解释的表型变异(R2)为5.85%~14.85%，前10个基因见表 6。QQ图是关联分析结果重要的质控图，λ能够判断群体是否分层，λ在1.05以上说明群体存在分层，本研究的λ值均小于1.05说明个体均匀分布，不存在群体分层现象(图 7B、7D、7F)。通过文献查阅和KOBAS对关联的基因进行KEGG富集分析(图 8)，在白羽和麻羽静原鸡中分别筛选出2个(TYR、TH)与羽色相关的基因，在黑羽静原鸡中筛选出两个基因(FZD6、PLCB2)与羽色相关，这些基因主要在黑色素合成和酪氨酸代谢中发挥重要的调控作用(表 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 |
|
|
表 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 |
|
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 |
|
|
表 7 GWAS筛选的候选基因 Table 7 Candidate genes screened by GWAS |
随着生产技术及生产力的发展，畜禽经过长期的自然选择、遗传漂变和人工选择，其遗传多样性随之发生改变[8-9]。静原鸡作为国家级地方保护资源，其遗传多样性同样受到了严重威胁。2013年，固原
市彭阳县建立了静原鸡原种保种场，对静原鸡进行系统保护和开发利用，但是在这期间由于保种知识的欠缺，导致系谱数据的大量缺失和错误，因此本研究通过分子保种技术保护静原鸡群体的遗传多样性并持续监测保种群的状态，评估和制定有效的保种技术手段，对静原鸡以及我国地方鸡种的保护与开发利用具有重要意义。Ho是反映等位基因频率的一个指标，其数值越接近，表明群体的分布越均匀[7]。Fis可以明确近交程度，对于选种选配防止近交衰退和品系繁育具有重要的指导意义[10]。Pi越高，说明物种的多样性越高[6]。本研究利用RAD-seq在静原鸡共检测出712 915个SNPs，3个类群的Ho、Pi和Fis分别在0.273 2~0.278 2、0.304 9~0.309 6和0.096 1~0.109 8之间，与韩威等[5]所研究的19个地方鸡种相比，静原鸡的Ho和Pi均高于19个地方鸡种，Fis低于瓢鸡、河南斗鸡、金湖乌凤鸡等8个地方鸡种，由此表明静原鸡遗传多样性丰富，品系选育良好，资源活力好，开发利用潜力较大[2]。LD是评价遗传多样性的一个重要指标[7]。本研究中，3个类群LD系数的衰减率差异较小，说明该群体没有经过较强的选择，受选择的强度低。主成分分析中除麻羽静原鸡144号和黑羽静原鸡180号两个个体之外，其他个体根据羽色聚在一起，这与进化树和群体结构分析的结果一致，说明静原鸡根据羽色分成了不同的类群，这为静原鸡的分群保种提供了理论依据。
3.2 选择性清除分析物种的进化离不开选择，在选择的过程中，群体内的有利突变发生后，这个突变基因的适合度越高，越容易被固定，与此基因座连锁的染色体区域，由于搭车效应也被固定下来，大片紧密连锁的染色体区域因此失去多态性，这种由于搭车效应引起多态性下降的现象，遗传上称为选择清除[11]。通过选择清除分析，可以鉴定出差异基因组区域，结合功能注释进一步预测表型相关的候选基因[12]。本研究通过选择性清除分别获得819、730、719个选择区域，这些区域内筛选到与静原鸡羽色相关的候选基因共11个(FZD4、WNT16、EDNRB、TYR、KRAS、CTNNB1、DDC、MC1R、CAMK2A、PRKCB、PRKCA)。TYR是一种酪氨酸酶，它的表达缺乏将导致白羽球茎黑色素生物合成不足，是白羽形成的直接原因[13]。李洪林[14]鉴定出控制兴义矮脚鸡隐性白羽的主效基因为TYR。EDNRB在色素细胞的增殖和分化中起重要作用，包括黑素母细胞和黑素细胞，研究表明EDNRB2在番鸭黑羽中的表达量高于白羽[15-16]。FZD4是一种Wnt蛋白的受体，通过Wnt信号途径调控黑色素生成[17]。WNT16是WNT信号通路中的重要基因，在神经嵴来源的黑素细胞发育和迁移等过程中发挥重要作用[18]。Wnt1或Wnt3a表达缺陷的小鼠表现出神经嵴黑色素细胞形成的丧失[19]。KRAS具有内在的GTPase活性，在调节细胞增殖中起重要作用，研究表明KRAS基因突变会造成黑色素细胞增生，从而导致黑色素瘤[20]。MITF是色素沉着的主要调节剂，是Wnt信号通路的靶标，降低CTNNB1的表达可能会下调MITF的表达，最终减少黑色素的合成[21]。在黑色素合成途径中，前体多巴(3, 4-二羟基苯丙氨酸)和多巴胺由TH和DDC合成[22]。MC1R编码一个七跨膜结构域G蛋白偶联受体，主要在发育中的羽毛和头发的黑色素细胞中表达[23]。Nam等[24]研究发现，位置和TYR表达水平的差异对韩国本土鸡品种的羽毛色素沉着的影响比MC1R表达水平的变化更大。CAMK2A是钙调蛋白依赖性蛋白激酶，研究表明，CAMK2A是调控中国五指山黑猪毛色的候选基因[25]。PRKCA和PRKCB是两种蛋白激酶，PRKCA和PRKCB的激活是TPA抑制黑色素合成的原因，研究已证实了PRKCA的表达与小鼠皮肤色素沉着大致相关[26]。综上所述，调控静原鸡羽色的基因主要与黑色素合成相关。
3.3 全基因组关联分析全基因组关联分析是利用自然群体探测物种的遗传变异，进而分析复杂性状遗传结构的有力工具[27]，近年来已在鸡、牛、猪、羊等物种中均有广泛研究[28-31]。本研究针对静原鸡的3种羽色性状开展全基因组关联分析，初步鉴定到影响静原鸡群体白、麻、黑3种羽色性状相关的致因突变位点和候选基因。其中鉴定到与白羽性状显著关联的位点7个，对这些基因进行功能注释，发现TYR富集在酪氨酸酶通路，该通路在合成黑色素和其他多酚化合物等色素时起重要作用[13]，关联分析结果表明，有2个显著关联位点189168614、189176762位于TYR基因的第3内含子上。该基因被选择性清除和全基因组关联分析同时筛选到，由此推测TYR是控制静原鸡白羽的关键候选基因，具体调控过程还需进一步验证。另外，对黑羽显著关联水平SNPs分析发现919883和129265651位点与黑色素形成有关，这两个位点分别位于PLCB2基因的第2内含子和FZD6的第2内含子上，PLCB2是一种磷酸肌醇磷脂酶，FZD6是一种Wnt蛋白的受体。Zhang等[32]通过逆转录聚合酶链式反应(RT-PCR)分析了PLCB2基因在各种癌症系中的mRNA表达情况，同时结合其他生物学手段发现PLCB2是影响黑色素瘤增殖和凋亡的关键因素。Korstanje等[33]基于生物信息学分析确定了两个导致色素沉着缺陷的寡聚体复合物FZD6和MC5R，研究发现FZD6敲除将导致色素沉着障碍，由此得出FZD6能够间接影响色素沉着。但有关PLCB2和FZD6对羽色的调控研究尚未见报道。除此之外与麻羽性状显著关联的位点有178个，其中与羽色相关的只有13923834位点，该位点位于TH的9号内含子上，TH是酪氨酸3-羟化酶，Fernstrom等[34]研究表明，TYR在TH酶的催化下，能生成DOPA，而DOPA是黑色素形成的关键物质。由此说明TH可能是调控羽色的基因。综上所述，静原鸡3种羽色形成的分子机制不同，遗传基础复杂，以现有的资源群体并不能准确关联到所有致因突变和潜在候选基因，在后续研究中将扩大资源群体规模、优化模型和方法，进一步对相关突变位点进行挖掘和验证。
4 结论本研究利用RAD-seq技术对白羽、麻羽和黑羽静原鸡群体进行了全基因组范围内的遗传变异检测，基于鉴定到的SNPs计算群体遗传学统计指标，系统评价了静原鸡群体的遗传多样性和遗传结构，通过选择性清除和全基因组关联分析鉴定了与羽色相关的FZD4、WNT16、EDNRB、TYR、KRAS、CTNNB1、DDC、MC1R、CAMK2A、PRKCB、PRKCA基因。综上所述，本研究结果将对静原鸡这一地方种质资源群体的保护奠定理论基础，同时为静原鸡不同羽色的品系化培育提供新的遗传标记和基因靶点。
参考文献
| [1] |
母童. 静原鸡ELOVL2和ELOVL5基因组织表达特性及生物信息学分析[D]. 银川: 宁夏大学, 2019. MU T. Tissue-specific expression and bioinformatics analysis of ELOVL2 and ELOVL5 gene strains in Jingyuan chicken[D]. Yinchuan: Ningxia University, 2019. (in Chinese) |
| [2] |
额尔和花, 丁伟, 杨文清, 等. 宁夏地方品种鸡的开发与利用[J]. 畜牧与饲料科学, 2013, 34(7-8): 106-107. E E H H, DING W, YANG W Q, et al. Development and utilization of local breed chicken in Ningxia[J]. Animal Husbandry and Feed Science, 2013, 34(7-8): 106-107. (in Chinese) |
| [3] |
OLIEHOEK P A, BIJMA P. Effects of pedigree errors on the efficiency of conservation decisions[J]. Genet Sel Evol, 2009, 41(1): 9. DOI:10.1186/1297-9686-41-9 |
| [4] |
李国辉, 魏清宇, 张会永, 等. 基于RAD-seq技术分析边鸡保种群体的遗传进化[J]. 农业生物技术学报, 2019, 27(12): 2198-2206. LI G H, WEI Q Y, ZHANG H Y, et al. The Analysis on genetic evolution of Bian chicken (Gallus gallus) population based on RAD-seq technology[J]. Journal of Agricultural Biotechnology, 2019, 27(12): 2198-2206. (in Chinese) |
| [5] |
韩威, 朱云芬, 殷建玫, 等. 基于RAD-seq简化基因组测序的19个地方鸡种遗传进化研究[J]. 畜牧兽医学报, 2020, 51(4): 670-678. HAN W, ZHU Y F, YIN J M, et al. Study on genetic evolution of 19 indigenous chicken breeds based on RAD-seq[J]. Acta Veterinaria et Zootechnica Sinica, 2020, 51(4): 670-678. (in Chinese) |
| [6] |
LIU C M, WANG S H, DONG X G, et al. Exploring the genomic resources and analysing the genetic diversity and population structure of Chinese indigenous rabbit breeds by RAD-seq[J]. BMC Genomics, 2021, 22(1): 573. DOI:10.1186/s12864-021-07833-6 |
| [7] |
LIU C M, CHEN H L, REN Z J, et al. Development of Genomic resources and identification of genetic diversity and genetic structure of the domestic Bactrian Camel in China by RAD sequencing[J]. Front Genet, 2020, 11: 797. DOI:10.3389/fgene.2020.00797 |
| [8] |
WANG M S, THAKUR M, PENG M S, et al. 863 genomes reveal the origin and domestication of chicken[J]. Cell Res, 2020, 30(8): 693-701. DOI:10.1038/s41422-020-0349-y |
| [9] |
ZHANG M M, HAN W, TANG H, et al. Genomic diversity dynamics in conserved chicken populations are revealed by genome-wide SNPs[J]. BMC Genomics, 2018, 19(1): 598. DOI:10.1186/s12864-018-4973-6 |
| [10] |
史良玉, 王立刚, 张鹏飞, 等. 不同来源大白猪总产仔数近交衰退评估[J]. 畜牧兽医学报, 2021, 52(10): 2772-2782. SHI L Y, WANG L G, ZHANG P F, et al. Evaluation of inbreeding depression on the total numbers of piglets born in different groups of large white pigs[J]. Acta Veterinaria et Zootechnica Sinica, 2021, 52(10): 2772-2782. DOI:10.11843/j.issn.0366-6964.2021.010.008 (in Chinese) |
| [11] |
STEPHAN W. Selective sweeps[J]. Genetics, 2019, 211(1): 5-13. DOI:10.1534/genetics.118.301319 |
| [12] |
吕世杰, 陈付英, 金磊, 等. 利用简化基因组测序筛选安格斯牛生长相关的受选择基因[J]. 畜牧兽医学报, 2020, 51(4): 713-721. LV S J, CHEN F Y, JIN L, et al. Identification of growth-related genes under selection in angus cattle using SLAF-seq[J]. Acta Veterinaria et Zootechnica Sinica, 2020, 51(4): 713-721. (in Chinese) |
| [13] |
YU S, LIAO J, TANG M, et al. A functional single nucleotide polymorphism in the tyrosinase gene promoter affects skin color and transcription activity in the black-boned chicken[J]. Poult Sci, 2017, 96(11): 4061-4067. DOI:10.3382/ps/pex217 |
| [14] |
李洪林. 贵州地方鸡种隐性白羽与慢羽分子标记研究与应用[D]. 贵州: 贵州大学, 2019. LI H L. Study and Application of recessive white feather and slow feather molecular markers in Guizhou local chicken breeds[D]. Guizhou: Guizhou University, 2019. (in Chinese) |
| [15] |
WU N, QIN H, WANG M, et al. Variations in endothelin receptor B subtype 2 (EDNRB2) coding sequences and mRNA expression levels in 4 Muscovy duck plumage colour phenotypes[J]. Br Poult Sci, 2017, 58(2): 116-121. DOI:10.1080/00071668.2016.1259531 |
| [16] |
SALDANA-CABOVERDE A, KOS L. Roles of endothelin signaling in melanocyte development and melanoma[J]. Pigment Cell Melanoma Res, 2010, 23(2): 160-170. DOI:10.1111/j.1755-148X.2010.00678.x |
| [17] |
DAVOODI P, EHSANI A, TORSHIZI R V, et al. New insights into genetics underlying of plumage color[J]. Anim Genet, 2022, 53(1): 80-93. DOI:10.1111/age.13156 |
| [18] |
BELLEI B, PITISCI A, CATRICALÀ C, et al. Wnt/β-catenin signaling is stimulated by α-melanocyte-stimulating hormone in melanoma and melanocyte cells: implication in cell differentiation[J]. Pigment Cell Melanoma Res, 2011, 24(2): 309-325. DOI:10.1111/j.1755-148X.2010.00800.x |
| [19] |
RABBANI P, TAKEO M, CHOU W, et al. Coordinated activation of Wnt in epithelial and melanocyte stem cells initiates pigmented hair regeneration[J]. Cell, 2011, 145(6): 941-955. DOI:10.1016/j.cell.2011.05.004 |
| [20] |
NAGARATHINAM S, BAALANN K P. Nevus pigmentosus et pilosus[J]. Pan Afr Med J, 2021, 40: 107. |
| [21] |
LEE J E, KIM S Y, JEONG Y M, et al. The regulatory mechanism of melanogenesis by FTY720, a sphingolipid analogue[J]. Exp Dermatol, 2011, 20(3): 237-241. DOI:10.1111/j.1600-0625.2010.01148.x |
| [22] |
XIAO D, CHEN X, TIAN R B, et al. Molecular and potential regulatory mechanisms of melanin synthesis in Harmonia axyridis[J]. Int J Mol Sci, 2020, 21(6): 2088. DOI:10.3390/ijms21062088 |
| [23] |
MUNDY N I. A window on the genetics of evolution: MC1R and plumage colouration in birds[J]. Proc Biol Sci, 2005, 272(1573): 1633-1640. |
| [24] |
NAM I S, OH M G, NAM M S, et al. Specific mutations in the genes of MC1R and TYR have an important influence on the determination of pheomelanin pigmentation in Korean native chickens[J]. J Adv Vet Anim Res, 2021, 8(2): 266-273. DOI:10.5455/javar.2021.h511 |
| [25] |
XU Q, LIU X M, CHAO Z, et al. Transcriptomic analysis of coding genes and non-coding RNAs reveals complex regulatory networks underlying the black back and white belly coat Phenotype in Chinese Wuzhishan pigs[J]. Genes (Basel), 2019, 10(3): 201. DOI:10.3390/genes10030201 |
| [26] |
LI L F, CHANG B D, RONINSON I B, et al. Alteration of skin protein kinase C α protein and mRNA levels during induced mouse hair growth[J]. J Dermatol, 1999, 26(4): 203-209. DOI:10.1111/j.1346-8138.1999.tb03457.x |
| [27] |
周庆红, 周灿, 郑伟, 等. 甘蓝型油菜角果长度全基因组关联分析[J]. 中国农业科学, 2017, 50(2): 228-239. ZHOU Q H, ZHOU C, ZHENG W, et al. Genome wide association analysis of silique length in Brassica napus L[J]. Scientia Agricultura Sinica, 2017, 50(2): 228-239. (in Chinese) |
| [28] |
LI Z Z, ZHENG M, ALI ABDALLA B, et al. Genome-wide association study of aggressive behaviour in chicken[J]. Sci Rep, 2016, 6(1): 30981. DOI:10.1038/srep30981 |
| [29] |
段星海, 安炳星, 杜丽丽, 等. 中国肉用西门塔尔牛生长曲线参数的全基因组关联分析[J]. 畜牧兽医学报, 2021, 52(5): 1267-1277. DUAN X H, AN B X, DU L L, et al. Genome-wide association study of growth curve parameters in Chinese Simmental beef cattle[J]. Acta Veterinaria et Zootechnica Sinica, 2021, 52(5): 1267-1277. (in Chinese) |
| [30] |
周李生, 杨杰, 刘先先, 等. 苏太猪宰后72 h pH和肉色性状的全基因组关联分析[J]. 中国农业科学, 2014, 47(3): 564-573. ZHOU L S, YANG J, LIU X X, et al. Genome-wide association analyses for musle pH 72 h value and meat color traits in Sutai Pigs[J]. Scientia Agricultura Sinica, 2014, 47(3): 564-573. DOI:10.3864/j.issn.0578-1752.2014.03.016 (in Chinese) |
| [31] |
王晓薇, 马青. 滩羊毛色的全基因组关联分析[J]. 浙江农业学报, 2020, 32(1): 28-34. WANG X W, MA Q. Genome-wide association studies for coat color in Tan sheep[J]. Acta Agriculturae Zhejiangensis, 2020, 32(1): 28-34. (in Chinese) |
| [32] |
ZHANG H H, XIE T, SHUI Y J, et al. Knockdown of PLCB2 expression reduces melanoma cell viability and promotes melanoma cell apoptosis by altering Ras/Raf/MAPK signals[J]. Mol Med Rep, 2019, 21(1): 420-428. |
| [33] |
KORSTANJE R, DESAI J, LAZAR G, et al. Quantitative trait loci affecting phenotypic variation in the vacuolated lens mouse mutant, a multigenic mouse model of neural tube defects[J]. Physiol Genomics, 2008, 35(3): 296-304. DOI:10.1152/physiolgenomics.90260.2008 |
| [34] |
FERNSTROM J D, FERNSTROM M H. Tyrosine, phenylalanine, and catecholamine synthesis and function in the brain[J]. J Nutr, 2007, 137(6 Suppl 1): 1539S-1547S. |
(编辑 郭云雁)