绵羊MTNR1A基因SNPs筛选及生物信息学分析

王明远，张宾，于刚，王雨曈，杨永林，余乾，赵宗胜，杨华

(1 石河子大学动物科技学院，新疆 石河子 832000；2 新疆农垦科学院省部共建绵羊遗传改良和健康养殖国家重点实验室，新疆 石河子 832000；3 北屯市阿山屯科肉用种羊繁育有限公司，新疆 北屯 836000；4 新疆生产建设兵团第十师181团农业发展服务中心，新疆 北屯 836001)

繁殖性状是绵羊重要的经济性状，繁殖力的高低直接影响养羊业经济效益。母羊的繁殖力主要表现为产羔数和季节发情特征，褪黑激素(melatonin，Mel)是一种重要的吲哚类内分泌激素，在哺乳动物中主要由视网膜和松果体分泌，在动物繁殖中具有重要的生物学作用。根据光照周期的表现，褪黑激素分泌呈现昼低夜高，其生物学功能主要通过与高亲和力的褪黑素受体(melatonin receptor，MTNR)结合对季节性发情动物的昼夜节律和生殖活动发挥调节作用[1]。

绵羊MTNR1A基因位于第26号染色体，包含2个外显子与1个内含子[2]。目前，研究人员已将MTNR1A基因作为调控绵羊繁殖性能的候选基因之一，其多态性位点大多数位于外显子2区域，并证明与产仔数、季节性发情有关。PELLETIER等[3]发现了绵羊MTNR1A基因外显子2的606和612两个沉默突变与繁殖相关。MATEESCU等[4]研究表明，多赛特母羊MTNR1A基因外显子2的612位点MM和Mm基因型相比于mm基因型具有更高的产羔数。MARTNEZ R等[5]证实在Rasa Aragonesa绵羊MTNR1A基因外显子2存在11个SNPs，启动子区存在17个SNPs，SNP606/RsaI等位基因与母羊常年发情相关，启动子的422和677位点和外显子2的612位点T、A和T等位基因的单倍型与常年发情相关。MUATEESCU等[6]研究表明，Sarda绵羊MTNR1A基因外显子2的612位点的MM基因型与603～606位点的RR基因型具有更高的的受孕率与产羔数。由此可见，研究人员已将MTNR1A基因作为绵羊繁殖力的关键基因来研究，并取得了较多的研究成果。我国拥有丰富的绵羊品种资源，其中湖羊是引入新疆的高繁殖力和常年发情绵羊品种，多胎萨福克羊是新疆农垦科学院培育的多胎新品系肉羊[7]，中国美利奴羊(新疆军垦型)是季节性繁殖的毛用绵羊品种。在这些不同繁殖特征的绵羊中，MTNR1A基因是否存在影响绵羊繁殖功能的新SNPs，SNPs是否引起MTNR1A基因功能改变值得进一步研究。

本试验以高繁殖力的湖羊和多胎萨福克羊，以及中国美利奴羊(新疆军垦型)单胎羊为研究对象，对MTNR1A基因具有高度多态性的外显子2区域进行克隆和测序，筛选重要的SNPs，并分析序列特征，预测基因功能，为进一步阐明MTNR1A基因在绵羊繁殖活动中的作用提供基础。

1 材料与方法

1.1 样品采集及DNA提取

试验样本选择周岁健康的201只母羊，包括45只湖羊，94只多胎萨福克羊和62只中国美利奴(新疆军垦型)单胎羊。采集绵羊颈静脉血，制备EDTA-Na2抗凝血，利用基因组DNA提取试剂盒(Omega公司)提取血液基因组DNA。

1.2 引物设计与合成

根据NCBI数据库(https://www.ncbi.nlm.nih.gov/)中公布的绵羊MTNR1A基因序列(GenBank登录号：NM_001009725.1)，与UCSC(http://genome.ucsc.edu/)中的绵羊参考基因组数据(Oar_v4.0，2015年11月)比对分析，获得MTNR1A基因的exon 2及其上下游序列，利用Primer Premier 5.0设计引物，引物MTA2扩增的序列覆盖MTNR1A基因exon2全长，引物信息见表1，引物由生工生物工程(上海)有限公司合成。

表1 扩增MTNR1A基因exon 2的引物信息

1.3 MTNR1A基因exon 2的克隆

随机选择3个品种(系)绵羊的基因组DNA各3只，应用PCR扩增MTNR1A基因exon 2，体系为50 μL：2×Es Taq MasterMix 25 μL，上下游引物各1 μL，基因组DNA模板2 μL，ddH2O 21 μL。MTNR1A基因exon 2的PCR扩增程序为94 ℃预变性5 min；94 ℃变性30 s，49 ℃ 30 s，72 ℃延伸30 s，共35个循环；最后72 ℃终延伸5 min；PCR扩增产物用1.5%琼脂糖凝胶电泳检测。

根据琼脂糖凝胶回收试剂盒说明书(天根生化科技(北京)有限公司)分别回收纯化目的PCR产物。经25 ℃ 20 min连接回收产物与pEASY-T1 Simple TA载体，连接产物转化到Trans1-T1感受态细胞中。筛选阳性重组子，将阳性克隆菌液、MTNR1A基因exon 2的 PCR产物送往生工生物工程(上海)有限公司测序。

1.4 MTNR1A基因exon 2的序列比对和多态性分析

根据测序获得的湖羊、多胎萨福克羊和中国美利奴羊MTNR1A基因exon 2序列，利用DNAMAN9.0对DNA序列进行比对分析，获得基因新的SNPs位点。以3个品种(系)201只绵羊为试验群体，应用Sequenom MassARRAY®SNP分型技术检测突变位点的真实性，并分析不同基因型在各绵羊群体中的分布，以及基因突变所引起的MTNR1A氨基酸变化。根据绵羊参考基因组数据获得突变位点上下游200 bp序列，MassARRAY引物由北京康普森生物技术有限公司设计并合成。

1.5 生物信息学分析

根据绵羊MTNR1A基因参考序列(GenBank登录号：NM_001009725.1)，利用DNAstar分析MTNR1A基因不同物种之间核苷酸序列的同源性；应用NovoPro在线转换工具获得氨基酸序列；通过ExPASy Protparam(https://web.expasy.org/protparam/)、SOPMA、PSORT Ⅱ prediction、TMHMM Server 2.0、NetPhos 3.1 Server(http://www.cbs.dtu.dk/services/NetPhos/)、SignalP(http://www.cbs.dtu.dk/services/SignalP/)、NetOGlyc 4.0 Server、NetNGlyc 1.0 Server、SWISS-MODEL(https://swissmodel.expasy.org/)在线软件分析MTNR1A蛋白的基本理化性质、二级结构、亚细胞定位、跨膜结构域、磷酸化位点、信号肽、O糖基化位点、N糖基化位点以及三级结构预测，通过Interpro(www.ebi.ac.uk/interpro/)在线分析网站对MTNR1A蛋白功能进行预测。

1.6 数据分析

计算绵羊群体中各位点的基因型频率和基因频率，利用Haploview 4.2软件进行哈代-温伯格平衡检验、群体杂合度(Het-erozygosity，Ho)和最小等位基因频率(MAF)的计算，利用PIC_CALC 0.6软件计算多态信息含量(PIC)。

2 结果与分析

2.1 MTNR1A基因的PCR扩增

MTNR1A基因exon2的PCR产物用1.5%琼脂糖凝胶电泳检测，PCR产物约为1 155 bp，扩增产物大小与预期结果一致，特异性好(图1)。

1-9为MTNR1A基因PCR产物，C为空白对照，M为DL2000 DNA Marker图1 MTNR1A基因PCR扩增产物的电泳检测

2.2 MTNR1A核苷酸序列比对分析

比对湖羊、多胎萨福克羊和中国美利奴羊MTNR1A基因exon 2序列，发现在绵羊群体中MTNR1A的exon 2存在8个突变位点。

其中，C378T、G405T、G564 A、A735G、A753G和T843C为同义突变，A658G和A845C为错义突变；其中A658G导致第220位氨基酸缬氨酸变成异亮氨酸(V220I)，A845C导致第282位氨基酸处天冬氨酸变成丙氨酸(D282A)。

2.3 MTNR1A基因单核苷酸多态性分析

在绵羊群体中应用Sequenom MassARRAY®SNP分型技术验证8个SNPs位点的真实性(图2)，进一步分析绵羊群体中MTNR1A基因8个SNPs的遗传多态性。MassARRAY系统具有灵敏度高，精确度高，检测速度快，成本低等优点并且可以根据需求灵活调节，本次试验8个突变位点在201个个体中检出率达到100%。表2可见，MTNR1A基因C378T位点在湖羊和多胎萨福克羊中只存在CC基因型和C等位基因，而在中国美利奴羊中存在CC、CT和TT3种基因型，并处于哈代温伯格不平衡状态(P<0.05)。其余7个位点在3种绵羊品种中均存在3种基因型，并处于哈代温伯格平衡状态。T405G位点在湖羊和多胎萨福克羊群体中G为优势等位基因，而在中国美利奴羊中T为优势等位基因；T843C位点在湖羊中C为优势等位基因，而在多胎萨福克羊和中国美利奴羊群体中T为优势等位基因；C378T、G564 A、A658G、A735G、A753G和A845C位点在3个绵羊品种中均有相同的优势等位基因。

图2 MTNR1A基因突变位点的MassARRAY分型检测结果

T405G、G564 A、A658G、A735G与A753G位点在湖羊，A845C位点在湖羊和中国美利奴羊中呈低度多态性(PIC<0.25)，这些位点在其余品种中均表现为中度多态性(0.25

表2 MTNR1A基因不同SNP位点在三个品种(系)绵羊中的群体遗传学分析

表2续

2.4 MTNR1A基因核苷酸序列同源性分析

经序列比对分析表明，3个绵羊品种MTNR1A基因的核苷酸序列与Genbank数据库中绵羊序列(登录号NM_001009725.1)的核苷酸序列同源性均达到99.27%以上，表明克隆的序列确实为绵羊MTNR1A基因869 bp的exon2区域。核苷酸序列与已发表的Hassan羊(Genbank登录号KC580735)、Malpura羊(Genbank登录号KC757710)、Kashmir羊(Genbank登录号HQ658144)、AvikaLIN羊(Genbank登录号JF901325)、Musimon羊(Genbank登录号FJ801038)MTNR1A基因同源性均在99.3%以上，与山羊(Genbank登录号XM_018041838.1)、牛(Genbank登录号XM_002698656.5)、人(Genbank登录号NM_005958.4)、小鼠(Genbank登录号NM_008639.3)同源性分别为98.2%、97.0%、84.6%、80.6%(图3)。

MTNR1A基因核苷酸序列同源性比较，1～12分别为多胎萨福克羊、湖羊、中国美利奴羊、Hassan羊、Malpura羊、Kashmir羊、AvikaLIN羊、Musimon羊、山羊、普通牛、人、小鼠。右上为核苷酸序列同源性结果，左下为核苷酸序列散度(分歧)结果图3 MTNR1A基因核苷酸序列同源性分析

2.5 MTNR1A基因编码蛋白质序列的基本信息和功能预测

利用ExPASy Protparam软件分析突变前后MTNR1A基因编码的蛋白质基本理化性质，结果显示突变前后MTNR1A的蛋白分子式均为C1874H2932N476O482S18，由366个氨基酸组成，总平均亲水性为0.614，表现为疏水性，蛋白质的理化性质均没有发生变化。

绵羊MTNR1A突变引起了7处蛋白二级结构的改变(图4)，突变前的α螺旋为37.43%，β转角为5.46%，无规则卷曲为36.34%和延伸链为20.77%。

突变后α螺旋为37.98%，β折叠为4.92%，无规则卷曲为36.07%和延伸链为21.04%。MTNR1A蛋白的亚细胞定位表明，MTNR1A主要分布在内质网，分布比例为55.6%；磷酸化位点预测结果显示MTNR1A蛋白存在20个丝氨酸残基磷酸化位点，6个苏氨酸残基磷酸化位点；糖基化位点预测结果显示MTNR1A蛋白含有6个O-糖基化潜在位点和2个N-糖基化潜在位点。基因8个位点的突变没有引起磷酸化、糖基化的变化。经蛋白功能预测，MTNR1A属于视紫红质样G蛋白偶联受体家族中褪黑素受体家族，MTNR1A在129-145含有G蛋白受体F1_2结构域蛋白。MTNR1A生物学过程参与G蛋白偶联受体信号通路 (GO：0007186)；分子功能为G蛋白偶联受体活性 (GO：0004930)、褪黑素受体活性 (GO：0008502)；细胞学组件均为膜整体组分 (GO：0016021)。

蓝色长竖线为α螺旋，红色中长线为延伸链，绿色中短线为β转角，紫色短线为无规则卷曲；箭头标识处为突变后二级结构有变化的位置图4 MTNR1A突变前后蛋白二级结构示意图

2.6 MTNR1A蛋白质三级结构分析

根据同源建模方法预测的MTNR1A蛋白质三级结构见图5。突变前序列与模板相似度61.36%，GMQE=0.72，QMEAN=-3.68，可信度较高；突变后序列与模板相似度61.36%，GMQE=0.71，QMEAN=-3.69，可信度较高。可见，2个错义突变引起的氨基酸改变没有导致蛋白质三级结构的显著变化。

图5 MTNR1A突变前后蛋白三级结构示意图

3 讨论

褪黑激素是一种在夜间由松果体分泌的多效性信号分子，光周期信号被褪黑素转化为神经内分泌信号后传递到下丘脑-垂体-性腺轴，能够调节多种生理功能。在调节动物繁殖过程中，光周期信号会通过视交叉上核(Suprachiasmatic nucleus，SCN)传递到松果体，使松果体分泌褪黑激素与垂体结节部褪黑激素受体结合产生促甲状腺激素，再与室管膜细胞上的促甲状腺激素受体结合引起室管膜细胞中脱碘酶2(DIO2)和脱碘酶3(DIO3)表达变化，而DIO2和DIO3能够调节T3/T4的水平，进而调节GnRH的释放，最终调控动物的季节性繁殖[8]。研究表明，褪黑素是通过与高亲和力的褪黑素受体结合发挥其生殖和生理效应，褪黑激素受体是一种经典的G蛋白，带有七个跨膜结构域，一个细胞外N-末端结构域和三个细胞外环的耦合受体[9]。

单核苷酸多态性(single nucleotide polymorphism，SNP)是在基因组水平上由单个核苷酸变异引起的DNA序列多态性，根据对氨基酸序列的影响分为同义突变、错义突变、无义突变以及延长突变，其中错义突变会导致蛋白质功能发生变化。基因突变所引起的任何氨基酸变化都可能影响受体与配体之间的相互作用。在MTNR1A基因多态性的研究报道中，MTNR1A基因存在较多SNPs位点，而且大多数都不相同，表现出品种差异性。程笃学等[10]以常年发情的小尾寒羊和季节性发情的多赛特羊共20只母羊为研究对象，在小尾寒羊的MTNR1A基因外显子2上发现C329T、G355T、C566T、C580 A和A675G5个核苷酸变化。CALVO等[11]证明Rasa-Aragonesa羊存在rs403212791和rs403212791突变位点。SAXENA等[8]以101只Chokla绵羊为研究对象，在MTNR1A基因ExonⅡ区域克隆测序发现10个突变，其中G706 A、C893 A和G931C为非同义突变，在印度热带干旱绵羊品种中发现C426T、G555 A和G706 A三个突变位点；LURIDIANA等[12]发现Sarda 母羊存在C606T和G612 A两个突变位点，与产羔数关联分析证明两位点与生殖活动相关联。本次试验以高繁殖力的湖羊和多胎萨福克羊，以及中国美利奴羊(新疆军垦型)单胎羊为研究对象，在MTNR1A基因exon2共发现8个SNPs(C378T、G405T、G564 A、G658 A、A735G、A753G、T843C、A845C)，均为未曾报道过的新发现的多态位点，其中A658G引起氨基酸变化(V220I)，该突变在第5跨膜结构域可能会导致蛋白结合能力改变。A845C引起氨基酸变化(D282 A)，该突变在第三个胞外环没有明确的生理功能。功能预测表明，A658G和A845C位点所引起的氨基酸的改变，导致蛋白质二级结构发生变化，但是没有引起翻译后修饰、拓扑结构以及三级结构的改变。王世银等[13]以季节性发情的阿勒泰羊和常年发情的湖羊为研究对象，通过PCR-SSCP检测g.15143697位点的T/C突变发现，湖羊以CC基因型为主,其基因型频率达到0.919，阿勒泰羊中以TT基因型为主,基因型频率为0.950,推测该位点可能与绵羊的季节性发情有关。与之相似的李广录[14]研究结果显示MTNR1A基因与绵羊的繁殖性状有一定联系，而王琼等[15]研究显示MTNR1A基因与策勒黑羊繁殖性能关系不大。本研究根据突变位点在2个高繁殖力绵羊品种(湖羊和多胎萨福克羊)与中国美利奴羊(新疆军垦型)单胎群体中的基因型分布，C378T位点在多胎的湖羊、多胎萨福克羊中只存在CC基因型，而在单胎中国美利奴羊中存在CC、CT和TT三种基因型，由此可见单胎中国美利奴羊与具有多胎性能的湖羊和多胎萨福克羊相比有特殊的T等位基因，初步认为C378T位点突变影响绵羊繁殖性能。然而，由于试验样本为周岁母羊，MTNR1A基因C378T突变位点引起的繁殖功能变化有待于试验母羊配种以后统计产羔情况后关联分析验证，或进一步扩大群体在经产母羊群体中验证功能。

蒋维东等[16]研究表明和田羊MTNR1A基因与已发表的绵羊属Chokla、Musimon和Patanwadi(登录号分别为KC727711、FJ801038、KJ865682)三个品种羊的同源性均高于98%；与牦牛、野牛、牛、山羊和藏羚羊(登录号分别为KF569810、XM010835821、U73327、KC865680、KJ005972610)的同源性为96%～98%。程笃学等[17]研究表明小尾寒羊MTNR1A基因外显子2氨基酸序列与反刍动物差异最小。赵艳红等[18]研究表明内蒙古绒山羊MNTR1A基因第2外显子氨基酸序列与已发表的山羊氨基酸序列的同源性为97.7%，与已发表的绵羊、牛、人和鼠的氨基酸同源性分别为94.1%、91.8%、82.2%、83.1%。本试验发现，绵羊品种间的MTNR1A核苷酸序列同源性达到99%以上，与牛科动物山羊、牛同源性较高也在98%、96%以上，与人、鼠同源性在80%以上。研究结果与蒋维东等、程笃学等与赵艳红等相似，说明MTNR1A基因的核苷酸序列具有较高的的保守性。生物信息学分析表明，MTNR1A是疏水性的非极性分子，MTNR1A蛋白存在20个丝氨酸残基磷酸化位点，6个苏氨酸残基磷酸化位点，这些修饰位点可能通过丝氨酸和苏氨酸磷酸化作用以激活MTNR1A蛋白质的活力。进一步基因功能预测也显示，MTNR1A的生物学作用主要富集在细胞表面受体信号传导途径与G -蛋白偶联受体信号通路，可以促进在相关三聚体G蛋白复合物的α亚基上GTP/GDP的交换，表明MTNR1A在绵羊繁殖调控信号传导中扮演了重要角色。

综上所述，在湖羊、多胎萨福克羊和中国美利奴羊中发现MTNR1A基因存在8个新的SNPs位点(C378T、G405T、G564 A、G658 A、A735G、A753G、T843C、A845C)，MTNR1A基因具有较高的遗传多态性。其中，C378T突变位点可能与绵羊繁殖性能密切相关。下一步将扩大试验绵羊群体，根据产羔数等表型数据进一步研究MTNR1A基因功能性SNPs，为阐明MTNR1A基因的繁殖学功能，以及开展高繁殖力绵羊的标记辅助育种提供基础。