海南黑山羊GH1基因nsSNPs功能性预测及生物信息学分析

管凇 周汉林 徐铁山 孙卫平 胡海超

（中国热带农业科学院热带作物品种资源研究所 571101）

生物信息学技术是现阶段筛选候选功能位点最有效的技术之一，因此实验研究中可利用该技术预测未知表型的nsSNPs的功能[1]。已知，现阶段研究人员利用各类生物信息分析工具对人类或动物表型和性状基因功能性突变进行了大量的分析和预测，且取得了突破性进展[2-6]。本试验主要通过利用生物信息学技术预测海南黑山羊GH1 基因编码区nsSNPs 的功能。旨在筛选可显著影响海南黑山羊生长的功能性突变位点，并分析该突变位点的作用机理，以期更好的理解黑山羊GH1 基因的作用机理，同时也对应用标记辅助选择改善山羊生长发育提供一定的参考。

1 材料和方法

1.1 样品采集和DNA 处理

海南黑山羊样品采自品资所山羊保种场山羊耳组织。对104 只成年山羊的体长、体高、体重进行测量。剪取耳组织，置于装有75%乙醇的离心管中，4℃保存。采用DNA 提取试剂盒提取耳组织总DNA，用紫外分光光度计检验每个样品DNA质量浓度。

1.2 引物设计

从NCBI 数据库中查询山羊GH1 和IGF2 基因DNA 序列（GenBank 登录号分别为：NC_030826 和NC_030836），根据单核苷酸多态性位点的序列信息，利用引物设计软件Primer Plex 2 设计引物，使用primer3plus 在线设计单碱基延伸引物，引物信息见表1。

表1 引物信息

1.3 SNP 位点的筛选和分型

对提取的DNA 基因组序列进行扩增，PCR 扩增反应总体积为10μL：5μL5×TaqPCRmastermix，上下游引物各0.5μL（20μmol/L），模板DNA1μL (20ng/ul），加ddH2O 至10μL。PCR 扩增程序：94℃预变性5min；45 个循环（94℃，20s；60℃，30s；72℃，30min）；72℃再延伸3min。将PCR 产物去除体系中游离的dNTPs，再进行单碱基延伸反应，反应体系总体积5μl。PCR 产物委托武汉天一辉远生物科技有限公司采用SnaPshot 方法，用Gene marker 分析得到的.fsa 结果，导出峰图及表格文件，并统计出每个样品的SNP 位点的基因型。

1.4 数据分析

用 ROVEN（Protein Variation Effect Analyzer）http://provean.jcvi.org/index.php 和SITF（Sorting Intolerant From Tolerant）https://sift.bii.a-star.edu.sg/和分析nsSNPs 对蛋白功能的影响，利用http://rna.tbi.univie.ac.at/cgi-bin/RNAfold.cgi 软件进行mRNA 二级结构预测[7]。利用Sopma（https://npsa-prabi.ib -cp.fr/cgi -bin/npsa_automat.pl?page=npsa_sop -ma.htm 和CPHmodels（http://www.cbs.dtu.dk/services/CPHmodels/）软件分别对突变体的二级结构和三级结构分别进行预测和分析，进而评估nsSNPs 可能造成的影响，三级结构3D 结果展示使用软件BIOVIA（https://www.3dsbiovia.com）。

2 结果与分析

2.1 单核苷酸多态性位点分析

对海南黑山羊GH1 基因的测序比对结果显示，在GH1 基因外显子1、2、3、4 区域检测到10 个SNP 位点，有4 个错义突变位点，4 个错义突变突变位点测序峰图见图1 和表2。在第一外显子350bp 处发生G→A 错义突变，突变位点导致氨基酸由原来的甘氨酸（Gly）变为丝氨酸（Ser），第二外显子739bp 处发生G→A 错义突变，突变位点导致氨基酸由原来甘氨酸（Gly）变为丝氨酸（Ser），第三外显子1101bp 处发生A→G 错义突变，突变位点导致氨基酸由原来天冬氨酸（Asp）变为甘氨酸（Gly），第三外显子1154bp 处发生C→T 错义突变，突变位点导致氨基酸由原来精氨酸（Arg）变为色氨酸（Trp），其余6 个同义突变位点没有导致氨基酸改变。GH1 基因Exon2-729G→A、Exon2-739G→A 两个位点在群体中分布的等位基因型频率突变前后相同，两个位点连锁，Exon3-1101A→G、xon3-1120G→A、Exon3-1154C→T3 个位个等位基因型频率突变前后相同，完全连锁。

图1 错义突变位点测序峰图

表2 突变位点信息

2.2 不同SNP 位点与性状关联分析

利用SPSS 软件统计各基因的不同SNP 对应的表型（体长，体高和体重）的平均值及标准误差，详细结果见下表4。利用SPSS 中的一般线性模型对各基因的SNPs 位点与表型进行关联分析，GH1 基因外显子1、2、3、4 检测到10 个SNP 位点，每个位点均只有2 种基因型。一般线性模型的主体间效应检测结果显示，所有SNP 位点基因型与表型（体长，体高和体重）均未有显著差异。

表4 突变位点与生长性状的关联性分析

2.3 SNP 位点的生物信息学分析

2.3.1 nsSNPs 对蛋白功能的影响

2.3.1.1 ROVEN 分析

用ROVEN 分析nsSNPs 对蛋白功能的影响，当SNP 的PROVEAN score 分数等于或低于-2.5 的变体被认为是“有害的”。当SNP 的PROVEAN score 分数高于-2.5 的变体被视为“中性”突变。分析结果显示，D129G（天冬氨酸变为甘氨酸）和R147W（精氨酸变为色氨酸）SNP 的PROVEAN score 分数低于-2.5，该突变位点属于有害突变，对蛋白功能有较大影响见表5。

表5 nsSNPs 的score 评判分值

2.3.1.2 SITF 分析

用SITF 分析nsSNPs 对蛋白功能的影响，评级为4 个标准：tolerated、low confidence tolerated、low confidence deleterious、deleterious。当SNP 的score 大于0.5 时，说明该SNP 位点对蛋白质功能没有影响或者影响较小；当SNP 的score 小于0.5 时，说明这个突变的SNP 位点对蛋白质功能有较大影响。一个SNP 的得分越低，危害性越大，分析结果显示，G85S（甘氨酸变为丝氨酸）、D129G（天冬氨酸变为甘氨酸）、R147W（精氨酸变为色氨酸）突变位点score＜0.5，该突变位点属于有害突变，对蛋白功能有较大影响，见表6。

表6 nsSNPs 的score 评判分值

2.3.2 mRNA 二级结构分析

对引起错义突变的GH1 基因Exon1-350G→A、Exon2-739G →A；Exon3-1101A →G、Exon3-1154C →T位点进行mRNA 二级结构预测，结果显示，4 个突变位点突变前后均引起mRNA 二级结构和自由能发生变化，Exon1-350G→A 位点的自由能和结构变化较明显。说明基因序列SNPs 突变前和突变后mRNA 自由能的改变、二级结构的改变均导致其结构的稳定性发生了改变，见图2 和表7。

表7 mRNA 二级结构突变前后自由能变化

图2 mRNA 二级结构突变前后的变化结果

2.3.3 蛋白质二级结构分析

对引起错义突变的GH1 基因Exon1-350G→A、Exon2-739G→A；Exon3-1120G→A、Exon3-1154C→T 位点进行进行蛋白质二级结构分析，结果表明：Exon1-350G→A（G31S）位点突变前后的无规则卷曲、α-螺旋、β-转角、扩展链均无变化，该突变位点对蛋白质二级结构没有发生影响；Exon2-739G→A（G85S）位点突变前后的无规则卷曲和扩展链发生改变，Exon3 -1120G →A（D129G）、Exon3 -1154C →T（R147W）位点突变前后的α-螺旋和无规则卷曲发生改变，这3 个突变位点对蛋白质二级结构发生影响。统计结果见表8。

表8 蛋白质突变前后二级结构预测分析

2.3.4 位点突变前后蛋白质三级结构分析

突变位点Exon1-350G →A（G31S）、Exon2-739G →A（G85S）、Exon3-1120G →A（D129G）、Exon3-1154C →T（R147W）突变前后蛋白质三级结构进行预测分析，结果三级结构与二级结构变化一致，Exon1-350G→A（G31S）位点突变前后三级结构无变化，该突变位点对蛋白质三级结构没有发生影响；Exon2-739G→A（G85S），Exon3-1120G→A（D129G）、Exon3-1154C→T（R147W）位点突变前后三级结构发生变化（见图3），预测这3 个突变位点对蛋白质的结构功能可能会造成影响。

图3 突变前后蛋白质三级结构

3 讨论

随着生物信息学的不断发展，生物信息学工具已广泛应用于人类医学和农业领域中对功能性SNP 研究中。该工具可以显著降低SNPs 的数目，且能筛选出最有可能的功能性nsSNPs，为后续功能试验验证提供一定的参考依据。因此，使用该工具可显著降低研究成本和工作量[7]。目前较常用的工具为SIFT和ROVEN，建议联合使用，以扩大数据库，提高检测的灵敏度和准确度。单独使用某一工具时预测结果不同，SIFT 预测结果假阳性率更低[8-9]。ROVEN 对高风险位点的预测结果更准确，但对中风险位点的预测准确度较低[10]，本研究同时采用SIFT 和ROVEN 两种方法对海南黑山羊GH1 基因的外显子区的功能性nsSNPs 进行了预测（表5-6）。预测G85S、D129G、R147W3 个位点可能是有害突变位点。D129G、R147W 两个位点突变前后二级结构的α-螺旋比列下降，无规则卷曲比例有上升，G85S 位点的扩展链发生变化。通过nsSNPs、mRNA二级结构、蛋白质二级结构和蛋白质三级结构分析结果，预测G85S、D129G、R147W3 个突变位点可能是导致海南黑山羊个体矮小，生长缓慢的功能性SNP。G85S、D129G、R147W3 个突变位点分别使氨基酸由甘氨酸（Gly）变为丝氨酸（Ser），天冬氨酸（Asp）变甘氨酸（Gly），精氨酸（Arg）变为色氨酸（Trp）。甘氨酸是在氨基酸系列中结构最为简单的非极性氨基酸，在中枢神经系统，尤其是在脊椎里，甘氨酸是抑制性神经递质；丝氨酸是一种非必需氨基酸，在脂肪和脂肪酸的新陈代谢及肌肉的生长中发挥着作用；天冬氨酸是一种α-氨基酸，对肝脏功能有增强作用，消除疲劳；精氨酸为碱性氨基酸，是维持婴幼儿生长发育必不可少的氨基酸；色氨酸对泌乳期的乳牛和母猪有促进泌乳作用，当畜禽缺乏色氨酸时，生长停滞，体重下降，脂肪积累降低。本文预测到G85S、D129G、R147W3 个突变位点属于有害突变，对蛋白功能有较大影响，可能是导致海南黑山羊个体矮小，生长缓慢的功能性SNP。