曹海月，譚雨革，董信陽，毛海光，馬有智，盧 磊，姜俊保，尹兆正*

（1.浙江大學(xué)動物科學(xué)學(xué)院，浙江杭州 310058；2.寧波市振寧牧業(yè)有限公司,浙江寧波 315600）

全基因組關(guān)聯(lián)分析（Genome-Wide Association Study,GWAS）是指在全基因組范圍內(nèi)找出存在的序列變異，即單核苷酸多態(tài)性（Single Nucleotide Polymorphisms,SNPs），并從中篩選出與目的性狀相關(guān)的SNPs。GWAS 研究最先被應(yīng)用于醫(yī)學(xué)，Klein 等[1]在《Science》雜志上首次發(fā)表了利用GWAS 研究人類視網(wǎng)膜黃斑病的文章，自此之后，GWAS 被廣泛應(yīng)用于人類疾病致病基因的篩選，并成功鑒定了糖尿病、抑郁癥等疾病的主要致病基因[2-3]。近年來，GWAS 開始逐漸應(yīng)用于畜禽數(shù)量性狀的研究，鑒定出大量與豬繁殖性狀、雞生長及蛋品質(zhì)等性狀相關(guān)的候選基因[4-6]。隨著高通量測序和高分辨代謝檢測技術(shù)不斷發(fā)展，以及多種生物信息學(xué)技術(shù)和統(tǒng)計學(xué)方法發(fā)展，使復(fù)雜性狀基因變異的定位更加準(zhǔn)確，為畜禽重要經(jīng)濟(jì)性狀主效基因和疾病性狀致病基因挖掘提供了更加有效的方法[7]。

寧海黃雞和廣西黃雞是我國南方地區(qū)的2 個優(yōu)質(zhì)地方雞種，相比而言，寧海黃雞產(chǎn)蛋率高，廣西黃雞生長速度快。由于表型差異通常由基因差異引起，因此探究影響2 個品種性狀特征的基因?qū)⒂欣诘胤狡贩N雞的遺傳改良和配套利用。前人研究表明，雞的大多數(shù)性狀受多個微效基因控制[8]。近年來，GWAS 被大量應(yīng)用于中低密度的單核苷酸多態(tài)性與雞生長、繁殖性狀相關(guān)性的研究[9-10]。高密度（600K）的基因芯片可以使目標(biāo)基因的篩選更加準(zhǔn)確有效[11]，但其目前主要應(yīng)用于雞體重、疾病和脂肪沉積等方面的研究[6,12-13]。為此，本研究采用雞600K SNP 芯片技術(shù)對寧海黃雞和廣西黃雞進(jìn)行基因分型和GWAS 分析，篩選鑒定與繁殖性狀相關(guān)的SNP 位點和候選基因。

1 材料與方法

1.1 實驗動物 寧海黃雞和廣西黃雞種母雞各59 只由浙江省寧波市振寧牧業(yè)有限公司提供，個體籠養(yǎng)，營養(yǎng)和環(huán)境條件一致，其中產(chǎn)蛋期光照時間每天14 h。

1.2 性狀測定 記錄每只母雞在300 日齡內(nèi)的產(chǎn)蛋數(shù)（Egg Number at the Age of 300 d，EN300）。28～29周齡為選留下一代種蛋期，依據(jù)種蛋大小、形狀和蛋殼厚度等選出合格種蛋，然后轉(zhuǎn)移到機(jī)器孵化器進(jìn)行系譜孵化。記錄每只母雞28～29 周齡合格種蛋的總數(shù),記為入孵蛋數(shù)（Setting Eggs，SE）。在孵化第5 天，照蛋除去未受精的蛋，記錄每只母雞剩余受精蛋的數(shù)量，并將其與對應(yīng)母雞產(chǎn)蛋數(shù)的比率計為受精率（Fertility Rate，F(xiàn)R）。孵化至第18 天時，把孵化器里的受精蛋移到出雛器里繼續(xù)孵化，統(tǒng)計移盤數(shù)并記為移盤數(shù)（Transferring Eggs，TE）。在出殼當(dāng)天，將每只母雞的雛雞數(shù)計為出雛數(shù)（Hatching Number，HN），將未孵化出雛雞的受精蛋的數(shù)量計為死胚蛋數(shù)（Addle Egg Number，AE）。將出雛數(shù)與入孵蛋數(shù)及受精蛋數(shù)的比例分別計為入孵蛋孵化率（Hatchability Of Setting Eggs，HSE）和受精蛋孵化率（Hatchability Of Fertile Eggs，HFE）。

1.3 基因分型和數(shù)據(jù)質(zhì)控 雞頸靜脈采血并儲存在-20℃。依據(jù)制造商說明，使用Genome DNA Extraction Kit（天根生化科技有限公司，中國北京）提取基因組DNA，并稀釋至50～100 ng/μL。SNP 芯片試驗的基因組DNA 應(yīng)達(dá)到以下條件：采用NanoDrop2000 紫外分光光度計測量OD 值，A260/A280在1.8～2.1；用1%瓊脂糖凝膠進(jìn)行電泳，條帶清晰，且長度大于10 kb。將合格的DNA 送紐勤生物科技（上海）有限公司進(jìn)行雞600K SNP 芯片檢測和基因分型。采用PLINK（v1.09）進(jìn)行基因型數(shù)據(jù)的質(zhì)控[14]。

1.4 統(tǒng)計分析 利用PLINK（v1.09）軟件對所有樣本的600K 芯片掃描數(shù)據(jù)進(jìn)行質(zhì)量控制，剔除SNPs 檢出率小于90% 且哈迪-溫伯格平衡檢驗P＜1×10-6（卡方檢驗）和最小等位基因頻率（MAF）≤0.05 的SNPs，剔除分型成功率小于90% 的樣品。使用SPSS 20.0 對表型數(shù)據(jù)進(jìn)行統(tǒng)計分析，使用PLINK v1.09 中的主成分分析法（Principal Component Analysis，PCA）對兩個群體進(jìn)行群體分層評估，使用PLINK（v1.09）的線性模型進(jìn)行全基因組關(guān)聯(lián)分析，所用模型：

其中，Y 是性狀表型值向量；β是表型均值、SNPs、群體結(jié)構(gòu)（校正群體結(jié)構(gòu)時，使用前5 個主成分作為協(xié)變量）3 個固定效應(yīng)向量；X 是β 的關(guān)聯(lián)矩陣；e 是殘差效應(yīng)向量，e～ N（0，Iσ2e），其中I 是單位陣，Iσ2e是隨機(jī)殘差方差，最終，每個SNP 位點都能得到一個關(guān)聯(lián)值。為減少多重檢驗帶來的假陽性率，以Bonferroni 校正法[15]對全基因組關(guān)聯(lián)分析結(jié)果的P值進(jìn)行校正。此處獨立檢驗數(shù)計算使用PLINK（v1.09）中的“indeppairwise 25 5 0.2”命令（即以25 個SNPs 為一個窗口，5 個SNP 為步移，r2閾值為0.2），估算出連鎖不平衡（LD）塊和單個獨立SNP 數(shù)目為149 653，因此Bonferroni校正的達(dá)5% 基因組水平顯著的P 值閾值為3.341E-07（0.05/149653），即P值低于此閾值的SNPs 則被認(rèn)為與繁殖性狀顯著關(guān)聯(lián)；達(dá)到基因組水平潛在關(guān)聯(lián)的P值閾值為6.682E-06（1/149653），即P值低于此閾值的SNPs 則被認(rèn)為與繁殖性狀潛在關(guān)聯(lián)。使用R 軟件作Quantile-Quantile（QQ）圖和曼哈頓圖。

2 結(jié)果

2.1 表型值的描述性統(tǒng)計 如表1 所示，寧海黃雞和廣西黃雞的EN300、SE、TE 和HN 存在極顯著差異。

2.2 基因分型和數(shù)據(jù)質(zhì)控 剔除SNP 檢出率小于90%的位點5 604 個、哈迪-溫伯格平衡檢驗P＜1×10-6（卡方檢驗）的位點2 734 個，剔除MAF ≤ 0.05 的SNP 位點78 879 個，剔除分型成功率小于90%的樣品0 個，剩余495 704 個SNPs 118 個樣本可用于后續(xù)與繁殖性狀的全基因組關(guān)聯(lián)分析。獨立SNPs 計算后，共得到獨立SNPs 標(biāo)記和Block 模塊149 653 個。通過統(tǒng)計基因組每100 K 范圍內(nèi)的SNP 數(shù)，得到SNPs 在基因各染色體上的分布情況（圖1），結(jié)果顯示基于高密度芯片技術(shù)檢測到的SNPs 在染色體上分布均勻，SNPs 數(shù)據(jù)可靠。

表1 寧海黃雞和廣西黃雞繁殖性狀的描述性統(tǒng)計

圖1 質(zhì)控后的SNPs 標(biāo)記在各染色體上的分布情況

2.3 群體層化分析 在GWAS 分析中，群體分層會導(dǎo)致分析的結(jié)果出現(xiàn)假陽性，但主成分分析結(jié)果顯示，在2個雞群體中未顯示明顯的亞群結(jié)構(gòu)（圖2），并且繁殖性狀的Quantile-Quantile 圖結(jié)果顯示，觀測值（縱坐標(biāo)）和期望值（橫坐標(biāo)）基本吻合（圖3），說明關(guān)聯(lián)分析不會因為群體分層而產(chǎn)生假陽性，基于線性模型的關(guān)聯(lián)分析結(jié)果可靠。

圖2 群體結(jié)構(gòu)主成分分析圖

2.4 全基因組關(guān)聯(lián)分析 繁殖性狀的曼哈頓圖見圖4，共有4 個SNPs 與繁殖性狀的相關(guān)性達(dá)到Bonferroni 校正的5%基因組顯著性水平，其中rs313221983 和rs314844182影響AE，rs14758703 影響FR，rs312721292 影響HSE，rs314844182 和rs312721292 影響HFE。在每個顯著的SNP 位點上、下游5 kb 內(nèi)尋找可能的候選基因，結(jié)果共找到3 個可能的候選基因，分別為唾液酸轉(zhuǎn)移酶1（ST8 alpha-N-acetyl-neuraminide alpha-2,8-sialyltransferase 1，ST8SIA1）基因，神經(jīng)外胚層皮質(zhì)1（Ectodermal-neural cortex 1，ENC1）基因和LOC101750905基因（表2）。

3 討 論

從本研究結(jié)果可知，寧海黃雞和廣西黃雞的EN300、SE、TE、HN 存在極顯著差異。此外，品種內(nèi)各個性狀都擁有較大的極值差距，顯示品種的遺傳多樣性得到了保持，這種多樣性將使GWAS 分析結(jié)果更加有效。本研究結(jié)果顯示，總共4 個SNPs 與繁殖性狀的關(guān)聯(lián)性達(dá)到了全基因組顯著性水平，其中rs313221983和rs314844182 位于1號染色體，rs14758703 位 于Z染色體，rs312721292 位于4 號染色體，值得注意的是，rs314844182 同時與AE、HFE 2 個性狀相關(guān)聯(lián)，rs312721292 同時與HSE、HFE 2 個性狀相關(guān)聯(lián)。之前的研究通過GWAS 已經(jīng)鑒定出了許多與雞繁殖性狀相關(guān)的QTL、候選區(qū)域和SNPs。對白色來航雞繁殖性狀進(jìn)行GWAS 分析發(fā)現(xiàn)，7 號和13 號染色體上各有1 個SNP 分別與21～56 周的產(chǎn)蛋數(shù)和初產(chǎn)日齡顯著相關(guān)[16]。另外，有研究發(fā)現(xiàn)5 號染色體上29.5～46.9 Mb 區(qū)域與300 日齡產(chǎn)蛋數(shù)和孵化率有關(guān)[10]，并確定了5 號染色體上與21～26 周產(chǎn)蛋數(shù)顯著相關(guān)的6 個基因座[17]。本研究結(jié)果所得到的與AE、FR、HSE 和HFE 顯著相關(guān)的SNPs 位于3 個候選基因的序列附近或內(nèi)部，分別是ST8SIA1、ENC1和LOC101750905。對ST8SIA1基因的報道主要在于其對乳腺癌細(xì)胞的作用，Bobowski 等[18]研究表明，在人乳腺癌細(xì)胞中，ST8SIA1的表達(dá)量受雌二醇抑制，其核心啟動子包含2 個推測的雌激素反應(yīng)元件，結(jié)合本研究結(jié)果表明，ST8SIA1基因是影響雞繁殖性狀的候選基因之一。ENC1主要在神經(jīng)系統(tǒng)中表達(dá)，編碼一種肌動蛋白結(jié)合蛋白，在神經(jīng)細(xì)胞分化和凸起生長中發(fā)揮重要作用[19]。Kim 等[20]研究表明，在排卵前顆粒細(xì)胞中，ENC1 被促黃體生成素/人絨毛膜促性腺激素（LH/ hCG）誘導(dǎo)與F-肌動蛋白物理結(jié)合，可能通過調(diào)節(jié)細(xì)胞骨架組織促進(jìn)顆粒細(xì)胞的分化，這與本研究結(jié)果不謀而合，ENC1基因可以作為影響雞繁殖性狀的候選基因之一。有關(guān)LOC101750905基因的報道目前較少，GWAS 研究表明，此基因下游的rs14490981位點對活體重、腳重、全凈膛重和半凈膛重4 個性狀顯著相關(guān)[21]。本研究表明，LOC101750905基因與雞繁殖性狀相關(guān)。今后尚需要進(jìn)一步研究ST8SIA1、ENC1和LOC101750905基因?qū)﹄uAE、FR、HSE 和HFE 的調(diào)節(jié)機(jī)制。

圖3 繁殖性狀的Quantile-Quantile 圖

表2 繁殖性狀達(dá)到5%基因組顯著水平的全基因組關(guān)聯(lián)分析結(jié)果

圖4 繁殖性狀的曼哈頓圖

4 結(jié) 論

本研究采用600K SNP 芯片檢測及GWAS，發(fā)現(xiàn)了4 個SNPs 對寧海黃雞和廣西黃雞繁殖性狀關(guān)聯(lián)性達(dá)到Bonferroni 校正5%的全基因組顯著性水平，并檢測到了可能與AE、FR、HSE 和HFE 相關(guān)的3 個近端基因（ST8SIA1、ENC1和LOC101750905）,為進(jìn)一步了解地方雞種繁殖性狀的遺傳基礎(chǔ)和基因組選擇提供了理論依據(jù)。