基于重測序數據的湖羊產羔數性狀全基因組關聯分析

謝銳，朱墨，童世鋒，趙福平，劉楊*

(1. 南京農業大學動物科技學院，江蘇 南京 210095； 2. 農業農村部都市農業重點實驗室，上海 200240；3. 中國農業科學院北京畜牧獸醫研究所，北京 100193)

繁殖性狀是羊的重要經濟性狀，繁殖性能的高低直接決定了整個養殖經濟效益的高低。如何提高羊產羔數一直是科學家們致力于研究解決的重要課題。提高母羊產羔數的前提是提高母羊的排卵數，而哺乳動物的排卵過程受到復雜的機制調控，而且不同品種羊之間還存在著顯著差異，如特塞爾綿羊和薩福克羊每個排卵周期只排1個卵，而高產的布魯拉美利奴羊每個排卵周期可以排10個卵[1]。除了國外的羊具有優秀繁殖性能外，我國本地綿羊品種湖羊也因具有“一年多胎，一胎多羔”的優秀繁殖性能而聞名世界。但影響湖羊產羔數的關鍵功能基因以及高繁殖力性狀的遺傳機理還有待進一步研究。

隨著二代基因組測序技術的快速發展，越來越多的研究利用基因組測序方法來尋找與羊繁殖性狀有關的候選基因，這為從基因組水平上研究目標經濟性狀提供了便利。全基因組關聯分析(GWAS)是在全基因組范圍內尋找與動植物重要經濟性狀相關聯的遺傳變異，是一種復雜性狀功能基因鑒定的分析策略[2]。GWAS概念是由Risch在1996年提出[3]，最早是被應用在醫學上研究復雜疾病的遺傳模式[4]。GWAS方法因具有速度快、通量高和高精度等顯著優點,其在畜牧學遺傳育種中的作用日益凸顯，目前已經成為挖掘與家畜重要經濟性狀相關候選基因最重要的方法之一[5]。狄江等[6]采用單標記回歸方法對中國美利奴(新疆型)周歲細毛羊進行全基因組關聯分析，檢測到3個與羊毛產量和2個與羊毛長度顯著相關的候選SNP位點，并根據這些SNP位點信息鑒定到3個與羊毛生長、發育相關的生物學過程相關基因，分別為IBIN、HSD17B11和PIAS1基因。王珊等[7]利用全基因組重測序數據對多浪羊產羔性狀進行了全基因組關聯分析，經過基因注釋和富集分析篩選出了2個與多浪羊高繁殖力相關的基因。張軍霞等[8]通過多態性位點與性狀的關聯分析發現小尾寒羊中高繁殖力群體產羔數與KITLG基因g.47468C>T位點顯著相關，可作為產羔性狀的候選基因。蘭蓉等[9]利用全基因組關聯分析方法對云南黑山羊產羔性狀進行研究，發現2個與山羊產羔數顯著相關的SNP位點。劉世佳等[10]利用關聯分析的方法發現BMP15基因B2突變(C718T)與產羔數顯著相關，可作為提高綿羊產羔數性狀的分子標記應用于育種。本研究旨在通過重測序數據，結合GWAS分析技術對影響湖羊產羔數性狀的候選基因進行挖掘，這有利于進一步揭示湖羊高繁殖性狀的遺傳基礎，縮短湖羊的育種年限，降低育種成本，提高經濟效益。

1 材料與方法

1.1 樣本采集及測序

從甘肅某羊場隨機挑選206只湖羊，統一采集耳組織樣，置于-20℃的干冰中冷凍保存，用于提取基因組DNA。樣品采集的DNA經過質控檢測合格之后，統一送到北京貝瑞和康生物技術有限公司進行雙端150 bp全基因組重測序，平均測序深度為6×。質控后共得到3 301 Gb的有效數據(clean reads)，用于后期的變異檢測。

1.2 變異檢測及SNP質控

使用綿羊4.0(Ovis_aries 4.0)參考基因組，并使用BWA[11](version 0.7.17)的mem模塊測序得到的clean reads 與參考基因組進行比對；使用GATK4[12](version 4.0.2.1)中的Picard 模塊標記PCR重復序列，變異檢測流程使用GATK4軟件中的Haplotype Caller模塊和Variant Filtration模塊進行SNP的鑒定和初步過濾，將結果輸出到VCF文件保存。最后使用VCFtools[13](version 0.1.17)對SNP進行質控，并保留符合以下標準的SNP：①所有SNP位點的平均測序深度大于2×；②只保留常染色體上的SNP；③SNP最大缺失率5%；④只保留二等位基因；⑤最小等位基因頻率大于0.05；⑥哈迪溫伯格平衡P值大于1×10-6。

1.3 GWAS分析

使用EMMAX[14](version beta-07)中的混合線性模型，利用湖羊產羔數表型性狀與全基因組上SNP來進行GWAS分析，混合線性模型如下：

y=μ+Xb+u+e，

其中，y表示湖羊產羔數性狀表型值；μ表示產羔數平均值；X表示固定效應矩陣，b為固定效應向量；u表示剩余多基因效應；e表示產羔數表型值的隨機殘差，u和e均服從正態分布。

本研究中對湖羊GWAS結果中的P值采用Bonferroni校正，并使用R軟件中的CMplot包進行曼哈頓圖和Q_Q圖(Quantile-Quantile Plot)繪制。

2 結果與分析

2.1 表型數據統計

選用的206只母羊均有第一胎產羔數數據，其中181只母羊有第二胎產羔數數據，其余25只母羊第二胎產羔數數據缺失。在湖羊第一胎產羔數數據中，產羔數最多的為4只，產羔數最少的為1只，平均產羔數為2.1只；在湖羊第二胎產羔數數據中，產羔數最多的為5只，產羔數最少的為1只，平均產羔數為2.36只；此外，還統計了湖羊2胎平均產羔數數據，2胎平均產羔數最多的為4只，產羔數最少的為1只。第二胎產羔數和2胎平均產羔數表型數據統計時已剔除數據缺失的個體。統計學結果顯示湖羊第一胎平均產羔數和第二胎平均產羔數差異顯著。

2.2 全基因組關聯分析

經SNP質控之后，在26條常染色體上一共得到1 604 526個SNP標記用于全基因組關聯分析研究。湖羊第一胎產羔數、第二胎產羔數和2胎平均產羔數GWAS結果分別如圖1、圖2和圖3所示。從曼哈頓圖結果來看，沒有SNP位點達到Bonferroni校正的閾值線，表明湖羊產羔數的全基因組關聯分析未發現顯著SNP位點。從第一胎產羔數、第二胎產羔數以及2胎平均產羔數的Q_Q圖結果來看，所有的P值的觀測值都沒有明顯超過期望值，也未找到與湖羊產羔數顯著相關的SNP位點，這與曼哈頓圖結果一致。

圖1 湖羊第一胎產羔數全基因組關聯分析曼哈頓圖(左)和Q_Q圖(右)

圖2 湖羊第二胎產羔數全基因組關聯分析曼哈頓圖(左)和Q_Q圖(右)

圖3 湖羊2胎平均產羔數全基因組關聯分析曼哈頓圖(左)和Q_Q圖(右)

3 討論

3.1 群體分層分析

GWAS研究是基于全基因組范圍內不同位點之間會出現連鎖不平衡的現象來確定影響某些表型性狀或數量性狀的基因，而影響連鎖不平衡的因素很多，如遺傳漂變、群體分層等，若群體出現分層現象會導致關聯分析的結果出現假陽性而不可靠[15-16]。根據本研究得到的Q_Q圖結果可以看出，P值觀察值和預測值基本相同。同時本研究還根據湖羊不同胎次產羔數GWAS結果的P值計算了實際的膨脹系數，其中，第一胎產羔數GWAS結果的實際膨脹系數為0.999；第二胎產羔數GWAS結果的實際膨脹系數為1.007；2胎平均產羔數GWAS結果的實際膨脹系數為0.999，實際膨脹系數與預期膨脹系數(1)很接近，說明湖羊樣本未出現群體分層的現象，本研究所得到的GWAS結果是可靠的，所使用的混合線性模型是合理的。

3.2 全基因組關聯分析

通過查閱現有的文獻得知，目前還未有與湖羊產羔數性狀的GWAS研究報道。GWAS已經是應用非常普遍的功能基因篩選方法，但仍存在一些缺陷，比如檢測結果中得不到顯著性的SNP位點，或者得到的假陽性的顯著性結果[17]。出現這樣的結果可能是由于以下因素導致：①GWAS研究群體的規模小，導致檢驗功效較低；②SNP標記的密度低或者標記分布不均勻；③選擇的統計模型不適合；④表型值檢測不準確；⑤目標性狀遺傳力低且受微效多基因控制。

本研究使用的是湖羊二代重測序數據，SNP密度很高，從曼哈頓圖中的SNP密度分布圖可以看出，鑒定到的所有SNP都均勻分布，未出現分布不均勻的情況。Q_Q圖結果也顯示選擇的混合線性模型是適合的。然而，湖羊產羔數的全基因組關聯分析的結果未鑒定到顯著性的SNP位點，推測可能由以下原因所導致：①湖羊的繁殖性狀的遺傳力低，只有0.125，并且繁殖性狀是一種數量性狀，受微效多基因控制[18]，導致單個基因的效應對湖羊整個產羔數貢獻不大，達不到GWAS研究檢測的顯著性閾值，從而不能通過GWAS的分析方法挖掘出來；②GWAS分析中使用的群體規模應該越大越好，群體大小是制約GWAS分析檢驗功效的第一要素，群體越大有助于陽性結果的檢出。