MET基因多态性与中国荷斯坦牛繁殖和泌乳性状的关联分析

许静漪，徐昊祺，胡丽蓉，张 帆，高 清，罗汉鹏，张海亮，师 睿，李 想，刘 林，郭 刚，王雅春*

(1.中国农业大学动物科学技术学院 农业农村部动物遗传育种与繁殖(家畜)重点实验室，畜禽育种国家工程实验室，北京100193； 2.新疆农业大学动物科学学院，乌鲁木齐 830052；3.北京奶牛中心，北京 100194； 4.北京首农畜牧发展有限公司，北京 100029)

奶业是健康中国、强壮民族不可或缺的产业，中国不仅是世界牛奶产量第三大经济体，而且也已经成长为世界第一大乳制品进口贸易国[1]。为了满足消费者对牛奶及乳制品的需求，在过去的几十年间，奶牛育种者们对奶牛的产奶性能进行了大幅度的选育，但也导致了繁殖力、抗病力的下降，严重阻碍了奶业的可持续发展。随着平衡育种理念广泛应用于实际选育，越来越多的性状，如繁殖、长寿等功能性状得到了育种工作者和牧场生产经营者的关注。

本课题组前期对中国荷斯坦牛的多项繁殖性状进行了全基因组关联分析(genome-wide association study，GWAS)，获得了一系列的候选基因。通过蛋白互作网络(protein-protein interaction networks，PPI)分析，发现MET(MET proto-oncogene, receptor tyrosine kinase)基因与多个候选基因存在互作关系，表明该基因可能是影响中国荷斯坦牛繁殖性状的重要基因。MET是一种原癌基因，又称c-MET。牛的MET基因包含21个外显子和20个内含子，其cDNA全长4.8 kb，包含4 152 bp的开放阅读框，编码由1 384个氨基酸组成的Met蛋白[2]。Met蛋白是一种由二硫键连接的异二聚体组成的酪氨酸激酶，由190 ku的前体裂解为成熟的胞外α链和较长的跨膜β链[3]。已有研究证明，肝细胞生长因子(hepatocyte growth factor，HGF)为Met的唯一配体[4-5]。Met与HGF结合能够激活多种信号转导通路，诱导多种生物反应，如增殖、凋亡抑制、运动、细胞外基质侵袭和血管生成的激活[6-10]。Parrott和Skinner[11]研究发现，MET与牛卵泡发育有关，HGF与Met在牛卵泡发育过程中高度表达，HGF由卵泡膜细胞表达，作用于卵泡膜细胞和颗粒细胞上的Met，雌激素和促黄体生成素(luteinizing hormone，LH)通过调节局部HGF的产生间接控制卵泡发育，促进颗粒细胞的增殖，从而促进卵泡生长。此外，HGF能够直接降低牛和大鼠颗粒细胞中芳香化酶活性的基础水平，并且降低促卵泡激素(follicle stimulating hormone，FSH)对芳香化酶活性的刺激，同时HGF通过抑制LH间接抑制孕酮等类固醇激素的合成，负向调节颗粒细胞的功能[11]。

众多研究表明，MET基因与奶牛的泌乳性能也具有潜在联系。Accornero等[12]通过细胞试验证明，MET的mRNA在牛乳腺上皮细胞系中高表达，且其表达量是同等条件下小鼠乳腺上皮细胞系的2倍。HGF与Met结合后，可激活丝裂原活化蛋白激酶(mitogen-activated protein kinase，MAPK)通路以促进细胞增殖通路和蛋白激酶B(protein kinase B，PKB or Akt)通路进而抑制细胞凋亡[13]，从而诱导牛乳腺上皮细胞的增殖[12]，在三维支架中还可观察到HGF与Met有助于乳腺导管的形成和分支[13]，进一步促进乳腺的发育和成熟。研究表明，在成年牛乳腺上皮细胞周期性的发育、退化和恢复阶段可够检测到Met的表达，而在泌乳阶段不表达；在乳腺肌上皮细胞中，整个生长周期均能检测到Met的表达[2]。此外，Met蛋白激活下游PI3K-AKT-mTOR通路控制脂肪酸和蛋白质的合成、葡萄糖代谢及摄取过程[14-15]，参与乳腺上皮细胞的功能分化和乳汁合成的代谢途径。因此，MET基因在奶牛泌乳系统的发育、乳腺组织周期性退化与恢复过程、乳成分合成中发挥着重要的作用。

综上，MET基因不仅在奶牛卵泡发育过程中发挥着重要的功能，还对奶牛乳腺组织的恢复和发育有着关键的作用。目前，大量对MET基因多态性的研究集中在其对人癌症细胞增殖和分化以及靶向治疗等方面[16-17]，在牛中仅有少量针对生理功能的研究，尚未见对基因多态性及其与生产、功能性状的关联进行挖掘。本研究旨在检测中国荷斯坦牛群体中MET基因的单核苷酸多态(single nucleotide polymorphisms，SNP)，并进一步筛选与繁殖及泌乳性状关联的重要SNPs，为奶牛个体选育策略的制定提供理论依据。

1 材料与方法

1.1 MET基因SNPs的筛查

1.1.1 DNA池的构建 选取70头无亲缘关系的中国荷斯坦牛公牛构建DNA池，冻精样品由北京奶牛中心提供。采用高盐法[18]提取基因组DNA，使用Nanodrop 2000检测浓度后调整为100 ng·μL-1，并随机将样品分成3个DNA池(20、20、30头)，-40 ℃保存备用。

1.1.2 多态位点的筛查 参考MET基因的DNA序列(Ensembl数据库：ENSBTAG00000006161)，共设计24对引物(部分引物序列见表1)，涵盖该基因全部外显子、部分内含子及上下游侧翼区各2 000 bp。

表1 MET基因部分引物信息

采用混池测序法对MET基因进行多态性筛选。PCR反应体系为25 μL：包含DNA混池样(100 ng·μL-1)1 μL，1×CataAmp Taq Plus PCR Mix(北京开拓赛思生物技术有限公司，北京，中国)22 μL，上、下游引物(10 μmol·μL-1)各1 μL。PCR反应条件：预变性(2 min，98 ℃)；变性(10 s，98 ℃)，退火(30 s，温度梯度55～58 ℃)，延伸(10 s，72 ℃)，循环35次；终延伸(1 min，72 ℃)。PCR产物经2%琼脂糖凝胶检测后，送至北京擎科生物科技有限公司测序。最后利用R语言中的“sangerseqR”包对测序结果进行可视化，并计算重叠峰中两个峰的比值，以比值0.1为阈值筛选潜在的多态位点，确定MET基因的19个SNPs。选取SNP中峰值比例高、物理位置间距远、位于基因特殊区域的7个SNPs用于后续大群分型。

1.2 MET基因多态性与繁殖和泌乳性状EBV的关联分析

1.2.1 DNA的提取及基因分型 根据70头公牛的后裔记录，挑选北京地区8个牛场(金银岛、半截河、南三、渠头、长三、长四、金星和中以)1 160头健康荷斯坦母牛为研究对象，使用一次性采血针和抗凝真空采血管采集试验牛尾根静脉血约10 mL，使用天根血液基因组DNA提取试剂盒(DP318，北京天根生化科技有限公司，北京，中国)提取DNA，随后利用Nanodrop 2000测定其浓度和纯度，并通过2%的琼脂糖凝胶电泳鉴定DNA片段完整性后，存放于-40 ℃冰箱中。

对7个候选SNPs设计检测引物和通用引物(表2)，采用竞争性等位基因特异性PCR(kompetitive allele-specific PCR，KASP)进行大群分型。本部分试验委托中玉金标记(北京)生物技术股份有限公司完成。

表2 多态位点竞争性等位基因特异性PCR(KASP)分型引物信息

1.2.2 繁殖与泌乳性状的估计育种值 课题组前期收集北京地区1998—2018年共4 272 839条奶牛生产性能测定(dairy herd improvement，DHI)日记录，利用单性状动物模型进行遗传评估，获得所需繁殖和泌乳性状的EBV[19-20]。繁殖性状包括初产日龄(age at the first calving，AFC)、初配日龄(age at the first service，AFS)、产犊至首次配种间隔(the interval from calving to the first insemination，ICF)、青年牛首末配种间隔(the interval from the first to last insemination in heifers，IFL_H)、经产牛首末配种间隔(the interval from the first to last insemination in cows，IFL_C)共5个性状；泌乳性状包括产奶量(milk yield，MY)、乳脂率(fat percentage，FP)、乳脂量(fat yield，FY)、乳蛋白率(protein percentage，PP)、乳蛋白量(protein yield，PY)和体细胞评分(somatic cell score，SCS)共6个性状，其中FY、PY、SCS 3个性状分别由FP、PY、体细胞数(somatic cell count，SCC，千·mL-1)3个性状转化得到，性状转化公式为：

1.2.3MET基因多态性位点及单倍型与繁殖和泌乳性状估计育种值关联分析 使用Haploview(version 4.2)软件对MET基因的7个SNPs进行连锁不平衡分析，以频率0.05为阈值筛选高度连锁的重要单倍型，获得单倍型块及对应单倍型。采用SAS 9.2软件的GLM过程对SNPs与5个繁殖和6个泌乳性状的EBV进行单位点和双倍型关联分析，模型为：

Yij=μ+Gi+eij

其中，Yij为个体繁殖性状(AFC、AFS、ICF、IFL_H、IFL_C)或泌乳性状(MY、PP、PY、FP、FY、SCS)的估计育种值(estimated breeding values，EBV)；μ为群体均值；Gi为基因型或单倍型效应；eij为随机残差效应。采用Bonferroni法进行多重比较，结果以“最小二乘均值±标准差”表示，P<0.05为显著，P<0.01为极显著。

1.3 生物信息学预测及MET基因表达量验证

1.3.1 生物信息学预测 使用“TFBSTools”R包[21]对MET基因上游突变位点附近潜在结合的转录因子进行预测，通过在线预测网站RNAfold[22-23]对位于外显子的突变位点进行mRNA二级结构预测。

1.3.2MET基因多态与基因表达量的关联分析MET基因的表达量数据来自转录组测序，为课题组前期所积累。采用Trizol法从血液样品中提取总RNA，使用NanoDrop 2000检测RNA的浓度和纯度，而RNA的完整性通过1%琼脂糖凝胶电泳进行检测。OD260 nm/OD280 nm比值为1.8～2.0，电泳条带无明显降解即符合送样条件。构建cDNA文库后采用Illumina Novaseq 6000测序平台进行双端测序。对RNA测序得到的数据进行分析，包括质量控制、比对参考基因组、表达量分析等，最终获得234个样品中基因的表达量数据。随后，匹配这234个样品对应的SNPs位点及双倍型。采用SAS 9.2软件的GLM过程对SNP与MET基因表达量进行单位点和双倍型关联分析，模型为：

Yijkm=μ+Gi+Yj+Pk+Dijkm+Lijkm+Mijkm+eijkm

其中，Yijkm为MET基因表达量；μ为群体均值；Gi为基因型或单倍型效应；Yj为年季效应；Pk为胎次效应；Dijkm为妊娠天数效应；Lijkm为泌乳天数效应；Mijkm为7日平均泌乳量效应；eijkm为随机残差效应。采用Bonferroni法进行多重比较，结果以“最小二乘均值±标准差”表示，P<0.05为显著，P<0.01为极显著。

2 结 果

2.1 中国荷斯坦牛MET基因SNPs分析

2%的琼脂糖凝胶电泳结果显示，24对引物的扩增条带大小均符合预期，其中MET-18、MET-19、MET-20的PCR产物电泳结果见图1。

M.DNA相对分子质量标准；18/19/20表示引物序号，每个引物同时扩增3份产物，其DNA模板分别混池1、混池2、混池3

根据测序峰图的重叠峰情况，本研究共筛选到了19个候选SNPs(表3)，均已被dbSNP数据库所收录，其中有16个SNPs位于内含子，1个SNP位于外显子，2个SNPs位于上游侧翼区。用于后续大群分型的7个SNPs的测序峰图见图2。

表3 MET基因候选多态性位点信息

箭头表示SNPs所在的位置

2.2 群体遗传学分析

通过KASP分型技术获得了7个SNPs在1 160个体中的基因型数据，群体遗传学分析显示(表4)，MET基因7个SNPs最小等位基因频率(minor allele frequency，MAF)的范围为0.184～0.419，其多态性信息含量(polymorphism information content，PIC)均位于0.25～0.5之间，表明这些位点为中度多态；此外，这些SNPs哈代温伯格平衡(H-D平衡)的P值均大于0.05，表明均处于哈代温伯格平衡状态；SNPs的遗传杂合度(heterozygosity，He)低于遗传纯合度(1-He)，说明其在测定群体中变异较小；有效等位基因数(Ne)除g.51736640A>C较小，其余SNPs均大于1.50，说明这些位点的等位基因在测验群体中分布均匀。

表4 中国荷斯坦牛MET基因7个SNPs位点遗传多样性

2.3 MET基因单标记与繁殖和泌乳性状EBV的关联分析

单标记与部分繁殖性状的关联分析及多重比较结果显示(表5)，MET基因的7个SNPs均与繁殖性状IFL_H_EBV关联显著(P<0.05)，位于上游调控区的g.51737889T>C、g.51736640A>C两个SNPs与AFC_EBV显著关联(P<0.05)，与IFL_C_EBV极显著关联(P<0.01)，其余位点关联不显著(P>0.05)；位于上游调控区的g.51737889T>C位点3种基因型AFC_EBV差异极显著(P<0.01)，总体呈现GA>AA>GG的趋势，表明具有等位基因A的个体初产日龄较小，为有利等位基因；同样位于上游调控区的g.51736640A>C位点3种基因型间IFL_C_EBV差异极显著(P<0.01)，总体趋势为TT>GT>GG，与该位点IFL_H_EBV趋势一致，GG基因型经产牛首末配种间隔更短，为有利基因型；位于外显子的g.51660569G>A位点的TT基因型IFL_H_EBV显著高于CC个体(P<0.05)，TC基因型型IFL_H_EBV最低，总体趋势为TT>CC>TC，该位点杂合子TC为有利基因型。

表5 MET基因SNPs多态与部分繁殖性状EBV的关联分析结果

MET基因单标记与部分泌乳性状的关联分析及多重比较结果显示(表6)，除g.51737889T>C外，其余6个SNPs位点均与PP_EBV显著相关(P<0.05)，g.51736640A>C与PP_EBV、FP_EBV、SCS_EBV均呈现显著相关(P<0.05)，g.51737889T>C与SCS显著相关，其余位点关联不显著(P>0.05)；g.51737889T>C位点的AA基因型SCS_EBV极显著高于GG基因型(P<0.01)，GA基因型SCS_EBV显著高于GG基因型(P<0.05)，表现为AA>GA>GG，说明等位基因A导致体细胞评分的升高，属于不利等位基因；g.51736640A>C中GT基因型与GG基因型FP_EBV差异显著(P<0.05)，GG基因型与GT基因型SCS_EBV差异极显著(P<0.01)，趋势表现为GT>TT>GG，GT基因型个体FP_EBV和PP_EBV较低而体细胞评分最高，该GT杂合基因型为不利基因型；g.51665308_51665316 dup中双插入(EE)基因型PP_EBV极显著低于双缺失(FF)基因型(P<0.01)，表现为FF>EF>EE，双插入纯合子(EE)基因型个体PP_EBV最低，为不利基因型；g.51660569G>A中CC基因型与TT、TC基因型PP_EBV差异极显著(P<0.01)，表现为TT>TC>CC，CC基因型个体SCS_EBV最低，但FP_EBV、PP_EBV也最低，为不利基因型。

表6 MET基因SNPs多态与部分泌乳性状EBV的关联分析结果

2.4 MET基因双倍型与繁殖和泌乳性状EBV的关联分析

连锁不平衡分析结果如图3所示，7个SNPs可形成2个单倍型块。第1个单倍型块(BLOCK1)包含5个SNPs，g.51666651G>C至g.51645814C>T具有较紧密的连锁关系，其中g.51666651G>C与g.51665308_51665316dup紧密连锁，D′=1.000，r2=1.000；g.51666651G>C与g.51640656C>A紧密连锁，D′=1.000，r2=1.000；g.51665308_51665316dup与g.51640656C>A紧密连锁，D′=1.000，r2=1.000；g.51660569G>A与g.51645814C>T紧密连锁，D′=1.000，r2=0.998，H1为优势单倍型，频率为0.689；第2个单倍型块(BLOCK2)包含2个SNP，g.51737889T>C与g.51736640A>C具有较紧密的连锁关系，D′=0.995，r2=0.318(r2<0.330)，H4为优势单倍型，频率为0.583。

方块颜色表示连锁程度，颜色越深，连锁程度越高；方块中数值为位点间相关性值的百分数，无数值代表两位点完全连锁

双倍型关联分析结果及多重比较结果显示(表7)，MET基因的BLOCK1仅与繁殖性状IFL_C_EBV关联显著(P<0.05)，与其余繁殖性状关联不显著(P>0.05)；BLOCK2与AFC_EBV关联显著(P<0.05)，与IFL_H_EBV和IFL_C_EBV关联极显著(P<0.001)，与AFS_EBV、ICF_EBV关联不显著(P>0.05)；BLOCK1中H1H1与H1H3两双倍型IFL_C_EBV之间存在显著差异(P<0.05)，总体上呈现H2H2>H3H3>H1H1>H1H2>H1H3>H2H3的趋势，H2H3双倍型个体的IFL_C_EBV最小，为优势双倍型；BLOCK2中IFL_C_EBV总体趋势为H6H6>H5H6>H4H6>H5H5>H4H4>H4H5，H4H4、H4H5与H4H6、H6H6双倍型IFL_H_EBV之间差异极显著(P<0.01)，趋势为H6H6>H4H6>H5H6>H5H5>H4H5>H4H4，H4H4双倍型个体AFC_EBV、IFL_H_EBV均为最低，为优势双倍型。

表7 MET基因单倍型块多态与部分繁殖性状的关联分析结果

不同双倍型的泌乳性状EBV多重比较结果显示(表8)，MET基因的BLOCK1仅与泌乳性状PP_EBV关联极显著(P<0.001)，与其余繁殖性状关联不显著(P>0.05)；BLOCK2与FY_EBV关联显著(P<0.05)，与PP_EBV和SCS_EBV关联极显著(P<0.01)，与MY_EBV、PY_EBV、FP_EBV关联不显著(P>0.05)；BLOCK1中H2H3双倍型PP_EBV与H1H1、H1H2差异极显著(P<0.01)，总体趋势为H1H1>H1H2>H1H3=H3H3>H2H2>H2H3；BLOCK2中H4H6与H5H5双倍型PP_EBV差异极显著(P<0.01)，H5H5与H4H4、H6H6双倍型PP_EBV之间存在显著差异(P<0.05)，总体趋势表现为H6H6>H4H6>H5H6>H4H4>H4H5>H5H5，FY_EBV趋势为H6H6>H5H6>H4H4>H4H6>H4H5>H5H5，H4H4与H4H6双倍型SCS_EBV差异极显著(P<0.01)，H4H4与H5H6双倍型SCS_EBV差异显著(P<0.05)，表现为H5H6>H4H6>H6H6>H5H5>H4H5>H4H4，H6H6双倍型个体具有PP_EBV、FY_EBV最高，SCS_EBV居中，为优势双倍型。

表8 MET基因单倍型块多态与部分泌乳性状EBV的关联分析结果

2.5 生物信息学分析与MET基因表达量验证

2.5.1 基序(motif)分析 g.51737889T>C与g.51736640A>C位点均位于MET基因转录起始位点(TSS)上游约2 000 bp左右，推测可能通过影响转录因子与MET上游序列结合进而影响MET基因的表达，所以进行了motif分析。预测结果表明(图4a和4b)，当g.51737889T>C位点的等位基因为A时，存在72个独有的转录因子，如BARX1和FOXO6等与间质细胞向上皮细胞转化调节的重要转录因子；当等位基因为G时，仅有4种独有的转录因子，如STAT家族和HIC2等与基因激活、新陈代谢相关的转录因子。当g.51736640A>C位点的等位基因为G时，存在9个独有的转录因子，如SIX1、VDR、RUNX3与细胞增殖、分化相关的重要转录因子；当等位基因为T时，有20种转录因子，如NR1I3等与葡萄糖代谢、脂质代谢等相关的重要转录因子。

2.5.2 mRNA结构预测 g.51660569G>A位点位于MET基因第6外显子，推测其可能是通过影响mRNA的二级结构，进而影响基因表达，所以进行了结构预测。预测结果表明(图4c)，当等位基因为C时，mRNA二级结构如左图所示，最小自由能(minimum free energy，MFE)为-165.69 kJ·mol-1，当等位基因为T时，mRNA二级结构如右图所示，MFE为-153.97 kJ·mol-1。因此，当g.51660569G>A位点等位基因为C时，MFE更低，形成的mRNA结构更稳定。

2.5.3MET基因多态与基因表达量的关联分析 关联分析结果及不同基因型的基因表达量多重比较结果显示(图4d)，位于上游调控区的g.51737889T>C、g.51736640A>C两个SNPs与基因表达量关联极显著(P<0.01)，位于内含子的g.51640656C>A与基因表达量关联显著(P<0.05)，其余4个SNPs与基因表达量关联不显著但具有显著趋势(0.05C位点GG与AA基因型个体的基因表达量差异极显著(P<0.01)，总体呈现GG>GA>AA的趋势，具有等位基因G的个体MET基因表达量高；同样位于上游调控区的g.51736640A>C位点GG与GT基因型的基因表达量差异显著(P<0.05)，总体趋势为GG>TT>GT，具有等位基因G的个体MET基因表达量高；位于外显子的g.51660569G>A位点多态与基因表达量存在显著趋势，趋势为TT>CC>TC；位于内含子的g.51640656C>A位点的3种基因型基因表达量差异不显著(P>0.05)，总体趋势为GG>TT>TG，GG基因型个体MET基因表达量最高。

关联分析结果及不同双倍型的基因表达量多重比较结果显示(图4e)，两个单倍型块与基因表达量关联均显著(P<0.05)，但各双倍型之间基因表达量无显著差异(P>0.05)，总体趋势为BLOCK1：H2H2>H1H1>H1H3>H1H2>H2H3，BLOCK2：H5H5>H4H4>H4H5>H4H6>H5H6>H6H6。BLOCK1中H2H2双倍型的个体MET基因表达量高；BLOCK2中H5H5双倍型个体MET基因表达量最高。

a. g.51737889T>C与g.51736640A>C不同等位基因结合转录因子预测结果；b. 部分转录因子示例；c. g.51660569G>A不同等位基因mRNA二级结构预测结果，圈处为多态位点；d. MET基因单位点多态与基因表达量的关联分析结果：*. P<0.05，**. P<0.01，NS. P>0.05；e. MET基因双倍型多态与基因表达量的关联分析结果

3 讨 论

3.1 MET基因SNPs与繁殖、泌乳性状EBV的关联性

本研究通过混池测序法，鉴定到MET基因在中国荷斯坦牛中共有19个SNPs，其中16个SNPs位于内含子，2个SNPs位于上游调控区，1个SNP位于外显子，表明该基因内含子区存在丰富的遗传变异。随后，筛选出7个SNPs进行基因分型，并与繁殖性状、泌乳性状和基因表达量进行关联分析，发现MET基因7个SNPs与中国荷斯坦牛繁殖和泌乳性状存在显著的相关性。研究发现，Met蛋白可以促进颗粒细胞的增殖，影响卵泡发育[11]。研究发现，Met蛋白可激活诱导牛乳腺上皮细胞和导管增殖的通路，促进乳腺的发育和成熟，间接调控葡萄糖摄取和代谢、蛋白质合成等生理过程[13-15]。本研究中，位于上游调控区的g.51737889T>C位点与AFC、IFL_H、IFL_C、SCS等4个性状的EBV具有显著或极显著的相关性，且GG基因型AFC、IFL_H、IFL_C、SCS等性状EBV最低；同样位于上游调控区的g.51736640A>C位点与AFC、IFL_H、IFL_C、PP、FP、SCS等6个性状的EBV具有显著或极显著的相关性，且GG基因型AFC、IFL_H、IFL_C、SCS等性状EBV最低，PP、FP等性状EBV最低；位于外显子的g.51660569G>A位点的TC基因型AFC较低，IFL_H、IFL_C等性状EBV最低。本研究所关注5个繁殖性状的EBV越小，配种时间越短，个体繁殖性能越好，能够保证奶牛养殖的经济效益[24]。对于产奶性状，除SCS外，其余5个性状的EBV越高，个体产奶性能越好。因此，g.51737889T>C位点GG基因型个体同时具有较好的繁殖性能和泌乳性能；g.51736640A>C位点GG基因型个体繁殖性能较好，但泌乳性能相对较差；g.51660569G>A位点TC基因型个体具有较好的繁殖性能。上述3个SNPs位点可作为中国荷斯坦牛繁殖和产奶性状的候选位点，应被重点关注。

此外，位于内含子的g.51640656C>A位点与IFL_H_EBV关联显著，与PP_EBV关联极显著，且与SCS_EBV存在显著关联趋势，该位点TG基因型IFL_H、IFL_C、FP最低，PP、SCS较低，具有该基因型个体的繁殖性能较好，但泌乳性能相对较差。

3.2 MET基因SNPs的生物学预测与表达量关联分析

g.51737889T>C和g.51736640A>C两位点均位于MET基因转录起始位点(transcription start site，TSS)上游约2 000 bp左右，研究表明，基因上游调控区多态性位点通过影响转录因子的结合调控基因表达[25]。在Campbell等[26]的研究中，在人MET基因启动子区(-20 bp)证明存在核心启动子区域，并且能够通过与SP-1、AP-2、PC-4的结合影响启动子活性。为进一步探究这两个突变位点的潜在作用机制，本研究进行了基序分析，结果表明，当g.51737889T>C位点的等位基因为G时，富集的多种转录因子与基因激活和新陈代谢相关；当此位点的等位基因为A时，富集的多种转录因子与间质上皮分化有关。推测，该SNP的GG基因型能够招募与基因激活相关的转录因子结合序列，促进MET基因的表达，最终导致该基因型个体繁殖和泌乳性能表现优异。该位点与MET基因表达量的关联分析结果证明，g.51737889T>C位点GG基因型个体的MET基因表达量最高，且与其他基因型存在显著差异。同样，当g.51736640A>C位点的等位基因为G时，富集的多种转录因子与细胞增殖、分化相关；当此位点的等位基因为T时，富集的多种转录因子与葡萄糖代谢、脂质代谢有关。推测，该SNP的GG基因型能够招募与细胞分化相关的转录因子结合序列，促进MET基因的表达，最终导致该基因型个体繁殖性能表现优异。后续的分析中，g.51736640A>C位点GG基因型个体的MET基因表达量最高，且与其他基因型存在显著差异，更进一步证明了此猜测。

g.51660569G>A为同义突变位点，参与编码Met蛋白第598位天冬氨酸，位于胞内近膜区。有研究表明，同义突变影响核糖体运输、蛋白质折叠以及细胞适应过程[27-28]，所以对该位点进行了mRNA结构预测。结果表明，g.51660569G>A位点不同基因型mRNA二级结构无明显差异，但等位基因为C时，mRNA的MFE更低，结构更加稳定。推测是该位点的同义突变通过改变mRNA的二级结构，进而影响mRNA稳定性，延长半衰期，提高基因表达量，使不同基因型个体MET基因表达量存在差异，造成了其在卵泡和乳腺发育过程的差异。但该位点多态与基因表达量的关联不显著，预测结果并不能完全解释杂合子基因型个体繁殖和泌乳性能表现优异。基因的表达存在时空特异性[29]，研究表明，Met的激活通过自分泌或旁分泌的方式启动[17]，且Met主要通过改变蛋白活性影响表型，而不是基因表达[30]。部分研究也表明，在乳腺发育阶段，Met蛋白高度表达，mRNA表达量很低，在哺乳期MetmRNA快速降解[2,31-32]。此次RNA提取样品来自于泌乳牛尾根血液，并非卵巢或乳腺组织，这可能是表达量检测结果与预测结果不一致的原因。此外，研究表明，内含子SNP能够通过影响可变剪切来影响基因表达[33-34]，或通过形成lncRNA调控基因表达[35]；g.51640656C>A位点多态与基因表达量关联显著，TG基因型表达量最低。但该SNP并非已证实的可变剪切位点，推测该位点对繁殖和产奶性状的影响可能与内含子形成的lncRNA有关，有待进一步验证。

3.3 MET基因单倍型与繁殖、泌乳性状EBV和基因表达量的关联分析

研究表明，基因对表型的影响可能受到多个突变位点的共同作用[29,36]。本研究中，在中国荷斯坦牛群体中MET基因的7个SNPs共形成2个单倍型块，g.51666651G>C、g.51665308_51665316 dup、g.51660569G>A、g.51640656C>A、g.51645814C>T共5个SNPs紧密连锁形成BLOCK1，g.51737889T>C与g.51736640A>C两个SNPs紧密连锁形成BLOCK2。BLOCK2与AFC、IFL_H、IFL_C、PP、FY、SCS多个繁殖、泌乳性状EBV关联显著，且与PY、FP两个泌乳性状EBV关联存在显著趋势。BLOCK2的2个SNPs g.51737889T>C与g.51736640A>C均位于上游调控区，这2个SNPs位点各自与多个繁殖、泌乳性状EBV显著相关。

另外，一个单倍型块上的多个突变会对性状产生更大的影响[37]。H4H4正是g.51737889T>C与g.51736640A>C两个SNPs的优势基因型GG基因型的组合，H4H5是g.51737889T>C的杂合GA基因型与g.51736640A>C的GG优势基因型的组合；H4H4与H4H5个体排在AFC、IFL_H、IFL_C 3种繁殖性状EBV最低的前三名；H4H4的PP_EBV第三高，FY_EBV第二高，SCS_EBV最低。H4H4双倍型个体具有较好的繁殖和泌乳性状。且与基因表达量关联分析结果也证明，H4H4双倍型个体基因表达量最高，充分验证了g.51737889T>C与g.51736640A>C两位点等位基因G对中国荷斯坦牛繁殖及泌乳性能的影响。

4 结 论

本研究获得了中国荷斯坦牛MET基因的多态性图谱，MET基因的7个SNPs及其形成的2个单倍型块与繁殖、泌乳性状存在显著或极显著关联。位点g.51737889T>C和g.51736640A>C可能通过改变转录因子结合影响基因的表达，对MET的表达量有显著影响；位点g.51660569G>A对mRNA的二级结构无显著影响但不同等位基因对最小自由能产生影响。上述3个位点可作为奶牛繁殖和泌乳性状的候选位点重点关注，研究结果为中国荷斯坦牛高产高效选育提供了可用的遗传标记。