王添禎，高 雪，宋文芹，姚大為，陳麗麗，陳成彬*，馬 毅*

(1. 天津市畜牧獸醫(yī)研究所，天津 300381；2. 南開大學(xué) 生命科學(xué)院，天津 300071;3. 中國(guó)農(nóng)業(yè)科學(xué)院 北京畜牧獸醫(yī)研究所，北京 100193)

伴隨全基因組關(guān)聯(lián)分析(Genome-wide Association Study, GWAS)在人類[1]和家畜[2]中的應(yīng)用，許多與復(fù)雜疾病和數(shù)量性狀相關(guān)的候選基因和因果變異位點(diǎn)被鑒定。該方法不僅有助于解析復(fù)雜表型的分子遺傳機(jī)制，同時(shí)有助于促進(jìn)基因組選擇的應(yīng)用及加快遺傳改良速度[3-4]。中國(guó)荷斯坦牛是我國(guó)主要奶牛品種，目前不少研究報(bào)道對(duì)中國(guó)荷斯坦牛的產(chǎn)奶量[5]、乳脂等[6]進(jìn)行了GWAS分析，確定了多個(gè)候選基因以及數(shù)量性狀位點(diǎn)(Quantitative trait locus, QTLs)。然而,其他重要的經(jīng)濟(jì)及健康性狀尚研究不足。

奶牛體型性狀包含肢蹄結(jié)構(gòu)(包括蹄角度、肢蹄評(píng)分以及整體的肢蹄結(jié)構(gòu)評(píng)分)和乳房形態(tài)(包括乳房深度、后乳頭位置、乳頭長(zhǎng)度、后乳房寬度、前乳房附著、后乳房高度以及整體的乳房形態(tài)評(píng)分)。具體的評(píng)分標(biāo)準(zhǔn)參照美國(guó)2018年荷斯坦牛的評(píng)分標(biāo)準(zhǔn)。肢蹄病是常見的疾病，可以造成生產(chǎn)性能下降；乳房形態(tài)(乳房長(zhǎng)、寬、基部圍等)則與產(chǎn)奶量呈中等以上正相關(guān)。在過去幾十年里，研究者主要關(guān)注產(chǎn)奶性狀，利用基因組選擇(Genetic Selection,GS)、調(diào)整牛場(chǎng)管理模式以及改良飼料配方等多種策略，使奶牛的產(chǎn)奶量實(shí)現(xiàn)快速增長(zhǎng)。當(dāng)前，奶牛育種更加關(guān)注平衡育種，體型性狀作為奶牛的重要性狀，影響奶牛的健康以及產(chǎn)奶量等重要經(jīng)濟(jì)指標(biāo)，因此體型性狀也逐漸被重視。然而體型性狀的遺傳基礎(chǔ)研究國(guó)內(nèi)外報(bào)道仍然較少。本研究旨在進(jìn)行中國(guó)荷斯坦牛肢蹄結(jié)構(gòu)與乳房形態(tài)兩個(gè)性狀的全基因組關(guān)聯(lián)分析，確定顯著的SNP候選位點(diǎn)，為進(jìn)一步開展奶牛分子育種提供參考。

1 材料與方法

1.1 樣本基因型和質(zhì)量控制

選取天津地區(qū)奶牛場(chǎng)中2014年出生的300頭中國(guó)荷斯坦母牛(使用8個(gè)父系凍精，從PCA的結(jié)果中也可以看出，樣本大致可分為3個(gè)家系)進(jìn)行血液樣本采集，并使用TiangenTIANamp Blood基因組提取試劑盒提取全基因組DNA，采用GeneSeek Genomic Profiler Bovine 50 K SNP chip芯片對(duì)所有個(gè)體進(jìn)行基因分型，SNP均勻分布在牛基因組上，平均標(biāo)記密度為59 kb。

用PLINK v1.9[7]軟件進(jìn)行基因型數(shù)據(jù)的質(zhì)量控制，包括個(gè)體水平和SNPs水平兩個(gè)方面：(1)刪除檢出率低于94%的個(gè)體數(shù)據(jù)；刪除所有基因型數(shù)據(jù)缺失超過10%的個(gè)體；去除缺失基因型超過1%、最小等位基因頻率小于5%和哈代-溫伯格平衡檢驗(yàn)P值小于1.0×10-7的SNP數(shù)據(jù)。經(jīng)過質(zhì)量控制后，共有274個(gè)樣本的40 501個(gè)高質(zhì)量SNPs用于GWAS分析，占總SNPs的84.7%，其所在牛染色體上的分布情況如表1所示。

1.2 預(yù)期傳遞力

根據(jù)美國(guó)農(nóng)業(yè)部奶牛育種委員會(huì)(The Council of Dairy Cattle Breeding, CDCB, https://www.uscdcb.com/)2018年更新的標(biāo)準(zhǔn)方法計(jì)算得到預(yù)期傳遞力(Predicted transfer Ability, PTA)作為表型[8]，對(duì)體型性狀中的肢蹄結(jié)構(gòu)、乳房形態(tài)性狀進(jìn)行全基因組關(guān)聯(lián)分析。

表1 質(zhì)控后SNPs的分布和相鄰SNPs的平均距離

1.3 統(tǒng)計(jì)分析

使用PLINK (v1.9)進(jìn)行群體主成分分析(Principal Components Aalysis,PCA)，對(duì)膨脹系數(shù)(λ)進(jìn)行評(píng)估。根據(jù)PCA結(jié)果對(duì)群體進(jìn)行校正，避免因群體分層而導(dǎo)致較高的假陽性[9]。采用EMMAX[10]軟件中基于混合線性模型[11]的單標(biāo)記回歸方法進(jìn)行全基因組關(guān)聯(lián)分析，建立模型如下：

y=Xβ+Zυ+e

1.4 關(guān)聯(lián)分析的顯著性檢驗(yàn)

采用Bonferroni檢驗(yàn)校正。單次檢驗(yàn)的顯著性水平為：染色體水平顯著(1/N)、全基因組水平顯著(0.05/N)，其中N為質(zhì)控后可用SNPs數(shù)為40 501。即染色體水平顯著閾值為(1/40 501=2.47×10-5)，全基因組水平顯著(0.05/40 501=1.23×10-6)。

1.5 群體分層檢驗(yàn)

利用R v3.5.1[12]繪制Q-Q圖(Quantile-Quantile plot)進(jìn)行群體分層檢驗(yàn)。以此來判斷樣本群體是否存在群體分層及分析性狀顯著相關(guān)的SNPs。

1.6 候選基因分析策略

基于芯片對(duì)應(yīng)的Bos_taurus genome UMD 3.1[4]在線數(shù)據(jù)庫(kù)(http://genome-asia.-ucsc. edu/)，在顯著SNPs的上、下游各50 kb區(qū)域內(nèi)進(jìn)行候選基因的篩查,并對(duì)基因功能進(jìn)行注釋確定影響奶牛肢蹄結(jié)構(gòu)和乳房形態(tài)性狀的候選基因。

2 結(jié) 果

2.1 全基因組關(guān)聯(lián)分析結(jié)果

通過對(duì)芯片數(shù)據(jù)進(jìn)行嚴(yán)格的質(zhì)量控制，共得到274頭母牛40 501個(gè)SNPs位點(diǎn)的基因型結(jié)構(gòu)。對(duì)中國(guó)荷斯坦牛肢蹄結(jié)構(gòu)、乳房形態(tài)這兩個(gè)性狀進(jìn)行全基因組關(guān)聯(lián)分析。利用R(3.5.1)繪制全基因組范圍SNPs的-log10(P-value)分布曼哈頓圖，如圖1所示。通過Bonferroni校正，對(duì)于肢蹄結(jié)構(gòu)共檢測(cè)到8個(gè)染色體水平顯著(P<1/40 501)SNPs位點(diǎn)，對(duì)于乳房形態(tài)共檢測(cè)到16個(gè)染色體水平顯著SNPs位點(diǎn)和1個(gè)全基因組水平顯著(P<0.05/40 501)位點(diǎn)。

2.2 群體分層檢驗(yàn)結(jié)果

主成分分析結(jié)果如圖2，從圖中可以看出，樣本群體可分為3個(gè)組，在全基因組關(guān)聯(lián)分析中，把群體親緣關(guān)系作為效應(yīng)考慮進(jìn)模型中，能夠避免因群體分層導(dǎo)致的假陽性結(jié)果。從Q-Q 圖結(jié)果(圖3)中顯示通過SNPs關(guān)聯(lián)分析計(jì)算的χ2統(tǒng)計(jì)量分布與假設(shè)檢驗(yàn)沒有過早偏離，說明采用的分析模型是合理的。而在圖形的右上角則是顯著較高的位點(diǎn)，是潛在與性狀相關(guān)的候選位點(diǎn)。

圖1 全基因組SNPs分布散點(diǎn)圖

圖2 主成分分析圖

2.3 顯著SNPs及相關(guān)基因

通過Bonferroni校正，共檢測(cè)出25個(gè)顯著SNPs，其中8個(gè)與肢蹄結(jié)構(gòu)相關(guān)，顯著SNP分別位于4(2個(gè))、6(1個(gè))、11(1個(gè))、14(2個(gè))、15(1個(gè))、26(1個(gè))等6條染色體上；17個(gè)與乳房形態(tài)相關(guān)，顯著SNP分別位于5(1個(gè))、6(1個(gè))、8(1個(gè))、9(1個(gè))、13(2個(gè))、14(5個(gè))、15(1個(gè))、16(1個(gè))、22(1個(gè))、23(3個(gè))等10條染色體上，其中2個(gè)位于基因內(nèi)部。對(duì)25個(gè)顯著SNPs位點(diǎn)進(jìn)行分析注釋，得到22個(gè)與其距離最近的候選基因。對(duì)候選基因進(jìn)行注釋，其中有3個(gè)編碼microRNA，19個(gè)蛋白質(zhì)編碼類基因(見表2)。

3 討 論

3.1 預(yù)期傳遞力

在GWAS中，表型校正的方式有所不同，有些研究采用EBV[12-14]; 然而美國(guó)研究機(jī)構(gòu)的大部分學(xué)者采用PTA[15-17]育種值作為表型進(jìn)行關(guān)聯(lián)分析。本研究采用PTA作為表型,可以有效排除其它環(huán)境效應(yīng)的干擾[18]，通過計(jì)算PTA估計(jì)基因組的遺傳效應(yīng)。出生奶牛的初始PTA為父母的均值，隨著后期成長(zhǎng)，綜合利用表型及多層面信息來校正PTA，PTA的參考群體每五年進(jìn)行更新以保證數(shù)據(jù)的可靠性。

圖3 EMMAX方法結(jié)果的QQ圖

表2 中國(guó)荷斯坦牛肢蹄結(jié)構(gòu)和乳房形態(tài)全基因組關(guān)聯(lián)分析顯著SNPs和相關(guān)基因

3.2 全基因組關(guān)聯(lián)分析

本研究采用GeneSeek Genomic Profiler Bovine 50 K SNP芯片對(duì)中國(guó)荷斯坦牛的肢蹄結(jié)構(gòu)和乳房形態(tài)兩個(gè)性狀進(jìn)行了全基因組關(guān)聯(lián)分析。利用PLINK軟件進(jìn)行主成分分析，前兩個(gè)主成分累計(jì)大于93%，基于前兩個(gè)組分進(jìn)行群體分層分析，研究發(fā)現(xiàn)樣本個(gè)體大致可分為3組。因此，為了提高檢測(cè)多效應(yīng)遺傳變異的效力及避免關(guān)聯(lián)分析的假陽性[19]，在GWAS分析中，本研究將群體分層和SNP效應(yīng)作為固定效應(yīng)，親緣關(guān)系矩陣作為隨機(jī)效應(yīng)[20-21]有效地提高了關(guān)聯(lián)分析的檢測(cè)效力。

3.3 候選基因功能注釋

通過數(shù)據(jù)庫(kù)對(duì)候選基因進(jìn)行功能注釋，發(fā)現(xiàn)多個(gè)涉及代謝通路調(diào)節(jié)的候選基因，例如：LEP、HS3ST1、RB1CC1、ASPH與生長(zhǎng)發(fā)育相關(guān)，OSTF1、FKBP5參與免疫調(diào)節(jié)，ZMYND8、PTK2與產(chǎn)奶性狀相關(guān)。

LEP(瘦素)是一種主要由脂肪細(xì)胞合成和分泌的肽激素，位于4號(hào)染色體上，在胎盤、胃和骨骼肌等組織中表達(dá)，參與多種信號(hào)通路，如脂質(zhì)代謝、葡萄糖轉(zhuǎn)運(yùn)、胰島素分泌[22]等，在調(diào)節(jié)能量穩(wěn)態(tài)中起著至關(guān)重要的作用。例如：LEP參與調(diào)節(jié)5'-AMP活化的蛋白激酶(AMPK)的信號(hào)傳導(dǎo)。AMPK作為能量傳感器，響應(yīng)AMP與ATP比率的增加而激活?；罨腁MPK通過調(diào)節(jié)脂肪酸生物合成-乙酰輔酶A羧化酶(ACC)的活性來調(diào)節(jié)脂肪酸生物合成。LEP表達(dá)量高的動(dòng)物，其日增重可能比其他飼喂量和營(yíng)養(yǎng)狀況相似的動(dòng)物要低，產(chǎn)犢期也可能更長(zhǎng)[23]，在牛皮膚組織中，LEP通過旁分泌參與控制表皮生長(zhǎng)和毛囊循環(huán)[24]。

OSTF1在細(xì)胞溶菌產(chǎn)物中表達(dá)，具有特殊的吸收功能，參與某些局部炎癥病灶的吸收。在荷斯坦奶牛中，OSTF1有5種基因型(AA、AB、BB、AC、CC)，其中BB型的牛奶中具有較低的體細(xì)胞，表明OSTF1的BB基因型可以作為乳房炎抗性遺傳標(biāo)記[25]。

ZMYND8是一個(gè)蛋白質(zhì)編碼基因，轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)的關(guān)鍵組成部分。在同源小鼠中，注射ZMYND8過表達(dá)的浸潤(rùn)性乳腺癌細(xì)胞，腫瘤體積和重量減小，與正常組織相比，乳腺癌組織中ZMYND8顯著下調(diào)[26]。Venturini等[27]通過GWAS發(fā)現(xiàn)有一個(gè)SNP位于ZMYND8基因內(nèi)，并與產(chǎn)奶量、乳脂量、乳蛋白量顯著相關(guān)。

PTK2作為整合素信號(hào)傳導(dǎo)的主要介質(zhì)，已被發(fā)現(xiàn)對(duì)乳腺上皮細(xì)胞的存活、增殖和分化具有重要作用[28]，可能參與幾種信號(hào)轉(zhuǎn)導(dǎo)途徑，如細(xì)胞運(yùn)動(dòng)[29]、微管穩(wěn)定[30]、細(xì)胞與細(xì)胞連接的調(diào)節(jié)[31]。此外，在維持乳腺發(fā)育和體內(nèi)功能方面發(fā)揮著重要作用[32]。Wang等[33]在GWAS分析中也鑒定出了PTK2與產(chǎn)奶性狀相關(guān)，并通過后續(xù)的試驗(yàn)驗(yàn)證了不同的突變體對(duì)PTK2的表達(dá)量有顯著相關(guān)，在所有檢測(cè)到的組織中均有表達(dá)，在乳腺，子宮和腎臟中表達(dá)水平較高,TT基因型的乳腺中表達(dá)量最高。

此外，Kim等[34]鑒定出HS3ST1與牛的胴體重量相關(guān)，Magalhaes等[35]鑒定出RB1CC1與生長(zhǎng)和肌肉發(fā)育相關(guān)。ASPH在牛中發(fā)現(xiàn)與肌肉肥大有關(guān)[36]，并且是胴體性狀的候選基因[37]。FKBP參與TNFα/NF-kB信號(hào)通路，是宿主對(duì)疾病和其它有害應(yīng)激物免疫反應(yīng)的主要途徑[14]。

4 結(jié) 論

本研究通過全基因組關(guān)聯(lián)分析，檢測(cè)了影響中國(guó)荷斯坦牛肢蹄結(jié)構(gòu)和乳房形態(tài)的候選基因，包括LEP、HS3ST1、RB1CC1、ASPH、OSTF1、FKBP5、ZMYND8、PTK2等，為揭示奶牛肢蹄結(jié)構(gòu)和乳房形態(tài)的分子遺傳基礎(chǔ)提供了重要參考。