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

(新疆農(nóng)業(yè)大學動物科學學院，烏魯木齊 830052)

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

1 材料與方法

1.1 樣品采集

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

1.2 DNA提取及重測序

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

1.3 SNP變異檢測

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

1.4 群體系統(tǒng)進化樹構(gòu)建

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

1.5 群體遺傳結(jié)構(gòu)分析

采用frappe進行群體結(jié)構(gòu)分析。首先創(chuàng)建PLINK的輸入文件-Ped文件，然后利用frappe軟件構(gòu)建群體遺傳結(jié)構(gòu)和群體世系信息。

1.6 選擇信號分析

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

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

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

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

2 結(jié) 果

2.1 全基因組重測序結(jié)果

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

表1 伊犁鵝與霍爾多巴吉鵝測序數(shù)據(jù)平均統(tǒng)計表

2.2 測序數(shù)據(jù)與參考基因組的比對分析

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

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

2.3 SNP變異檢測結(jié)果

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

表3 SNPs分類結(jié)果

2.4 群體系統(tǒng)進化樹分析

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

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

2.5 群體遺傳結(jié)構(gòu)分析

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

K=2代表有2個不同的亞群，不同顏色代表不同的亞群，如某個樣品由多種顏色構(gòu)成，則對應(yīng)著縱坐標可看出每個亞群所占的比例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 伊犁鵝與霍爾多巴吉鵝遺傳結(jié)構(gòu)圖Fig.2 Genetic structure diagram of Yili and Hortobgy geese

2.6 群體分化結(jié)果

由圖3可知，以top 5%為閾值作為篩選標準，即top 5%的Fst值為0.228的窗口被視為受選擇的窗口，共篩選到1 231個受選擇的候選區(qū)域，其中有5個選擇信號窗口的Fst值>0.5，位于5和22號染色體，F(xiàn)st值最高的選擇信號窗口位于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聯(lián)合分析結(jié)果

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

2.8 基因富集分析

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

橫坐標為Pi值，縱坐標為Fst值，分別對應(yīng)上面和右側(cè)的頻率分布圖。中部的點代表不同窗口內(nèi)的相應(yīng)的Fst和Pi值，其中藍色和綠色區(qū)域為Fst和Pi選擇出來的top 5%區(qū)域，紅色區(qū)域為兩者共同所選擇的top 5%區(qū)域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通路中的基因富集情況

續(xù)表

3 討 論

3.1 測序數(shù)據(jù)分析

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

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

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

3.2 候選基因推測分析

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

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

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

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

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

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

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

4 結(jié) 論

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