范晨宇,单艳菊,章 明,姬改革,巨晓军,屠云洁,贺 喜,束婧婷,刘一帆,3*,张海涵*
(1.湖南农业大学动物科学技术学院,湖南家禽安全生产工程技术研究中心,长沙 410128;2.江苏省家禽科学研究所,江苏省家禽遗传育种重点实验室,扬州 225125;3.江苏省家禽科学研究所科技创新有限公司,扬州 211412)
黄羽肉鸡产业是我国畜牧业的支柱产业之一,出栏数约占肉鸡总数的一半,为满足城乡居民肉食品供应和养殖户增收做出了重要的贡献。根据生长速度和出栏日龄,黄羽肉鸡一般分为快速型(49～70 d)、中速型(71～90 d)和慢速型(91～180 d)三类[1]。快速型黄羽肉鸡生产时间短,饲料报酬高,具有最低的生产成本,并且肉质优于白羽肉鸡,因此越来越受到消费者的喜爱。随着肉鸡消费逐步转向冰鲜上市,快速型黄羽肉鸡能够更好地适应屠宰需求,市场份额逐渐增加,其选育提高也越来越受到育种人员的关注。
生长速度和肉品质是影响快速型肉鸡经济效益最重要的两类性状。虽然肉鸡体重在过去几十年中得到了显着提高,但由于瓶颈效应的存在,进一步的提升存在明显困难[2-3]。在实际选育工作中,肉色、剪切力等肉质性状由于遗传力低、测量难度大等原因,常规选择遗传进展明显滞后[4-5]。因此,深入鉴定与生长和肉品质性状相关的分子标记,并基于此开展分子选育,对于加快肉鸡品种选育进展有着重要的意义。
近年来,全基因组关联分析(genome-wide associated study,GWAS)已成为研究人类疾病和畜禽经济性状遗传调控机制的重要手段[6-7]。在肉鸡上,已有许多学者针对重要经济性状开展了相关研究,包括屠宰性能[8]、体重[9-10]、抗病性状等[11-12],并报道了大量的候选SNPs和基因。然而,大部分研究集中于国外品种和部分地方品种,在快速型黄羽肉鸡品种上的研究较为有限。
本研究选择本团队与江苏立华牧业股份有限公司合作培育的立华麻黄鸡终端父系种鸡作为研究对象,该品系具有生长速度快、优质、早熟的特点,60日龄(配套系上市日龄)体重可达到2.5 kg。利用重测序获得393只公鸡的基因组遗传变异信息,通过GWAS方法获得与体重和肉品质性状显着关联的候选基因和位点,为快速型黄羽肉鸡选育提高提供理论基础。
1 材料与方法
1.1 试验动物
本研究选择393只立华麻黄终端父系公鸡作为试验材料,饲养于江苏立华牧业股份有限公司金坛种鸡场。试验鸡群均同批饲养于同一鸡舍,采用公司标准的种鸡饲粮进行饲喂,试验期自由饮水和采食。采用常规光照和免疫程序。对所有试验鸡进行翅静脉采血保存备用。
1.2 体重和肉品质测定
鸡群饲养至60日龄(上市日龄),禁食12 h后对所有个体进行称重。随后进行屠宰,统一采集左侧胸肌样品进行剪切力、系水力、pH和肉色测定。肉品质测定方法参考《NY/T 1333—2007畜禽肉质的测定》进行。剪切力测定:使用C-LM4型电脑测控式肌肉嫩度仪进行测定,胸肌样品修剪成长条状(厚0.5 cm,宽1.0 cm),测定3次取平均值。系水力测定:使用YYW-2型应变控制式无侧限压力仪对样品进行加压测定,取新鲜胸肌肉样1 g左右(w1)放置于上下压台间,并于样品上下垫8层滤纸,施压35 kg,保持5 min后,去掉压力和滤纸,测定肉样重量(w2),计算系水力(w1-w2)/w1。pH测定:使用胴体肌肉pH测定仪PH-STAR,屠宰后45 min测定胸肌pH,选取3个不同位置取平均值。肉色测定:使用NR10QC-3nh色差仪进行肉色测定,在相同环境条件下,选取3个不同位置进行3次测定,分别记录亮度值(L*)、红度值(a*)、黄度值(b*)。
1.3 全基因组重测序
使用TIANGEN基因组试剂盒提取和纯化高质量的基因组DNA,通过凝胶迁移(1%琼脂糖凝胶)测定DNA降解和污染程度,DNA纯度采用Nano Photometer分光光度计测定,DNA浓度通过NanoDrop荧光仪进行定量。将合格样品送至北京博致格雅生物科技有限公司建库,采用MGISEQ-2000平台进行全基因组重测序,测序深度为10×。获得原始测序数据后,使用FASTP(v.0.20.0)[13]进行过滤,去除低质量reads,并对得到的Clean Reads进行质量控制,后用BWA(v0.7.17)[14]软件将质控后得到的reads比对至参考基因组GRcg7b。用Picard软件将重复比对的reads移除。通过GATK4[15]软件分析获得高质量的SNP位点。
1.4 全基因组关联分析
使用Plink v1.9[16]软件对所得SNP数据进行质量控制,合格的个体和位点用于GWAS分析。筛选标准为剔除基因频率缺失高于2%的SNP(geno<0.02),个体基因频率缺失超过4%的个体(mind<0.04),次等位基因频率小于0.05(MAF<0.05),哈代温伯格频率小于0.000 1(hwe<1×10-4),杂合率平均值中偏离3倍标准差的个体(mean±3SD)。利用GEMMA(v0.98.5)[17]软件进行关联分析,模型如下:
y=Wα+Xβ+μ+e;
其中,y是表型向量;W是协方差矩阵;α是以第一、二、三主成分作为的协变量;X是固定效应矩阵;β是有效SNP位点的数量;μ代表随机效应向量;e代表残差。
使用Plink(v1.90p)软件中的参数(-indep-pairwise 50 10 0.2)来推断独立检验的有效SNP数量,最终推断出有效独立检验的SNPs为432 529个。多重检验后,采用Bonferroni校正法来设定显着阈值,全基因组水平显着阈值和潜在显着阈值分别为1.15×10-7(0.05/432 529)和2.31×10-6(1/432 529)。使用Cmplot程序包对GWAS结果进行可视化。通过Bedtools (v2.30.0)[18]软件和参考基因组GRcg7b对显着位点的上、下游100 kb进行基因注释。
2 结 果
2.1 表型描述统计
本研究使用的393只立华父系麻黄鸡7种表型的基本统计描述如表1所示。体重、剪切力、系水力、pH、肉色L值、a值和b值的变异系数分别为8.18%、22.17%、10.45%、5.57%、4.91%、69.57%、22.58%。其中a值变异达到高度水平,剪切力和b值为中等变异性状,表明这些性状具有较大的遗传选择潜力。本群体中体重的变异系数较低,可能由于该品系经过系统的体重选择导致。
2.2 全基因组SNP质量控制
本研究获得的原始下机数据共包含16 411 878个SNPs和393个个体,首先剔除SNP缺失率大于2%的位点1 767 209个;2个个体因缺失率大于2%而被剔除;以MAF<0.05和HWE<1×10-4为标准,剔除位点106 298个;根据杂合度剔除个体3个。经过质控后,剩余388个个体和10 266 594个SNPs可用于后续的关联分析。结合异常表型值的剔除,最终用于各性状GWAS分析个体数为:体重382个,剪切力361个,系水力354个,pH 360个,L值365个,a值362个,b值364个。
2.3 体重性状全基因组关联分析
利用GEMMA软件构建混合线性模型对382个个体的体重性状进行关联分析。如图1和表2所示,本研究鉴定到2个与体重显着关联的SNPs位点,分别位于4号和28号染色体上,通过对SNP上、下游基因进行搜索注释,发现了11个基因位于与体重显着关联的SNPs位点上、下游100 kb的区域,包括FSTL3、DOCK11、IL13RA1、STK26等基因。而与体重潜在显着关联的位点有52个,主要分布于4、5、19、28号染色体上,位于NPFF、PFKM、MBNL3等基因附近。
图1 体重性状的全基因组关联分析Fig.1 Genome-wide association study of body weight trait
表2 关联SNPs注释结果Table 2 Associated SNPs annotation results
2.4 肉品质性状全基因组关联分析
肉品质关联分析结果如图2-7和表2所示。结果发现位于4号染色体的1个位点与剪切力性状显着相关;还发现15个位点与剪切力性状存在潜在显着关联,主要分布于4 号和19号染色体上,位于SLC13A5等基因附近。与系水力显着相关的位点有1个,位于17号染色体PTGS1基因附近;与系水力存在潜在关联的位点有25个,主要分布在1号、14号和17号染色体,位于MDFIC等基因附近。共发现与pH性状潜在关联的位点有36个,主要分布在2号和6号染色体上,位于FGD3、AQP1等基因附近。位于2号染色体的1个SNP与肉色L值存在显着关联,与a值存在显着关联的位点有5个,主要分布在1号和4号染色体上。与L*、a*、b*值存在潜在关联的位点分别有124个、83个和13个。
图2 剪切力性状的全基因组关联分析Fig.2 Genome-wide association study of shear force trait
图3 系水力性状的全基因组关联分析Fig.3 Genome-wide association study of water binding capacity trait
图4 pH的全基因组关联分析Fig.4 Genome-wide association study of pH
图5 L*值的全基因组关联分析Fig.5 Genome-wide association study of L* value
图6 a*值的全基因组关联分析Fig.6 Genome-wide association study of a* value
图7 b*值的全基因组关联分析Fig.7 Genome-wide association study of b* value
本试验QQ-plot分析结果见图8,各性状膨胀系数介于1.008～1.037之间,表明本研究试验群体中不存在明显的种群分层现象。
3 讨 论
由于传统的饮食习惯,我国肉鸡生产呈现出独特的产业结构。随着活禽交易的逐步减少,消费者开始接受冰鲜鸡肉产品,使得育种企业和研究人员对兼顾成本和品质的快速型黄羽肉鸡品种的关注日益增加。利用GWAS方法挖掘与重要性状相关的分子标志并用于分子育种,是加快畜禽品种选育进展的有效措施。目前许多研究已经应用GWAS对京海黄鸡[19]、汶上芦花鸡[20]、茶花鸡、狼山鸡、竹丝鸡[21]等品种的体重性状开展了研究。鸡肉品质方面,许多学者已经开展了肌肉糖原与血糖水平[22-23]、肌内脂肪[24]、pH[25]、胸肌脂肪酸组成[26]、胸肌游离氨基酸和核苷酸组成[27]等性状的关联研究,为相关性状的选育提高积累了重要的基础数据。
本研究利用全基因组重测序,对393只立华父系麻黄鸡的60日龄体重和肉品质性状进行了关联分析。结果共发现了10个与目标性状显着关联的SNPs位点,338个潜在关联位点。这些与体重相关的SNPs主要位于4号、5号、19号和28号染色体上,其中以4号染色体上的位点最为丰富。先前的研究也发现,与体重有关的位点集中在1号和4号染色体上,张涛等[28]基于简化基因组测序对京海黄鸡上市体重进行研究时,筛选出了9个位于4号染色体75.641～79.592 Mb区域的关联位点,这说明本研究构建的模型相对准确。
在本研究中,发现一个与体重显着关联的位点位于28号染色体,其临近基因为FSTL3。FSTL3是一种分泌型糖蛋白,属于卵泡抑素模块蛋白家族,与肥胖水平和脂肪细胞炎症发生相关。血清中的FSTL3水平可反映脂肪中FSTL3表达水平及肥胖相关指标,被认为是内脏肥胖的生物标志[29]。此外,FSTL3的过表达减少了肥胖发展过程中的脂肪积累,并可能通过抑制肌肉生长抑制素等因子来控制体重[30]。在潜在关联位点中,rs43415532位于4号染色体上,其临近基因为MBNL3。MBNL3属于肌盲样蛋白家族成员,与肌肉发育和RNA的剪接调控密切相关。研究发现,MBNL3能够通过调控Myod的转录参与肌细胞的分化过程,抑制肌管的形成,并拮抗肌原素和肌球蛋白重链的表达。MBNL3过表达也会导致肌肉发育和收缩相关的基因表达下调[31]。一项研究报道了MBNL3敲除的小鼠出现肌肉再生功能的明显受损[32]。另外,位于34号染色体附近的一个体重相关SNP注释到基因NPFF和PFKM。NPFF是一种由两个氨基酸(Phe-Leu)组成的神经肽,在中枢神经系统中发挥多种作用。研究表明,NPFF可能通过调控相关激素的分泌来控制动物的摄食行为[33]。PFKM是磷酸化果糖激酶,在糖酵解途径中起重要作用,并在肌肉组织中广泛表达。Tan等[34]报道该基因与鸡和猪的肌肉发育有关。综上所述,本研究推测FSTL3、MBNL3、NPFF、PFKM可能通过控制脂肪沉积、肌肉发育和摄食行为参与快速型黄羽肉鸡的生长发育,因此可作为候选基因进行进一步研究。
目前,用于评价鸡肉品质的常见指标包括剪切力、系水力、pH和肉色。肉色通常被认为是评价肉质的首要指标,直接影响消费者的购买选择,一般用L(亮度)、a(红度)和b(黄度)来衡量。本研究共鉴定到220个与肉色相关联的位点,位点数量多于前人报道的结果[35-36],可能是由于本研究选择重测序作为分型工具,可用于关联的位点更多。本研究发现一个与胸肌黄度值相关的位点位于PRKG1基因附近。Zheng等[37]研究发现,PRKG1可以通过调控MYOD、MYOG和MEF2C的表达,从而参与肌肉品质性状的调控。剪切力是另一项重要的肉质指标,可直接反映肉的鲜嫩程度。本研究发现一个与剪切力显着关联的位点位于SLC13A5基因附近,该基因是溶质载体家族成员,在脂肪酸合成中起重要作用[38]。此外,Luo等[39]报道SLC13A5与肌内脂肪含量(IMF)显着相关,而IMF是影响肉类多汁性、嫩度和风味的关键因素之一[40]。系水力即肌肉的保水能力[41],它决定了肉的多汁性、营养组成、肉色外观、酸碱度和其他感官特性[42]。此外,系水力与pH存在联系。本研究发现了一个与系水力相关位点附近的基因为MDFIC。MDFIC属于一个包含富含半胱氨酸的c端结构域的小家族蛋白,研究发现该基因可能通过PI3K/AKT等通路调控肌细胞增殖分化,从而参与肌肉功能的调控[43-44]。本研究还发现AQP1基因旁的一个SNP与鸡肉的pH潜在相关。AQP编码一类水通道蛋白,参与水的跨膜运输。在一项关于猪肉品质的GWAS研究中也报道了AQP1与滴水损失显着关联[45]。另一个与pH潜在关联的位点位于FGD3基因附近。Marchesi等[46]在白条纹鸡肉上的研究中发现,FGD3基因在患病和正常个体的胸肌组织中表达存在显着差异,提示该基因与这种鸡肉问题的发生有关。本试验中,pH性状关联信号较强,与之相关的位点集中于1、2、5、6、12和23号染色体。葡萄糖的代谢途径之一即为合成肝糖原与肌糖原,肌糖原分解为乳酸影响pH,已有研究发现鸡肌肉糖原和血糖水平关联位点位于1、2、3、6和12号染色体[22-23],与本试验结果相似,提示相关位点在肉质调控中发挥关键作用。根据以上结果,认为PRKG1、SLC13A5、MDFIC、AQP1、FGD3等基因可能是快速型黄羽肉鸡肉品质性状的候选基因。然而,具体的作用机制仍需进一步进行试验验证。
4 结 论
本研究基于立华麻黄鸡全基因组重测序数据进行了GWAS分析,在与体重性状关联的结果中,筛选到2个显着关联,52个潜在关联的位点;肉品质性状结果中,筛选到8个显着关联,296个潜在关联的位点;并注释到99个与体重、225个与肉品质相关的功能基因。其中FSTL3、MBNL3、NPFF、PFKM、PRKG1、SLC13A5、MDFIC、FGD3、AQP1等基因可作为肉鸡体重和肉质性状的关键候选基因。本研究为快速型黄羽肉鸡重要经济性状提供了重要的遗传变异位点和候选基因,可以为相关品种开展分子育种提供理论基础。
