東湖雜交羊產(chǎn)羔性狀的GWAS分析及候選基因GRID2的驗(yàn)證

趙雪洋,李丹妮,王鈺晨,郭 磊,王 麗,黃 杰,焦儀強(qiáng),安小鵬,張希云, 張 磊*,宋宇軒*

(1.西北農(nóng)林科技大學(xué)動物科技學(xué)院,楊凌 712100;2.永昌縣畜牧獸醫(yī)站,金昌 737200; 3.甘肅元生農(nóng)牧科技有限公司,金昌 737100)

綿羊產(chǎn)羔性狀是指母羊每胎次的產(chǎn)羔數(shù),它是衡量綿羊繁殖力的重要指標(biāo)[1]。提高綿羊產(chǎn)羔數(shù)可降低飼養(yǎng)成本,提高生產(chǎn)效益,促進(jìn)綿羊養(yǎng)殖業(yè)的高效發(fā)展。綿羊產(chǎn)羔性狀受多種因素的影響,包括遺傳因素和環(huán)境因素[2-3]。遺傳因素主要是指不同品種間的產(chǎn)羔數(shù)差異,環(huán)境因素主要是指飼養(yǎng)管理、配種控制、疾病防治等方面[4-7]。為了提高綿羊產(chǎn)羔性狀,需要采取多種措施,包括優(yōu)化飼料營養(yǎng)、調(diào)整配種時間、加強(qiáng)孕期護(hù)理等[8]。此外,利用分子遺傳技術(shù)對綿羊產(chǎn)羔數(shù)相關(guān)基因進(jìn)行挖掘和標(biāo)記輔助選擇,也是一種有效的方法[9-10]。

全基因組關(guān)聯(lián)分析(genome wide association study,GWAS)是一種常用的分子遺傳學(xué)方法,旨在鑒定復(fù)雜性狀的遺傳基礎(chǔ)。通過GWAS分析,研究者可以深入了解綿羊的遺傳特征和種群結(jié)構(gòu),發(fā)現(xiàn)與性狀相關(guān)的基因和變異,為綿羊育種和遺傳改良提供重要的信息和指導(dǎo)[11-14]。近年來,GWAS已成為提高綿羊繁殖性能的有效方法之一[15]?？梢栽谌蚪M范圍內(nèi)鑒定與繁殖性能相關(guān)的基因,進(jìn)而深入了解繁殖性狀的遺傳機(jī)制和調(diào)控網(wǎng)絡(luò)[10,16-18]。這些發(fā)現(xiàn)有望為改良和優(yōu)化綿羊繁殖性能提供理論指導(dǎo)和科學(xué)依據(jù)。同時,GWAS還可以幫助篩選和鑒定出對繁殖性能有積極影響的遺傳變異,進(jìn)而為選育高產(chǎn)優(yōu)質(zhì)綿羊提供有力的選育手段[19]?？偟膩碚f,GWAS已經(jīng)成為綿羊育種中不可或缺的工具之一。隨著技術(shù)的不斷進(jìn)步和方法的不斷創(chuàng)新,GWAS在綿羊育種中的應(yīng)用前景將更加廣闊。

現(xiàn)代數(shù)量遺傳學(xué)在動物育種實(shí)踐中的應(yīng)用取得了巨大成功,如今它對一些遺傳力高且呈連續(xù)性正態(tài)分布的數(shù)量性狀仍然是必不可少的選擇方法。但對于低遺傳力性狀,如母畜的產(chǎn)仔數(shù)等,選擇反應(yīng)并不理想。因此,可以從分子遺傳水平上找到性狀的遺傳差異或與數(shù)量性狀連鎖的遺傳標(biāo)記,從而實(shí)現(xiàn)真正的基因型選擇[20-22]。

本試驗(yàn)以東湖雜交奶綿羊?yàn)檠芯繉ο?對168只東湖雜交奶綿羊進(jìn)行高深度全基因組重測序,再根據(jù)全基因組關(guān)聯(lián)分析技術(shù)確定與產(chǎn)羔性狀相關(guān)的關(guān)鍵基因,最后挑選關(guān)鍵基因的變異位點(diǎn)采用KASP分型技術(shù)進(jìn)行SNP分型檢測,并將其與產(chǎn)羔性狀進(jìn)行關(guān)聯(lián)分析,探究其潛在的育種價(jià)值。

1 材料與方法

1.1 試驗(yàn)材料

試驗(yàn)動物來自西北農(nóng)林科技大學(xué)金昌奶綿羊試驗(yàn)示范基地。選擇體況良好、平均年齡一歲半的東湖級雜二代母羊748只,頸靜脈采集5 mL血液于抗凝管中,并于-80 ℃保存。其中168只進(jìn)行全基因組重測序,580只進(jìn)行SNP分型檢測并統(tǒng)計(jì)試驗(yàn)羊群不同胎次的產(chǎn)羔數(shù)據(jù)。

1.2 DNA提取與建庫測序

本研究通過標(biāo)準(zhǔn)苯-酚氯仿法提取基因組DNA。將檢測合格的DNA樣品用超聲波破碎儀打斷至350 bp大小,再經(jīng)末端修復(fù)、加A尾、加測序接頭、純化和PCR擴(kuò)增等步驟進(jìn)行文庫的制備。文庫檢測合格后,使用Illumina高通量測序平臺Nova-Seq 6000進(jìn)行測序,生成150 bp的成對末端reads,共獲得47 616 147個原始位點(diǎn)。

1.3 質(zhì)量控制與變異檢測

本研究利用fastp軟件對原始的fastq格式文件進(jìn)行過濾及質(zhì)控。質(zhì)控后的數(shù)據(jù)使用BWA-MEN對比到Oar_rambouillet_v1.0上,重復(fù)reads使用Picard軟件進(jìn)行排序。采用genome analysis toolkit (GATK, version 3.8)提取基因組高質(zhì)量變異,獲取變異數(shù)據(jù)的VCF集合。使用GATK進(jìn)行過濾(參數(shù):QD<2.0, FS>60.0, MQ<40.0, MQRankSum<-12.5, ReadPosRankSum<-8.0 and SOR>3.0)。使用GATK軟件的SelectVariants模塊分別提取SNP和InDel數(shù)據(jù)集合(參數(shù):-restrict-alleles-to BIALLELIC,--remove-unused-alternates)。為了降低結(jié)果的假陽性,本研究將PLINK 1.9軟件[23]獲取的主成分分析結(jié)果作為協(xié)變量加入,此外還根據(jù)PLINK--missing 0.01參數(shù)和—maf 0.01參數(shù)進(jìn)行了SNP位點(diǎn)的過濾,共得到27 750 039個位點(diǎn)。

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

先使用GEMMA的-gk 2參數(shù)生成親緣關(guān)系矩陣,然后結(jié)合生成的親緣關(guān)系矩陣使用混合線性模型(MLM)進(jìn)行GWAS分析(參數(shù):gemma-bfile-k-lmm 1-p)。MLM模型結(jié)合了固定效應(yīng)(配種方式和產(chǎn)羔季節(jié))和隨機(jī)效應(yīng)。校正閾值線的確定為0.05/n(總位點(diǎn)數(shù))。根據(jù)輸出結(jié)果,后續(xù)使用R語言包CMplot進(jìn)行曼哈頓圖和QQ(Quantiles-Quantiles)圖的繪制。

1.5 顯著位點(diǎn)的注釋

使用SnpEff軟件構(gòu)建綿羊參考基因組的本地注釋數(shù)據(jù)庫,并使用該軟件對本研究構(gòu)建的VCF集合進(jìn)行注釋,使用自己編寫的perl腳本從注釋后的VCF集合中提取注釋信息。

1.6 注釋基因的富集分析

本研究中使用北京大學(xué)構(gòu)建的KOBAS 3.0數(shù)據(jù)庫(http://bioinfo.org/kobas/retrieve/?taskid=5f7a3a18b1de4e53a517ffc7c7c86cfb)對注釋得到的基因進(jìn)行富集分析,以探究GWAS顯著信號區(qū)域的基因在生物學(xué)網(wǎng)絡(luò)中的作用和相互作用。

1.7 連鎖不平衡分析

使用LDBlockShow軟件進(jìn)行連鎖不平衡熱圖的繪制,突出顯著性最強(qiáng)SNPs周邊的連鎖區(qū)域。

1.8 基因型檢測

從NCBI(national center for biotechnology information)獲取GRID2的序列信息(NC_040257.1),根據(jù)g.35685218A>G,和g.35685499C>T兩個突變的位置,分別設(shè)計(jì)兩條不同的上游引物和一條相同的下游引物。使用KASP進(jìn)行特異性識別SNP基因型。

1.9 數(shù)據(jù)統(tǒng)計(jì)與分析

根據(jù)等位基因頻率、基因型頻率、純合度和雜合度等參數(shù)的計(jì)算原理,自行編寫Excel函數(shù)進(jìn)行群體遺傳學(xué)參數(shù)的計(jì)算,同時檢驗(yàn)群體是否符合哈代-溫伯格平衡原理。使用SPSS 20.0軟件分析基因型與產(chǎn)羔性狀之間的關(guān)系。線性分析模型為Y=μ+a+b+c+e,其中Y為性狀的觀察值,μ為均值,a為年季節(jié)效應(yīng),b為基因型固定效應(yīng),c為配種方式固定效應(yīng),e為隨機(jī)誤差。

2 結(jié) 果

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

使用GEMMA軟件的混合線性模型進(jìn)行全基因組SNP的關(guān)聯(lián)分析,設(shè)置顯著性水平的校正閾值線為1.80×10-9,共檢測到1 163個達(dá)到全基因組顯著性水平的SNPs位點(diǎn),所有達(dá)到顯著性水平的位點(diǎn)均分布于綿羊6號染色體上,如圖1A所示。通過SnpEff對所有顯著位點(diǎn)進(jìn)行注釋,其中位于外顯子區(qū)域4個,基因間區(qū)708個,內(nèi)含子區(qū)域445個,其他區(qū)域6個。注釋到與產(chǎn)羔性狀相關(guān)的基因有BMPR1B、PDLIM5、PDHA2、UNC5C、GRID2、NFKB1、TSPAN5、STPG2、PPP3CA等。其中處于顯著性前十的SNPs注釋到了BMPR1B、GRID2、PDLIM5、PDHA2基因。

基于GATK提取的InDel集合利用GEMMA混合線性模型進(jìn)行全基因組InDel關(guān)聯(lián)分析,與SNP的結(jié)果相比,信號雜亂(圖2),不具有較好參考價(jià)值,這也符合前人研究中多胎性狀的主效位點(diǎn)主要為SNP的結(jié)果。

圖1 SNP全基因組關(guān)聯(lián)分析的曼哈頓圖(A)和QQ-plot圖(B)Fig.1 Manhattan map (A)and QQ-plot(B) of SNP genome-wide association study

2.2 KEGG富集分析

通過KOBAS對GWAS顯著位點(diǎn)的注釋基因進(jìn)行KEGG富集分析后,共得到了17個校正P值達(dá)到顯著性水平(correctedP<0.05)的通路,如圖3所示,表1列出了KEGG富集分析的前20個通路。其中有多個信號通路已被證明與生殖過程有關(guān)。如Wnt信號通路、axon guidance、Th1和Th2細(xì)胞分化通路等。

圖3 KEGG富集分析氣泡圖Fig.3 Bubble diagram of KEGG enrichment analysis

2.3 連鎖不平衡熱圖

根據(jù)BMPR1B基因在綿羊參考基因組的位置,選擇信號最高點(diǎn)6∶34198746上下游10 kb區(qū)域、FecB(BMPR1B基因上A746G突變)突變位置上下游10 kb區(qū)域進(jìn)行三角熱圖的繪制,如圖4A、4B所示,并沒有發(fā)現(xiàn)顯著強(qiáng)連鎖的區(qū)域。此外,為了研究候選基因GRID2兩個SNPs位點(diǎn)和周邊區(qū)域的連鎖情況,選擇g.35685218A>G和g.35685499C>T兩個位置上下游10 kb區(qū)域進(jìn)行三角熱圖繪制,發(fā)現(xiàn)周邊SNP連鎖情況相對較強(qiáng),如圖4C所示。

2.4 GRID2基因g.35685218A>G和g.35685499C>T的多態(tài)位點(diǎn)檢測

通過KASP基因分型技術(shù)對GRID2基因的兩個多態(tài)位點(diǎn)在560只東湖雜交奶綿羊群體中進(jìn)行分型檢測(圖5),其中g(shù).35685218A>G位點(diǎn)共檢測到野生型(AA)的個體有211個,雜合型(AG)266個,突變純合型(GG)83個。g.35685499C>T位點(diǎn)共檢測到野生型(CC)的個體有220個,雜合型(CT)256個,突變純合型(TT)84個。

2.5 GRID2基因g.35685218A>G和g.35685499C>T位點(diǎn)遺傳多樣性參數(shù)分析

本研究對560只東湖雜交奶綿羊GRID2的兩個多態(tài)位點(diǎn)進(jìn)行遺傳多樣性參數(shù)的計(jì)算,兩個位點(diǎn)的等位基因頻率和和基因型頻率均十分相似(表2)。g.35685218A>G野生型(AA)基因型頻率為0.377,雜合型(AG)基因型頻率為0.475,突變純合型(GG)基因型頻率為0.148,A基因頻率為0.614,G基因頻率為0.386。g.35685499C>T野生型(CC)基因型頻率為0.393,雜合型(CT)基因型頻率為0.457,突變純合型(TT)基因型頻率為0.150,C基因頻率為0.621,T基因頻率為0.379。本研究中,g.35685218A>G位點(diǎn)的實(shí)際數(shù)值接近預(yù)測值(χ2<0.05),而g.35685499C>T位點(diǎn)的實(shí)際數(shù)值偏離預(yù)測值(χ2>0.05)。兩個位點(diǎn)的多態(tài)信息含量均為中度多態(tài)性(0.250.05的哈代-溫伯格平衡狀態(tài),表明這兩個位點(diǎn)均沒有受到過強(qiáng)選擇。

2.6 GRID2基因g.35685218A>G和g.35685499C>T與產(chǎn)羔數(shù)關(guān)聯(lián)分析

為了系統(tǒng)地研究GRID2基因g.35685218A>G和g.35685499C>T兩個位點(diǎn)對產(chǎn)羔數(shù)的影響,本研究將不同胎次產(chǎn)羔數(shù)與基因型進(jìn)行了關(guān)聯(lián)分析,如表3所示。g.35685218A>G位點(diǎn)第一胎的野生型相較于突變雜合型平均產(chǎn)羔數(shù)增加了0.34只(P<0.05),相較于突變純合型增加了0.64只(P<0.05)。第二胎的野生型相較于突變雜合型平均產(chǎn)羔數(shù)平均增加了0.57只(P<0.05),相較于突變純合型增加了1.11只(P<0.05)。第三胎的野生型相較于突變雜合型平均產(chǎn)羔數(shù)平均增加了0.20只(P<0.05),相較于突變純合型增加了1.06只(P<0.05)。

g.35685499C>T位點(diǎn)第一胎的野生型相較于突變雜合型平均產(chǎn)羔數(shù)增加了0.14只(P<0.05),相較于突變純合型增加了0.50只(P<0.05)。第二胎的野生型相較于突變雜合型平均產(chǎn)羔數(shù)增加了0.41只(P<0.05),相較于突變純合型增加了1.0只(P<0.05)。第三胎的野生型相較于突變雜合型平均產(chǎn)羔數(shù)增加了0.20只(P<0.05),相較于突變純合型增加了1.15只(P<0.05)。

表1 KEGG富集通路Table 1 KEGG enrichment pathway

顏色從白色到紅色,代表連鎖程度從低到高。A.信號最高點(diǎn)6:34198746上、下游10 kb區(qū)域;B.FecB突變位置上、下游10 kb區(qū)域;C.g.35685218A>G和g.35685499C>T兩個位置上、下游10 kb區(qū)域 The color ranges from white to red represent a low to high degree of linkage.A.The 10 kb region upstream and downstream of signal peak 6:34198746; B. The 10 kb region upstream and downstream of the FecB mutation site;C.The 10 kb region upstream and downstream of g.35685218A>G and g.35685499C>T圖4 SNP連鎖不平衡熱圖Fig.4 SNP linkage disequilibrium heat map

圖5 GRID2基因g.35685218A>G(A)和g.35685499C>T(B)的KASP基因分型圖Fig.5 KASP genotyping map of g.35685218A> G (A) and g.35685499C> T (B) in GRID2 gene

3 討 論

本研究對168只東湖雜交奶綿羊的重測序數(shù)據(jù)與產(chǎn)羔性狀進(jìn)行了GWAS分析,共獲得了1 163個達(dá)到顯著性水平的SNPs位點(diǎn)。根據(jù)注釋結(jié)果,發(fā)現(xiàn)了與產(chǎn)羔性狀有關(guān)的多個基因和信號通路。其中,BMPR1B基因中FecB(A746G)錯義突變導(dǎo)致氨基酸變化(Q249R),從而增加羊的排卵率[24],此突變對羊排卵具有加性效應(yīng)[25-27],本課題組前期已經(jīng)在奶綿羊群體中大量檢測了FecB突變,研究結(jié)果顯示,FecB突變顯著影響東湖雜交二代羊的產(chǎn)羔數(shù)[28]。

根據(jù)KEGG富集分析的結(jié)果,同樣有一些和繁殖過程有關(guān)的通路。Wnt信號通路是一個重要的調(diào)控機(jī)制,在細(xì)胞分化、增殖、凋亡和胚胎發(fā)育過程中發(fā)揮著關(guān)鍵作用。研究表明,Wnt信號通路在小鼠和人類卵巢中的活性會隨著動情周期的變化而變化,表明其在生殖周期中起著重要作用[29]。另外一個與繁殖過程相關(guān)的通路是軸突導(dǎo)向通路(axon guidance),這個通路在胚胎發(fā)育期間,對于神經(jīng)系統(tǒng)的形成和連接至關(guān)重要。研究發(fā)現(xiàn),軸突導(dǎo)向通路的一些關(guān)鍵基因在胚胎期間還會表達(dá)于生殖細(xì)胞系中,這些基因?qū)τ谏臣?xì)胞系的發(fā)育也有重要的調(diào)控作用[30]。Th1和Th2細(xì)胞分化通路在免疫調(diào)節(jié)過程中發(fā)揮著關(guān)鍵作用,研究表明,妊娠期間,母體的免疫系統(tǒng)處于一種Th2傾向狀態(tài),Th1/Th2平衡的變化可能與早期流產(chǎn)、胎兒生長受限和先兆子癇等疾病有關(guān)[31]。此外,HIF-1還被認(rèn)為在雌性哺乳動物中調(diào)節(jié)子宮內(nèi)膜生長和修復(fù),從而有助于維持正常的妊娠[32]。

表2 GRID2基因g.35685218A>G和g.35685499C>T的遺傳多態(tài)性參數(shù)Table 2 Genetic polymorphism parameters for the g.35685218A>G and g.35685499C>T loci in GRID2

同時,本研究也發(fā)現(xiàn)與東湖雜交奶綿羊產(chǎn)羔數(shù)相關(guān)的顯著性SNPs位點(diǎn)中有幾個位于GRID2基因,此基因在大腦皮層和小腦中廣泛表達(dá)[33]。GRID2基因位于染色體4q22-q23上,包含約118 000個堿基對,并編碼約3 930個氨基酸,編碼蛋白質(zhì)為Delta-2,在中樞神經(jīng)系統(tǒng)中發(fā)揮重要作用,參與了許多神經(jīng)功能的調(diào)節(jié),包括學(xué)習(xí)、記憶、情緒和感知[34]。有報(bào)道顯示,GRID2可能與動物繁殖性狀有關(guān),一項(xiàng)對31只大足黑山羊產(chǎn)羔數(shù)的全基因組選擇信號研究發(fā)現(xiàn),GRID2可能是影響產(chǎn)羔數(shù)的候選基因,LD分析表明,高產(chǎn)大足黑山羊表現(xiàn)出與低產(chǎn)山羊顯著不同的連鎖不平衡模式[35]。根據(jù)Chen等[36]對兩個不同F(xiàn)ecB基因型小尾寒羊不同生殖周期的垂體組織轉(zhuǎn)錄組研究發(fā)現(xiàn),GRID2在不同基因型羊的垂體組織間差異表達(dá),推測該基因可能介導(dǎo)了一些生殖相關(guān)激素的釋放。此外,GRID2基因也被證明和豬卵母細(xì)胞的增殖和凋亡過程相關(guān)[37]。另外,6號染色體上存在大量的連鎖區(qū)域,可能還存在一些其他類型的變異,例如結(jié)構(gòu)變異(structural variation, SV)、拷貝數(shù)變異(copy number variation, CNV)等,受二代測序數(shù)據(jù)的限制,很難檢測到這些變異,需要進(jìn)一步的挖掘。并且,統(tǒng)計(jì)學(xué)分析發(fā)現(xiàn)GRID2基因g.35685218A>G和g.35685499C>T兩個位點(diǎn)的野生型均處于優(yōu)勢狀態(tài),且兩個位點(diǎn)在胎次增加時均值也增加。

4 結(jié) 論

本研究結(jié)果表明,通過全基因組關(guān)聯(lián)分析確定了6號染色體上1 163個達(dá)到顯著性水平的SNPs位點(diǎn),注釋到候選基因有BMPR1B、PDLIM5、PDHA2、UNC5C、GRID2、NFKB1、TSPAN5、STPG2、PPP3CA。其中GRID2基因2個SNPs(g.35685218A>G和g.35685499C>T)對東湖雜交奶綿羊產(chǎn)羔數(shù)有負(fù)效應(yīng),可以作為分子標(biāo)記輔助選擇的候選位點(diǎn)。