MET基因多態(tài)性與中國荷斯坦牛繁殖和泌乳性狀的關(guān)聯(lián)分析

許靜漪，徐昊祺，胡麗蓉，張 帆，高 清，羅漢鵬，張海亮，師 睿，李 想，劉 林，郭 剛，王雅春*

(1.中國農(nóng)業(yè)大學(xué)動物科學(xué)技術(shù)學(xué)院 農(nóng)業(yè)農(nóng)村部動物遺傳育種與繁殖(家畜)重點(diǎn)實(shí)驗(yàn)室，畜禽育種國家工程實(shí)驗(yàn)室，北京100193； 2.新疆農(nóng)業(yè)大學(xué)動物科學(xué)學(xué)院，烏魯木齊 830052；3.北京奶牛中心，北京 100194； 4.北京首農(nóng)畜牧發(fā)展有限公司，北京 100029)

奶業(yè)是健康中國、強(qiáng)壯民族不可或缺的產(chǎn)業(yè)，中國不僅是世界牛奶產(chǎn)量第三大經(jīng)濟(jì)體，而且也已經(jīng)成長為世界第一大乳制品進(jìn)口貿(mào)易國[1]。為了滿足消費(fèi)者對牛奶及乳制品的需求，在過去的幾十年間，奶牛育種者們對奶牛的產(chǎn)奶性能進(jìn)行了大幅度的選育，但也導(dǎo)致了繁殖力、抗病力的下降，嚴(yán)重阻礙了奶業(yè)的可持續(xù)發(fā)展。隨著平衡育種理念廣泛應(yīng)用于實(shí)際選育，越來越多的性狀，如繁殖、長壽等功能性狀得到了育種工作者和牧場生產(chǎn)經(jīng)營者的關(guān)注。

本課題組前期對中國荷斯坦牛的多項(xiàng)繁殖性狀進(jìn)行了全基因組關(guān)聯(lián)分析(genome-wide association study，GWAS)，獲得了一系列的候選基因。通過蛋白互作網(wǎng)絡(luò)(protein-protein interaction networks，PPI)分析，發(fā)現(xiàn)MET(MET proto-oncogene, receptor tyrosine kinase)基因與多個(gè)候選基因存在互作關(guān)系，表明該基因可能是影響中國荷斯坦牛繁殖性狀的重要基因。MET是一種原癌基因，又稱c-MET。牛的MET基因包含21個(gè)外顯子和20個(gè)內(nèi)含子，其cDNA全長4.8 kb，包含4 152 bp的開放閱讀框，編碼由1 384個(gè)氨基酸組成的Met蛋白[2]。Met蛋白是一種由二硫鍵連接的異二聚體組成的酪氨酸激酶，由190 ku的前體裂解為成熟的胞外α鏈和較長的跨膜β鏈[3]。已有研究證明，肝細(xì)胞生長因子(hepatocyte growth factor，HGF)為Met的唯一配體[4-5]。Met與HGF結(jié)合能夠激活多種信號轉(zhuǎn)導(dǎo)通路，誘導(dǎo)多種生物反應(yīng)，如增殖、凋亡抑制、運(yùn)動、細(xì)胞外基質(zhì)侵襲和血管生成的激活[6-10]。Parrott和Skinner[11]研究發(fā)現(xiàn)，MET與牛卵泡發(fā)育有關(guān)，HGF與Met在牛卵泡發(fā)育過程中高度表達(dá)，HGF由卵泡膜細(xì)胞表達(dá)，作用于卵泡膜細(xì)胞和顆粒細(xì)胞上的Met，雌激素和促黃體生成素(luteinizing hormone，LH)通過調(diào)節(jié)局部HGF的產(chǎn)生間接控制卵泡發(fā)育，促進(jìn)顆粒細(xì)胞的增殖，從而促進(jìn)卵泡生長。此外，HGF能夠直接降低牛和大鼠顆粒細(xì)胞中芳香化酶活性的基礎(chǔ)水平，并且降低促卵泡激素(follicle stimulating hormone，F(xiàn)SH)對芳香化酶活性的刺激，同時(shí)HGF通過抑制LH間接抑制孕酮等類固醇激素的合成，負(fù)向調(diào)節(jié)顆粒細(xì)胞的功能[11]。

眾多研究表明，MET基因與奶牛的泌乳性能也具有潛在聯(lián)系。Accornero等[12]通過細(xì)胞試驗(yàn)證明，MET的mRNA在牛乳腺上皮細(xì)胞系中高表達(dá)，且其表達(dá)量是同等條件下小鼠乳腺上皮細(xì)胞系的2倍。HGF與Met結(jié)合后，可激活絲裂原活化蛋白激酶(mitogen-activated protein kinase，MAPK)通路以促進(jìn)細(xì)胞增殖通路和蛋白激酶B(protein kinase B，PKB or Akt)通路進(jìn)而抑制細(xì)胞凋亡[13]，從而誘導(dǎo)牛乳腺上皮細(xì)胞的增殖[12]，在三維支架中還可觀察到HGF與Met有助于乳腺導(dǎo)管的形成和分支[13]，進(jìn)一步促進(jìn)乳腺的發(fā)育和成熟。研究表明，在成年牛乳腺上皮細(xì)胞周期性的發(fā)育、退化和恢復(fù)階段可夠檢測到Met的表達(dá)，而在泌乳階段不表達(dá)；在乳腺肌上皮細(xì)胞中，整個(gè)生長周期均能檢測到Met的表達(dá)[2]。此外，Met蛋白激活下游PI3K-AKT-mTOR通路控制脂肪酸和蛋白質(zhì)的合成、葡萄糖代謝及攝取過程[14-15]，參與乳腺上皮細(xì)胞的功能分化和乳汁合成的代謝途徑。因此，MET基因在奶牛泌乳系統(tǒng)的發(fā)育、乳腺組織周期性退化與恢復(fù)過程、乳成分合成中發(fā)揮著重要的作用。

綜上，MET基因不僅在奶牛卵泡發(fā)育過程中發(fā)揮著重要的功能，還對奶牛乳腺組織的恢復(fù)和發(fā)育有著關(guān)鍵的作用。目前，大量對MET基因多態(tài)性的研究集中在其對人癌癥細(xì)胞增殖和分化以及靶向治療等方面[16-17]，在牛中僅有少量針對生理功能的研究，尚未見對基因多態(tài)性及其與生產(chǎn)、功能性狀的關(guān)聯(lián)進(jìn)行挖掘。本研究旨在檢測中國荷斯坦牛群體中MET基因的單核苷酸多態(tài)(single nucleotide polymorphisms，SNP)，并進(jìn)一步篩選與繁殖及泌乳性狀關(guān)聯(lián)的重要SNPs，為奶牛個(gè)體選育策略的制定提供理論依據(jù)。

1 材料與方法

1.1 MET基因SNPs的篩查

1.1.1 DNA池的構(gòu)建 選取70頭無親緣關(guān)系的中國荷斯坦牛公牛構(gòu)建DNA池，凍精樣品由北京奶牛中心提供。采用高鹽法[18]提取基因組DNA，使用Nanodrop 2000檢測濃度后調(diào)整為100 ng·μL-1，并隨機(jī)將樣品分成3個(gè)DNA池(20、20、30頭)，-40 ℃保存?zhèn)溆谩?/p>

1.1.2 多態(tài)位點(diǎn)的篩查 參考MET基因的DNA序列(Ensembl數(shù)據(jù)庫：ENSBTAG00000006161)，共設(shè)計(jì)24對引物(部分引物序列見表1)，涵蓋該基因全部外顯子、部分內(nèi)含子及上下游側(cè)翼區(qū)各2 000 bp。

表1 MET基因部分引物信息

采用混池測序法對MET基因進(jìn)行多態(tài)性篩選。PCR反應(yīng)體系為25 μL：包含DNA混池樣(100 ng·μL-1)1 μL，1×CataAmp Taq Plus PCR Mix(北京開拓賽思生物技術(shù)有限公司，北京，中國)22 μL，上、下游引物(10 μmol·μL-1)各1 μL。PCR反應(yīng)條件：預(yù)變性(2 min，98 ℃)；變性(10 s，98 ℃)，退火(30 s，溫度梯度55～58 ℃)，延伸(10 s，72 ℃)，循環(huán)35次；終延伸(1 min，72 ℃)。PCR產(chǎn)物經(jīng)2%瓊脂糖凝膠檢測后，送至北京擎科生物科技有限公司測序。最后利用R語言中的“sangerseqR”包對測序結(jié)果進(jìn)行可視化，并計(jì)算重疊峰中兩個(gè)峰的比值，以比值0.1為閾值篩選潛在的多態(tài)位點(diǎn)，確定MET基因的19個(gè)SNPs。選取SNP中峰值比例高、物理位置間距遠(yuǎn)、位于基因特殊區(qū)域的7個(gè)SNPs用于后續(xù)大群分型。

1.2 MET基因多態(tài)性與繁殖和泌乳性狀EBV的關(guān)聯(lián)分析

1.2.1 DNA的提取及基因分型 根據(jù)70頭公牛的后裔記錄，挑選北京地區(qū)8個(gè)牛場(金銀島、半截河、南三、渠頭、長三、長四、金星和中以)1 160頭健康荷斯坦母牛為研究對象，使用一次性采血針和抗凝真空采血管采集試驗(yàn)牛尾根靜脈血約10 mL，使用天根血液基因組DNA提取試劑盒(DP318，北京天根生化科技有限公司，北京，中國)提取DNA，隨后利用Nanodrop 2000測定其濃度和純度，并通過2%的瓊脂糖凝膠電泳鑒定DNA片段完整性后，存放于-40 ℃冰箱中。

對7個(gè)候選SNPs設(shè)計(jì)檢測引物和通用引物(表2)，采用競爭性等位基因特異性PCR(kompetitive allele-specific PCR，KASP)進(jìn)行大群分型。本部分試驗(yàn)委托中玉金標(biāo)記(北京)生物技術(shù)股份有限公司完成。

表2 多態(tài)位點(diǎn)競爭性等位基因特異性PCR(KASP)分型引物信息

1.2.2 繁殖與泌乳性狀的估計(jì)育種值 課題組前期收集北京地區(qū)1998—2018年共4 272 839條奶牛生產(chǎn)性能測定(dairy herd improvement，DHI)日記錄，利用單性狀動物模型進(jìn)行遺傳評估，獲得所需繁殖和泌乳性狀的EBV[19-20]。繁殖性狀包括初產(chǎn)日齡(age at the first calving，AFC)、初配日齡(age at the first service，AFS)、產(chǎn)犢至首次配種間隔(the interval from calving to the first insemination，ICF)、青年牛首末配種間隔(the interval from the first to last insemination in heifers，IFL_H)、經(jīng)產(chǎn)牛首末配種間隔(the interval from the first to last insemination in cows，IFL_C)共5個(gè)性狀；泌乳性狀包括產(chǎn)奶量(milk yield，MY)、乳脂率(fat percentage，F(xiàn)P)、乳脂量(fat yield，F(xiàn)Y)、乳蛋白率(protein percentage，PP)、乳蛋白量(protein yield，PY)和體細(xì)胞評分(somatic cell score，SCS)共6個(gè)性狀，其中FY、PY、SCS 3個(gè)性狀分別由FP、PY、體細(xì)胞數(shù)(somatic cell count，SCC，千·mL-1)3個(gè)性狀轉(zhuǎn)化得到，性狀轉(zhuǎn)化公式為：

1.2.3MET基因多態(tài)性位點(diǎn)及單倍型與繁殖和泌乳性狀估計(jì)育種值關(guān)聯(lián)分析 使用Haploview(version 4.2)軟件對MET基因的7個(gè)SNPs進(jìn)行連鎖不平衡分析，以頻率0.05為閾值篩選高度連鎖的重要單倍型，獲得單倍型塊及對應(yīng)單倍型。采用SAS 9.2軟件的GLM過程對SNPs與5個(gè)繁殖和6個(gè)泌乳性狀的EBV進(jìn)行單位點(diǎn)和雙倍型關(guān)聯(lián)分析，模型為：

Yij=μ+Gi+eij

其中，Yij為個(gè)體繁殖性狀(AFC、AFS、ICF、IFL_H、IFL_C)或泌乳性狀(MY、PP、PY、FP、FY、SCS)的估計(jì)育種值(estimated breeding values，EBV)；μ為群體均值；Gi為基因型或單倍型效應(yīng)；eij為隨機(jī)殘差效應(yīng)。采用Bonferroni法進(jìn)行多重比較，結(jié)果以“最小二乘均值±標(biāo)準(zhǔn)差”表示，P<0.05為顯著，P<0.01為極顯著。

1.3 生物信息學(xué)預(yù)測及MET基因表達(dá)量驗(yàn)證

1.3.1 生物信息學(xué)預(yù)測 使用“TFBSTools”R包[21]對MET基因上游突變位點(diǎn)附近潛在結(jié)合的轉(zhuǎn)錄因子進(jìn)行預(yù)測，通過在線預(yù)測網(wǎng)站RNAfold[22-23]對位于外顯子的突變位點(diǎn)進(jìn)行mRNA二級結(jié)構(gòu)預(yù)測。

1.3.2MET基因多態(tài)與基因表達(dá)量的關(guān)聯(lián)分析MET基因的表達(dá)量數(shù)據(jù)來自轉(zhuǎn)錄組測序，為課題組前期所積累。采用Trizol法從血液樣品中提取總RNA，使用NanoDrop 2000檢測RNA的濃度和純度，而RNA的完整性通過1%瓊脂糖凝膠電泳進(jìn)行檢測。OD260 nm/OD280 nm比值為1.8～2.0，電泳條帶無明顯降解即符合送樣條件。構(gòu)建cDNA文庫后采用Illumina Novaseq 6000測序平臺進(jìn)行雙端測序。對RNA測序得到的數(shù)據(jù)進(jìn)行分析，包括質(zhì)量控制、比對參考基因組、表達(dá)量分析等，最終獲得234個(gè)樣品中基因的表達(dá)量數(shù)據(jù)。隨后，匹配這234個(gè)樣品對應(yīng)的SNPs位點(diǎn)及雙倍型。采用SAS 9.2軟件的GLM過程對SNP與MET基因表達(dá)量進(jìn)行單位點(diǎn)和雙倍型關(guān)聯(lián)分析，模型為：

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

其中，Yijkm為MET基因表達(dá)量；μ為群體均值；Gi為基因型或單倍型效應(yīng)；Yj為年季效應(yīng)；Pk為胎次效應(yīng)；Dijkm為妊娠天數(shù)效應(yīng)；Lijkm為泌乳天數(shù)效應(yīng)；Mijkm為7日平均泌乳量效應(yīng)；eijkm為隨機(jī)殘差效應(yīng)。采用Bonferroni法進(jìn)行多重比較，結(jié)果以“最小二乘均值±標(biāo)準(zhǔn)差”表示，P<0.05為顯著，P<0.01為極顯著。

2 結(jié) 果

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

2%的瓊脂糖凝膠電泳結(jié)果顯示，24對引物的擴(kuò)增條帶大小均符合預(yù)期，其中MET-18、MET-19、MET-20的PCR產(chǎn)物電泳結(jié)果見圖1。

M.DNA相對分子質(zhì)量標(biāo)準(zhǔn)；18/19/20表示引物序號，每個(gè)引物同時(shí)擴(kuò)增3份產(chǎn)物，其DNA模板分別混池1、混池2、混池3

根據(jù)測序峰圖的重疊峰情況，本研究共篩選到了19個(gè)候選SNPs(表3)，均已被dbSNP數(shù)據(jù)庫所收錄，其中有16個(gè)SNPs位于內(nèi)含子，1個(gè)SNP位于外顯子，2個(gè)SNPs位于上游側(cè)翼區(qū)。用于后續(xù)大群分型的7個(gè)SNPs的測序峰圖見圖2。

表3 MET基因候選多態(tài)性位點(diǎn)信息

箭頭表示SNPs所在的位置

2.2 群體遺傳學(xué)分析

通過KASP分型技術(shù)獲得了7個(gè)SNPs在1 160個(gè)體中的基因型數(shù)據(jù)，群體遺傳學(xué)分析顯示(表4)，MET基因7個(gè)SNPs最小等位基因頻率(minor allele frequency，MAF)的范圍為0.184～0.419，其多態(tài)性信息含量(polymorphism information content，PIC)均位于0.25～0.5之間，表明這些位點(diǎn)為中度多態(tài)；此外，這些SNPs哈代溫伯格平衡(H-D平衡)的P值均大于0.05，表明均處于哈代溫伯格平衡狀態(tài)；SNPs的遺傳雜合度(heterozygosity，He)低于遺傳純合度(1-He)，說明其在測定群體中變異較?。挥行У任换驍?shù)(Ne)除g.51736640A>C較小，其余SNPs均大于1.50，說明這些位點(diǎn)的等位基因在測驗(yàn)群體中分布均勻。

表4 中國荷斯坦牛MET基因7個(gè)SNPs位點(diǎn)遺傳多樣性

2.3 MET基因單標(biāo)記與繁殖和泌乳性狀EBV的關(guān)聯(lián)分析

單標(biāo)記與部分繁殖性狀的關(guān)聯(lián)分析及多重比較結(jié)果顯示(表5)，MET基因的7個(gè)SNPs均與繁殖性狀I(lǐng)FL_H_EBV關(guān)聯(lián)顯著(P<0.05)，位于上游調(diào)控區(qū)的g.51737889T>C、g.51736640A>C兩個(gè)SNPs與AFC_EBV顯著關(guān)聯(lián)(P<0.05)，與IFL_C_EBV極顯著關(guān)聯(lián)(P<0.01)，其余位點(diǎn)關(guān)聯(lián)不顯著(P>0.05)；位于上游調(diào)控區(qū)的g.51737889T>C位點(diǎn)3種基因型AFC_EBV差異極顯著(P<0.01)，總體呈現(xiàn)GA>AA>GG的趨勢，表明具有等位基因A的個(gè)體初產(chǎn)日齡較小，為有利等位基因；同樣位于上游調(diào)控區(qū)的g.51736640A>C位點(diǎn)3種基因型間IFL_C_EBV差異極顯著(P<0.01)，總體趨勢為TT>GT>GG，與該位點(diǎn)IFL_H_EBV趨勢一致，GG基因型經(jīng)產(chǎn)牛首末配種間隔更短，為有利基因型；位于外顯子的g.51660569G>A位點(diǎn)的TT基因型IFL_H_EBV顯著高于CC個(gè)體(P<0.05)，TC基因型型IFL_H_EBV最低，總體趨勢為TT>CC>TC，該位點(diǎn)雜合子TC為有利基因型。

表5 MET基因SNPs多態(tài)與部分繁殖性狀EBV的關(guān)聯(lián)分析結(jié)果

MET基因單標(biāo)記與部分泌乳性狀的關(guān)聯(lián)分析及多重比較結(jié)果顯示(表6)，除g.51737889T>C外，其余6個(gè)SNPs位點(diǎn)均與PP_EBV顯著相關(guān)(P<0.05)，g.51736640A>C與PP_EBV、FP_EBV、SCS_EBV均呈現(xiàn)顯著相關(guān)(P<0.05)，g.51737889T>C與SCS顯著相關(guān)，其余位點(diǎn)關(guān)聯(lián)不顯著(P>0.05)；g.51737889T>C位點(diǎn)的AA基因型SCS_EBV極顯著高于GG基因型(P<0.01)，GA基因型SCS_EBV顯著高于GG基因型(P<0.05)，表現(xiàn)為AA>GA>GG，說明等位基因A導(dǎo)致體細(xì)胞評分的升高，屬于不利等位基因；g.51736640A>C中GT基因型與GG基因型FP_EBV差異顯著(P<0.05)，GG基因型與GT基因型SCS_EBV差異極顯著(P<0.01)，趨勢表現(xiàn)為GT>TT>GG，GT基因型個(gè)體FP_EBV和PP_EBV較低而體細(xì)胞評分最高，該GT雜合基因型為不利基因型；g.51665308_51665316 dup中雙插入(EE)基因型PP_EBV極顯著低于雙缺失(FF)基因型(P<0.01)，表現(xiàn)為FF>EF>EE，雙插入純合子(EE)基因型個(gè)體PP_EBV最低，為不利基因型；g.51660569G>A中CC基因型與TT、TC基因型PP_EBV差異極顯著(P<0.01)，表現(xiàn)為TT>TC>CC，CC基因型個(gè)體SCS_EBV最低，但FP_EBV、PP_EBV也最低，為不利基因型。

表6 MET基因SNPs多態(tài)與部分泌乳性狀EBV的關(guān)聯(lián)分析結(jié)果

2.4 MET基因雙倍型與繁殖和泌乳性狀EBV的關(guān)聯(lián)分析

連鎖不平衡分析結(jié)果如圖3所示，7個(gè)SNPs可形成2個(gè)單倍型塊。第1個(gè)單倍型塊(BLOCK1)包含5個(gè)SNPs，g.51666651G>C至g.51645814C>T具有較緊密的連鎖關(guān)系，其中g(shù).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為優(yōu)勢單倍型，頻率為0.689；第2個(gè)單倍型塊(BLOCK2)包含2個(gè)SNP，g.51737889T>C與g.51736640A>C具有較緊密的連鎖關(guān)系，D′=0.995，r2=0.318(r2<0.330)，H4為優(yōu)勢單倍型，頻率為0.583。

方塊顏色表示連鎖程度，顏色越深，連鎖程度越高；方塊中數(shù)值為位點(diǎn)間相關(guān)性值的百分?jǐn)?shù)，無數(shù)值代表兩位點(diǎn)完全連鎖

雙倍型關(guān)聯(lián)分析結(jié)果及多重比較結(jié)果顯示(表7)，MET基因的BLOCK1僅與繁殖性狀I(lǐng)FL_C_EBV關(guān)聯(lián)顯著(P<0.05)，與其余繁殖性狀關(guān)聯(lián)不顯著(P>0.05)；BLOCK2與AFC_EBV關(guān)聯(lián)顯著(P<0.05)，與IFL_H_EBV和IFL_C_EBV關(guān)聯(lián)極顯著(P<0.001)，與AFS_EBV、ICF_EBV關(guān)聯(lián)不顯著(P>0.05)；BLOCK1中H1H1與H1H3兩雙倍型IFL_C_EBV之間存在顯著差異(P<0.05)，總體上呈現(xiàn)H2H2>H3H3>H1H1>H1H2>H1H3>H2H3的趨勢，H2H3雙倍型個(gè)體的IFL_C_EBV最小，為優(yōu)勢雙倍型；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雙倍型個(gè)體AFC_EBV、IFL_H_EBV均為最低，為優(yōu)勢雙倍型。

表7 MET基因單倍型塊多態(tài)與部分繁殖性狀的關(guān)聯(lián)分析結(jié)果

不同雙倍型的泌乳性狀EBV多重比較結(jié)果顯示(表8)，MET基因的BLOCK1僅與泌乳性狀PP_EBV關(guān)聯(lián)極顯著(P<0.001)，與其余繁殖性狀關(guān)聯(lián)不顯著(P>0.05)；BLOCK2與FY_EBV關(guān)聯(lián)顯著(P<0.05)，與PP_EBV和SCS_EBV關(guān)聯(lián)極顯著(P<0.01)，與MY_EBV、PY_EBV、FP_EBV關(guān)聯(lián)不顯著(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)，總體趨勢表現(xiàn)為H6H6>H4H6>H5H6>H4H4>H4H5>H5H5，F(xiàn)Y_EBV趨勢為H6H6>H5H6>H4H4>H4H6>H4H5>H5H5，H4H4與H4H6雙倍型SCS_EBV差異極顯著(P<0.01)，H4H4與H5H6雙倍型SCS_EBV差異顯著(P<0.05)，表現(xiàn)為H5H6>H4H6>H6H6>H5H5>H4H5>H4H4，H6H6雙倍型個(gè)體具有PP_EBV、FY_EBV最高，SCS_EBV居中，為優(yōu)勢雙倍型。

表8 MET基因單倍型塊多態(tài)與部分泌乳性狀EBV的關(guān)聯(lián)分析結(jié)果

2.5 生物信息學(xué)分析與MET基因表達(dá)量驗(yàn)證

2.5.1 基序(motif)分析 g.51737889T>C與g.51736640A>C位點(diǎn)均位于MET基因轉(zhuǎn)錄起始位點(diǎn)(TSS)上游約2 000 bp左右，推測可能通過影響轉(zhuǎn)錄因子與MET上游序列結(jié)合進(jìn)而影響MET基因的表達(dá)，所以進(jìn)行了motif分析。預(yù)測結(jié)果表明(圖4a和4b)，當(dāng)g.51737889T>C位點(diǎn)的等位基因?yàn)锳時(shí)，存在72個(gè)獨(dú)有的轉(zhuǎn)錄因子，如BARX1和FOXO6等與間質(zhì)細(xì)胞向上皮細(xì)胞轉(zhuǎn)化調(diào)節(jié)的重要轉(zhuǎn)錄因子；當(dāng)?shù)任换驗(yàn)镚時(shí)，僅有4種獨(dú)有的轉(zhuǎn)錄因子，如STAT家族和HIC2等與基因激活、新陳代謝相關(guān)的轉(zhuǎn)錄因子。當(dāng)g.51736640A>C位點(diǎn)的等位基因?yàn)镚時(shí)，存在9個(gè)獨(dú)有的轉(zhuǎn)錄因子，如SIX1、VDR、RUNX3與細(xì)胞增殖、分化相關(guān)的重要轉(zhuǎn)錄因子；當(dāng)?shù)任换驗(yàn)門時(shí)，有20種轉(zhuǎn)錄因子，如NR1I3等與葡萄糖代謝、脂質(zhì)代謝等相關(guān)的重要轉(zhuǎn)錄因子。

2.5.2 mRNA結(jié)構(gòu)預(yù)測 g.51660569G>A位點(diǎn)位于MET基因第6外顯子，推測其可能是通過影響mRNA的二級結(jié)構(gòu)，進(jìn)而影響基因表達(dá)，所以進(jìn)行了結(jié)構(gòu)預(yù)測。預(yù)測結(jié)果表明(圖4c)，當(dāng)?shù)任换驗(yàn)镃時(shí)，mRNA二級結(jié)構(gòu)如左圖所示，最小自由能(minimum free energy，MFE)為-165.69 kJ·mol-1，當(dāng)?shù)任换驗(yàn)門時(shí)，mRNA二級結(jié)構(gòu)如右圖所示，MFE為-153.97 kJ·mol-1。因此，當(dāng)g.51660569G>A位點(diǎn)等位基因?yàn)镃時(shí)，MFE更低，形成的mRNA結(jié)構(gòu)更穩(wěn)定。

2.5.3MET基因多態(tài)與基因表達(dá)量的關(guān)聯(lián)分析 關(guān)聯(lián)分析結(jié)果及不同基因型的基因表達(dá)量多重比較結(jié)果顯示(圖4d)，位于上游調(diào)控區(qū)的g.51737889T>C、g.51736640A>C兩個(gè)SNPs與基因表達(dá)量關(guān)聯(lián)極顯著(P<0.01)，位于內(nèi)含子的g.51640656C>A與基因表達(dá)量關(guān)聯(lián)顯著(P<0.05)，其余4個(gè)SNPs與基因表達(dá)量關(guān)聯(lián)不顯著但具有顯著趨勢(0.05C位點(diǎn)GG與AA基因型個(gè)體的基因表達(dá)量差異極顯著(P<0.01)，總體呈現(xiàn)GG>GA>AA的趨勢，具有等位基因G的個(gè)體MET基因表達(dá)量高；同樣位于上游調(diào)控區(qū)的g.51736640A>C位點(diǎn)GG與GT基因型的基因表達(dá)量差異顯著(P<0.05)，總體趨勢為GG>TT>GT，具有等位基因G的個(gè)體MET基因表達(dá)量高；位于外顯子的g.51660569G>A位點(diǎn)多態(tài)與基因表達(dá)量存在顯著趨勢，趨勢為TT>CC>TC；位于內(nèi)含子的g.51640656C>A位點(diǎn)的3種基因型基因表達(dá)量差異不顯著(P>0.05)，總體趨勢為GG>TT>TG，GG基因型個(gè)體MET基因表達(dá)量最高。

關(guān)聯(lián)分析結(jié)果及不同雙倍型的基因表達(dá)量多重比較結(jié)果顯示(圖4e)，兩個(gè)單倍型塊與基因表達(dá)量關(guān)聯(lián)均顯著(P<0.05)，但各雙倍型之間基因表達(dá)量無顯著差異(P>0.05)，總體趨勢為BLOCK1：H2H2>H1H1>H1H3>H1H2>H2H3，BLOCK2：H5H5>H4H4>H4H5>H4H6>H5H6>H6H6。BLOCK1中H2H2雙倍型的個(gè)體MET基因表達(dá)量高；BLOCK2中H5H5雙倍型個(gè)體MET基因表達(dá)量最高。

a. g.51737889T>C與g.51736640A>C不同等位基因結(jié)合轉(zhuǎn)錄因子預(yù)測結(jié)果；b. 部分轉(zhuǎn)錄因子示例；c. g.51660569G>A不同等位基因mRNA二級結(jié)構(gòu)預(yù)測結(jié)果，圈處為多態(tài)位點(diǎn)；d. MET基因單位點(diǎn)多態(tài)與基因表達(dá)量的關(guān)聯(lián)分析結(jié)果：*. P<0.05，**. P<0.01，NS. P>0.05；e. MET基因雙倍型多態(tài)與基因表達(dá)量的關(guān)聯(lián)分析結(jié)果

3 討 論

3.1 MET基因SNPs與繁殖、泌乳性狀EBV的關(guān)聯(lián)性

本研究通過混池測序法，鑒定到MET基因在中國荷斯坦牛中共有19個(gè)SNPs，其中16個(gè)SNPs位于內(nèi)含子，2個(gè)SNPs位于上游調(diào)控區(qū)，1個(gè)SNP位于外顯子，表明該基因內(nèi)含子區(qū)存在豐富的遺傳變異。隨后，篩選出7個(gè)SNPs進(jìn)行基因分型，并與繁殖性狀、泌乳性狀和基因表達(dá)量進(jìn)行關(guān)聯(lián)分析，發(fā)現(xiàn)MET基因7個(gè)SNPs與中國荷斯坦牛繁殖和泌乳性狀存在顯著的相關(guān)性。研究發(fā)現(xiàn)，Met蛋白可以促進(jìn)顆粒細(xì)胞的增殖，影響卵泡發(fā)育[11]。研究發(fā)現(xiàn)，Met蛋白可激活誘導(dǎo)牛乳腺上皮細(xì)胞和導(dǎo)管增殖的通路，促進(jìn)乳腺的發(fā)育和成熟，間接調(diào)控葡萄糖攝取和代謝、蛋白質(zhì)合成等生理過程[13-15]。本研究中，位于上游調(diào)控區(qū)的g.51737889T>C位點(diǎn)與AFC、IFL_H、IFL_C、SCS等4個(gè)性狀的EBV具有顯著或極顯著的相關(guān)性，且GG基因型AFC、IFL_H、IFL_C、SCS等性狀EBV最低；同樣位于上游調(diào)控區(qū)的g.51736640A>C位點(diǎn)與AFC、IFL_H、IFL_C、PP、FP、SCS等6個(gè)性狀的EBV具有顯著或極顯著的相關(guān)性，且GG基因型AFC、IFL_H、IFL_C、SCS等性狀EBV最低，PP、FP等性狀EBV最低；位于外顯子的g.51660569G>A位點(diǎn)的TC基因型AFC較低，IFL_H、IFL_C等性狀EBV最低。本研究所關(guān)注5個(gè)繁殖性狀的EBV越小，配種時(shí)間越短，個(gè)體繁殖性能越好，能夠保證奶牛養(yǎng)殖的經(jīng)濟(jì)效益[24]。對于產(chǎn)奶性狀，除SCS外，其余5個(gè)性狀的EBV越高，個(gè)體產(chǎn)奶性能越好。因此，g.51737889T>C位點(diǎn)GG基因型個(gè)體同時(shí)具有較好的繁殖性能和泌乳性能；g.51736640A>C位點(diǎn)GG基因型個(gè)體繁殖性能較好，但泌乳性能相對較差；g.51660569G>A位點(diǎn)TC基因型個(gè)體具有較好的繁殖性能。上述3個(gè)SNPs位點(diǎn)可作為中國荷斯坦牛繁殖和產(chǎn)奶性狀的候選位點(diǎn)，應(yīng)被重點(diǎn)關(guān)注。

此外，位于內(nèi)含子的g.51640656C>A位點(diǎn)與IFL_H_EBV關(guān)聯(lián)顯著，與PP_EBV關(guān)聯(lián)極顯著，且與SCS_EBV存在顯著關(guān)聯(lián)趨勢，該位點(diǎn)TG基因型IFL_H、IFL_C、FP最低，PP、SCS較低，具有該基因型個(gè)體的繁殖性能較好，但泌乳性能相對較差。

3.2 MET基因SNPs的生物學(xué)預(yù)測與表達(dá)量關(guān)聯(lián)分析

g.51737889T>C和g.51736640A>C兩位點(diǎn)均位于MET基因轉(zhuǎn)錄起始位點(diǎn)(transcription start site，TSS)上游約2 000 bp左右，研究表明，基因上游調(diào)控區(qū)多態(tài)性位點(diǎn)通過影響轉(zhuǎn)錄因子的結(jié)合調(diào)控基因表達(dá)[25]。在Campbell等[26]的研究中，在人MET基因啟動子區(qū)(-20 bp)證明存在核心啟動子區(qū)域，并且能夠通過與SP-1、AP-2、PC-4的結(jié)合影響啟動子活性。為進(jìn)一步探究這兩個(gè)突變位點(diǎn)的潛在作用機(jī)制，本研究進(jìn)行了基序分析，結(jié)果表明，當(dāng)g.51737889T>C位點(diǎn)的等位基因?yàn)镚時(shí)，富集的多種轉(zhuǎn)錄因子與基因激活和新陳代謝相關(guān)；當(dāng)此位點(diǎn)的等位基因?yàn)锳時(shí)，富集的多種轉(zhuǎn)錄因子與間質(zhì)上皮分化有關(guān)。推測，該SNP的GG基因型能夠招募與基因激活相關(guān)的轉(zhuǎn)錄因子結(jié)合序列，促進(jìn)MET基因的表達(dá)，最終導(dǎo)致該基因型個(gè)體繁殖和泌乳性能表現(xiàn)優(yōu)異。該位點(diǎn)與MET基因表達(dá)量的關(guān)聯(lián)分析結(jié)果證明，g.51737889T>C位點(diǎn)GG基因型個(gè)體的MET基因表達(dá)量最高，且與其他基因型存在顯著差異。同樣，當(dāng)g.51736640A>C位點(diǎn)的等位基因?yàn)镚時(shí)，富集的多種轉(zhuǎn)錄因子與細(xì)胞增殖、分化相關(guān)；當(dāng)此位點(diǎn)的等位基因?yàn)門時(shí)，富集的多種轉(zhuǎn)錄因子與葡萄糖代謝、脂質(zhì)代謝有關(guān)。推測，該SNP的GG基因型能夠招募與細(xì)胞分化相關(guān)的轉(zhuǎn)錄因子結(jié)合序列，促進(jìn)MET基因的表達(dá)，最終導(dǎo)致該基因型個(gè)體繁殖性能表現(xiàn)優(yōu)異。后續(xù)的分析中，g.51736640A>C位點(diǎn)GG基因型個(gè)體的MET基因表達(dá)量最高，且與其他基因型存在顯著差異，更進(jìn)一步證明了此猜測。

g.51660569G>A為同義突變位點(diǎn)，參與編碼Met蛋白第598位天冬氨酸，位于胞內(nèi)近膜區(qū)。有研究表明，同義突變影響核糖體運(yùn)輸、蛋白質(zhì)折疊以及細(xì)胞適應(yīng)過程[27-28]，所以對該位點(diǎn)進(jìn)行了mRNA結(jié)構(gòu)預(yù)測。結(jié)果表明，g.51660569G>A位點(diǎn)不同基因型mRNA二級結(jié)構(gòu)無明顯差異，但等位基因?yàn)镃時(shí)，mRNA的MFE更低，結(jié)構(gòu)更加穩(wěn)定。推測是該位點(diǎn)的同義突變通過改變mRNA的二級結(jié)構(gòu)，進(jìn)而影響mRNA穩(wěn)定性，延長半衰期，提高基因表達(dá)量，使不同基因型個(gè)體MET基因表達(dá)量存在差異，造成了其在卵泡和乳腺發(fā)育過程的差異。但該位點(diǎn)多態(tài)與基因表達(dá)量的關(guān)聯(lián)不顯著，預(yù)測結(jié)果并不能完全解釋雜合子基因型個(gè)體繁殖和泌乳性能表現(xiàn)優(yōu)異?；虻谋磉_(dá)存在時(shí)空特異性[29]，研究表明，Met的激活通過自分泌或旁分泌的方式啟動[17]，且Met主要通過改變蛋白活性影響表型，而不是基因表達(dá)[30]。部分研究也表明，在乳腺發(fā)育階段，Met蛋白高度表達(dá)，mRNA表達(dá)量很低，在哺乳期MetmRNA快速降解[2,31-32]。此次RNA提取樣品來自于泌乳牛尾根血液，并非卵巢或乳腺組織，這可能是表達(dá)量檢測結(jié)果與預(yù)測結(jié)果不一致的原因。此外，研究表明，內(nèi)含子SNP能夠通過影響可變剪切來影響基因表達(dá)[33-34]，或通過形成lncRNA調(diào)控基因表達(dá)[35]；g.51640656C>A位點(diǎn)多態(tài)與基因表達(dá)量關(guān)聯(lián)顯著，TG基因型表達(dá)量最低。但該SNP并非已證實(shí)的可變剪切位點(diǎn)，推測該位點(diǎn)對繁殖和產(chǎn)奶性狀的影響可能與內(nèi)含子形成的lncRNA有關(guān)，有待進(jìn)一步驗(yàn)證。

3.3 MET基因單倍型與繁殖、泌乳性狀EBV和基因表達(dá)量的關(guān)聯(lián)分析

研究表明，基因?qū)Ρ硇偷挠绊懣赡苁艿蕉鄠€(gè)突變位點(diǎn)的共同作用[29,36]。本研究中，在中國荷斯坦牛群體中MET基因的7個(gè)SNPs共形成2個(gè)單倍型塊，g.51666651G>C、g.51665308_51665316 dup、g.51660569G>A、g.51640656C>A、g.51645814C>T共5個(gè)SNPs緊密連鎖形成BLOCK1，g.51737889T>C與g.51736640A>C兩個(gè)SNPs緊密連鎖形成BLOCK2。BLOCK2與AFC、IFL_H、IFL_C、PP、FY、SCS多個(gè)繁殖、泌乳性狀EBV關(guān)聯(lián)顯著，且與PY、FP兩個(gè)泌乳性狀EBV關(guān)聯(lián)存在顯著趨勢。BLOCK2的2個(gè)SNPs g.51737889T>C與g.51736640A>C均位于上游調(diào)控區(qū)，這2個(gè)SNPs位點(diǎn)各自與多個(gè)繁殖、泌乳性狀EBV顯著相關(guān)。

另外，一個(gè)單倍型塊上的多個(gè)突變會對性狀產(chǎn)生更大的影響[37]。H4H4正是g.51737889T>C與g.51736640A>C兩個(gè)SNPs的優(yōu)勢基因型GG基因型的組合，H4H5是g.51737889T>C的雜合GA基因型與g.51736640A>C的GG優(yōu)勢基因型的組合；H4H4與H4H5個(gè)體排在AFC、IFL_H、IFL_C 3種繁殖性狀EBV最低的前三名；H4H4的PP_EBV第三高，F(xiàn)Y_EBV第二高，SCS_EBV最低。H4H4雙倍型個(gè)體具有較好的繁殖和泌乳性狀。且與基因表達(dá)量關(guān)聯(lián)分析結(jié)果也證明，H4H4雙倍型個(gè)體基因表達(dá)量最高，充分驗(yàn)證了g.51737889T>C與g.51736640A>C兩位點(diǎn)等位基因G對中國荷斯坦牛繁殖及泌乳性能的影響。

4 結(jié) 論

本研究獲得了中國荷斯坦牛MET基因的多態(tài)性圖譜，MET基因的7個(gè)SNPs及其形成的2個(gè)單倍型塊與繁殖、泌乳性狀存在顯著或極顯著關(guān)聯(lián)。位點(diǎn)g.51737889T>C和g.51736640A>C可能通過改變轉(zhuǎn)錄因子結(jié)合影響基因的表達(dá)，對MET的表達(dá)量有顯著影響；位點(diǎn)g.51660569G>A對mRNA的二級結(jié)構(gòu)無顯著影響但不同等位基因?qū)ψ钚∽杂赡墚a(chǎn)生影響。上述3個(gè)位點(diǎn)可作為奶牛繁殖和泌乳性狀的候選位點(diǎn)重點(diǎn)關(guān)注，研究結(jié)果為中國荷斯坦牛高產(chǎn)高效選育提供了可用的遺傳標(biāo)記。