全基因组关联分析黄羽肉鸡专门化品系体斜长性状相关的分子标记

Identification of molecular markers associated with body oblique length through genome-wide association study on a specialized strain of yellow feathered broilers

摘要

为找到适应于黄羽肉鸡专门化品系体斜长的SNP分子标记，加速体斜长世代选育进展，本研究对快速型黄羽肉鸡（60日龄）的终端父系400只种公鸡进行全基因组重测序，最终366只鸡被纳入GWAS分析挖掘快速型黄羽肉鸡专门化品系体斜长发育相关的候选基因和分子标记。结果表明：1）11个SNPs与体斜长潜在显著相关，位于8个注释基因，包括VPS36、COG6、LHFPL1、KIAA1210、DARS、THSD7B、TSHZ3和MVP。2）SNP3、SNP4、SNP10、SNP11 4个位点偏离H-W平衡状态（P<0.05），其他7个SNPs位点均处于H-W平衡状态（P>0.05）。3）单标记与体斜长性状关联分析表明，SNP1（rs312577341）AA、SNP2（rs313388213）AA和AG、SNP3（rs730967108）GG、SNP4（rs315520013）AA和AG、SNP5（rs317946632）CC、SNP6（rs315619998）CC、SNP7（rs317342517）CC、SNP8（rs313349990）CC、SNP9（rs731600098）AA和AG、SNP10（rs13792170）CC和CT、SNP11（rs7076259500）GG和AG分别为较大体斜长表型相关的基因型。4）SNP5和SNP6之间的强连锁产生了3种单倍型：H1（TG）、H2（CC）和H3（CG），其中，H2H2单倍型（CCCC）为较大体斜长表型相关的基因型。综上，筛选到与快速型黄羽肉鸡专门化品系体斜长相关的11个SNPs分子标记和8个重要注释基因，这些发现为进一步理解体斜长发育的潜在机制和更有效地培育体型更大的快速型黄羽肉鸡新品种（配套系）提供科学理论支撑。

Abstract

To identify SNP molecular markers associated with the body oblique length of specialized yellow-feathered broiler lines and accelerate generational selection progress， herein， a total of 400 terminal male chickens of the fast-type yellow feathered broiler （60 day old） were selected as the research subjects， and their genomic genetic variation information was obtained using the whole genome resequencing and 366 chickens were ultimately included in genome-wide association study （GWAS） to identify candidate genes and molecular markers associated with body length development. The results were as follows： 1） Eleven SNPs were potentially significantly correlated with body oblique length， located at 8 annotated genes， including VPS36，COG6，LHFPL1，KIAA1210，DARS，THSD7B，TSHZ3，and MVP. 2） The four loci of SNP3， SNP4，SNP10， and SNP11 deviated from the Hardy-Weinberg equilibrium state （P<0.05）， and the other 7 SNPs loci were all in the Hardy-Weinberg equilibrium state （P>0.05）. 3） The association analysis between single SNP marker and body oblique length traits showed that AA at SNP1 （rs312577341）， AA and AG at SNP2 （rs313388213）， GG at SNP3 （rs730967108）， AA and AG at SNP4 （rs315520013）， CC at SNP5 （rs317946632）， CC at SNP6 （rs315619998）， CC at SNP7 （rs317342517）， CC at SNP8 （rs313349990）， AA and AG at SNP9 （rs731600098）， CC and CT at SNP10 （rs13792170） ， and GG and AG at SNP11 （rs7076259500） were all the genotypes associated with a larger body oblique length. 4） The strong linkage between SNP5 and SNP6 resulted in three haplotypes H1 （TG）， H2 （CC） and H3 （CG）. Among them， the H2H2 haplotype （CCCC） was a genotype associated with a larger body oblique length. In summary， 11 SNPs molecular markers and 8 important annotation genes related to the body oblique length in specialized strains of fast-type yellow feathered broilers were screened. These findings provide important references for further studying the regulatory mechanisms of body oblique length development and cultivating new high-quality broiler breeds （crossed lines） with larger body sizes.

Graphical abstract

关键词

全基因组关联分析 / 鸡 / 体斜长 / 肉鸡育种 / 体尺

Key words

GWAS / chicken / body oblique length / broiler breeding / body size

鸡品种受到长期的自然选择或/和人工选择，因而表现出很大的体型差异^［1］。鸡的体型大小反映了其身体发育情况，尤其是骨骼和肌肉的生长情况^［2］。体尺性状是直接反映体型大小和生长发育程度的度量指标，与生产性能密切相关，对畜牧业的经济效益有重大影响。体斜长、龙骨长、胸宽、胸深等体尺性状是反应体型的重要指标，其中体斜长作为一个重要的体尺性状，可以反映鸡骨骼的生长发育状况^［3］。体斜长与体重、屠体重、半净膛重、全净膛重等经济性状显著相关，可以间接用于体重、屠宰性能等性状的选择，也是鸡早期生长性状选择的重要候选性状之一^［4-5］。黄羽肉鸡配套系选育过程中都要对其专门化品系（特别是终端父系）进行体重和体型的选育，以生产生长速度更快、体型较大的商品代鸡。目前黄羽肉鸡育种过程中主要依靠简单测量体斜长性状进行体型选育，遗传进展较慢，筛选鉴定体斜长发育相关的候选基因和分子标记应用于体型选育，以更快提高体斜长选育进展具有十分重要的意义。随着高通量测序技术的出现，能更准确地检测出体斜长性状基因组中的重要候选基因。Wang等^［6］通过重测序和全基因组关联研究（GWAS）发现SLIT2和LCORL为18周龄汶上芦花鸡体斜长重要候选基因，Li等^［7］特利用异性位点扩增片段测序技术和GWAS筛选到ZCCHC7、PAX5和MELK为汶上芦花鸡56日龄体斜长的重要候选基因。骆娜等^［3］通过“京芯一号”55 K芯片和GWAS筛选到KCNIP4、LDB2、TAPT1、NRXN3 可以作为113日龄文昌鸡体斜长性状潜在的候选基因。综合多项研究表明，不同品种和不同鸡群间体斜长性状候选基因不同，表明体斜长性状是由多个基因和信号通路控制， 其遗传机制在不同品种和不同群体中差异明显。

黄羽肉鸡（地方鸡种以及由地方鸡种杂交的黄羽肉鸡配套系）抗病力强、肉品质优良、遗传和表型多样性丰富，具有极高的经济价值和文化价值。关于黄羽肉鸡专门化品系体斜长性状相关的候选基因和分子标记研究相对较少，了解鸡体斜长变异背后的遗传机制仍然不足。本研究以快速型黄羽肉鸡（60日龄上市）的终端父系400只种公鸡为研究对象，利用重测序获得其基因组遗传变异信息，通过GWAS挖掘适应于快速型黄羽肉鸡专门化品系体斜长发育相关的候选基因和分子标记，以更快地提高鸡体斜长的遗传进展，为更有效地培育体型更大、体型均匀度更好的黄羽肉鸡新品种（配套系）提供科学理论支撑。

1　材料与方法

1.1　实验动物

400只实验鸡来自江苏立华食品集团有限公司黄羽肉鸡终端父系A系第12世代种公鸡，参与配套形成的商品鸡为快速型肉鸡品种，60 d上市体重2.5 kg以上。实验鸡群饲养于江苏立华育种有限公司，所有试验鸡均为同一批次孵化并同批饲养于同一鸡舍内。采用公司标准的种鸡饲粮进行饲喂，免疫遵循公司制定的标准程序。试验期鸡群自由饮水。

1.2　体斜长的测定

实验鸡饲养至60日龄时，禁食12 h后进行称重，记录脚号。体斜长测定方法按照农业行业标准NY/T 823—2020 《家禽生产性能名词术语和度量计算方法》执行^［8］，具体测定方法如下：用皮尺沿体表测量肩关节至坐骨结节间的距离即为鸡体斜长，所有鸡体斜长的测定均由同一人测量，所测部位一致。

1.3　全基因组重测序

60日龄时，对所有测定过体斜长性状的群体进行翅下静脉采血1.5 mL，EDTA抗凝。采用常规苯酚氯仿方法提取DNA^［9］。

全基因组重测序：通过质检的DNA样品，统一交至深圳华大进行全基因组重测序，包括制备基因组文库、MGISEQ-2000平台测序，根据标准BGISEQ程序进行全基因组测序后，366个个体产生了4.25 Tb的原始数据，平均深度为10.17×。过滤过程产生了4.24 Tb的高质量测序数据，与参考鸡基因组GRcg7b比对时，平均作图率为98.73%。GC的平均含量为45.65%。

质量控制：对测序的基因组原始reads进行过滤，使用fastp（v.0.20.0）预处理器去除低质量reads、接头以及包含poly-N的reads^［10］。检查污染情况时，随机选择10万个读数，并与NT（核苷酸序列数据库）库中的序列进行比对，然后用BWA（v0.7.17）软件将质控后得到的reads比对至参考基因组GRcg7b^［11］。用Picard3.1和SAMtools 1.18进行格式转换和序列排序，将重复比对的reads移除。通过GATK 4.0软件检测SNPs，获得高质量的SNPs位点，共检测到16 411 878个SNPs位点^［12］。使用PLINK（v1.9）软件对群体遗传变异的SNPs进行质控^［13］，确定以下位点和个体质控标准：SNPs缺失频率<2%、次等位基因频率maf>0.05、哈迪温伯格P<1×10^-4，同时剔除基因型缺失率高于4%的个体。最终获得366个个体和10 266 594个SNPs位点用于关联分析。

1.4　全基因组关联研究

通过GEMMA（v0.985）软件的线性混合模型（Linear mixed model，LMM）进行体斜长性状的GWAS分析^［14］。所使用的线性混合模型如下：

y＝Wα+xβ+u+ϵu~MVNn（0，λτ^-1K），

ϵ~MVNn（0，λτ^-1In）

式中：y是表型向量；W是固定效应矩阵（PCA在内的前3个主成分）；α是固定效应向量；x为基因型向量，β为SNP的效应向量；u为随机效应向量，u是一个具有协方差结构的随机效应向量，指的是n维多元正态分布u∼MVNn（0，λτ^-1K）。MVNn为n维多元正态分布，λ为遗传方差与残差方差（随机效应向量u与残差向量ϵ）的比值，τ^-1为残差方差。K为使用GEMMA基于SNP计算的亲缘关系矩阵，In为单位矩阵。ϵ代表残差向量。

使用PLINK软件通过将参数设置为50个SNPs窗口大小、11个SNPs步长和0.2的r²阈值进行连锁不平衡（LD）修剪，得到432 529个有效的独立检验SNP标记。多重检验后，采用5%的Bonferroni校正法来设定显著阈值，本研究中全基因组水平显著阈值0.05/N=0.05/432 529≈ 1.16×10⁻⁷（N为经LD修剪后的有效独立SNP标记数）［-lgP］=6.94；潜在显著阈值参考Lander等^［15］提出的suggestive significance概念，即在一次全基因组扫描中预期产生一个假阳性结果的P值水平，在Bonferroni校正框架下近似计算为1/N= 1/432 529≈2.31×10⁻⁶［-lgP］=5.64。R语言环境下CMplot软件包（https：//github.com/YinLiLin/CMplot）绘制曼哈顿图和Q-Qplot^［16］。通过Bedtools（v2.30.0）软件对显著位点上下游100 kb内的临近基因进行注释^［13］。质检后，对366只鸡体斜长表型值进行全基因组关联分析。

1.5　数据统计分析

使用SPSS 18.0软件中的一般线性模型中的单变量方差分析进行基因型与体斜长关联分析，找出优势基因型，结果用平均值±标准差表示。SNP标记的不同基因型为固定因子，体斜长为因变量，采用Tukey HSD法比较不同标记基因型之间的体斜长的差异显著性，P<0.05表示差异显著。利用Haploview4.1对SNP位点进行连锁不平衡分析，统计各个SNP位点基因型频率、基因频率、期望杂合度（He）和观察杂合度（Ho）。

2　结果与分析

2.1　体斜长性状的测定

60日龄366只公鸡实验鸡平均体斜长197.26 mm，最大值和最小值分别为220.91和163.87 mm，体斜长变异系数为5.05%。

2.2　体斜长GWAS分析

在过滤SNP后，对10 266 594个SNPs进行了GWAS，以鉴定与体斜长相关的基因组区域。体斜长的曼哈顿图和Q-Q图见图1和图2。由图1和表1可知，在1、4、7、11、28号染色体鉴定到与体斜长潜在显著相关的11个SNPs（P<2.31×10^-6），这些位点分别位于VPS36、COG6、LHFPL1、KIAA1210、DARS、THSD7B、TSHZ3和MVP基因的内含子、上游或下游。Q-Q结果如图2所示，体斜长的膨胀系数（λ）为0.958（λ≈1），表明模型选择合理，筛选到与目标性状关联性较强位点，分析结果可靠。在GWAS分析中，QQ图 的λ是衡量群体分层或表型分布异常的关键指标，用于检测数据是否存在系统性偏差^［17］。λ>1.2时，QQ图可能出现系统性偏移，需排查样本是否包含不同群体，或重新校正群体分层^［18-19］。本试验模型中λ≈1，说明，QQ图进一步验证GWAS结果可靠。

2.3　SNPs位点遗传多态性分析

由表2可知，11个SNPs位点中除了SNP8具有2种基因型，其他10个SNPs位点均具有3种基因型。经H-W平衡检测，11个SNPs位点，SNP3、SNP4、SNP10和SNP11 4个位点偏离H-W平衡状态（P<0.05），其他7个SNPs位点均处于H-W平衡状态（P>0.05）。3个SNPs位点（SNP8、SNP9和SNP10）遗传多样性较低（He<0.25），其他8个SNPs位点均为中度遗传多样性（0.25<He<0.50）^［20］。

2.4　SNPs位点与体斜长相关性分析

首先分析单个SNP标记与体斜长的相关性，11个SNPs位点不同基因型间体斜长均差异显著，如表3所示。SNP1（rs312577341）AA基因型个体的体斜长显著大于AG和GG基因型（P<0.05）；SNP2（rs313388213）AA和AG基因型个体的体斜长显著大于GG基因型（P<0.05）；SNP3（rs730967108）GG基因型个体的体斜长显著大于AA和AG基因型（P<0.05）；SNP4（rs315520013）AA和AG基因型个体的体斜长显著大于GG基因型（P<0.05）；SNP5（rs317946632）CC基因型个体的体斜长显著大于CT和TT基因型（P<0.05）；SNP6（rs315619998）CC基因型个体的体斜长显著大于CG和GG基因型（P<0.05）；SNP7（rs317342517）CC基因型个体的体斜长显著大于CT和TT基因型（P<0.05）；SNP8（rs313349990）具有2种基因型CC和CT，其中CC基因型个体的体斜长显著大于CT基因型（P<0.05）；SNP9（rs731600098）AA和AG基因型个体的体斜长显著大于GG基因型（P<0.05）；SNP10（rs13792170）CC和CT基因型个体的体斜长显著大于TT基因型（P<0.05）；SNP11（rs7076259500）GG和AG基因型个体的体斜长显著大于AA基因型（P<0.05）。

与体斜长潜在显著相关的11个SNPs位点进行连锁不平衡分析（LD）结果分别见图3，方框中的数值为D'值乘以100后得到。由图3可知，SNP5和SNP6位点处于强连锁状态，SNP5和SNP6之间的连锁产生了3种单倍型：H1（TG）、H2（CC）和H3（CG），单倍型频率分别为0.795、0.177和0.019（图3）。其他 9个SNPs位点不处于强连锁状态，所以只分析SNP5和SNP6的合并基因型与体斜长的相关性，不再分析其他9个SNPs分子标记的合并基因型与体斜长的相关性。SNP5和SNP6连锁产生的3种单倍型组合后，产生了5种基因型（表4）。由表4可知，H2H2单倍型（CCCC基因型）体斜长显著高于其他单倍型组合（P<0.05）。

2.5　GEMMA线性混合模型潜在显著相关的11个SNP位点效应值分析

GEMMA线性混合模型潜在显著相关的11个SNP位点效应值参数见表5。由表5可知，11个位点中，SNP1~SNP7和SNP11的β值为正，表示其minor allele（allele1）与体斜长增加相关；SNP8~SNP10的β值为负，表示其minor allele与体斜长降低相关。所有位点的β方向与表3中各基因型均值排序完全一致，验证了GWAS结果与单标记关联分析的内在逻辑一致性。所有位点的95%置信区间均不包含0，效应方向估计具有统计稳健性。

3　讨 论

本研究通过重测序和GWAS首次挖掘到与黄羽肉鸡专门化品系体斜长潜在显著相关的11个SNPs位点和一些重要注释基因，与前人挖掘到汶上芦花鸡和文昌鸡体斜长SNP分子标记和候选基因不同^{［3，6-7］}，这可能由于体斜长是一个复杂的数量性状，由多个基因和信号通路决定，在不同品种不同日龄鸡群中遗传机制不同。筛选到与体斜长潜在显著相关的11个SNPs位点，有7个位于非编码区内含子，尽管编码区序列SNPs可能影响蛋白质的氨基酸序列，但位于非编码区SNPs会影响RNA加工过程，进而可能改变蛋白质的序列或产量^［21］，从而影响相关表型性状。

H-W平衡假设的前提条件包括群体足够大、随机交配、没有突变和迁移、没有基因型选择优势等^［22］。人工配种、群体数量不够大等因素会导致一些基因位点偏离H-W平衡状态^［23］。本研究实验群体为快速型黄羽肉鸡终端父系A系第12世代种公鸡，经历了长期高强度的人工选择，由于采用人工配种而非随机交配，有效群体含量相对有限，经历了12个世代的持续定向选择（主要选育目标为体重、体型等生长性能），这些因素均可使一些基因位点偏离H-W平衡状态，这在专门化品系中是常见的现象，也是选育成效的反映。另外，为排除基因分型技术偏差，对这4个偏离H-W平衡的SNP位点进行了分型质量重检，各等位基因reads覆盖深度均匀，未发现特定等位基因在低深度样本中系统性丢失的现象；检查了B-allele frequency分布，未发现异常的模糊分型模式，低质量变异位点已被排除，综合以上检查，这4个位点的H-W偏离更可能源于生物学因素（人工选择效应）而非技术误差。GEMMA模型通过亲缘关系矩阵K校正了群体结构和亲缘关系，在一定程度上可以缓解H-W偏离对关联分析的影响。此外，本研究的PLINK质控中H-W平衡过滤阈值设定为P<1×10⁻⁴（仅剔除极度偏离的位点），这4个位点的偏离水平（P值为0.007~0.047）远高于该阈值，因此被保留在分析中，这也是畜禽GWAS中的常规做法。考虑到SNP3、SNP4、SNP10和SNP11偏离H-W平衡，其关联信号和效应估计可能受到一定影响，需在独立群体中进一步验证。

本研究筛选到的与体斜长发育潜在相关的重要候选基因主要集中参与骨骼发育、钙信号稳态、免疫和细胞生长、增殖、存活和细胞骨架结构调控等各种生物过程。如血小板反应蛋白1型结构域7B（Thrombospondin type 1 domain containing protein 7B， THSD7B）编码含凝血酶敏感蛋白1型结构域的蛋白7B，在肌动蛋白细胞骨架重组中起关键作用，运动小鼠的THSD7B基因表达水平显著上调^［24］，运动会激活AMP活化蛋白激酶（AMPK）和过氧化物酶体增殖物激活受体γ共激活因子1-α（PGC-1α），它们是肿瘤细胞表达的核心，这些分子可能会影响控制THSD7B表达的转录因子和辅助调节因子，THSD7B基因突变可能抑制与细胞死亡相关的信号通路，同时增强与侵袭、转移和下调免疫反应相关的通路^［24-26］。由于鸡在平养过程中经常运动，THSD7B在鸡生长发育过程中可能发挥着重要的作用。THSD7B参与成骨和软骨发育过程，是肉鸡股骨头分离症状发生发展的重要候选基因^［27］。THSD7B基因与牦牛体长、体高、胸围和体重呈极显著相关^［28］，本研究中THSD7B基因SNP9（rs731600098，A→G）位点与体斜长显著相关，具有AA和AG基因型个体的体斜长显著大于GG基因型个体，由于GG基因型样本量只有3个，无法进行有效的正态性检验，该比较结果仅作为参考，后续将扩大群体进一步验证。

空泡蛋白分选蛋白家族成员36（Vacuolar protein sorting 36，VPS36），是一种泛素结合蛋白，主要参与内体分选转运复合体（ESCRT）途径，在泛素化内吞受体进入溶酶体降解的通路中发挥重要作用^［29］，对于形成多泡体和空泡生物发生泛素化质膜蛋白的降解至关重要^［30-31］。也有研究表明，VPS36在鸡卵巢卵泡发育过程中发挥重要作用，VPS36 mRNA随雌激素浓度升高而逐渐下降，并受鸡卵泡中的促卵泡激素（FSH）和雌二醇调控^［32-33］。Shen等^［34］通过600芯片和GWAS发现VPS36为鸡冠发育的重要候选基因。本研究发现VPS36基因上的SNP1（rs312577341， G→A）突变与体斜长显著相关，且突变后体斜长性状显著增加，表明VPS36可能在鸡体斜长发育中发挥重要作用。

KIAA1210 是日本研究所“Kazusa Institute of Advanced Biosciences”对大量人类 cDNA 进行测序和功能注释研究发现，并以项目名称KIAA命名系列基因中一个基因。KIAA1210蛋白是一种必需的细胞连接蛋白，可能涉及细胞骨架结构的调控，主要在睾丸组织中表达水平较高，参与人精子发生过程^［35-37］，KIAA1210基因的突变或异常表达会导致精子畸形^［38］。本研究首次发现KIAA1210基因与鸡体斜长相关，KIAA1210基因3个SNPs（SNP5 rs317946632，T→C；SNP6 rs315619998，G→C；SNP7 rs317342517，T→C）位点，SNP5和SNP6强连锁后CCCC基因型和SNP7位点CC基因型均为较大体斜长表型相关的基因型，说明该基因的3个SNPs突变均为有益突变。SNP5、SNP6和SNP7均位于KIAA1210基因内含子区域，间距仅约30 kb。图4已显示SNP5和SNP6处于强连锁状态（D'>0.9）。从效应值来看，SNP5（β=+5.04）>SNP6（β=+4.74）>SNP7（β=+4.40），三者高度接近且方向一致，说明3个SNP的关联信号很可能来源于同一个功能性单倍型或因果变异，三者作为标记SNP共同标记了该区域的同一QTL信号。本研究的分辨率无法区分哪个SNP更接近因果变异，精细定位需要更高密度的分型数据或功能实验，后续将进一步研究。

GEMMA模型本身采用的是加性效应模型，表5中β系数即为每增加一个allele1拷贝的加性效应估计。以效应最大的SNP5为例进行验证：加性模型预测CC-TT差值应为2×5.04=10.08 mm，而表3实测CC-TT=210.22-195.61=14.61 mm，实测差值大于加性预期，提示该位点可能存在一定的显性效应（CC纯合子的增益超过加性预期）。但考虑到CC基因型仅n=12，均值估计精度有限，该超额可能部分来源于抽样波动。对于其他位点（如SNP1），加性模型预测AA-GG差值为2×3.67=7.34 mm，表3实测AA-GG=201.96-194.57=7.39 mm，两者高度吻合，支持加性效应模型在该位点的适用性。

保守的寡聚高尔基复合体亚基6（Conserved oligomeric Golgi complex subunit 6， COG6），在高尔基体内的囊泡逆向运输中发挥着关键作用，COG6基因突变影响了高尔基体形态以及蛋白的定位，在拟南芥花粉管快速生长过程中发挥着十分关键的作用^［39］，COG6是对甲型流感病毒（IAV）复制至关重要的新型宿主因子^［40］，COG6基因变异会引起先天性糖基化障碍，主要表现为生长迟缓、发育迟缓、小头畸形、肝脏和胃肠道疾病、关节挛缩等^［41-42］。COG6基因SNP2（rs313388213，A→G）位点与体斜长潜在显著相关，具体作用机制需要进一步验证。

天冬氨酰-tRNA合成酶（Aspartyl-tRNA synthetase，DARS）是氨酰-tRNA合成酶（Aminoacyl-tRNA Synthetases，ARS）家族中的重要成员之一，将天冬氨酸与其相应的tRNA配对，编码的胞质天冬氨酰-tRNA合成酶是所有细胞中蛋白质合成所必需的一种蛋白质。DARS基因突变可以导致人髓鞘形成不足，脑干和脊髓受累，腿部痉挛^［43-44］。DARS参与调节炎症、血管生成和细胞凋亡等，其突变可能接受错误的tRNA或氨基酸，从而导致翻译错误，进而影响ERK通路上游的钙信号稳态^［45-46］。DARS基因 SNP8（rs313349990，C→T）位点CC基因型为较大体斜长表型相关的基因型，其具体作用机制需要进一步研究。

脂肪瘤HMGIC融合伴侣样1（Lipoma HMGIC fusion partner-like 1，LHFPL1）在人骨骼肌、肺、胸腺、结肠和卵巢中表达量较高，在六种肝肿瘤细胞系中均存在转录^［47］，在家禽方面未见该基因的研究报道。本研究发现SNP3（rs730967108，G→A）和SNP4（rs315520013，G→A）均位于LHFPL1基因下游，SNP3位点具有GG基因型的个体体斜长显著高于具有GA和AA基因型的个体，而SNP4位点具有AA和GA个体的体斜长显著高于具有GG基因型的个体，并且这2个SNP位点不处于强连锁状态，说明这2个SNP位点可能对体斜长的作用不同。

锌指同源盒3（Teashirt zinc finger homeobox 3，TSHZ3）基因编码一种锌指转录因子，是鸡苗勒管形成的重要候选调节因子^［48］，TSHZ3蛋白可能与鸡雌性胚胎右侧性腺的退化有关^［49］。另外，TSHZ3基因表达异常参与多种肿瘤发生、发展过程。主要拱顶蛋白（Major vault protein， MVP）是真核细胞中最大的核糖核蛋白穹窿体的主要组成部分^［50］。MVP是磷酸化蛋白，可以被蛋白激酶C、酪蛋白激酶和Src激酶磷酸化，是细胞信号通路的重要调控者，MVP与参与不同细胞信号通路的蛋白酪氨酸激酶相互作用^［51］，MVP 参与MAPK通路^［52］。SNP10（rs13792170，C→T）和SNP11（rs7076259500，A→G）分别位于TSHZ3和MVP基因上游，SNP10位点CC和CT为较大体斜长表型相关的基因型，SNP11位点GG和AG为体斜长的优势基因型，这2个位点可能受到TSHZ3和MVP基因调控，进而影响体斜长的发育，其具体作用机制有待于进一步研究。

体斜长等体型发育是由骨骼发育决定，骨骼除了提供生物力学支持和调节矿物质平衡的典型作用外，还可作为一个内分泌器官发挥着额外的生物学功能，骨骼对全身能量代谢的变化做出反应，例如骨钙素在骨骼重塑活动和全身能量需求之间提供了联系，骨髓内分泌蛋白在组织生长、免疫调节和能量调节方面发挥重要作用^［53-55］。本研究筛选到一些SNP位点的注释基因具有免疫调节、能量代谢、组织生长等多种生物学效应。上述基因的功能角色目前仅基于物理位置关系和其他物种/组织的文献推测，本研究未进行eQTL分析、转录组验证或基因功能试验等直接功能验证，这些基因与体斜长发育的因果关系尚需通过基因表达分析、基因敲除/敲降等试验进一步确认。

本研究存在以下局限性：①样本量为366只，属于畜禽GWAS中的中等规模，对中等效应QTL位点的检测效力有限，11个SNP位点经Bonferroni校正后均达到潜在显著阈值但未达到全基因组严格显著阈值，其关联信号的真实性需在更大样本量的群体中验证；②研究对象仅为单一专门化品系的种公鸡，筛选到的关联位点在其他黄羽肉鸡品系、其他日龄阶段和母鸡群体中的有效性尚不确定；③缺乏独立群体的外部验证，筛选到的关联位点和优势基因型仅在发现群体内进行了初步验证；④未进行eQTL分析或转录组验证，注释基因的功能推断基于物理位置关系和文献推测，其与体斜长发育的因果关系尚待确认；⑤未评估这些SNP位点与体重、料肉比、腿病发生率等其他重要经济性状的遗传关联，无法排除拮抗多效性的风险。后续研究将扩大样本规模，在独立群体中验证上述关联位点的效应和注释基因的生物学作用。

4　结 论

本研究筛选到与体斜长潜在显著相关的11个SNPs分子标记和8个注释基因，并分析影响较大体斜长表型相关的基因型，THSD7B、VPS36、COG6和KIAA1210基因可能在体斜长发育中发挥重要作用，这些发现为进一步研究体斜长发育的调控机制和体斜长分子标记辅助育种提供参考。

参考文献

[1] Lyu S J， Arends D， Nassar M K， Weigend A， Weigend S， Wang E Y， Brockmann G A. High-density genotyping reveals candidate genomic regions for chicken body size in breeds of Asian origin［J］. Poultry Science， 2023， 102（1）：102303

[2] Geng A L， Zhang Y， Zhang J， Zeng L C， Chang C， Wang H H， Yan Z X， Chu Q， Liu H G. Effects of light regime on the hatching performance， body development and serum biochemical indexes in Beijing You Chicken［J］. Poultry Science， 2021， 100（8）：101270

[3] 骆娜，安炳星，魏立民，文杰，赵桂苹.全基因组关联分析筛选文昌鸡体尺性状相关分子标记，中国农业科学 2024，57（10）：2046-2060

[4] Luo N， An B X， Wei L M， Wen J， Zhao G P. Identification of molecular markers associated with body size traits through genome-wide association analysis in Wenchang chickens［J］. Scientia Agricultura Sinica， 2024， 57（10）：2046 - 2060 （in Chinese）

[5] 巨玉鑫，史兆国，郭志明，徐巍芹.康县太平鸡体尺性状与屠宰性状的测定及相关性分析［J］.畜牧兽医杂志，2024，43（1）：7-9，12

[6] Ju Y X， Shi Z G， Guo Z M， Xu W Q. Determination and correlation analysis of body size traits and slaughter traits of taiping chicken in Kang County［J］. Journal of Animal Science and Veterinary Medicine， 2024， 43（1）：7-9，12 （in Chinese）

[7] 蒋发俊，许文坤，范新阳，苗永旺. 武定鸡、大围山微型鸡及尼西鸡体尺对体重影响的比较研究［J］.中国畜牧杂志，2025，61（2）：129-137

[8] Jiang F J， Xu W K， Fan X Y， Miao Y W. Comparative study on the influence of body size on body weight of Wuding chicken， Daweishan miniature chicken and Nixi chicken［J］. Chinese Journal of Animal Science， 2025， 61（2）： 129-137 （in Chinese）

[9] Wang J， Liu J， Lei Q X， Liu Z H， Han H X， Zhang S E， Qi C， Liu W， Li D P， Li F W， Cao D G， Zhou Y. Elucidation of the genetic determination of body weight and size in Chinese local chicken breeds by large-scale genomic analyses［J］. BMC Genomics， 2024， 25（1）：296

[10] Li F W， Liu J， Liu W， Gao J B， Lei Q X， Han H X， Yang J C， Li H M， Cao D G， Zhou Y. Genome-wide association study of body size traits in Wenshang Barred chickens based on the specific-locus amplified fragment sequencing technology ［J］. Animal Science Journal， 2021， 92（1）：e13506

[11] NY/T 823—2004家禽生产性能名词术语和度量统计方法 . ［S］. 北京：中国农业出版社，2004

[12] NY/T 823—2004. Performance ferms and measurement for poultry ［S］. Beijing： China Agriculture Press， 2004 （in Chinese）

[13] 萨姆布鲁克，弗里奇， 曼尼阿蒂斯著，金冬雁，黎孟枫 译.分子克隆实验指南［M］.第2版. 北京： 科学出版社， 1992： 477-479

[14] Sambrook J， Fritsch E F， Manniadis T. Translated by Jin D Y， Li M F， Hou Y D .Molecular Cloning［M］. 2nd ed.Beijing： Science Press， 1992： 477-479 （in Chinese）

[15] Bush S J. Read trimming has minimal effect on bacterial SNP-calling accuracy［J］. Microbial Genomics， 2020， 6（12）：mgen000434

[16] Li H， Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform［J］. Bioinformatics， 2009， 25（14）：1754-1760

[17] Niaré K， Greenhouse B， Bailey J A. An optimized GATK4 pipeline for Plasmodium falciparum whole genome sequencing variant calling and analysis［J］. Malaria Journal，2023， 22（1）： 207

[18] Purcell S， Neale B， Todd-Brown K， Thomas L， Ferreira M A R， Bender D， Maller J， Sklar P， de Bakker P I W， Daly M J， Sham P C. PLINK： a tool set for whole-genome association and population-based linkage analyses ［J］. The American Journal of Human Genetics， 2007， 81（3）： 559-575

[19] Fan C Y， Shan Y J， Zhang M， Tu Y J， Ji G G， Ju X J， Shu J T， Liu Y F. Genome-wide association study of myofiber type composition traits in a yellow-feather broiler population ［J］. Poultry Science， 2025， 104（1）： 104634

[20] Lander E， Kruglyak L. Genetic dissection of complex traits： Guidelines for interpreting and reporting linkage results［J］. Nature Genetics， 1995， 11（3）： 241-247

[21] Yin L L， Zhang H H， Tang Z S， Xu J Y， Yin D， Zhang Z W， Yuan X H， Zhu M J， Zhao S H， Li X Y， Liu X L. rMVP： A memory-efficient， visualization-enhanced， and parallel-accelerated tool for genome-wide association study ［J］. Genomics， Proteomics & Bioinformatics， 2021， 19（4）： 619-628

[22] Quinlan A R. BEDTools： The Swiss-army tool for genome feature analysis ［J］. Current Protocols in Bioinformatics， 2014， 47： 11.12.1-11.1234

[23] Guintivano J， Shabalin A A， Chan R F， Rubinow D R， Sullivan P F， Meltzer-Brody S， Aberg K A， van den Oord E J C G. Test-statistic inflation in methylome-wide association studies ［J］. Epigenetics， 2020， 15（11）： 1163-1166

[24] Georgiopoulos G， Evangelou E. Power considerations for λ inflation factor in meta-analyses of genome-wide association studies . Genetics Research， 2016， 98： e9

[25] Zanetti E， De Marchi M， Abbadi M， Cassandro M. Variation of genetic diversity over time in local Italian chicken breeds undergoing in situ conservation［J］. Poultry Science， 2011， 90（10）： 2195-2201

[26] Rohmah L， Darwati S， Ulupi N， Khaerunnisa I， Sumantri C. Polymorphism of prolactin （PRL） gene exon 5 and its association with egg production in IPB-D1 chickens ［J］. Archives Animal Breeding， 2022， 65（4）： 449-455

[27] Guang Y Z， Donner A. The merits of testing hardy-Weinberg equilibrium in the analysis of unmatched case-control data： A cautionary note ［J］. Annals of Human Genetics， 2006， 70（6）： 923-933

[28] 李青，卢曾奎，金美林，权凯，储明星，魏彩虹. 绵羊BMP2基因型与尾性状的关联分析［J］. 农业生物技术学报，2019，27（11）：1985-1995

[29] Li Q， Lu Z K， Jin M L， Quan K， Chu M X， Wei C H. Association analysis of the BMP2 genotype with sheep （Ovis aries） tail traits ［J］. Journal of Agricultural Biotechnology， 2019， 27（11）： 1985-1995 （in Chinese）

[30] Luo Z W， Zhu J G， Fang Z Y， Xu R， Wan R W， He Y W， Chen Y S， Chen S， Wang Q， Liu Q Z， Chen S Y. Exercise-augmented THSD7B exhibited a positive prognostic implication and tumor-suppressed functionality in pan-cancer ［J］. Frontiers in Immunology， 2024，15： 1440226

[31] Zhang Q Q， Zhao Y W， Song Z N， Zhang Q， Tian C H， Li R R， Zheng J， Yan L L， Gu M L， Jia X D， Li M J. Identification of THSD7B and PRMT9 mutations as risk factors for familial lung adenocarcinoma： A case report［J］. Medicine， 2023， 102（6）： e32872

[32] Yao Z， Lin A， Yi Y， Shen W， Zhang J， Luo P. THSD7B mutation induces platinum resistance in small cell lung cancer patients［J］. Drug Design， Development and Therapy， 2022， 16： 1679-1695

[33] Ibelli A M G， de Oliveira Peixoto J， Zanella R， de Simoni Gouveia J J， Cantão M E， Coutinho L L， Marchesi J A P， Pizzol M S D， Marcelino D E P， Ledur M C. Downregulation of growth plate genes involved with the onset of femoral head separation in young broilers［J］. Frontiers in Physiology， 2022， 13： 941134

[34] 任稳稳，马晓明，余道宁，欧杰次仁，喇永富，郭宪，褚敏，阎萍，梁春年. 牦牛THSD7B和COL26A1基因多态性与生长性状的关联分析［J］.中国草食动物科学， 2025，45（6）：105-110

[35] Ren W W， Ma X M， Yu D N， Ou J， La Y F， Guo X， Chu M， Yan P， Liang C N. Association analysis of THSD7B and COL26A1 genes polymorphisms with growth traits in yaks［J］. China Herbivore Science， 2025， 45（6）： 105-110 （in Chinese）

[36] Raiborg C， Stenmark H. The ESCRT machinery in endosomal sorting of ubiquitylated membrane proteins［J］. Nature， 2009， 458（7237）： 445-452

[37] Lin S J， Lin M C， Liu T J， Tsai Y T， Tsai M T， Lee F S. Endosomal Arl4A attenuates EGFR degradation by binding to the ESCRT-II component VPS36［J］. Nature Communications， 2023， 14： 7859

[38] Wang H J， Hsu Y W， Guo C L， Jane W N， Wang H， Jiang L W， Jauh G Y. VPS36-dependent multivesicular bodies are critical for plasma membrane protein turnover and vacuolar biogenesis［J］. Plant Physiology， 2017， 173（1）： 566-581

[39] Cui X X， Yang C H， Kang L， Zhu G Y， Wei Q Q， Jiang Y L. Expression pattern and regulation of head-to-head genes Vps36 and Ckap2 during chicken follicle development ［J］. Frontiers of Agricultural Science and Engineering， 2014， 1（2）： 130

[40] Kang L， Zhang Y J， Zhang N B， Zang L， Wang M， Cui X X， Jiang Y L. Identification of differentially expressed genes in ovaries of chicken attaining sexual maturity at different ages［J］. Molecular Biology Reports， 2012， 39（3）： 3037-3045

[41] Shen M M， Qu L， Ma M， Dou T C， Lu J， Guo J， Hu Y P， Yuan J， Sun C J， Wang K H， Yang N. Genome-wide association studies for comb traits in chickens［J］. PLoS One， 2016， 11（7）： e0159081

[42] Iwamori T， Iwamori N， Matsumoto M， Ono E， Matzuk M M. Identification of KIAA1210 as a novel X-chromosome-linked protein that localizes to the acrosome and associates with the ectoplasmic specialization in testes ［J］. Biology of Reproduction， 2017， 96（2）： 469-477

[43] Iwamori T， Iwamori N， Matsumoto M， Imai H， Ono E. Novel localizations and interactions of intercellular bridge proteins revealed by proteomic profiling［J］. Biology of Reproduction， 2020， 102（5）： 1134-1144

[44] Malcher A， Stokowy T， Berman A， Olszewska M， Jedrzejczak P， Sielski D， Nowakowski A， Rozwadowska N， Yatsenko A N， Kurpisz M K. Whole-genome sequencing identifies new candidate genes for nonobstructive azoospermia［J］. Andrology， 2022， 10（8）： 1605-1624

[45] Arora M， Mehta P， Sethi S， Anifandis G， Samara M， Singh R. Genetic etiological spectrum of sperm morphological abnormalities［J］. Journal of Assisted Reproduction and Genetics， 2024， 41（11）： 2877-2929

[46] Rui Q C， Wang J X， Li Y B， Tan X Y， Bao Y Q. Arabidopsis COG6 is essential for pollen tube growth and Golgi structure maintenance［J］. Biochemical and Biophysical Research Communications， 2020， 528（3）： 447-452

[47] Feng D B， Guo J H， Yan J J， Chen J， Ding L F， Zhu X Y， Chen Z L， Hu Y Y， Zhang M M， Liu J， Zhu C S， Liu M B， Zhao C， Zhang X Y， Xu J Q. COG6 is an essential host factor for influenza A virus infection［J］. Microbiology Spectrum， 2025， 13（10）： e01362-e0136225

[48] Pakanová Z， Krchňák M， Nemčovič M， Kodríková R， Ondrušková N， Štufková H， Giertlová M， Okáľová K， Stretavská P， Martineková S， Mezenská R Z， Urminská M， Škopková M， Andrésová A， Lysinová M， Belujská L， Šalingová A， Beke G， Račková L， Honzík T， Hansíková H， Baráth P. Insights into the pathological glycosylation associated with COG6-CDG［J］. Human Mutation， 2025， 2025（1）： 7948771

[49] Ververi A， Stathopoulou T， Kontou A， Farini M， Vlahou G， Demiris N， Sarafidis K. Lethal COG6-CDG in neonatal patient with arachnodactyly， joint contractures， and skin manifestations： Founder mutation in the Southeastern European population ［J］. Pediatric Dermatology， 2022， 39（2）： 314-315

[50] Taft R J， Vanderver A， Leventer R J， Damiani S A， Simons C， Grimmond S M， Miller D， Schmidt J， Lockhart P J， Pope K， Ru K L， Crawford J， Rosser T， de Coo I F M， Juneja M， Verma I C， Prabhakar P， Blaser S， Raiman J， Pouwels P J W， Bevova M R， Abbink T E M， van der Knaap M S， Wolf N I. Mutations in DARS cause hypomyelination with brain stem and spinal cord involvement and leg spasticity［J］. The American Journal of Human Genetics， 2013， 92（5）： 774-780

[51] Fröhlich D， Suchowerska A K， Spencer Z H T， von Jonquieres G， Klugmann C B， Bongers A， Delerue F， Stefen H， Ittner L M， Fath T， Housley G D， Klugmann M. In vivo characterization of the aspartyl-tRNA synthetase DARS： Homing in on the leukodystrophy HBSL ［J］. Neurobiology of Disease， 2017， 97： 24-35

[52] Zhang M， Zhou Z， Liu Z， Liu F， Zhao C. Exploring the potential biomarkers for prognosis of glioblastoma via weighted gene co-expression network analysis ［J］. Peer J， 2022， 10： e12768

[53] 刘宾. 基于慢病毒shRNA文库的功能基因组学筛选胃癌治疗靶基因及DARS在胃癌中的生物学功能研究［D］.长春：吉林大学，2022

[54] Liu B. Identifying therapeutic target genes for gastric cancer using functional genomics screening based on lentiviral shRNA library and studies on the biological function of DARS in gastric cancer ［D］. Changchun： Jilin University， 2022 （in Chinese）

[55] Huang C Q， Guo J H， Liu S， Shan Y X， Wu S L， Cai Y P， Yu L. Isolation， tissue distribution and prokaryotic expression of a novel human X-linked gene LHFPL1 ［J］. DNA Sequence， 2004， 15（4）：299-302

[56] Roly Z Y， Godini R， Estermann M A， Major A T， Pocock R， Smith C A. Transcriptional landscape of the embryonic chicken Mullerian duct［J］. BMC Genomics， 2020， 21（1）： 688

[57] Qin H M， Wang J Y， Jia X X， Zhi Y F， Sun L L， Zhang J L， Wang J L， Lu Y Q. Quantitative proteomics analysis of chicken embryos reveals key proteins that affect right gonadal degeneration in females［J］. Proteomics， 2025， 25（13）： e2200428

[58] Sutovsky P， Manandhar G， Laurincik J， Letko J， Caamaño J N， Day B N， Lai L X， Prather R S， Sharpe-Timms K L， Zimmer R， Sutovsky M. Expression and proteasomal degradation of the major vault protein （MVP） in mammalian oocytes and zygotes ［J］. Reproduction， 2005， 129（3）： 269-282

[59] Kim E， Lee S， Mian M F， Yun S U， Song M， Yi K S， Ryu S H， Suh P G. Crosstalk between Src and major vault protein in epidermal growth factor-dependent cell signaling［J］. The FEBS Journal， 2006， 273（4）： 793-804

[60] Kolli S， Zito C I， Mossink M H， Wiemer E A C， Bennett A M. The major vault protein is a novel substrate for the tyrosine phosphatase SHP-2 and scaffold protein in epidermal growth factor signaling［J］. The Journal of Biological Chemistry， 2004， 279（28）： 29374-29385

[61] Ducy P. The role of osteocalcin in the endocrine cross-talk between bone remodelling and energy metabolism［J］. Diabetologia， 2011， 54（6）：1291-1297

[62] Karsenty G. Update on the biology of osteocalcin［J］. Endocrine Practice， 2017， 23（10）： 1270-1274

[63] Zhang W J， Bai Y J， Hao L L， Zhao Y Q， Zhang L J， Ding W Q， Qi Y P， Xu Q. One-carbon metabolism supports S-adenosylmethionine and m6A methylation to control the osteogenesis of bone marrow stem cells and bone formation ［J］. Journal of Bone and Mineral Research， 2024，39（9）：1356-1370

基金资助

江苏省种业振兴揭榜挂帅项目(JBGS［2021］107)

国家肉鸡产业技术体系项目(CARS-41)

生物育种—国家科技重大专项(2023ZD04064)