鵝產蛋性狀全基因組選擇信號分析

唐碧徽，潘 璐，李海英，趙曉鈺，吳盈萍，蔣 騰，丁雅文，吳利楠，曹 妍，梅志勇

(新疆農業大學動物科學學院，烏魯木齊 830052)

隨著人們對鵝產品需求的不斷增長，低繁殖率嚴重阻礙了鵝產業的規模化發展[1]，鵝的生產和繁殖能力取決于遺傳因素和許多非遺傳因素(飼養密度、公母比例、性成熟時體重、飼料組成和光照方案等)，其中遺傳因素更為重要，當非遺傳因素得到充分解決，再加上適當的基因型選擇與人工選擇措施，便可以起到改善鵝的繁殖性能的作用[2]。全基因組重測序是指在已知某一物種基因組序列的基礎上對該物種不同個體的整個基因組序列進行測序[3]。通過比對測序獲得的序列和已知該物種的參考序列可以得到大量的基因組遺傳信息，挖掘不同個體或群體基因組間的結構差異，探討功能基因或區域定位[4]。全基因組重測序數據分析內容中的高級數據分析包含選擇信號分析[5]。劉刁等[6]對深縣豬和梅山豬進行全基因組重測序，采用滑動窗口方法進行群體選擇信號分析，結果顯示，在受選擇區域篩選后得到受選擇位點580個，定位于23個基因，通過GO功能和KEGG通路富集分析顯示，篩選到3個與深縣豬繁殖相關的候選基因。呂世杰等[7]對郟縣紅牛和中國荷斯坦奶牛進行選擇信號檢測，篩選種間差異的基因組區域，通過在動物數量性狀座位(QTL)數據庫中的比對，尋找到3個與牛繁殖性狀相關的候選基因。Liu等[8]對獅頭鵝、浙東白鵝、太湖鵝和籽鵝4個品種共35個個體運用全基因組重測序方法進行選擇信號分析，在4個鵝品種基因組中篩選出許多候選基因，發現了可能與獅頭鵝產肉性狀和籽鵝產蛋性狀的相關通路。伊犁鵝是新疆特有的優良地方禽種，耐粗放飼養、肉質優良，但產蛋性能較低，且就巢性強，成年母鵝年均產蛋僅為8～12枚[9]。霍爾多巴吉鵝是由歐洲引進的一種絨蛋兼用型優良品種，其繁殖能力優于伊犁鵝，成年母鵝平均年產蛋為40～50枚[10]。伊犁鵝與霍爾多巴吉鵝在產蛋數量上差異較大，適合作為研究鵝產蛋性狀的素材。本研究利用全基因組重測序技術，通過選擇信號分析，篩選與鵝產蛋相關的候選基因，以期為伊犁鵝和霍爾多巴吉鵝的群體遺傳分析和產蛋性能的研究提供新的思路。

1 材料與方法

1.1 樣品采集

隨機選取健康狀況良好、飼養管理水平一致的2歲伊犁鵝母鵝和霍爾多巴吉鵝母鵝各24只(伊犁鵝(Y1～Y24)；霍爾多巴吉鵝(H1～H24))，翅下靜脈采血5 mL于EDTA抗凝管中，置于-20 ℃保存，用于DNA提取。所用試驗動物均由新疆額敏縣恒鑫實業有限公司國家級伊犁鵝保種場提供。

1.2 DNA提取及重測序

采用磁珠法從全血樣本中提取基因組DNA，通過1.0%瓊脂糖凝膠電泳檢測DNA完整性，并使用分光光度計(NanoDrop 2000)測定DNA濃度。檢測合格的DNA樣品通過Covaris超聲波破碎儀隨機打斷，經末端修復、加A尾、加測序接頭、純化和PCR擴增等步驟完成整個文庫制備。庫檢合格后，將不同文庫按照有效濃度及目標下機數據量的需求進行Illumina Hiseq PE150測序，經過對測序數據的嚴格過濾，得到高質量的clean base。

1.3 SNP變異檢測

使用BWA軟件(參數：mem-t4-k32-M)將讀取到的高質量數據比對到天府鵝的參考基因組(https:∥ftp.cngb.org/pub/gigadb/pub/10.5524/100001_101000/100789/)，采用SAMTOOLS等軟件對樣本進行群體SNP檢測。經過條件為dp 7(保留樣本深度不低于7×的數據)、Miss 0.1(保留樣本缺失率不高于0.1的數據)、maf 0.05(保留樣本最小等位基因頻率不小于0.05的數據)過濾后，獲得高質量的SNP位點用于后續分析。

1.4 群體系統進化樹構建

經過SNP檢測后，得到的個體SNP可用于計算種群之間的距離。運用TreeBest(http:∥treesoft.sourceforge.net/treebest.shtml)軟件計算距離矩陣，通過鄰接法(Neighbor-Joining method)構建系統進化樹。

1.5 群體遺傳結構分析

采用frappe進行群體結構分析。首先創建PLINK的輸入文件-Ped文件，然后利用frappe軟件構建群體遺傳結構和群體世系信息。

1.6 選擇信號分析

1.6.1 群體分化指數 利用工具vcftools計算伊犁鵝和霍爾多巴吉鵝2個群體之間的群體分化指數(Fst)值，使用窗口大小為40 kb、滑動步長為20 kb的滑動窗口，計算每個窗口內所有標記的Fst平均值，將其作為此窗口的Fst值。

1.6.2Fst和核苷酸多樣性的聯合分析 基于Fst值再分別計算伊犁鵝和霍爾多巴吉鵝的核苷酸多樣性(Pi)值，使用窗口大小為40 kb、滑動步長為20 kb的滑動窗口，計算每個滑動窗口SNP的Pi平均值作為此窗口的Pi值，進一步與Fst值聯合進行窗口的篩選。

1.7 候選基因注釋及基因功能富集分析

從參考基因組對應的基因注釋文件中，將篩選出的受選擇位點定位到基因。使用pfam和KOBAS 2.0數據庫進一步對基因進行GO功能和KEGG通路富集分析。

2 結 果

2.1 全基因組重測序結果

伊犁鵝和霍爾多巴吉鵝2個群體的全基因組重測序結果見表1，測序共獲得1 066.665 Gb的原始數據，各樣本的原始數據在12 911 727 038～31 532 656 750 bp之間，經過濾處理后得到的有效數據在12 493 713 313～31 423 739 525 bp之間，平均有效數據率達96.76%，Q20平均達到96.47%，Q30平均達到91.80%，GC含量在44.32%～45.90%之間。表明本次測序數據質量較高，可滿足后續試驗要求。

表1 伊犁鵝與霍爾多巴吉鵝測序數據平均統計表

2.2 測序數據與參考基因組的比對分析

質控過濾后的測序數據通過BWA軟件比對到鵝參考基因組。由表2可知，本試驗測序數據與參考基因組的比對率在97.31%～97.44%之間，測序深度在9.49×～21.06×之間，1×覆蓋率在98.18%以上，4×覆蓋率在83.90%以上。表明本試驗的測序物種與參考基因組相似度高，測序數據準確，可進行后續分析。

表2 伊犁鵝與霍爾多巴吉鵝平均測序深度及覆蓋度統計表

2.3 SNP變異檢測結果

本試驗對2種鵝群體SNPs進行變異檢測，結果見表3。由表3可知，共檢測到15 406 790個SNPs位點，對檢測到的SNPs位點通過ANNOVAR軟件進行注釋發現，SNPs在基因間、內含子和外顯子區域的比例分別為65.35%、29.49%和1.77%。其中外顯子區域中，少量的SNPs獲得終止密碼子的變異(2 003個)或失去終止密碼子的變異(282個)，但大多數SNPs導致同義突變(154 517個)或非同義突變(117 659個)，轉換(Ts)為10 838 735個，顛換(Tv)為4 568 055個，Ts/Tv為2.372。

表3 SNPs分類結果

2.4 群體系統進化樹分析

伊犁鵝群體和霍爾多巴吉鵝群體的系統進化樹結果見圖1。由圖1可知，兩種鵝群體明顯分為兩大支，伊犁鵝群體和霍爾多巴吉鵝群體各自聚為一支，說明兩種鵝群體之間的親緣關系存在一定的距離。

圖1 伊犁鵝與霍爾多巴吉鵝群體系統進化樹Fig.1 Phylogenetic tree of Yili and Hortobgy geese populations

2.5 群體遺傳結構分析

由圖2可知，K=2時，霍爾多巴吉鵝群體全部為一個顏色，說明霍爾多巴吉鵝群體全部在一個亞群群體中；而伊犁鵝群體出現了兩種顏色，說明其可能來自于2個亞群群體。

K=2代表有2個不同的亞群，不同顏色代表不同的亞群，如某個樣品由多種顏色構成，則對應著縱坐標可看出每個亞群所占的比例K=2 represents the two of different subgroups,and different colors represent different subgroups，if a sample is composed of multiple colors,the proportion of each subgroup can be seen according to the ordinate圖2 伊犁鵝與霍爾多巴吉鵝遺傳結構圖Fig.2 Genetic structure diagram of Yili and Hortobgy geese

2.6 群體分化結果

由圖3可知，以top 5%為閾值作為篩選標準，即top 5%的Fst值為0.228的窗口被視為受選擇的窗口，共篩選到1 231個受選擇的候選區域，其中有5個選擇信號窗口的Fst值>0.5，位于5和22號染色體，Fst值最高的選擇信號窗口位于22號染色體。

虛線代表選擇閾值(top 5%)The dotted line represents the selection threshold (top 5%)圖3 伊犁鵝與霍爾多巴吉鵝Fst圖Fig.3 Fst diagram of Yili and Hortobgy geese

2.7 Fst & Pi聯合分析結果

本試驗中Fst&Pi選取閾值為top 5%，關聯Fst&Pi提取相應候選區域，并提取相應區域內的變異位點信息，采用選擇性清除方法檢測到相應的選擇信號區域，伊犁鵝和霍爾多巴吉鵝基于Fst&Pi的選擇消除分析結果見圖4。由圖4可知，伊犁鵝群體與霍爾多巴吉鵝群體之間的Fst值為0.198，Pi值分別為-0.203 93和2.050 78，在伊犁鵝群體中Fst和Pi的選擇性消除分析中有10個受到選擇的區域，其中得到注釋的基因有5個；在霍爾多巴吉鵝群體中Fst和Pi選擇消除分析中有353個受到選擇的區域，其中得到注釋的基因有263個。

2.8 基因富集分析

對伊犁鵝和霍爾多巴吉鵝2個群體中受到選擇的基因分別進行GO功能和KEGG通路富集分析。在GO功能富集分析中，2個鵝群體中受到選擇區域中的基因主要富集在生物過程，如DNA重組(DNA recombination)、內穩定過程(homeostatic process)、Wnt受體信號通路(Wnt receptor signaling pathway)和細胞的自動調節(cellular homeostasis)等過程。在KEGG通路富集分析中，伊犁鵝群體顯著富集到細胞黏附分子(cell adhesion molecules)信號通路，霍爾多巴吉鵝群體主要富集到了Wnt信號通路(Wnt signaling pathway)、TGF-β信號通路(TGF-beta signaling pathway)、MAPK信號通路(MAPK signaling pathway)、胞質DNA感應途徑(cytosolic DNA-sensing pathway)和細胞因子-細胞因子受體相互作用(cytokine-cytokine receptor interaction)等通路。本試驗從這些通路中篩選到了6個與產蛋性狀相關的候選基因：BMP2、BMP6、ENO1、MIS、LIF和EP300(表4)，還篩選到了1個與免疫相關的基因IL-18。

橫坐標為Pi值，縱坐標為Fst值，分別對應上面和右側的頻率分布圖。中部的點代表不同窗口內的相應的Fst和Pi值，其中藍色和綠色區域為Fst和Pi選擇出來的top 5%區域，紅色區域為兩者共同所選擇的top 5%區域The horizontal coordinates is Pi values，and the vertical coordinates are Fst values,which correspond to the frequency distribution chart above and the frequency distribution chart on the right,respectively.The dots in the middle represent the corresponding Fst and Pi values in different windows,where the blue and green areas are the top 5% areas selected by Fst and Pi,and the red area is the top 5% area selected by both together圖4 伊犁鵝(A)和霍爾多巴吉鵝(B)Fst & Pi選擇消除分析圖Fig.4 Fst & Pi selection elimination analysis chart of Yili (A) and Hortobágy (B) geese

表4 伊犁鵝和霍爾多巴吉鵝群體在KEGG通路中的基因富集情況

續表

3 討 論

3.1 測序數據分析

本研究對24只伊犁鵝和24只霍爾多巴吉鵝進行全基因組重測序，測序共獲得1 066.665 Gb的原始數據，其中Q30平均達到91.80%，群體樣本的平均比對率達97.31%，平均測序深度為15.28×，4×覆蓋度在83.90%以上，與參考基因組相似性極高，表明測序數據均一性較好，結果可靠。在此基礎上，通過注釋分析得到SNPs位點共有15 406 790個，大部分(65.35%)SNPs位點分布在基因間區中。SNP的質量一般通過計算Ts/Tv來估計。在人、豬和山羊中，該比值為2.0～2.4[11-13]。本研究中Ts/Tv為2.372，表明潛在隨機測序的誤差較低，與前人研究結果一致，說明本研究中識別的變異具有較高的準確性。

群體進化樹用來表示群體間的進化關系，根據群體的物理或遺傳學特征等方面的共同點或差異可以推斷出它們的親緣關系遠近，即群體個體間由于共同祖先而產生的相互關系，親緣關系較近的物種進化分支往往聚成一簇。本試驗中2個鵝群體的系統發育進化樹表現出較明顯的遺傳分化。群體遺傳結構指遺傳變異在物種或群體中的一種非隨機分布，按照地理分布或其他標準可將一個群體分為若干亞群，處于同一亞群內的不同個體親緣關系較高，而亞群與亞群之間則親緣關系稍遠，不同顏色代表不同的亞群體。群體結構分析發現，在K=2時，伊犁鵝群體表現出內部分化，個別個體與其他個體間存在交叉，可能是基因交流的結果，也可能與所選樣本有關[14]，但整體來看，2個群體之間是可以明顯區分開的，與系統發育樹結果一致，二者相互印證。從SNP層次來看，2種鵝存在較大的差異性，可能是主導其不同性狀的潛在SNP位點造成的。

選擇信號是動物經歷人工選擇或自然選擇的過程中在基因組上留下了一些遺傳特征的現象[15]。吳旭東[16]研究表明，Pi值越高說明相同基因組區域物種的核苷酸多樣性越低。本研究中，霍爾多巴吉鵝作為從歐洲引進的培育品種，經歷了很多人工馴化培育的過程，而伊犁鵝相對來說，并未受到過強烈的人工選育過程，其自身還保留著一定的野性，伊犁鵝的Pi值(-0.20393)低于霍爾多巴吉鵝的Pi值(2.05078)，說明伊犁鵝的核苷酸多樣性較高，推測在選擇消除分析中伊犁鵝的受選擇區域比霍爾多巴吉鵝少與其自身的遺傳信息有關。

3.2 候選基因推測分析

本研究對伊犁鵝和霍爾多巴吉鵝群體進行選擇信號分析，篩選出6個與產蛋性狀相關的候選基因(BMP2、BMP6、MIS、ENO1、LIF和EP300)，以及1個與禽類免疫相關的IL-18基因。

BMP2和BMP6是骨形態蛋白(BMP)家族的成員。BMP2是來源于顆粒細胞的生長因子，在卵泡發育和排卵過程中扮演著重要的角色。研究發現，BMP2及其受體在卵巢內的顆粒細胞、黃體和卵泡液內都被檢測到[17-19]，說明BMP2可能在調節卵泡發育和黃體功能方面起著重要的作用。Haugen等[20]研究表明，在卵泡發育啟動之前，BMP2能夠維持雞顆粒細胞的未分化狀態，同時降低雞顆粒細胞的卵泡刺激素(FSH)敏感性。Selvaraju等[21]研究發現，BMP2可以促進雞卵巢顆粒細胞內雌二醇的合成，并抑制孕酮的生成。魏強[22]對與豬繁殖相關的基因進行基因分型，結果發現，BMP2基因不同基因型在豬的總產仔數、產活仔數和斷奶頭數上差異極顯著，說明BMP2基因與母豬繁殖力顯著相關。BMP6在哺乳動物生殖中也有著重要的作用，BMP6來源于卵母細胞，有刺激顆粒細胞增殖的作用[23]。研究表明，BMP6直接參與了調控繁殖相關激素的合成與分泌，并影響繁殖器官的發育和形成[24]。Sugiura等[25]對敲除BMP6基因的小鼠生殖性能進行測量發現，雌鼠每窩的產仔數只有自然小鼠的78%，推測BMP6可能通過適當增強黃體生成素(LH)信號水平和維持卵母細胞正常功能來增強雌性小鼠的生育能力。姚國瑜[26]研究發現，金寨黑雞BMP6基因存在A76150G位點突變，該位點與金寨黑雞體重和340日齡總產蛋數顯著相關。本研究中，BMP2和BMP6基因富集到TGF-β信號通路，推測BMP2和BMP6基因可能通過調節禽類顆粒細胞的增殖和分化來影響鵝的卵泡發育情況。

ENO1是糖類代謝中的關鍵酶之一，在細胞的新陳代謝中發揮著重要作用，廣泛存在于動物機體的各個組織和器官之中[27]。張偉宏等[28]研究發現，在體外培養的籽鵝顆粒細胞中，ENO1基因過表達會促進孕酮的分泌，推測ENO1基因過表達可能會促進禽類的卵泡排卵。計紅等[29]研究發現，ENO1基因過表達時，可能通過提高籽鵝卵泡顆粒細胞的糖酵解水平從而促進細胞增殖，抑制凋亡，進而改變細胞周期的時相性。本研究發現，ENO1基因在細胞黏附分子信號通路中顯著富集，推測ENO1基因可能與鵝的繁殖性能有關。

抗苗勒氏管激素(AMH)又稱為苗勒管抵制物(MIS)[30]。研究發現，AMH在雌性動物中主要由卵泡顆粒細胞合成和分泌，是目前發現的體內唯一可以調控卵泡募集的細胞因子，對卵泡的生長發育及優勢卵泡的選擇起重要作用[31]。吳玉麟等[32]研究發現，在AMH基因第2外顯子中檢測出的4個SNPs位點與京海黃雞的繁殖性狀均呈顯著相關。陳蓉等[33]采用實時熒光定量PCR技術檢測了連城白鴨不同等級卵泡中AMH等基因的表達情況，結果發現，AMH基因的表達量隨著卵泡的發育而逐漸降低，認為在鴨等級前卵泡中AMH基因的高表達量有助于維持顆粒細胞未分化的狀態，抑制原始卵泡的初始募集，進而維持卵泡的發育潛能。葛麗巖等[34]研究發現，AMH基因中存在的g.1731 G>A和g.2044 T>C位點與黑羽番鴨產蛋性狀存在相關。本研究中，MIS基因富集到TGF-β信號通路，推測其與鵝產蛋性狀相關。

白血病抑制因子(LIF)是一種分泌型糖蛋白，是白細胞介素-6(IL-6)細胞因子家族中的一員，是一個在各種組織中都能發揮生物學作用的多效性細胞因子[35]。研究發現，LIF參與了大鼠原始卵泡的活化[36]；LIF可作為顆粒細胞的增殖因子，促進小鼠培養的竇前卵泡生長[37]；LIF基因與豬的繁殖性能有關[38-40]。本研究中，LIF基因富集到細胞因子-細胞因子受體相互作用的通路，可能促進原始卵泡的激活和早期卵泡的生長發育，推測LIF可能與鵝繁殖性能有關。

E1A結合蛋白P300(EP300)是一種蛋白質編碼基因，該基因編碼腺病毒E1A相關細胞p300轉錄共激活蛋白。EP300作為組蛋白乙酰轉移酶，通過染色質重塑調節轉錄，在細胞增殖和分化過程中起重要作用。陳宇[41]對豬的全基因組SNP篩選時發現，EP300可能是與豬繁殖性狀相關的候選基因。本研究中，EP300雖然富集到了Wnt信號通路，但對于EP300是否和鵝產蛋相關還需要進一步的驗證。

IL-18編碼一種促炎性細胞因子，可增強脾細胞的自然殺傷細胞活性，刺激輔助性T細胞Ⅰ型細胞產生干擾素，是近年發現的一種重要的免疫調節因子[42]。研究發現，IL-18可顯著增強禽流感重組禽痘病毒疫苗的保護效力，在雞中是一種安全的免疫刺激劑[43]。本研究發現，IL-18基因在胞質DNA感應途徑和細胞因子-細胞因子受體相互作用通路中都有富集，推測IL-18基因和禽類免疫有關。

4 結 論

本試驗采用全基因組重測序的方法收集了伊犁鵝和霍爾多巴吉鵝的全基因組SNP變異情況，通過選擇信號分析篩選出了一些受選擇的區域，篩選了6個可能與鵝產蛋性狀相關的候選基因(BMP2、BMP6、MIS、ENO1、LIF、EP300)，以及1個與禽類免疫相關的IL-18基因，為今后在基因組層面選擇鵝的產蛋性狀提供了有價值的信息。