戴瑞奶綿羊產奶性狀的全基因組關聯分析

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

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

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

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

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

1 材料與方法

1.1 試驗動物來源與產奶性狀

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

表1 試驗動物精補料營養水平

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

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

1.3 重測序基因分型和質量控制

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

1.4 樣本群體結構分析

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

1.5 全基因組關聯分析及統計推斷

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

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

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

1.6 候選基因注釋

在ENSEMBL網站下載綿羊基因組序列版本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注釋到其對應的基因上。

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

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

2 結 果

2.1 表型數據分析

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

表2 產奶性狀的表型描述性統計

2.2 基因分型和數據質量控制

通過基因分型共獲得22 683 429個SNPs。在進行質量控制時，SNP缺失大于10%的位點，剔除696 899個；最小等位基因頻率小于1%的位點，剔除0個；哈迪-溫伯格平衡檢驗P<10-6的位點，剔除69 658個。經過質量控制后，有135個樣本和21 020 909個 SNPs用于后續關聯分析。基因分型率達到98.54%。

2.3 樣本群體結構分析

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

圖1 親緣關系矩陣熱圖

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

圖2 PCA群體結構圖

2.4 全基因組關聯分析

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

A.90天日均產奶量性狀；B.150天日均產奶量性狀；C.泌乳周期性狀

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

表3 與產奶性狀關聯的SNPs

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

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

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

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

3 討 論

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

3.1 90天日均產奶量相關的候選基因分析

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

3.2 150天日均產奶量相關的候選基因分析

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

3.3 泌乳周期相關的候選基因分析

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

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

4 結 論

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