绵羊HERC1 和ALDH6A1 多态性与胴体性状的关联分析

郭思武，刘玉芳※，储明星

（1.中国农业科学院北京畜牧兽医研究所，农业农村部动物遗传育种与繁殖重点实验室，北京 100193；2.河北工程大学生命科学与食品工程学院，河北 邯郸 056001）

胴体屠宰率是绵羊生产中的重要经济性状。胸椎数和胴体长是提高胴体屠宰率的重要指标，但其遗传机制尚不明晰。因此，挖掘相关候选基因多态性位点对提高绵羊胴体产量的分子标记辅助选育具有重要意义。HERC 蛋白家族由4 个蛋白组成，分为2 个亚家族（大型HERC 和小型HERC），其中大型HERCs（包括HERC1 和HERC2 两种蛋白）是长度近5 000 个氨基酸残基的巨型蛋白质，人中位于15 号染色体，绵羊中位于7 号染色体[1]。而小型HERC（包括HERC3 和HERC4 两种蛋白）的大小不到大型HERCs 的1/4[2]。HERC1（HECT and RLD domain containing E3 ubiquitin protein ligase family member 1）是该家族中最大的成员（532 kDa），包含许多保守的序列特征，这些特征被认为在蛋白质的整体功能中起着多种不同的作用[3]。

醛脱氢酶（ALDHs）是一个具有类似主要结构的NAD(P)+-依赖性酶的超家族，催化氧化广泛的内源性和外源性脂肪族和芳香族醛类[4]。到目前为止，在人类基因组中已经发现了16 个ALDH 家族基因，分布在染色体的不同位置，其中，ALDH3A2、ALDH4A1、ALDH5A1 和ALDH6A1 的多态性与通常以神经系统并发症为特征的代谢性疾病相关[4]。ALDH 基因超家族包括原核生物和真核生物物种中的250 多个ALDH基因[5]。迄今为止，真核生物ALDH 基因超家族至少由18 个家族组成，包括37 个亚家族[6]。ALDH6A1 参与缬氨酸和嘧啶的降解途径，分别催化丙二酸和丙二酸甲酯不可逆氧化脱羧为乙酰和丙酰辅酶A[7]。在大鼠的研究中发现，ALDH6A1 在肝脏、肾脏、心脏、肌肉和大脑中均有表达[7]。ALDH6A1 缺陷与发育迟缓有关，但目前尚不清楚这种缺陷的分子基础[8]。

迄今为止，未有HERC1 和ALDH6A1 突变位点与家畜胸椎数性状相关的文献报道，相关功能多集中在骨髓增生肿瘤[9]、乳腺癌[10]和脑部生长发育[11,12]等疾病研究领域。因此，本研究针对绵羊胸椎数变异相关候选基因HERC1 和ALDH6A1 SNPs 位点与胸椎数等性状进行关联分析，并通过对比突变前后蛋白理化性质和空间结构的差异，探讨HERC1 和ALDH6A1基因多态位点与胸椎数性状的关联性，为进一步拓展绵羊胸椎数变异的研究和分子标记辅助育种提供理论基础。

1 材料与方法

1.1 样品采集

实验选取天津市畜牧兽医研究所畜禽繁育基地饲养的383 只绵羊。其中苏尼特羊（SNT）来自内蒙古巴彦淖尔市乌拉特中旗，小尾寒羊（STH）来自山东郓城。胸椎数13 根（T13）样本259 只（SNT 122 只，STH 137只），胸椎数14 根（T14）样本124 只（SNT 66 只，STH 58 只）。绵羊健康状况良好，饲养环境一致。绵羊屠宰后，用2 mL RNase-Free 离心管快速采集新鲜背最长肌肌肉组织，用液氮冷冻运输带回实验室转移至-80℃超低温冰箱保存备用，记录小尾寒羊胴体长度数据。

1.2 DNA提取

使用DNA 提取试剂盒（天根生物科技有限公司，北京）对组织样品进行DNA 提取。提取后样品DNA用NanoDrop2000 和1.5%琼脂糖凝胶电泳进行浓度和质量评估，合格的DNA 样品供后续使用。

1.3 基因分型

根据绵羊HERC1 和ALDH6A1 SNPs 位点序列信息，使用Assay design3.1（Sequenom 公司）引物设计软件设计PCR 反应和单碱基扩展引物合成，引物信息见表1。采用Sequenom MassARRAY SNP[13]分型技术对HERC1 c.8950G ＞A（rs428482815）、g.43407058T ＞G、c.11902A ＞G（rs590493234）和ALDH6A1 c.1126A＞G（rs429743276）进行检测分型，分型样品为DNA，每个反应需要20L，DNA 浓度为40~80 ng/L。

表1 HERC1 与ALDH6A1 分型引物序列Table 1 Primer sequences for HERC1 and ALDH6A1 typing

1.4 数据统计分析

利用Excel 计算位点基因频率、基因型频率、有效等位基因数（Ne）、多态信息含量（PIC）、杂合度（He）及Hardy-Weinberg 平衡检验。计算PIC、He、有效等位基因数和Hardy-Weinberg 平衡检验等群体遗传学参数方法详见参考文献[14]。利用SPSS19.0 运用一般线性模型单变量分析与单因素方差分析LSD 对HERC1和ALDH6A1 不同基因型绵羊群体进行胸椎数关联分析。模型采用以下公式：

1.5 生物信息学分析

采用ProtParam（https://web.expasy.org/protparam/）分析目的蛋白理化性质；ProtScale（https://web.expasy.org/protscale/）分析目的蛋白疏水性；TMHMM Server v.2.0（http://www.cbs.dtu.dk/services/TMHMM/）分析目的蛋白跨膜区；运用SOPMA（https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl? page=npsa_sopma.html）分析目的蛋白二级结构，Phyre2（https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl? page=npsa_sopma.html）分析目的蛋白三级结构；STRING 11.0（https://string-db.org/cgi/input.pl）对目的蛋白进行蛋白互作分析。

2 结果与分析

2.1 HERC1 和ALDH6A1 多态性分析

383 只绵羊群体基因型分型（图1）分别是：HERC1 c.8950G ＞A 位点GG 型365 只、GA 型15 只，g.43407058T ＞G 位点TT 型376 只、GG 型1 只、GT型4 只，c.11902A ＞G 位点AA 型364 只、GG 型3只、GT 型13 只；ALDH6A1 c.1126A ＞G 位点TT 型364 只、CT 型18 只、CC 型1 只。群体遗传学分析表明（表2 和表3），STH 群体中HERC1 c.8950G ＞A、g.43407058T ＞G、c.11902A ＞G 位点的多态信息含量分别是0.035、0.030 和0.021，属低度多态（PIC＜0.25）；SNT 群体中HERC1 基因c.8950G ＞A、c.11902A ＞G位点的多态信息含量分别是0.041、0.075，属低度多态（PIC ＜0.25），SNT 群体中无g.43407058T ＞G 突变位点。STH 与SNT 群体中ALDH6A1 基因c.1126A ＞G位点的多态信息含量分别是0.084 及0.011，均属低度多态（PIC ＜0.25）。群体中SNPs 位点卡方检验值表明：STH 群体中HERC1 c.8950G ＞A、g.43407058T ＞G、c.11902A ＞G 与ALDH6A1 c.1126A ＞G 分别是0.798、0、0.885 和0.342；SNT群体中HERC1c.8950G＞A、c.11902A ＞G 与ALDH6A1 c.1126A ＞G 分别是0.764、0 和0.941，SNT 群体中无g.43407058T ＞G 突变位点，不计算其群体遗传学信息。结果表明，HERC1 基因g.43407058T ＞G 在STH 群体中不符合Hardy-Weinberg 平衡状态，且c.11902A ＞G 位点在SNT 群体中也不符合Hardy-Weinberg 平衡状态。其余SNPs位点在两种绵羊群体中均处于Hardy-Weinberg平衡状态（P ＞0.05）。

表2 HERC1 群体遗传学分析Table 2 Population genetic analysis of HERC1

表3 ALDH6A1 群体遗传学分析Table 3 Population genetic analysis of ALDH6A1

2.2 HERC1 和ALDH6A1 SNPs与胸椎数的关联分析

为探究HERC1 和ALDH6A1 SNPs 对SNT 和STH 胸椎数的影响，利用SPSS19.0 单因素方差分析分别在不同群体中分析SNPs 对胸椎数的影响。结果发现，HERC1 c.8950G ＞A、g.43407058T ＞G、c.11902A＞G 在SNT 与STH 中均与胸椎数呈不显著相关（P ＞0.05）；ALDH6A1 c.1126A > G 在SNT 与STH 中与胸椎数呈不显著相关（P ＞0.05）（表4 和表5）。采用一般线性模型单变量分析品种和基因型对绵羊不同胸椎数的影响。结果表明，在考虑品种效应时，基因位点基因型对于绵羊胸椎数的影响均不显著（P ＞0.05）（表6）。

表4 SNPs 位点基因型在不同品种中与胸椎数的关联分析结果Table 4 Results of association analysis between genotypes of SNPs loci in different varieties and the number of thoracic vertebrae

表5 SNPs 位点基因型小尾寒羊和苏尼特羊胸椎数最小二乘均值及标准误Table 5 Least squares means and standard errors for the number of thoracic vertebrae of genotypes of SNPs loci in Small Tail Han sheep and Sunit sheep

表6 SNPs 位点基因型、品种与胸椎数关联分析结果Table 6 Results of association analysis between genotype,variety and number of thoracic vertebrae at SNPs loci

2.3 HERC1 和ALDH6A1 SNPs位点与STH 胴体长的关联分析

在对SNT 和STH 群体HERC1 和ALDH6A1 突变位点多态性检测的基础上，通过单因素方差分析对记录胴体长数据的90只STH 进行SNPs 位点与胴体长的关联分析（表7）。结果显示，HERC1c.8950G ＞A 对STH 胴体 长 有 显 著 影 响（P ＜ 0.05），HERC1 g.43407058T ＞G、c.11902A ＞G 与ALDH6A1 c.1126A ＞G 在STH 胴体长性状中均不存在显著相关（P ＞0.05）。

表7 SNPs 位点基因型与STH 胴体体长的关联分析（平均值±标准误）Table 7 Association analysis of genotypes of SNPs loci with STH carcass length (mean±standard error )

2.4 HERC1生物信息学分析

2.4.1 HERC1 蛋白理化性质分析 针对HERC1 c.8950G ＞A，经ProtParam 在线软件预测得知（表8），HERC1 蛋白共4 882 个氨基酸，分子量由突变前的534 217.22 增加到534 231.25，等电点5.7，其中亮氨酸（Leu）所占比例最高（11.50%），色氨酸（Trp）所占比例最低（1.60%），不稳定性指数为47.42，属于不稳定蛋白。

表8 氨基酸数量与比例Table 8 Number and proportion of amino acids

2.4.2 HERC1 蛋白亲/疏水性分析 运用ExPASy 服务器上的ProtScale 程序预测HERC1 蛋白突变前后的疏水性变化，结果显示，HERC1 蛋白在1 766 位苏氨酸（T）亲水分数最大，达3.056，在2422位赖氨酸（K）亲水分数最小，仅为-3.756，总平均亲水性为-0.250，属于亲水蛋白（图2），突变前后并未引起明显的疏水性变化。

图2 亲水性分析Fig.2 Hydrophilic analysis

2.4.3 HERC1 蛋白跨膜结构预测与分析 运用TMHMM Serverv.2.0 对HERC1 蛋白跨膜结构进行分析（图3），结果显示，HERC1 蛋白突变前后均不存在跨膜螺旋区。

图3 蛋白质跨膜结构预测（突变前后一致）Fig.3 Protein transmembrane structure prediction(consistent before and after mutation)

2.4.4 HERC1 蛋白二级结构预测 通过ExPASy 的SOPMA 程序，针对HERC1 c.8950G ＞A 突变位置，预测其蛋白二级结构（图4），结果表明，突变前至突变后HERC1 蛋白-螺旋区域由1 720 个增加至1 721个，占比由37.22%增加至37.24%；-转角区域由327个增加至328 个，占比由7.08% 增加至7.10%；无规卷曲区域由1 780 个减少至1 778 个，占比由38.52%减少至38.48%；延伸链区域并无变化，共794 个，占17.18%。

图4 HERC1 蛋白质二级结构预测（截取部分，其他区域完全一致）Fig.4 Predicted secondary structure of HERC1 protein (intercepted part,other regions are identical)

2.4.5 HERC1 蛋白三级结构预测 采用Phyre2 软件对HERC1 蛋白突变前后进行同源建模获得三级结构模型。结果显示，HERC1 蛋白三级结构与二级结构一致（图5~6）。

图5 HERC1 蛋白质三级结构预测Fig.5 Predicted tertiary structure of HERC1 protein

图6 HERC1 蛋白质三级结构预测（氨基酸变化）Fig.6 Predicted tertiary structure of HERC1 protein (amino acid changes)

2.4.6 HERC1蛋白质的互作分析 采用STRING 11.0软件对HERC1 蛋白进行蛋白质互作分析，结果如图7 所示。

图7 蛋白质互作分析Fig.7 Protein interactions analysis

3 讨论

对绵羊等脊椎动物来说，脊椎的发育是从原肠胚中胚层开始的，其控制脊椎生长发育的初级神经胚主要受到骨形成蛋白调控[15]。其中，各阶段发育极为复杂，出现任何错误都可能造成脊椎分化异常[16]。由于猪繁殖快、产仔多的特点，使其在多脊椎性状相关候选基因研究中较其他家畜成果更多、进展更快。因此，可以将家猪多脊椎性状的研究作为参照，从遗传育种的角度发掘更多对绵羊胸椎数有显著影响的相关候选基因。目前，在绵羊多脊椎变异的研究中，已发现多个能够引起脊椎数变异的有效基因。刘建明等[17]在STH脊椎数变异的关联研究中指出，ActRⅡB 基因第85 位点发现的C ＞T 单核苷酸突变引发酶切位点消失，导致STH 第1 腰椎变异为第14 胸椎起到一定的作用，而该基因内含子上的1 个SNP 也对小尾寒羊脊椎发育存在影响。张立岭等[18]研究发现多脊椎畜禽在脊椎数目上比正常个体多1~2 个甚至更多，其中以蒙古羊和乌珠穆沁羊为首，多胸、腰椎个体的比例明显高于其他品种。Krumlauf等[19]也指出动物脊椎数的变化是由调控脊椎发育的同源框基因（简写Hox）发生突变导致的。虽然已存在许多绵羊胸椎数变异等性状的报道，但胸椎数变异受到多个QTLs的共同作用影响，其控制多脊椎性状的分子机制仍不明晰，从遗传育种角度有必要进一步发掘更多的对绵羊胸椎数有显著影响的候选基因。

HERC1 c.8950G ＞A 和c.11902A ＞G 在两种绵羊群体中表现为低度多态（PIC ＜0.25），表明这两个位点存在较强的选择潜力，这与Ensembl 数据库发现世界山羊群体中HERC1 c.8950G ＞A 位点G 等位基因与c.11902A ＞G 位点A 等位基因为优势等位基因（G 等位基因频率90%，A等位基因频率97%）相一致。有趣的是，HERC1 g.43407058T ＞G 是一个新位点，在SNT 群体中并未发现该突变位点，且本研究发现，该位点的突变杂合型（GT）与突变纯和型（GG）均高于野生型（GG），说明该基因突变位点在不同绵羊品种中可能存在选择性，具有成为影响胸椎数性状的候选基因潜力，这仍需进一步验证研究。另外，STH 群体HERC1 g.43407058T ＞G 位点以及SNT群体HERC1c.11902A ＞G 位点处于Hardy-Weinberg不平衡状态，可能是因为自然或人工选择干预，也可能与样品数量不足有关。关联分析结果显示，ALDH6A1 c.1126A>G 位点不同基因型与STH 及SNT 胸椎数和胴体长均无显著关联（P ＞0.05），但在胴体长关联分析中，该位点突变杂合型（CT）胴体长长度均值比野生型（TT）多2 cm 左右，说明该位点仍有可能影响胴体长性状，其中影响机制尚待研究发掘。而STH 群体中HERC1 c.8950G ＞A 位点在胴体长性状上GG 型显著高于GA 型，推测该位点突变在一定程度上降低了小尾寒羊的胴体产量。

通过HERC1 蛋白二级结构预测发现该蛋白没有跨膜结构域，这与HERC1 定位在高尔基体和胞浆等细胞内膜中[20]的结果相一致，表明该蛋白在细胞内膜转运中具有重要功能。而HERC1 c.8950G ＞A 位点的改变影响了HERC1 蛋白的二级与三级结构，使得突变前蛋白的两个无规则卷曲区域转变成一个-螺旋区域和一个-转角区域，且使蛋白氨基酸序列中突变位点的缬氨酸转变为突变后的异亮氨酸，推测其蛋白功能发生变化。同时，在群体遗传学分析中也发现该突变与胴体长显著相关。通过进一步对HERC1 蛋白进行蛋白质的互作分析发现，与其关联的SOX10 在神经发育及幼儿生长发育方面起着重要作用[21,22]。推测HERC1 可能与脊椎发育的起始原肠胚时期中胚层的初级神经胚的发育有关，这需要进一步探索研究。

4 结论

综上所述，HERC1 c.8950G ＞A、g.43407058T ＞G、c.11902A ＞G 和ALDH6A1 c.1126A>G 在小尾寒羊和苏尼特羊中均与胸椎数不存在显著相关（P ＞0.05）。小尾寒羊HERC1 c.8950G ＞A 位点与胴体长显著相关（P ＜0.05）。绵羊HERC1 为亲水蛋白，蛋白结构不稳定，无跨膜结构域，可能不参与跨膜间生化反应。HERC1 c.8950G ＞A 位点可作为潜在影响绵羊胴体长度的分子辅助选择标记。