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

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

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

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

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

1 材料與方法

1.1 試驗材料

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

1.2 血液DNA提取

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

1.3 全基因組重測序

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

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

表1 預測堿基判別發生錯誤概率結果

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

1.4 Small Indel的檢測與注釋

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

1.5 DNA水平的變異基因挖掘

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

2 結 果

2.1 數據質量控制

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

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

通過對獅頭鵝和吉林白鵝全基因組重測序，原始數據過濾處理后共得到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)，結果表明，測序質量良好，可進行后續分析。

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

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

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

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

表2 測序數據質量表

表3 參考基因組信息

表4 測序深度及覆蓋度統計

2.3 Indel序列長度及分布分析

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

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

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

表5 Indel 位點的檢測及注釋

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

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

10只試驗鵝的Indel位點在CDS區的分布趨勢和數量也與參考基因組相似，分別有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變異在全基因組和編碼區的統計情況

2.4 變異基因注釋分析

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

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

表7 變異基因功能富集分析的相關注釋

(轉下頁 Carried forward)

3 討 論

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

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

表8 產蛋相關通路及變異基因

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

4 結 論

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