儋州鸡体尺性状全基因组关联分析

姜宏正，杨德智，马中华，荀文娟，施力光，侯冠彧

(1.海南大学动物科技学院，海口 570228；2.中国热带农业科学院热带作物品种资源研究所，海口 571101)

儋州鸡又名儋州小种鸡、北岸小种鸡、石鸡等，原产于海南省儋州市，由当地人民经过长期的驯养而来，具有耐粗饲、抗病力强、肉质鲜嫩、营养价值高等种质特性。然而作为新发现的地方遗传资源，缺乏对儋州鸡系统地选育。全基因组关联分析(genome wide association study，GWAS)作为后基因组时代有效的遗传变异(标记)多态性检测方法，在畜禽育种工作中应取得显著的效果。Huang等[1]利用600K芯片对中国地方鸡的胴体性状进行了GWAS分析，发现了91个与鸡的胴体性状关联潜在显著的SNPs位点，筛选出细胞角化包膜前体(sciellin，SCEL)和MYC结合蛋白2(MYC binding protein 2，E3 ubiquitin protein ligase，MYCBP2)等6个可能的候选基因，丰富了鸡胴体性状遗传机理的认知。Bai等[2]利用鸡600K高密度SNP芯片对喙畸性状进行GWAS分析，结果发现8个与喙畸形性状关联的位点，类似转录延伸因子B多肽3(transcription elongation factor B polypeptide 3-like，TCEB3)、Tudor结构域蛋白3(tudor domain containing 3，TDRD3)、RET原癌基因(ret proto-oncogene，RET)和微管解聚蛋白1(stathmin 1，STMN1) 4个基因可能是鸡喙畸形性状研究的重要候选基因，得到钙离子信号调控通路在内的6条与喙畸形性状显著关联的通路。Liu等[3]利用SNP芯片对产蛋后期的蛋品质进行了GWAS分析，发现13号染色体约有0.36 Mb(8.95～9.31 Mb)区域与蛋白高度和哈氏单位显著相关，Ras同源基因家族成员A(ras homolog family member A，RHOA)、基质细胞衍生因子4 (stromal cell derived factor 4，SDF4)和受体超家族成员4 (TNF receptor superfamily member 4，TNF4)这3个基因可能是与蛋壳颜色相关的候选基因。体尺性状作为评价家禽体型外貌的主要指标，是衡量其生长发育和品种鉴定的重要特征。由于家禽的体尺性状与其产肉或产蛋性能密切相关，通过培育代次的体尺测量和选种来进行经济性状表型选择，可有效缩短育种进程。张建丰等[4]研究了中国地方鸡的体尺性状对繁殖性状的影响，发现胫长越长平均蛋重越大,体斜长越长年产蛋数越多,龙骨长越长鸡的平均开产日龄越晚。李万贵等[5]对兴义鸭体尺性状和屠宰性能进行相关性分析，结果发现，胸宽与全净膛重、半净膛重、肌胃重呈显著正相关。徐明明等[6]探讨了雪峰乌骨鸡体尺性状与产肉性能的相关关系，发现可通过胸深这一体尺指标预测雪峰乌骨鸡公鸡的屠宰性能。目前，利用GWAS对地方鸡体尺性状的研究尚不全面。本试验利用GWAS技术对70日龄儋州鸡体尺性状进行分析，旨在通过GWAS识别影响儋州鸡体尺性状的分子标记，探讨遗传标记的分布规律，寻找影响儋州鸡体尺性状的候选基因，为儋州鸡早期选种育种及开发利用奠定理论基础。

1 材料与方法

1.1 试验动物与饲养管理

随机选取200只同一批次孵化的1日龄儋州鸡(中国热带农业科学院热带作物品种资源研究所提供)，公母各半，按照《NY/T 33-2004鸡饲养标准》[7]进行饲养。

1.2 主要试剂和仪器

动物血液/细胞/组织基因组DNA提取试剂盒(DP304)购自天根生化科技(北京)有限公司。紫外分光光度计(NanoDrop 2000)购自Thermo Fisher Scientific公司。

1.3 样品采集和分析

1.3.1 血样采集及DNA提取 70日龄时，翅静脉采血，EDTA抗凝处理，使用动物血液/细胞/组织基因组DNA提取试剂盒提取基因组DNA，利用紫外分光光度计对DNA样品的浓度和纯度进行测定，测定标准为A260 nm/A280 nm值为1.8～2.0。质量合格的样品送杭州联川生物技术股份有限公司进行10×深度的全基因组重测序。

1.3.2 体尺性状测定 儋州鸡饲养至70日龄时，按照《NY/T 823-2004家禽生产性能名词术语和度量统计方法》[8]进行体尺性状测定，包括体斜长、龙骨长、胫长、胫围、胸宽、胸深等。体斜长测定肩关节至坐骨结节间的距离；龙骨长测定龙骨突前端到龙骨末端的距离；胫长测定胫部上关节到第3、4趾间的直线距离；胫围测定胫骨中部的周长；胸深测量第1胸椎到龙骨前缘的距离。测得数据经过Excel整理，进行GWAS分析。

1.4 测序数据的质控

使用FastP(0.19.3)对原始数据(Raw data)进行预处理。主要步骤为：去除接头(Adapter)；去除含有N(N表示无法确定碱基信息)的比例>5%的reads；去除低质量reads；利用BWA将Clean data比对到参考基因组上，并使用GATK软件进行单核苷酸多态性(single nucleotide polymorphysm，SNP)的检测，并对检测到的变异位点进行质量过滤。SNP的质量过滤标准：QD>2.0、FS>60.0、MQ>40.0、MQRankSum>-12.5、ReadPosRankSum>-8.0、SOR>3.0。

1.5 GWAS分析

1.5.1 关联分析模型 使用EMMAX[9](versionbeta-07)中的混合线性模型，对儋州鸡体尺性状与全基因组上SNP进行GWAS分析，混合线性模型为：

y=μ+Xb+u+e

其中，y表示儋州鸡体尺性状表型值；μ表示性状的平均值；X表示固定效应矩阵，其中固定效应包括养殖环境、育种方式和性别；b为固定效应向量；u表示剩余多基因效应；e表示表型值的随机残差，u和e均服从正态分布。

1.5.2 遗传进化与群体结构分析 基于SNP，通过Admixture(1.3.0)软件，分析所有样本的群体结构，分别假设分群数(K值)为1～10，进行聚类，预测样本之间遗传距离，对关联分析进行辅助。使用SPAGeDi(1.4)对所有样本两两之间的亲缘关系进行计算。

1.5.3 多重检验 利用Bonferroni方法对SNPs的显著性进行判定。若SNP的值<0.01/N的值(N表示质控后的SNP数目)，则SNP与性状显著关联；若SNP的值<0.1/N的值，则表示SNP与性状潜在显著关联。

1.6 基因注释及候选基因筛选

利用Ensembl (http:∥oct2018.archive.ensembl.org/Gallus-gallus/Info/Index)和NCBI(http:∥www.ncbi.nlm.nih.gov/)搜索显著及潜在关联SNP位点前后各100 kb区域的上、下游基因，将所得到的基因与GO、KEGG等数据库进行比对，根据基因注释情况分析可能的候选基因。

2 结 果

2.1 儋州鸡体尺性状测定结果

儋州鸡体尺性状统计结果见表1。70日龄儋州鸡平均胫长66.55 mm，胫围20.37 mm，体斜长89.06 mm，胸宽27.81 mm，髋骨宽43.23 mm，胸深44.57 mm，龙骨长63.62 mm。

表1 儋州鸡体尺性状描述性统计(n=200)

2.2 儋州鸡遗传进化与群体结构分析

利用Admixture软件评估群体分层，根据ΔK的分析结果确定最优分群数为1(图1)，表明试验群体未发生亚群分化； 亲缘关系热图如图2所示， 发现群体未出现家系分化，适宜进行后续的GWAS分析。

A中每种颜色代表一个群，每行代表一个分群值的情况

图2 亲缘关系热图

2.3 儋州鸡体尺性状GWAS结果

采用EMMAX模型对儋州鸡7个体尺指标性状(胫长、体斜长、胫围、胸宽、胸深、髋骨宽、龙骨长)GWAS分析结果如图3～9所示。胫长和胫围性状达到基因组显著水平的SNPs位点分别有12和8个。其中与胫长性状相关的SNPs分别定位于1、2、4和8号染色体上；与胫围性状相关的SNPs定位于2、4、8和13号染色体上。未发现与其他5个性状显著关联的SNP位点。

图3 胫长的全基因组关联分析曼哈顿图

图4 胫围的全基因组关联分析曼哈顿图

图5 体斜长的全基因组关联分析曼哈顿图

图6 胸宽的全基因组关联分析曼哈顿图

图7 髋骨宽的全基因组关联分析曼哈顿图

图8 胸深的全基因组关联分析曼哈顿图

图9 龙骨长的全基因组关联分析曼哈顿图

2.4 候选基因功能注释

各显著SNP上、下游预测的基因如表2、3所示，发现8个基因可作为胫长和胫围性状的候选基因，分别是KCNA1、TPK1、EZH2、FSTL5、AMY2A、TGFBI、LECT2和IL-9，其中FSTL5、AMY2A和TPK1基因可作为同时影响胫长和胫围性状的候选基因。使用GO和KEGG数据库对其进行分析，8个基因参与钾离子跨膜转运、硫胺素新陈代谢、细胞增殖、钙离子结合、骨骼肌卫星细胞维持与骨骼肌再生、细胞受体相互作用、生长因子活性等生物学进程(表4)。

表2 胫长性状显著相关的SNPs位点

表3 胫围性状显著相关的SNPs位点

表4 影响儋州鸡胫长和胫围的关键候选基因

3 讨 论

体尺性状是评价畜禽生长性能的重要指标。Emrani等[10]研究发现，影响鸡不同周龄胫骨发育的显著SNP位点分布于1、4、7、9、21和Z号染色体上。孙艳发等[11]使用60K SNP芯片对F2资源群体的胫长和胫围研究发现，与胫长相关的SNPs位点主要分布于1、9、13和27号染色体上，4号染色体上71.01—85.94 Mb区域是与鸡胫围相关的重要候选区域。本试验定位到与儋州鸡胫长胫围性状相关的SNPs位于1、2、4和8号染色体上，大部分位点是未报道的新位点，这可能与儋州鸡地方品种的遗传多样性及所采用的EMMAX分析方法差异所致。本试验未预测到与胸深、胸宽等性状显著相关的SNP，可能是由于样本量偏小导致检验功效较低造成的。但是基于GWAS分析自身对微效多基因控制的数量性状检测能力不足的局限性[12]，许多贡献度偏小的性状基因在GWAS分析中往往因计算模型的选择差异而达不到显著水平，进而不能准确挖掘显著性SNP，或者得到假阳性的显著性结果。复杂性状受到基因-基因、基因-环境相互作用的影响，GWAS分析很难将上述重要的影响因素考虑在内。

Zhang等[13]在基于单倍型的GWAS分析发现，8号染色体上10.13-12.97 Mb区域是与鸡胫围相关的重要候选区域。本试验预测到在儋州鸡8号染色体0.41 Mb附近的AMY2A基因与胫围性状显著相关。AMY2A是α-淀粉酶家族的一员，编码胰腺产生的淀粉酶同工酶。Endo等[14]研究也发现，抗AMY2A自身抗体可以作为自身免疫性胰腺炎和暴发性Ⅰ型糖尿病标志物。本试验中还检测到在儋州鸡13号染色体中有4个与胫围性状相关的SNPs位点达到基因组显著水平，定位到IL9、TGFBI和LECT2共3个候选基因。TGF-β超家族信号通路参与了广泛的生物学过程，对调控早期胚胎发育、细胞生长、干细胞自我更新、肿瘤发生发展等具有十分重要的调控作用。Lorda-Diez等[15]研究发现，TGFBI基因可促进TGF-β信号转导的促纤维化作用，诱导合成的细胞外基质蛋白还能中和缺氧诱导因子对软骨细胞的影响，参与细胞-细胞及细胞-基质的黏附。LECT2基因能够通过介导Wnt/β-Catenin信号通路促进间质干细胞的成骨分化[16]。IL-9是白细胞介素家族的一员，作为机体重要的免疫因子，可由Th9、Th17等多种细胞分泌,通过与受体结合发挥生物学作用。Biscarini等[17]研究发现，IL-9基因与鸡的羽毛状况相关，揭示了免疫系统和啄羽行为之间的关系。 以上结果说明,AMY2A、IL-9、TGFBI和LECT2基因的相关研究多集中在机体免疫和细胞分化方面,在动物尤其是家禽方面的研究较少,因此本试验所预测到的基因对体尺性状的影响需要进一步进行功能验证筛选。

本试验在儋州鸡1、2、4号染色体找到显著的SNP位点,定位到EZH2、TPK1、FSTL5和KCNA1 4个候选基因。 Akizu等[18]研究发现，EZH2基因对于维持脊椎动物早期发育的神经管内环境平衡是必不可少的，也与体温调节系统的发育和功能有关[19]。FSTL5是卵泡抑素样家族的一个成员，编码一种分泌性糖蛋白，在细胞生长和分化方面具有重要作用，其表达量的升高会促进细胞的凋亡[20]。Zhang等[21]研究证实，FSTL5基因在肝细胞癌中起抑制作用，表明FSTL5基因表达下调可以激活Wnt/β-Catenin信号通路，进而提高肝细胞癌细胞的生长和存活能力。KCNA1是钾离子通道蛋白家族成员之一，调节脑和肾中枢神经系统的钾离子通道，介导跨膜钾离子在兴奋性膜的传导，以及调节膜电位和神经信号，以防止神经元过度兴奋[22]。研究还发现，KCNA1基因与肿瘤细胞的细胞增殖、凋亡等过程显著相关[23-24]。TPK1是硫胺素代谢过程中的关键酶之一，催化硫胺素磷酸化为二磷酸硫胺和三磷酸硫胺。Fradin等[25]研究结果表明，胎儿和母体TPK1基因的基因组变异可能导致正常人出生体重的变化。这4个候选基因与细胞增殖、分化和生长调控相关,由于体尺性状是复杂性状,受到基因-基因、基因-环境的相互作用的影响，仍需在其他品种及更大的群体中进行后续的功能验证。本研究得到的候选基因和位点多次为首次发现，这可能是品种及分析方法不同所致，后续仍有待进行更深入的研究。

4 结 论

本研究利用10×重测序对儋州鸡体尺性状进行GWAS分析，发现了20个SNPs位点与体尺性状显著相关，并筛选到KCNA1、EZH2、TGFBI、LECT2IL-9、TPK1、FSTL5和AMY2A共8个目标性状候选基因,为儋州鸡的育种工作提供候选的分子标记，为地方品种鸡标记辅助选择提供了新的思路。