米布農(nóng)，張立果，烏日漢，郭玉林，王春偉，徐全忠，馮 爽，李光鵬，蘇小虎*，張 立*

(1.內(nèi)蒙古大學(xué)生命科學(xué)學(xué)院 省部共建草原家畜生殖調(diào)控與繁育國家重點實驗室，呼和浩特 010020；2.烏蘭察布市農(nóng)牧局 烏蘭察布市家畜畜牧工作站，烏蘭察布 012000)

奶綿羊具有優(yōu)良的產(chǎn)奶性能，且課題組前期的研究表明，綿羊奶的乳蛋白和乳脂顯著優(yōu)于山羊奶和牛奶[1]。戴瑞羊(Dairy Meade)是以新西蘭本地柯泊華茲羊為母本、東佛里生羊為父本，經(jīng)過20多年培育形成的新型乳用綿羊品種，于2016年被新西蘭綿羊育種聯(lián)合會認(rèn)定為新西蘭首個乳用綿羊品種，具有體型大、泌乳量高、繁殖性能好等優(yōu)良生產(chǎn)性狀[2]。戴瑞奶綿羊的泌乳期正常為180～230 d，初產(chǎn)平均產(chǎn)奶量為150～250 kg，經(jīng)產(chǎn)平均產(chǎn)奶量為350～600 kg。我國歷史上沒有專用的奶綿羊品種。為填補這一空白，內(nèi)蒙古大學(xué)于2016年引入了戴瑞奶綿羊胚胎，通過胚胎移植獲得了國內(nèi)純種戴瑞奶綿羊。隨后以本地小尾寒羊為母本，以戴瑞羊為父本開展級進(jìn)雜交獲得了一批新型奶綿羊繁育群體。

傳統(tǒng)的育種方式主要是利用表型信息和系譜信息進(jìn)行雜交選育，對遺傳力低和測量復(fù)雜的性狀，存在測定成本高、時間長、育種過程復(fù)雜的問題[3]。隨著高通量測序技術(shù)、生物信息學(xué)及統(tǒng)計學(xué)方法的快速發(fā)展，全基因組關(guān)聯(lián)分析(genome wide association study,GWAS)已成為禽畜重要性狀主效基因挖掘與鑒定的主要手段[4-5]。目前，GWAS已在牛、羊等產(chǎn)奶性狀的候選基因篩選方面取得了大量的研究成果。Liu等[6]使用90K Affymetrix Buffalo SNP陣列對意大利地中海水牛進(jìn)行GWAS，發(fā)現(xiàn)與產(chǎn)奶性能相關(guān)的4個SNPs，2個基因組區(qū)域和5個候選基因。Deng等[7]對地中海水牛產(chǎn)奶量全基因組關(guān)聯(lián)分析，發(fā)現(xiàn)了9個候選基因。Silva等[8]使用加權(quán)單步GWAS(weighted single-step GWAS，wssGWAS)對葡萄牙荷斯坦牛產(chǎn)奶性能全基因組關(guān)聯(lián)分析，發(fā)現(xiàn)與產(chǎn)奶量相關(guān)的QTL有51個，乳脂含量相關(guān)的QTL有5個，乳蛋白含量相關(guān)的QTL有24個。一項使用90 K Axiom?水牛陣列對埃及水牛產(chǎn)奶性狀進(jìn)行全基因組關(guān)聯(lián)分析的研究中，發(fā)現(xiàn)了47個顯著SNPs和11個基因座[9]。García-Gámez等[10]使用Illumina OvineSNP50 BeadChip對西班牙Churra羊產(chǎn)奶性能進(jìn)行全基因組關(guān)聯(lián)分析，發(fā)現(xiàn)14個染色體水平顯著的QTL和1個候選基因。Moioli等[11]使用Illumina SNP50K Bead chip對意大利Altamurana綿羊的產(chǎn)奶性能進(jìn)行了基因組掃描和全基因組關(guān)聯(lián)分析，發(fā)現(xiàn)2個候選基因。Li等[12]使用Illumina OvineSNP50對東弗里生和拉科訥雜交綿羊產(chǎn)奶性能全基因組關(guān)聯(lián)分析，發(fā)現(xiàn)與180天產(chǎn)奶量相關(guān)的SNP有115個，180天乳脂含量相關(guān)的SNP有288個，180天乳蛋白含量相關(guān)的SNP有74個。Sutera等[13]使用Illumina OvineSNP50K BeadChip對意大利Valle del Belice綿羊產(chǎn)奶性能進(jìn)行全基因組關(guān)聯(lián)分析，發(fā)現(xiàn)與產(chǎn)奶性能相關(guān)的9個SNPs。

然而針對奶綿羊的研究較少，一定程度限制了奶綿羊育種的發(fā)展。本研究基于前期已形成的戴瑞羊品種群體，通過對其產(chǎn)奶性狀測定和基因組重測序，利用GWAS方法挖掘其產(chǎn)奶性狀相關(guān)的候選基因及關(guān)鍵位點，為奶綿羊新品種的培育提供理論技術(shù)支持，為我國綿羊奶產(chǎn)業(yè)發(fā)展提供參考。

1 材料與方法

1.1 試驗動物來源與產(chǎn)奶性狀

本試驗所用的135只戴瑞奶綿羊繁育群體均來源于內(nèi)蒙古自治區(qū)烏蘭察布市蒙天然牧業(yè)科技發(fā)展有限公司。所選試驗羊均為初產(chǎn)，15～18月齡。于相同的環(huán)境、飼料、攝食比例條件下飼喂直至泌乳期結(jié)束。每天06：00及17：00飼喂試驗動物，試驗期間試驗動物可自由飲水。粗飼料是羊草，試驗日糧飼料組成粗精比為6∶4，母羊泌乳精補料營養(yǎng)水平見表1。產(chǎn)奶性狀為90天日均產(chǎn)奶量，150天日均產(chǎn)奶量和泌乳周期。所研究性狀均為分級性狀。根據(jù)群體泌乳量信息，以27%標(biāo)準(zhǔn)進(jìn)行產(chǎn)奶性狀等級劃分。產(chǎn)奶量以前后連續(xù)5 d測定的平均值為準(zhǔn)，90天日均產(chǎn)奶量性狀根據(jù)低于0.5 kg和高于1 kg的標(biāo)準(zhǔn)劃分為低產(chǎn)組和高產(chǎn)組，150天日均產(chǎn)奶量性狀以0.5 kg劃分低產(chǎn)組和高產(chǎn)組，泌乳周期性狀根據(jù)低于125 d和高于180 d的標(biāo)準(zhǔn)劃分為短周期組和長周期組。

表1 試驗動物精補料營養(yǎng)水平

1.2 血液基因組DNA提取與檢測

將采集的全血樣本通過康為世紀(jì)CWE9600DNA提取試劑盒V2提取基因組DNA。瓊脂糖凝膠電泳分析DNA降解程度以及是否有RNA污染；Nanodrop檢測DNA的純度(OD260 nm/OD280 nm比值)；Qubit對DNA濃度進(jìn)行精確定量。其中OD值在1.8～2.0之間，含量在1.5 μg以上的DNA樣品被用來建庫。

1.3 重測序基因分型和質(zhì)量控制

為了獲得核苷酸多態(tài)性信息，用illumina Novaseq6000平臺對奶綿羊基因組進(jìn)行重測序。有效測序數(shù)據(jù)通過BWA 0.7.17軟件[14]比對到參考基因組，比對結(jié)果經(jīng)SAMTOOLS 1.7軟件[15]去除重復(fù)；并將其映射到Oar_rambouillet_v1.0參考序列；映射的讀段通過SAMTOOLS分析獲取變異位點的分型；使用Plink1.90 軟件[16]進(jìn)行質(zhì)控，刪除 SNP缺失大于10%的位點、最小等位基因頻率小于1%的位點、哈迪-溫伯格平衡檢驗P<10-6的位點。

1.4 樣本群體結(jié)構(gòu)分析

由于群體分層會使全基因組關(guān)聯(lián)分析的結(jié)果出現(xiàn)假陽性[17]，因此采用Plink1.90對SNP分型數(shù)據(jù)計算IBS(Identity-by-state)遺傳距離矩陣和主成分分析，以研究奶綿羊群體的親緣關(guān)系以及群體分布。通過Plink軟件構(gòu)建親緣關(guān)系矩陣(G陣)，之后采用R語言的gg plots包的heatmap.2函數(shù)繪制熱圖，采用解釋方差最大的前三個主成分進(jìn)行散點圖繪制以分析群體結(jié)構(gòu)。

1.5 全基因組關(guān)聯(lián)分析及統(tǒng)計推斷

采用GEMMA v0.98.1軟件[18]進(jìn)行混合線性模型[19]分析：

y=Xβ+Ζκγκ+ξ+e

式中，y為表型向量，Xβ為PCA固定效應(yīng)，Ζκγκ為待檢驗標(biāo)記效應(yīng)，ξ～N(0,Kφ2)為多基因效應(yīng)，e～N(0,1σ2)為殘差效應(yīng)。多基因效應(yīng)中的K為標(biāo)記推斷的親緣關(guān)系矩陣。最終，每個SNP位點對應(yīng)一個關(guān)聯(lián)值。采用Bonferroni校正的閾值為全基因組顯著閾值設(shè)置為0.05/21 020 909，基因組水平潛在顯著閾值設(shè)置為1/21 020 909。曼哈頓圖以閾值線將關(guān)聯(lián)信息與顯著SNP進(jìn)行了可視化。曼哈頓圖和分位數(shù)圖由GEMMA繪制。

1.6 候選基因注釋

在ENSEMBL網(wǎng)站下載綿羊基因組序列版本Oar_rambouillet_v1.0(ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/002/742/125/GCF_002742125.1_Oar_rambouillet_v1.0/GCF_002742125.1Oar_rambouillet_v1.0_genomic.fna.g)，使用ANNOVAR軟件[20]將顯著SNP注釋到其對應(yīng)的基因上。

1.7 基因本體論和代謝途徑分析

對于基因本體(GO)和代謝途徑富集分析，使用cluster Profiler包[21]將注釋到的候選基因根據(jù)GO和KEGG數(shù)據(jù)庫進(jìn)行基因功能富集分析。利用R語言的enrich plot包進(jìn)行GO/KEGG的可視化畫圖。通過文獻(xiàn)查閱基因生物學(xué)功能，綜合分析確定影響奶綿羊產(chǎn)奶性狀的候選基因。

2 結(jié) 果

2.1 表型數(shù)據(jù)分析

表2展示了產(chǎn)奶性狀的平均值、標(biāo)準(zhǔn)差和變異系數(shù)。從表2中可以看出，90天日均產(chǎn)奶量的平均值是0.86 kg，150天日均產(chǎn)奶量的平均值是0.27 kg，泌乳周期的平均值是148 d，90天日均產(chǎn)奶量性狀的變異系數(shù)是86%，150天日均產(chǎn)奶量性狀的變異系數(shù)是140%，泌乳周期的變異系數(shù)是129%。

表2 產(chǎn)奶性狀的表型描述性統(tǒng)計

2.2 基因分型和數(shù)據(jù)質(zhì)量控制

通過基因分型共獲得22 683 429個SNPs。在進(jìn)行質(zhì)量控制時，SNP缺失大于10%的位點，剔除696 899個；最小等位基因頻率小于1%的位點，剔除0個；哈迪-溫伯格平衡檢驗P<10-6的位點，剔除69 658個。經(jīng)過質(zhì)量控制后，有135個樣本和21 020 909個 SNPs用于后續(xù)關(guān)聯(lián)分析?；蚍中吐蔬_(dá)到98.54%。

2.3 樣本群體結(jié)構(gòu)分析

圖1為135只奶綿羊的親緣關(guān)系矩陣熱圖。親緣關(guān)系矩陣熱圖是用顏色的深淺不同表示每只羊之間的親緣關(guān)系遠(yuǎn)近。顏色越淺，數(shù)值越小，親緣關(guān)系越遠(yuǎn)。因此，從圖1中可以看出，3個顏色深的地方代表3個家系。

圖1 親緣關(guān)系矩陣熱圖

圖2顯示了該奶綿羊群體結(jié)構(gòu)的分布，從圖中可以看出試驗群體分成了3個部分，說明該群體存在人工選育所引起的群體結(jié)構(gòu)差異。因此在后續(xù)全基因組關(guān)聯(lián)分析過程中需要考慮群體分層問題，并對其進(jìn)行校正。

圖2 PCA群體結(jié)構(gòu)圖

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

圖3為產(chǎn)奶性狀的分位數(shù)-分位數(shù)(Quantile-Quantile)圖，展示了90天日均產(chǎn)奶量、150天日均產(chǎn)奶量和泌乳周期的GWAS觀測值和預(yù)期值。虛線表示在SNP與所研究性狀不相關(guān)的原假設(shè)下SNP的分布。3個Q-Q圖的觀測值與預(yù)期P值的強烈偏差表明，與性狀顯著相關(guān)的SNP比預(yù)期的要多。

A.90天日均產(chǎn)奶量性狀；B.150天日均產(chǎn)奶量性狀；C.泌乳周期性狀

圖4為產(chǎn)奶性狀的曼哈頓圖，對分布在26條染色體上的21 020 909個SNPs進(jìn)行了關(guān)聯(lián)分析。Bonferroni校正后的顯著水平閾值為2.38×10-9，基因組水平潛在顯著閾值為4.76× 10-8。表3展示了與產(chǎn)奶性狀顯著相關(guān)的SNPs名稱、物理位置以及鄰近基因等。

表3 與產(chǎn)奶性狀關(guān)聯(lián)的SNPs

黑色水平線代表全基因組水平顯著閾值線，紅色水平線為全基因組潛在顯著關(guān)聯(lián)閾值線。A.90天日均產(chǎn)奶量性狀；B.150天日均產(chǎn)奶量性狀；C.泌乳周期性狀

與90天日均產(chǎn)奶量潛在顯著關(guān)聯(lián)的SNP有8個，分別位于2、13、18號染色體上。與90天日均產(chǎn)奶量顯著關(guān)聯(lián)的SNP有1個，位于1號染色體上。1號染色體上的1個顯著SNP(849495)位于TRNAQ-CUG-2(未表征的基因)下游102 468 bp處和LOC114117240(未表征的基因)上游151 509 bp處；2號染色體3個潛在顯著SNPs位于基因間區(qū)，鄰近ACADL(?；o酶A脫氫酶長鏈)和MYL1(肌球蛋白輕鏈1)3 kb左右；13號染色體上4個潛在顯著性SNPs位于CHD6(染色體域解旋酶DNA結(jié)合蛋白6)的內(nèi)含子區(qū)域；18號染色體上1個潛在顯著SNP(17666953)位于SLCO3A1(溶質(zhì)載體有機陰離子轉(zhuǎn)運蛋白家族成員3A1)的內(nèi)含子區(qū)域。

與150天日均產(chǎn)奶量潛在顯著關(guān)聯(lián)的SNP共2個，分別位于1和16號染色體上。1號染色體上1個潛在顯著SNP(727824)位于PRMT6(蛋白質(zhì)精氨酸甲基轉(zhuǎn)移酶6)上游499 695 bp處；16號染色體上1個潛在顯著性SNP(16293028)位于RNF180(無名指蛋白180)的內(nèi)含子區(qū)域。

與泌乳周期潛在顯著關(guān)聯(lián)的SNP共 2個，分別位于 1和6號染色體上。1號染色體上1個潛在顯著SNP(727824)位于基因PRMT6(蛋白質(zhì)精氨酸甲基轉(zhuǎn)移酶6)上游499 695 bp處；6號染色體上1個 潛在顯著SNP(9582478)位于TRNAW-CCA-68(未表征的基因)下游1 049 756 bp處和TRNAS-GGA-61(未表征的基因)上游198 531 bp處。

3 討 論

在禽畜育種中，GWAS分析已成為識別表型變異的有力策略[22]。目前，尚無關(guān)于戴瑞奶綿羊全基因組關(guān)聯(lián)分析對產(chǎn)奶性狀影響的研究報道。本研究對奶綿羊群體產(chǎn)奶性狀進(jìn)行全基因組關(guān)聯(lián)分析，共得到13個相關(guān)的SNPs。對于90天日均產(chǎn)奶量性狀，共檢測到9個顯著關(guān)聯(lián)的SNPs，1號染色體上有1個SNP達(dá)到基因組顯著水平，達(dá)到潛在顯著水平的8個SNPs分別位于2、13和18號染色體上；對于150天日均產(chǎn)奶量性狀，有2個SNPs達(dá)到潛在顯著水平，分別位于1和16號染色體上；對于泌乳周期性狀，達(dá)到潛在顯著水平的2個SNPs分別位于1和6號染色體上。

3.1 90天日均產(chǎn)奶量相關(guān)的候選基因分析

與90天日均產(chǎn)奶量關(guān)聯(lián)的9個SNPs，snp849495鄰近TRNAQ-CUG-2、LOC114117240基因，其功能尚未有報道；2號染色體46 800～43 690 bp區(qū)間內(nèi)存在3個SNPs與90天日均產(chǎn)奶量潛在關(guān)聯(lián)，鄰近ACADL和MYL1基因。ACADL是長鏈脂肪酸β-氧化起始步驟的催化酶，在長鏈脂肪酸β-氧化中發(fā)揮著重要作用，與脂肪代謝有著密切的關(guān)系[23]?；虮倔w分析表明，ACADL主要在脂肪酸分解的生物過程(GO:0009062)、氧化還原酶活性的分子功能(GO:0016627)和線粒體的細(xì)胞成分(GO:0005759)中起作用，并參與脂肪酸降解(KEGG：ssc00071)、脂肪酸代謝(KEGG：ssc01212)、PPAR信號通路(KEGG：ssc03320)。Zhang等[24]的研究表明，ACADL催化長鏈脂肪酸(C13-C22)的β-氧化，并調(diào)節(jié)線粒體中脂肪酸的代謝途徑，抑制ACADL會導(dǎo)致脂肪酸蓄積，從而激活過氧化物酶體增殖物激活受體(PPAR)，PPARα激活劑通過拮抗NF-κB的轉(zhuǎn)錄活性而在幾種細(xì)胞中顯示出抗炎活性。有研究發(fā)現(xiàn)，16碳及以下脂肪酸合成增加會促使產(chǎn)奶量升高，而乳中16碳以上脂肪酸比例的增加使產(chǎn)奶量降低[25]。ACADL可能通過降解長鏈脂肪酸進(jìn)而影響產(chǎn)奶量，因此，該基因可能是與90天日均產(chǎn)奶量相關(guān)的重要候選基因。MYL1基因是肌球蛋白輕鏈1，在骨骼肌不同發(fā)育時期發(fā)揮著重要作用[26]。對于MYL1基因與產(chǎn)奶相關(guān)的報道較少，是因為該品種也具有良好的產(chǎn)肉性能，所以才鑒定到MYL1基因。13號染色體上有4個SNPs位于CHD6內(nèi)。CHD6是一種蛋白質(zhì)編碼基因。Lutz等[27]研究發(fā)現(xiàn)，CHD6和轉(zhuǎn)錄因子(p300,SRC1)結(jié)合形成PRIC復(fù)合物，該復(fù)合物可以促使PPAR激動劑和PPARα結(jié)合，PPARα參與脂質(zhì)代謝和抗氧化過程。張雪瑩[28]的研究表明，在促進(jìn)脂肪酸分解代謝的組織中PPARα表達(dá)量相對較高。有文獻(xiàn)報道過氧化物酶體增殖物激活受體γ(PPARγ)可能是乳脂合成調(diào)控核心[29]。PPARα和PPARγ屬于同一家族。CHD6可能通過與PPARα作用參與脂質(zhì)代謝調(diào)控。因此，CHD6可能也是與90天日均產(chǎn)奶量相關(guān)的候選基因。snp17666953位于SLCO3A1基因內(nèi)。SLCO3A1屬于溶質(zhì)載體有機陰離子轉(zhuǎn)運蛋白家族成員，是一種蛋白質(zhì)編碼基因。該基因與維生素、核苷、葡萄糖、膽汁鹽和金屬離子等的運輸有關(guān)。SLCO3A1主要涉及的生物學(xué)過程有類花生酸轉(zhuǎn)運(GO:0071715)和脂肪酸衍生物轉(zhuǎn)運(GO:1901571)。與該基因相關(guān)的分子功能是有機陰離子跨膜轉(zhuǎn)運蛋白活性(GO:0008514)。先前的研究中，Ibeagha-Awemu等[30]通過全基因組關(guān)聯(lián)分析篩選到SLCO3A1是影響奶牛產(chǎn)奶性狀和乳腺功能的新候選基因之一。Liu等[31]使用Illumina BovineSNP150 BeadChip對中國荷斯坦奶牛產(chǎn)奶性能進(jìn)行全基因組關(guān)聯(lián)分析，發(fā)現(xiàn)了10個數(shù)量性狀基因座和8個候選基因,其中4個與脂肪和蛋白質(zhì)顯著相關(guān)的SNPs位于EPHA6、SLCO1A2、DGAT1和EP400內(nèi)部。SLCO1A2 是溶質(zhì)載體有機陰離子轉(zhuǎn)運蛋白家族成員1A2,此基因和SLCO3A1基因都有轉(zhuǎn)運有機陰離子的功能。故推測SLCO3A1可能是與90天日均產(chǎn)奶量相關(guān)的關(guān)鍵候選基因。

3.2 150天日均產(chǎn)奶量相關(guān)的候選基因分析

與150天日均產(chǎn)奶量存在潛在顯著相關(guān)的2個SNPs，1號染色體上的snp727824，其下游499 695 bp處是PRMT6基因，是蛋白質(zhì)精氨酸甲基轉(zhuǎn)移酶6，編碼的蛋白質(zhì)屬于精氨酸N-甲基轉(zhuǎn)移酶家族。PRMTs主要參與轉(zhuǎn)錄、RNA剪接、DNA損傷修復(fù)、細(xì)胞信號轉(zhuǎn)導(dǎo)、蛋白質(zhì)運輸?shù)戎匾飳W(xué)過程[32-36]。該酶可催化甲基從S-腺苷甲硫氨酸(SAM)轉(zhuǎn)移至組蛋白和其他蛋白質(zhì)上的精氨酸殘基。這種甲基化的失調(diào)在某些癌癥的發(fā)展中至關(guān)重要。Chan等[37]研究發(fā)現(xiàn)，PRMT6在乳腺癌中異常表達(dá)。有研究發(fā)現(xiàn)同屬于蛋白質(zhì)精氨酸甲基轉(zhuǎn)移酶家族的PRMT2基因表達(dá)的多少及部位與乳腺癌細(xì)胞的表達(dá)密切相關(guān)[38]。因此，該基因可能是影響150天日均產(chǎn)奶量的候選基因。snp16293028位于RNF180基因內(nèi)。RNF180是E3泛素連接酶家族的重要成員，通過調(diào)節(jié)靶蛋白底物的泛素化而影響細(xì)胞的生長和發(fā)育，與各種生物過程如細(xì)胞生殖、分化、增殖、凋亡和腫瘤的發(fā)生等高度相關(guān)[39]。GO注釋表明，RNF180主要在兒茶酚胺代謝的生物過程(GO:0006584)、與泛素結(jié)合酶結(jié)合的分子功能(GO:0031624)和內(nèi)質(zhì)網(wǎng)膜的細(xì)胞成分(GO:0031227)中起作用。有研究表明，轉(zhuǎn)運膜蛋白的泛素化利于乳腺上皮細(xì)胞外脂肪酸的轉(zhuǎn)入，使細(xì)胞利用更多的脂肪酸合成乳脂[40]。因此，RNF180可能是影響150天日均產(chǎn)奶量的重要候選基因。

3.3 泌乳周期相關(guān)的候選基因分析

snp9582478位于6號染色體，鄰近TRNAW-CCA-68和TRNAS-GGA-61基因。TRNAW-CCA-68和TRNAS-GGA-61的基因功能尚未有報道。snp727824位于1號染色體，鄰近PRMT6基因。

通過對以上候選基因的功能分析，ACADL、SLCO3A1可作為影響90天日均產(chǎn)奶量的候選基因，PRMT6可作為影響150天日均產(chǎn)奶量和泌乳周期的候選基因，后期仍需要進(jìn)一步對挖掘到的候選基因進(jìn)行功能驗證。

4 結(jié) 論

本研究檢測到13個與90天、150天日均產(chǎn)奶量和泌乳周期性狀相關(guān)的SNPs，發(fā)現(xiàn)了10個候選基因，基因功能分析發(fā)現(xiàn)ACADL、SLCO3A1有可能是影響產(chǎn)奶性狀的候選基因。本研究結(jié)果為深入揭示奶綿羊產(chǎn)奶性狀的遺傳效應(yīng)和標(biāo)記輔助選種提供理論技術(shù)支持。