基于全基因組重測序技術(shù)的獅頭鵝Indel標記分析

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

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

中國是世界上家鵝品種最多、馴養(yǎng)歷史最悠久的國家之一。由于我國幅員遼闊，各地自然生態(tài)條件復雜多樣，不同時期的經(jīng)濟文化背景不同，對鵝的選擇利用目的不同，逐步形成了具有不同遺傳特性和生產(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]。獅頭鵝起源于鴻雁，是我國第一批受保護的畜禽品種資源之一，因長期受人類和自然環(huán)境的影響，有著高度的遺傳多樣性，其具有耐粗飼、生長速度快、抗逆性強、耗料少、肉質(zhì)好等優(yōu)良特性，被譽為“世界鵝王”[11-12]。獅頭鵝成年種公鵝體重為10～12 kg，母鵝為9～10 kg，母鵝年產(chǎn)蛋25～32枚[13-14]。吉林白鵝是由吉林農(nóng)業(yè)大學科技工作者歷經(jīng)20余年精心培育而成，其成年公鵝體重4.1 kg，母鵝3.2 kg，年產(chǎn)蛋數(shù)達100枚以上，具有產(chǎn)蛋多、耐粗飼、抗病力強等特點[15]。獅頭鵝與吉林白鵝產(chǎn)蛋量相差懸殊，適合作為探索鵝繁殖性狀相關(guān)研究的素材。

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

1.3 全基因組重測序

1.3.1 文庫構(gòu)建與測序 樣品基因組DNA檢測合格后，用機械打斷的方法(超聲波)將DNA片段化，然后進行片段純化，末端修復，3′端加A以及連接測序接頭。通過瓊脂糖凝膠電泳選擇片段的大小，PCR擴增目的片段形成測序文庫，對質(zhì)檢合格的文庫進行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進行過濾得到Clean Reads后，方可用于后續(xù)的信息分析。數(shù)據(jù)過濾的主要步驟：1) 去除帶接頭(adapter)的reads；2) 過濾N含量超過10%的reads；3) 去除質(zhì)量值低于10的堿基超過50%的reads。

1.4 Small Indel的檢測與注釋

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

1.5 DNA水平的變異基因挖掘

比較篩選獅頭鵝與吉林白鵝兩個群體特有的差異Indel位點，通過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ì)量分布的分析來進行質(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平均值分別達到96.74%、92.19%。GC含量在43.67%～45.64%之間(表2)，結(jié)果表明，測序質(zhì)量良好，可進行后續(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é)果表明，測序物種同參考基因組相似度高，證明測序準確，可進行后續(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序列的分布位置進行注釋及分析(表5)。

通過SnpEff軟件注釋分析Indel位點，分別從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位點(平均有748 821.2個)，其數(shù)量和分布與參考基因組相似。在這10只試驗鵝中，Indel位點主要分布在內(nèi)含子，吉林白鵝與獅頭鵝Indel位點在內(nèi)含子中分別占有54.18%和52.12%；其次是基因間區(qū)，兩種鵝分別占有36.05%與33.91%(圖3)。大部分的Indel位點分布在內(nèi)含子與基因間區(qū)中，CDS區(qū)中分布的Indel位點最少。

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

表5 Indel 位點的檢測及注釋

在這10只試驗鵝全基因組中，Indel位點的長度分布趨勢也基本一致。對群體Indel突變位點的堿基長度進行統(tǒng)計并分析，得到圖4，發(fā)現(xiàn)全基因組中長度為1 bp的Indel位點數(shù)量超過310 000個(約占全基因組所有Indel的43%)；長度為2～4 bp的Indel位點數(shù)量為2 529 498個；長度為5～9 bp的Indel位點數(shù)量為860 657個；長度大于等于10 bp 的Indel位點數(shù)量共808 872個。并且在編碼區(qū)中，長度為1 bp的Indel位點數(shù)量共22 754個(約占編碼區(qū)所有Indel的33%)；長度為2～4 bp的Indel位點數(shù)量為22 426個；長度為5～9 bp的Indel位點數(shù)量為10 080個；長度大于等于10 bp的Indel位點數(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位點在CDS區(qū)的分布趨勢和數(shù)量也與參考基因組相似，分別有7 090、7 026、 7 022、7 060、6 696、6 649、6 651、6 862、6 868和6 827個Indel位點(表6)。其中共有11 737個Indel位點導致了插入或刪除。

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

2.4 變異基因注釋分析

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

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

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

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

3 討 論

獅頭鵝是我國乃至世界著名的大型水禽品種，也是我國南方養(yǎng)殖范圍較廣的經(jīng)濟家禽，其產(chǎn)蛋期一般為9月至次年4月[24-25]。吉林白鵝是我國北方養(yǎng)殖較廣的經(jīng)濟家禽，其產(chǎn)蛋旺季在2～6月，研究表明，兩種鵝日照長短不一樣[26]。甄霆等[27]研究發(fā)現(xiàn)，在浙東白鵝的繁殖季節(jié)通過調(diào)節(jié)光照可影響FSH和INH的表達，從而延長其產(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)蛋性狀提供了新的思路。由于獅頭鵝較強的就巢性和產(chǎn)蛋季節(jié)性導致其繁殖力較低，從而造成獅頭鵝產(chǎn)蛋性能低下。提高獅頭鵝的產(chǎn)蛋性能一直是研究的主要目標之一。

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

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

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

4 結(jié) 論

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