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

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

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

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

1 材料與方法

1.1 樣本采集及測序

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

1.2 變異檢測及SNP質(zhì)控

使用綿羊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的鑒定和初步過濾，將結(jié)果輸出到VCF文件保存。最后使用VCFtools[13](version 0.1.17)對SNP進行質(zhì)控，并保留符合以下標準的SNP：①所有SNP位點的平均測序深度大于2×；②只保留常染色體上的SNP；③SNP最大缺失率5%；④只保留二等位基因；⑤最小等位基因頻率大于0.05；⑥哈迪溫伯格平衡P值大于1×10-6。

1.3 GWAS分析

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

y=μ+Xb+u+e，

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

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

2 結(jié)果與分析

2.1 表型數(shù)據(jù)統(tǒng)計

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

2.2 全基因組關聯(lián)分析

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

圖1 湖羊第一胎產(chǎn)羔數(shù)全基因組關聯(lián)分析曼哈頓圖(左)和Q_Q圖(右)

圖2 湖羊第二胎產(chǎn)羔數(shù)全基因組關聯(lián)分析曼哈頓圖(左)和Q_Q圖(右)

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

3 討論

3.1 群體分層分析

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

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

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

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