儋州雞體尺性狀全基因組關(guān)聯(lián)分析

姜宏正，楊德智，馬中華，荀文娟，施力光，侯冠彧

(1.海南大學(xué)動(dòng)物科技學(xué)院，海口 570228；2.中國(guó)熱帶農(nóng)業(yè)科學(xué)院熱帶作物品種資源研究所，?？?571101)

儋州雞又名儋州小種雞、北岸小種雞、石雞等，原產(chǎn)于海南省儋州市，由當(dāng)?shù)厝嗣窠?jīng)過(guò)長(zhǎng)期的馴養(yǎng)而來(lái)，具有耐粗飼、抗病力強(qiáng)、肉質(zhì)鮮嫩、營(yíng)養(yǎng)價(jià)值高等種質(zhì)特性。然而作為新發(fā)現(xiàn)的地方遺傳資源，缺乏對(duì)儋州雞系統(tǒng)地選育。全基因組關(guān)聯(lián)分析(genome wide association study，GWAS)作為后基因組時(shí)代有效的遺傳變異(標(biāo)記)多態(tài)性檢測(cè)方法，在畜禽育種工作中應(yīng)取得顯著的效果。Huang等[1]利用600K芯片對(duì)中國(guó)地方雞的胴體性狀進(jìn)行了GWAS分析，發(fā)現(xiàn)了91個(gè)與雞的胴體性狀關(guān)聯(lián)潛在顯著的SNPs位點(diǎn)，篩選出細(xì)胞角化包膜前體(sciellin，SCEL)和MYC結(jié)合蛋白2(MYC binding protein 2，E3 ubiquitin protein ligase，MYCBP2)等6個(gè)可能的候選基因，豐富了雞胴體性狀遺傳機(jī)理的認(rèn)知。Bai等[2]利用雞600K高密度SNP芯片對(duì)喙畸性狀進(jìn)行GWAS分析，結(jié)果發(fā)現(xiàn)8個(gè)與喙畸形性狀關(guān)聯(lián)的位點(diǎn)，類(lèi)似轉(zhuǎn)錄延伸因子B多肽3(transcription elongation factor B polypeptide 3-like，TCEB3)、Tudor結(jié)構(gòu)域蛋白3(tudor domain containing 3，TDRD3)、RET原癌基因(ret proto-oncogene，RET)和微管解聚蛋白1(stathmin 1，STMN1) 4個(gè)基因可能是雞喙畸形性狀研究的重要候選基因，得到鈣離子信號(hào)調(diào)控通路在內(nèi)的6條與喙畸形性狀顯著關(guān)聯(lián)的通路。Liu等[3]利用SNP芯片對(duì)產(chǎn)蛋后期的蛋品質(zhì)進(jìn)行了GWAS分析，發(fā)現(xiàn)13號(hào)染色體約有0.36 Mb(8.95～9.31 Mb)區(qū)域與蛋白高度和哈氏單位顯著相關(guān)，Ras同源基因家族成員A(ras homolog family member A，RHOA)、基質(zhì)細(xì)胞衍生因子4 (stromal cell derived factor 4，SDF4)和受體超家族成員4 (TNF receptor superfamily member 4，TNF4)這3個(gè)基因可能是與蛋殼顏色相關(guān)的候選基因。體尺性狀作為評(píng)價(jià)家禽體型外貌的主要指標(biāo)，是衡量其生長(zhǎng)發(fā)育和品種鑒定的重要特征。由于家禽的體尺性狀與其產(chǎn)肉或產(chǎn)蛋性能密切相關(guān)，通過(guò)培育代次的體尺測(cè)量和選種來(lái)進(jìn)行經(jīng)濟(jì)性狀表型選擇，可有效縮短育種進(jìn)程。張建豐等[4]研究了中國(guó)地方雞的體尺性狀對(duì)繁殖性狀的影響，發(fā)現(xiàn)脛長(zhǎng)越長(zhǎng)平均蛋重越大,體斜長(zhǎng)越長(zhǎng)年產(chǎn)蛋數(shù)越多,龍骨長(zhǎng)越長(zhǎng)雞的平均開(kāi)產(chǎn)日齡越晚。李萬(wàn)貴等[5]對(duì)興義鴨體尺性狀和屠宰性能進(jìn)行相關(guān)性分析，結(jié)果發(fā)現(xiàn)，胸寬與全凈膛重、半凈膛重、肌胃重呈顯著正相關(guān)。徐明明等[6]探討了雪峰烏骨雞體尺性狀與產(chǎn)肉性能的相關(guān)關(guān)系，發(fā)現(xiàn)可通過(guò)胸深這一體尺指標(biāo)預(yù)測(cè)雪峰烏骨雞公雞的屠宰性能。目前，利用GWAS對(duì)地方雞體尺性狀的研究尚不全面。本試驗(yàn)利用GWAS技術(shù)對(duì)70日齡儋州雞體尺性狀進(jìn)行分析，旨在通過(guò)GWAS識(shí)別影響儋州雞體尺性狀的分子標(biāo)記，探討遺傳標(biāo)記的分布規(guī)律，尋找影響儋州雞體尺性狀的候選基因，為儋州雞早期選種育種及開(kāi)發(fā)利用奠定理論基礎(chǔ)。

1 材料與方法

1.1 試驗(yàn)動(dòng)物與飼養(yǎng)管理

隨機(jī)選取200只同一批次孵化的1日齡儋州雞(中國(guó)熱帶農(nóng)業(yè)科學(xué)院熱帶作物品種資源研究所提供)，公母各半，按照《NY/T 33-2004雞飼養(yǎng)標(biāo)準(zhǔn)》[7]進(jìn)行飼養(yǎng)。

1.2 主要試劑和儀器

動(dòng)物血液/細(xì)胞/組織基因組DNA提取試劑盒(DP304)購(gòu)自天根生化科技(北京)有限公司。紫外分光光度計(jì)(NanoDrop 2000)購(gòu)自Thermo Fisher Scientific公司。

1.3 樣品采集和分析

1.3.1 血樣采集及DNA提取 70日齡時(shí)，翅靜脈采血，EDTA抗凝處理，使用動(dòng)物血液/細(xì)胞/組織基因組DNA提取試劑盒提取基因組DNA，利用紫外分光光度計(jì)對(duì)DNA樣品的濃度和純度進(jìn)行測(cè)定，測(cè)定標(biāo)準(zhǔn)為A260 nm/A280 nm值為1.8～2.0。質(zhì)量合格的樣品送杭州聯(lián)川生物技術(shù)股份有限公司進(jìn)行10×深度的全基因組重測(cè)序。

1.3.2 體尺性狀測(cè)定 儋州雞飼養(yǎng)至70日齡時(shí)，按照《NY/T 823-2004家禽生產(chǎn)性能名詞術(shù)語(yǔ)和度量統(tǒng)計(jì)方法》[8]進(jìn)行體尺性狀測(cè)定，包括體斜長(zhǎng)、龍骨長(zhǎng)、脛長(zhǎng)、脛圍、胸寬、胸深等。體斜長(zhǎng)測(cè)定肩關(guān)節(jié)至坐骨結(jié)節(jié)間的距離；龍骨長(zhǎng)測(cè)定龍骨突前端到龍骨末端的距離；脛長(zhǎng)測(cè)定脛部上關(guān)節(jié)到第3、4趾間的直線距離；脛圍測(cè)定脛骨中部的周長(zhǎng)；胸深測(cè)量第1胸椎到龍骨前緣的距離。測(cè)得數(shù)據(jù)經(jīng)過(guò)Excel整理，進(jìn)行GWAS分析。

1.4 測(cè)序數(shù)據(jù)的質(zhì)控

使用FastP(0.19.3)對(duì)原始數(shù)據(jù)(Raw data)進(jìn)行預(yù)處理。主要步驟為：去除接頭(Adapter)；去除含有N(N表示無(wú)法確定堿基信息)的比例>5%的reads；去除低質(zhì)量reads；利用BWA將Clean data比對(duì)到參考基因組上，并使用GATK軟件進(jìn)行單核苷酸多態(tài)性(single nucleotide polymorphysm，SNP)的檢測(cè)，并對(duì)檢測(cè)到的變異位點(diǎn)進(jìn)行質(zhì)量過(guò)濾。SNP的質(zhì)量過(guò)濾標(biāo)準(zhǔn)：QD>2.0、FS>60.0、MQ>40.0、MQRankSum>-12.5、ReadPosRankSum>-8.0、SOR>3.0。

1.5 GWAS分析

1.5.1 關(guān)聯(lián)分析模型 使用EMMAX[9](versionbeta-07)中的混合線性模型，對(duì)儋州雞體尺性狀與全基因組上SNP進(jìn)行GWAS分析，混合線性模型為：

y=μ+Xb+u+e

其中，y表示儋州雞體尺性狀表型值；μ表示性狀的平均值；X表示固定效應(yīng)矩陣，其中固定效應(yīng)包括養(yǎng)殖環(huán)境、育種方式和性別；b為固定效應(yīng)向量；u表示剩余多基因效應(yīng)；e表示表型值的隨機(jī)殘差，u和e均服從正態(tài)分布。

1.5.2 遺傳進(jìn)化與群體結(jié)構(gòu)分析 基于SNP，通過(guò)Admixture(1.3.0)軟件，分析所有樣本的群體結(jié)構(gòu)，分別假設(shè)分群數(shù)(K值)為1～10，進(jìn)行聚類(lèi)，預(yù)測(cè)樣本之間遺傳距離，對(duì)關(guān)聯(lián)分析進(jìn)行輔助。使用SPAGeDi(1.4)對(duì)所有樣本兩兩之間的親緣關(guān)系進(jìn)行計(jì)算。

1.5.3 多重檢驗(yàn) 利用Bonferroni方法對(duì)SNPs的顯著性進(jìn)行判定。若SNP的值<0.01/N的值(N表示質(zhì)控后的SNP數(shù)目)，則SNP與性狀顯著關(guān)聯(lián)；若SNP的值<0.1/N的值，則表示SNP與性狀潛在顯著關(guān)聯(lián)。

1.6 基因注釋及候選基因篩選

利用Ensembl (http:∥oct2018.archive.ensembl.org/Gallus-gallus/Info/Index)和NCBI(http:∥www.ncbi.nlm.nih.gov/)搜索顯著及潛在關(guān)聯(lián)SNP位點(diǎn)前后各100 kb區(qū)域的上、下游基因，將所得到的基因與GO、KEGG等數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，根據(jù)基因注釋情況分析可能的候選基因。

2 結(jié) 果

2.1 儋州雞體尺性狀測(cè)定結(jié)果

儋州雞體尺性狀統(tǒng)計(jì)結(jié)果見(jiàn)表1。70日齡儋州雞平均脛長(zhǎng)66.55 mm，脛圍20.37 mm，體斜長(zhǎng)89.06 mm，胸寬27.81 mm，髖骨寬43.23 mm，胸深44.57 mm，龍骨長(zhǎng)63.62 mm。

表1 儋州雞體尺性狀描述性統(tǒng)計(jì)(n=200)

2.2 儋州雞遺傳進(jìn)化與群體結(jié)構(gòu)分析

利用Admixture軟件評(píng)估群體分層，根據(jù)ΔK的分析結(jié)果確定最優(yōu)分群數(shù)為1(圖1)，表明試驗(yàn)群體未發(fā)生亞群分化； 親緣關(guān)系熱圖如圖2所示， 發(fā)現(xiàn)群體未出現(xiàn)家系分化，適宜進(jìn)行后續(xù)的GWAS分析。

A中每種顏色代表一個(gè)群，每行代表一個(gè)分群值的情況

圖2 親緣關(guān)系熱圖

2.3 儋州雞體尺性狀GWAS結(jié)果

采用EMMAX模型對(duì)儋州雞7個(gè)體尺指標(biāo)性狀(脛長(zhǎng)、體斜長(zhǎng)、脛圍、胸寬、胸深、髖骨寬、龍骨長(zhǎng))GWAS分析結(jié)果如圖3～9所示。脛長(zhǎng)和脛圍性狀達(dá)到基因組顯著水平的SNPs位點(diǎn)分別有12和8個(gè)。其中與脛長(zhǎng)性狀相關(guān)的SNPs分別定位于1、2、4和8號(hào)染色體上；與脛圍性狀相關(guān)的SNPs定位于2、4、8和13號(hào)染色體上。未發(fā)現(xiàn)與其他5個(gè)性狀顯著關(guān)聯(lián)的SNP位點(diǎn)。

圖3 脛長(zhǎng)的全基因組關(guān)聯(lián)分析曼哈頓圖

圖4 脛圍的全基因組關(guān)聯(lián)分析曼哈頓圖

圖5 體斜長(zhǎng)的全基因組關(guān)聯(lián)分析曼哈頓圖

圖6 胸寬的全基因組關(guān)聯(lián)分析曼哈頓圖

圖7 髖骨寬的全基因組關(guān)聯(lián)分析曼哈頓圖

圖8 胸深的全基因組關(guān)聯(lián)分析曼哈頓圖

圖9 龍骨長(zhǎng)的全基因組關(guān)聯(lián)分析曼哈頓圖

2.4 候選基因功能注釋

各顯著SNP上、下游預(yù)測(cè)的基因如表2、3所示，發(fā)現(xiàn)8個(gè)基因可作為脛長(zhǎng)和脛圍性狀的候選基因，分別是KCNA1、TPK1、EZH2、FSTL5、AMY2A、TGFBI、LECT2和IL-9，其中FSTL5、AMY2A和TPK1基因可作為同時(shí)影響脛長(zhǎng)和脛圍性狀的候選基因。使用GO和KEGG數(shù)據(jù)庫(kù)對(duì)其進(jìn)行分析，8個(gè)基因參與鉀離子跨膜轉(zhuǎn)運(yùn)、硫胺素新陳代謝、細(xì)胞增殖、鈣離子結(jié)合、骨骼肌衛(wèi)星細(xì)胞維持與骨骼肌再生、細(xì)胞受體相互作用、生長(zhǎng)因子活性等生物學(xué)進(jìn)程(表4)。

表2 脛長(zhǎng)性狀顯著相關(guān)的SNPs位點(diǎn)

表3 脛圍性狀顯著相關(guān)的SNPs位點(diǎn)

表4 影響儋州雞脛長(zhǎng)和脛圍的關(guān)鍵候選基因

3 討 論

體尺性狀是評(píng)價(jià)畜禽生長(zhǎng)性能的重要指標(biāo)。Emrani等[10]研究發(fā)現(xiàn)，影響雞不同周齡脛骨發(fā)育的顯著SNP位點(diǎn)分布于1、4、7、9、21和Z號(hào)染色體上。孫艷發(fā)等[11]使用60K SNP芯片對(duì)F2資源群體的脛長(zhǎng)和脛圍研究發(fā)現(xiàn)，與脛長(zhǎng)相關(guān)的SNPs位點(diǎn)主要分布于1、9、13和27號(hào)染色體上，4號(hào)染色體上71.01—85.94 Mb區(qū)域是與雞脛圍相關(guān)的重要候選區(qū)域。本試驗(yàn)定位到與儋州雞脛長(zhǎng)脛圍性狀相關(guān)的SNPs位于1、2、4和8號(hào)染色體上，大部分位點(diǎn)是未報(bào)道的新位點(diǎn)，這可能與儋州雞地方品種的遺傳多樣性及所采用的EMMAX分析方法差異所致。本試驗(yàn)未預(yù)測(cè)到與胸深、胸寬等性狀顯著相關(guān)的SNP，可能是由于樣本量偏小導(dǎo)致檢驗(yàn)功效較低造成的。但是基于GWAS分析自身對(duì)微效多基因控制的數(shù)量性狀檢測(cè)能力不足的局限性[12]，許多貢獻(xiàn)度偏小的性狀基因在GWAS分析中往往因計(jì)算模型的選擇差異而達(dá)不到顯著水平，進(jìn)而不能準(zhǔn)確挖掘顯著性SNP，或者得到假陽(yáng)性的顯著性結(jié)果。復(fù)雜性狀受到基因-基因、基因-環(huán)境相互作用的影響，GWAS分析很難將上述重要的影響因素考慮在內(nèi)。

Zhang等[13]在基于單倍型的GWAS分析發(fā)現(xiàn)，8號(hào)染色體上10.13-12.97 Mb區(qū)域是與雞脛圍相關(guān)的重要候選區(qū)域。本試驗(yàn)預(yù)測(cè)到在儋州雞8號(hào)染色體0.41 Mb附近的AMY2A基因與脛圍性狀顯著相關(guān)。AMY2A是α-淀粉酶家族的一員，編碼胰腺產(chǎn)生的淀粉酶同工酶。Endo等[14]研究也發(fā)現(xiàn)，抗AMY2A自身抗體可以作為自身免疫性胰腺炎和暴發(fā)性Ⅰ型糖尿病標(biāo)志物。本試驗(yàn)中還檢測(cè)到在儋州雞13號(hào)染色體中有4個(gè)與脛圍性狀相關(guān)的SNPs位點(diǎn)達(dá)到基因組顯著水平，定位到IL9、TGFBI和LECT2共3個(gè)候選基因。TGF-β超家族信號(hào)通路參與了廣泛的生物學(xué)過(guò)程，對(duì)調(diào)控早期胚胎發(fā)育、細(xì)胞生長(zhǎng)、干細(xì)胞自我更新、腫瘤發(fā)生發(fā)展等具有十分重要的調(diào)控作用。Lorda-Diez等[15]研究發(fā)現(xiàn)，TGFBI基因可促進(jìn)TGF-β信號(hào)轉(zhuǎn)導(dǎo)的促纖維化作用，誘導(dǎo)合成的細(xì)胞外基質(zhì)蛋白還能中和缺氧誘導(dǎo)因子對(duì)軟骨細(xì)胞的影響，參與細(xì)胞-細(xì)胞及細(xì)胞-基質(zhì)的黏附。LECT2基因能夠通過(guò)介導(dǎo)Wnt/β-Catenin信號(hào)通路促進(jìn)間質(zhì)干細(xì)胞的成骨分化[16]。IL-9是白細(xì)胞介素家族的一員，作為機(jī)體重要的免疫因子，可由Th9、Th17等多種細(xì)胞分泌,通過(guò)與受體結(jié)合發(fā)揮生物學(xué)作用。Biscarini等[17]研究發(fā)現(xiàn)，IL-9基因與雞的羽毛狀況相關(guān)，揭示了免疫系統(tǒng)和啄羽行為之間的關(guān)系。 以上結(jié)果說(shuō)明,AMY2A、IL-9、TGFBI和LECT2基因的相關(guān)研究多集中在機(jī)體免疫和細(xì)胞分化方面,在動(dòng)物尤其是家禽方面的研究較少,因此本試驗(yàn)所預(yù)測(cè)到的基因?qū)w尺性狀的影響需要進(jìn)一步進(jìn)行功能驗(yàn)證篩選。

本試驗(yàn)在儋州雞1、2、4號(hào)染色體找到顯著的SNP位點(diǎn),定位到EZH2、TPK1、FSTL5和KCNA1 4個(gè)候選基因。 Akizu等[18]研究發(fā)現(xiàn)，EZH2基因?qū)τ诰S持脊椎動(dòng)物早期發(fā)育的神經(jīng)管內(nèi)環(huán)境平衡是必不可少的，也與體溫調(diào)節(jié)系統(tǒng)的發(fā)育和功能有關(guān)[19]。FSTL5是卵泡抑素樣家族的一個(gè)成員，編碼一種分泌性糖蛋白，在細(xì)胞生長(zhǎng)和分化方面具有重要作用，其表達(dá)量的升高會(huì)促進(jìn)細(xì)胞的凋亡[20]。Zhang等[21]研究證實(shí)，F(xiàn)STL5基因在肝細(xì)胞癌中起抑制作用，表明FSTL5基因表達(dá)下調(diào)可以激活Wnt/β-Catenin信號(hào)通路，進(jìn)而提高肝細(xì)胞癌細(xì)胞的生長(zhǎng)和存活能力。KCNA1是鉀離子通道蛋白家族成員之一，調(diào)節(jié)腦和腎中樞神經(jīng)系統(tǒng)的鉀離子通道，介導(dǎo)跨膜鉀離子在興奮性膜的傳導(dǎo)，以及調(diào)節(jié)膜電位和神經(jīng)信號(hào)，以防止神經(jīng)元過(guò)度興奮[22]。研究還發(fā)現(xiàn)，KCNA1基因與腫瘤細(xì)胞的細(xì)胞增殖、凋亡等過(guò)程顯著相關(guān)[23-24]。TPK1是硫胺素代謝過(guò)程中的關(guān)鍵酶之一，催化硫胺素磷酸化為二磷酸硫胺和三磷酸硫胺。Fradin等[25]研究結(jié)果表明，胎兒和母體TPK1基因的基因組變異可能導(dǎo)致正常人出生體重的變化。這4個(gè)候選基因與細(xì)胞增殖、分化和生長(zhǎng)調(diào)控相關(guān),由于體尺性狀是復(fù)雜性狀,受到基因-基因、基因-環(huán)境的相互作用的影響，仍需在其他品種及更大的群體中進(jìn)行后續(xù)的功能驗(yàn)證。本研究得到的候選基因和位點(diǎn)多次為首次發(fā)現(xiàn)，這可能是品種及分析方法不同所致，后續(xù)仍有待進(jìn)行更深入的研究。

4 結(jié) 論

本研究利用10×重測(cè)序?qū)僦蓦u體尺性狀進(jìn)行GWAS分析，發(fā)現(xiàn)了20個(gè)SNPs位點(diǎn)與體尺性狀顯著相關(guān)，并篩選到KCNA1、EZH2、TGFBI、LECT2IL-9、TPK1、FSTL5和AMY2A共8個(gè)目標(biāo)性狀候選基因,為儋州雞的育種工作提供候選的分子標(biāo)記，為地方品種雞標(biāo)記輔助選擇提供了新的思路。