王慧芳，周光現(xiàn)，孫永峰，陳新企，雷栩斌，張銳毅，賈汝敏*，趙志輝*

(1. 廣東海洋大學(xué)濱海農(nóng)業(yè)學(xué)院，湛江 524088； 2. 吉林農(nóng)業(yè)大學(xué)動物科學(xué)技術(shù)學(xué)院，長春 130118)

中國是世界上家鵝品種最多、馴養(yǎng)歷史最悠久的國家之一。由于我國幅員遼闊，各地自然生態(tài)條件復(fù)雜多樣，不同時期的經(jīng)濟(jì)文化背景不同，對鵝的選擇利用目的不同，逐步形成了具有不同遺傳特性和生產(chǎn)性能的地方品種，分布遍及全國[1-4]。目前，隨著鵝相關(guān)產(chǎn)品需求水平的增長，低繁殖力成為了制約鵝產(chǎn)業(yè)發(fā)展的主要因素[5-6]。影響鵝產(chǎn)蛋性能的因素是多面的，如遺傳基因、生理狀況、環(huán)境條件、飼料營養(yǎng)和使用年限等，各種應(yīng)激(如停電停水、噪音等)都可顯著降低鵝的產(chǎn)蛋量[7-10]。獅頭鵝起源于鴻雁，是我國第一批受保護(hù)的畜禽品種資源之一，因長期受人類和自然環(huán)境的影響，有著高度的遺傳多樣性，其具有耐粗飼、生長速度快、抗逆性強(qiáng)、耗料少、肉質(zhì)好等優(yōu)良特性，被譽(yù)為“世界鵝王”[11-12]。獅頭鵝成年種公鵝體重為10～12 kg，母鵝為9～10 kg，母鵝年產(chǎn)蛋25～32枚[13-14]。吉林白鵝是由吉林農(nóng)業(yè)大學(xué)科技工作者歷經(jīng)20余年精心培育而成，其成年公鵝體重4.1 kg，母鵝3.2 kg，年產(chǎn)蛋數(shù)達(dá)100枚以上，具有產(chǎn)蛋多、耐粗飼、抗病力強(qiáng)等特點(diǎn)[15]。獅頭鵝與吉林白鵝產(chǎn)蛋量相差懸殊，適合作為探索鵝繁殖性狀相關(guān)研究的素材。

隨著基因組學(xué)技術(shù)的發(fā)展，使得大規(guī)模、高通量檢測獅頭鵝基因組內(nèi)的變異位點(diǎn)，篩選產(chǎn)蛋相關(guān)通路以及相關(guān)產(chǎn)蛋基因成為可能[16-18]。Indel (insertion deletion) 是影響基因組中生產(chǎn)繁殖性能的關(guān)鍵遺傳標(biāo)記，研究發(fā)現(xiàn)，生物學(xué)基因組中Indel起著關(guān)鍵性的作用[18-20]。Indel變異一般比SNP變異少，但同樣反映樣品與參考基因組之間的差異，并且編碼區(qū)的Indel會引起移碼突變，導(dǎo)致基因功能上的變化[21-22]。本研究采用第三代重測序技術(shù)，獲得了獅頭鵝插入/缺失(Indel)位點(diǎn)，經(jīng)過數(shù)據(jù)的比較分析發(fā)現(xiàn)，大量基因與獅頭鵝的產(chǎn)蛋性能相關(guān)。經(jīng)過數(shù)據(jù)分析發(fā)現(xiàn)，GnRH信號通路、雌激素信號通路、催產(chǎn)素信號通路、催乳素信號通路以及孕酮介導(dǎo)的卵母細(xì)胞成熟通路與獅頭鵝產(chǎn)蛋性狀相關(guān)。利用比較基因組學(xué)篩選鑒定兩種鵝群體中特異性Indel位點(diǎn)[23]，為獅頭鵝的群體遺傳分析提供基礎(chǔ)，并進(jìn)一步篩選產(chǎn)蛋相關(guān)基因，為獅頭鵝生產(chǎn)研究提供新的思路。

1 材料與方法

1.1 試驗材料

本研究選取健康成年獅頭鵝母鵝(年均產(chǎn)蛋25～32枚) 和吉林白鵝母鵝(年均產(chǎn)蛋≥100枚)各5只(R1～R5：獅頭鵝；R6～R10：吉林白鵝)，靜脈翅下采血3 mL，加入抗凝劑，置于-20 ℃保存?zhèn)溆谩?/p>

1.2 血液DNA提取

血液樣品DNA使用美基生物公司的HiPure Blood DNA Mini Kit試劑盒進(jìn)行提取。經(jīng)1%凝膠電泳檢測DNA完整性，并使用Nanodrop檢測其純度。合格樣品DNA送北京百邁客公司測序。

1.3 全基因組重測序

1.3.1 文庫構(gòu)建與測序 樣品基因組DNA檢測合格后，用機(jī)械打斷的方法(超聲波)將DNA片段化，然后進(jìn)行片段純化，末端修復(fù)，3′端加A以及連接測序接頭。通過瓊脂糖凝膠電泳選擇片段的大小，PCR擴(kuò)增目的片段形成測序文庫，對質(zhì)檢合格的文庫進(jìn)行Illumina測序。

1.3.2 測序數(shù)據(jù)質(zhì)控 測序數(shù)據(jù)通過公式Qphred =-10lg(e)(Qphred為堿基質(zhì)量值；e為測序錯誤率)轉(zhuǎn)化可以得到堿基測序的錯誤率。在堿基識別過程中計算Phred數(shù)值(一種預(yù)測堿基判別發(fā)生錯誤的概率模型)，其對應(yīng)關(guān)系如表1所示。

表1 預(yù)測堿基判別發(fā)生錯誤概率結(jié)果

1.3.3 低質(zhì)量數(shù)據(jù)過濾 原始測序序列(Raw Reads)含有帶接頭的以及低質(zhì)量的Reads，需對Raw Reads進(jìn)行過濾得到Clean Reads后，方可用于后續(xù)的信息分析。數(shù)據(jù)過濾的主要步驟：1) 去除帶接頭(adapter)的reads；2) 過濾N含量超過10%的reads；3) 去除質(zhì)量值低于10的堿基超過50%的reads。

1.4 Small Indel的檢測與注釋

將在參考基因組上定位后的Clean Reads結(jié)果進(jìn)行分析，檢測是否存在Small Indel變異。使用GATK軟件工具包檢測Small Indel插入與缺失；使用Picard(http://sourceforge.net/projects/ picard/)過濾冗余reads(MarkDuplicates)，以保證檢測結(jié)果的準(zhǔn)確性；然后使用GATK的HaplotypeCaller(局部單體型組裝)算法進(jìn)行Indel的變異檢測；通過SnpEff軟件實現(xiàn)Small Indel的注釋。

1.5 DNA水平的變異基因挖掘

比較篩選獅頭鵝與吉林白鵝兩個群體特有的差異Indel位點(diǎn)，通過BLAST變異基因與GO、COG、KEGG等功能數(shù)據(jù)庫比對，篩選產(chǎn)蛋相關(guān)通路。通過KEGG注釋篩選產(chǎn)蛋相關(guān)候選基因，用以分析基因功能。

2 結(jié) 果

2.1 數(shù)據(jù)質(zhì)量控制

由電泳檢測結(jié)果(圖1)可以看出，DNA純度以及完整性良好。另外檢測得出OD260 nm/OD280 nm值均在1.88左右，其濃度、體積和質(zhì)量平均值分別為：173.4 ng·μL-1、99.1 μL和18.7 μg。

通過對原始數(shù)據(jù)各堿基比例分布和堿基質(zhì)量分布的分析來進(jìn)行質(zhì)量控制，得到圖2。從圖中可以看出，AT曲線以及GC曲線均重合，證明AT堿基以及GC堿基比例均相等，無明顯偏向性，曲線較平緩，未發(fā)生分離，測序結(jié)果正常，該數(shù)據(jù)可用于下一步分析。

通過對獅頭鵝和吉林白鵝全基因組重測序，原始數(shù)據(jù)過濾處理后共得到350.35 Gb的Clean Data，各樣本的Raw data在71 553 908～208 410 322 Mb之間， Q20與Q30平均值分別達(dá)到96.74%、92.19%。GC含量在43.67%～45.64%之間(表2)，結(jié)果表明，測序質(zhì)量良好，可進(jìn)行后續(xù)分析。

R1～R5. 獅頭鵝；R6～R10. 吉林白鵝。下同 R1-R5. Shitou goose；R6-R10. Jilin white goose. The same as below圖1 瓊脂糖凝膠電泳檢測結(jié)果Fig.1 The detection results of agarose gel electrophoresis

2.2 與參考基因組的比對分析

參考基因組為Anser cygnoides (Swan goose)，其版本為AnsCyg_PRJNA183603_v1.0；下載網(wǎng)址：https://www.ncbi.nlm.nih.gov/genome/?term=Anser+cygnoides，參考基因組的大小為1 134 Mb(表3)。由表4可知，所有樣本數(shù)據(jù)與參考基因組的比對率在93.18%～94.96%之間，平均比對率為94.15%，平均覆蓋深度為26×，1×覆蓋度(至少有1個堿基的覆蓋)在99.47%以上，5×覆蓋度(至少有5個堿基的覆蓋)在96.84%以上。結(jié)果表明，測序物種同參考基因組相似度高，證明測序準(zhǔn)確，可進(jìn)行后續(xù)的相關(guān)分析。

N. 測序中識別不出的堿基 N. The unrecognized bases in sequencing圖2 樣品各堿基比例分布Fig.2 The proportion of each base of the samples

表2 測序數(shù)據(jù)質(zhì)量表

表3 參考基因組信息

表4 測序深度及覆蓋度統(tǒng)計

2.3 Indel序列長度及分布分析

利用GATK從10只試驗鵝與參考基因組比對結(jié)果中檢測Indel變異(小于50 bp堿基片段的插入與缺失)。然后對Indel序列的分布位置進(jìn)行注釋及分析(表5)。

通過SnpEff軟件注釋分析Indel位點(diǎn)，分別從10只試驗鵝中發(fā)現(xiàn)了786 695、782 331、766 853、775 202、718 792、717 274、708 825、738 934、753 315、 739 991個Indel位點(diǎn)(平均有748 821.2個)，其數(shù)量和分布與參考基因組相似。在這10只試驗鵝中，Indel位點(diǎn)主要分布在內(nèi)含子，吉林白鵝與獅頭鵝Indel位點(diǎn)在內(nèi)含子中分別占有54.18%和52.12%；其次是基因間區(qū)，兩種鵝分別占有36.05%與33.91%(圖3)。大部分的Indel位點(diǎn)分布在內(nèi)含子與基因間區(qū)中，CDS區(qū)中分布的Indel位點(diǎn)最少。

圖3 吉林白鵝和獅頭鵝相同Indels在基因組上的分布Fig.3 Genomic distribution of the same Indels of Jilin white goose and Shitou goose

表5 Indel 位點(diǎn)的檢測及注釋

在這10只試驗鵝全基因組中，Indel位點(diǎn)的長度分布趨勢也基本一致。對群體Indel突變位點(diǎn)的堿基長度進(jìn)行統(tǒng)計并分析，得到圖4，發(fā)現(xiàn)全基因組中長度為1 bp的Indel位點(diǎn)數(shù)量超過310 000個(約占全基因組所有Indel的43%)；長度為2～4 bp的Indel位點(diǎn)數(shù)量為2 529 498個；長度為5～9 bp的Indel位點(diǎn)數(shù)量為860 657個；長度大于等于10 bp 的Indel位點(diǎn)數(shù)量共808 872個。并且在編碼區(qū)中，長度為1 bp的Indel位點(diǎn)數(shù)量共22 754個(約占編碼區(qū)所有Indel的33%)；長度為2～4 bp的Indel位點(diǎn)數(shù)量為22 426個；長度為5～9 bp的Indel位點(diǎn)數(shù)量為10 080個；長度大于等于10 bp的Indel位點(diǎn)數(shù)量共13 491個。另外，圖4說明了在片段變異10 bp長度的范圍內(nèi)，隨著堿基片段長度增加，發(fā)生Indel變異的數(shù)量隨之減少，且數(shù)量波動較大。

圖4 全基因組和編碼區(qū)的Indel長度分布情況Fig.4 Distribution map of Indel length in the population’s entire genome and coding region

10只試驗鵝的Indel位點(diǎn)在CDS區(qū)的分布趨勢和數(shù)量也與參考基因組相似，分別有7 090、7 026、 7 022、7 060、6 696、6 649、6 651、6 862、6 868和6 827個Indel位點(diǎn)(表6)。其中共有11 737個Indel位點(diǎn)導(dǎo)致了插入或刪除。

表6 Indel變異在全基因組和編碼區(qū)的統(tǒng)計情況

2.4 變異基因注釋分析

通過GO注釋分析富集了到細(xì)胞組分(cellular component)17個、分子功能(molecular function)20個和生物過程(biological process)23個，詳見表7。

通過KEGG注釋分析富集到253條通路，其中篩選出與產(chǎn)蛋相關(guān)的通路，分別為GnRH信號通路、雌激素信號通路、催乳素信號通路、催產(chǎn)素信號通路以及孕酮介導(dǎo)的卵母細(xì)胞成熟通路。通過比對測序得到的群體基因注釋列表發(fā)現(xiàn)了36個產(chǎn)蛋相關(guān)變異基因，詳見表8。這為獅頭鵝后期高產(chǎn)基因的篩選奠定了基礎(chǔ)。

表7 變異基因功能富集分析的相關(guān)注釋

(轉(zhuǎn)下頁 Carried forward)

3 討 論

獅頭鵝是我國乃至世界著名的大型水禽品種，也是我國南方養(yǎng)殖范圍較廣的經(jīng)濟(jì)家禽，其產(chǎn)蛋期一般為9月至次年4月[24-25]。吉林白鵝是我國北方養(yǎng)殖較廣的經(jīng)濟(jì)家禽，其產(chǎn)蛋旺季在2～6月，研究表明，兩種鵝日照長短不一樣[26]。甄霆等[27]研究發(fā)現(xiàn)，在浙東白鵝的繁殖季節(jié)通過調(diào)節(jié)光照可影響FSH和INH的表達(dá)，從而延長其產(chǎn)蛋時間，縮短抱窩和就巢時間，增加鵝的產(chǎn)蛋量，說明鵝的繁殖性能受光照的影響，合理地控制光照可延長鵝的產(chǎn)蛋期，縮短休產(chǎn)期，提高鵝的產(chǎn)蛋性能。姚穎等[28-31]研究發(fā)現(xiàn)，浙東白鵝約83.04%的產(chǎn)蛋是在有光照的情況下，僅有16.96%的產(chǎn)蛋是在無光照的情況下，說明產(chǎn)蛋行為與光照密切相關(guān)[28]。產(chǎn)蛋性狀的遺傳力較低，其由微效多基因控制，三代測序技術(shù)的發(fā)展為研究獅頭鵝產(chǎn)蛋性狀提供了新的思路。由于獅頭鵝較強(qiáng)的就巢性和產(chǎn)蛋季節(jié)性導(dǎo)致其繁殖力較低，從而造成獅頭鵝產(chǎn)蛋性能低下。提高獅頭鵝的產(chǎn)蛋性能一直是研究的主要目標(biāo)之一。

本研究利用重測序技術(shù)對5只獅頭鵝和5只吉林白鵝全基因組測序后共得到350.35 Gb數(shù)據(jù)，其中Q20與Q30平均值分別達(dá)到96.74%、92.19%，測序數(shù)據(jù)結(jié)果可靠。在此基礎(chǔ)上，通過注釋分析得到Indel位點(diǎn)平均有748 821.2個變異。大部分(53%左右)的Indel位點(diǎn)分布在內(nèi)含子中，蛋白質(zhì)編碼區(qū)中分布的Indel位點(diǎn)最少(11 737個)。研究還發(fā)現(xiàn)，Indel位點(diǎn)的長度分布趨勢也基本一致，10只 試驗鵝的Indel位點(diǎn)在CDS區(qū)的分布趨勢和數(shù)量也與參考基因組相似。在此基礎(chǔ)上，通過plink軟件對SNPs進(jìn)行主成分分析(PCA)發(fā)現(xiàn)(數(shù)據(jù)未顯示)，吉林白鵝與獅頭鵝兩個群體能夠有效地分離。從SNP層次觀測，兩種鵝具有較大的差異性，這種差異性包含潛在主導(dǎo)其不同性狀的SNP位點(diǎn)。吉林白鵝體型小，年產(chǎn)蛋量高，而獅頭鵝體型大，年產(chǎn)蛋量低。張澤賓[32]對22只野鴨和56只家鴨全基因組測序得到310萬個Indel變異位點(diǎn)，對Indel注釋發(fā)現(xiàn)，絕大多數(shù)的變異位于基因間區(qū)及內(nèi)含子等非編碼區(qū)域，此研究發(fā)現(xiàn)，編碼區(qū)中的Indel突變比率與其它區(qū)間相比較低。劉杰[33]分別選取北京油雞和科寶肉雞中高采食量個體48只和低采食量個體48只進(jìn)行全基因組重測序，得到的Indel變異位點(diǎn)在內(nèi)含子和基因間區(qū)最多，編碼區(qū)較少。章雙杰等[34]對4種鴨群進(jìn)行全基因組重測序分析，通過檢測到的Indel位點(diǎn)發(fā)現(xiàn)，其Indel變異大部分存在于基因間或是基因內(nèi)的內(nèi)含子區(qū)，證明編碼區(qū)的變異相對較少。以上研究結(jié)果均與本試驗Indel變異分布結(jié)果一致，非編碼區(qū)的Indel變異多于編碼區(qū)，這一結(jié)果符合生物學(xué)特性。對差異分析篩選的候選Indel標(biāo)記的作用可進(jìn)一步研究，為鵝生長性狀遺傳機(jī)理的相關(guān)研究奠定理論基礎(chǔ)。

表8 產(chǎn)蛋相關(guān)通路及變異基因

在此基礎(chǔ)上，通過數(shù)據(jù)庫比對分析進(jìn)一步篩選產(chǎn)蛋相關(guān)通路，發(fā)現(xiàn)在促性腺激素釋放激素(GnRH)信號通路、催乳素(PRL)信號通路、催產(chǎn)素(OT)信號通路和孕酮介導(dǎo)的卵母細(xì)胞成熟通路均有Indel位點(diǎn)變異。GnRH作為HPG軸的關(guān)鍵信號分子與繁殖周期的變化呈相關(guān)性。產(chǎn)蛋性狀是由多個相關(guān)基因共同表達(dá)和相互作用來調(diào)控的，并受多種外源因素影響[35-37]。常見的產(chǎn)蛋相關(guān)主要基因有促性腺激素釋放激素(GnRH)[38]、促性腺激素抑制激素(GnIH)[39]、催乳素(PRL)[40]、促卵泡激素(FSH)[41]、雌激素受體(ESR)[42]等。張克山等[43]研究表明，在產(chǎn)蛋前期或高峰期，四川白鵝的GnRHmRNA表達(dá)量高于休產(chǎn)期，并指出在調(diào)控鵝產(chǎn)蛋性狀中GnRH信號通路發(fā)揮重要作用。催產(chǎn)素通路在鵝的繁殖、就巢和產(chǎn)蛋等方面有重要作用。本課題組成功克隆了獅頭鵝PRL基因并進(jìn)行了原核表達(dá)，研究發(fā)現(xiàn)，獅頭鵝與其他禽類PRL具有高度保守性[44]。催乳素通路是細(xì)胞生長和分化過程中調(diào)控鵝性成熟的主要通路，同時也是鵝產(chǎn)蛋性狀的候選遺傳標(biāo)記基因之一[45]。PRL基因的多態(tài)性主要與產(chǎn)卵性狀相關(guān)。催乳素通路在鵝卵巢發(fā)育和產(chǎn)卵過程中也發(fā)揮重要作用。孕酮介導(dǎo)的卵母細(xì)胞成熟通路對卵泡發(fā)育起調(diào)節(jié)作用。卵泡和排卵取決于促性腺激素的釋放。在家禽生產(chǎn)中，孕酮的分泌與卵巢的生長和狀態(tài)緊密相關(guān)。在排卵過程中，孕酮作為正反饋信號在排卵前觸發(fā)GnRH/LH通路[46]?？傊?，本研究篩選出的產(chǎn)蛋相關(guān)變異通路和調(diào)節(jié)途徑，以及進(jìn)一步篩選出的與產(chǎn)蛋相關(guān)的候選基因，為研究獅頭鵝繁殖性能提供有效基礎(chǔ)和依據(jù)。

4 結(jié) 論

通過全基因組測序后注釋分析Indel變異位點(diǎn)，發(fā)現(xiàn)大部分的Indel位點(diǎn)分布在內(nèi)含子中，蛋白質(zhì)編碼區(qū)中分布的Indel位點(diǎn)最少；另外發(fā)現(xiàn)，Indel位點(diǎn)長度分布趨勢基本一致，在CDS區(qū)的分布趨勢和數(shù)量與參考基因組相似，非編碼區(qū)的Indel變異數(shù)量多于編碼區(qū)，這一結(jié)果符合生物學(xué)特性。比對測序得到的群體基因注釋列表篩選出了產(chǎn)蛋相關(guān)變異基因，本研究突破了傳統(tǒng)育種的局限性，有效地為低遺傳力產(chǎn)蛋性狀的遺傳進(jìn)展加快了步伐，提供研究鵝繁殖性能重要的分子標(biāo)記理論數(shù)據(jù)，可以在實際選育工作中提高鵝養(yǎng)殖產(chǎn)業(yè)的經(jīng)濟(jì)效益。