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

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

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

1 材料與方法

1.1 試驗動物與飼養(yǎng)管理

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

1.2 主要試劑和儀器

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

1.3 樣品采集和分析

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

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

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

使用FastP(0.19.3)對原始數(shù)據(jù)(Raw data)進行預處理。主要步驟為：去除接頭(Adapter)；去除含有N(N表示無法確定堿基信息)的比例>5%的reads；去除低質(zhì)量reads；利用BWA將Clean data比對到參考基因組上，并使用GATK軟件進行單核苷酸多態(tài)性(single nucleotide polymorphysm，SNP)的檢測，并對檢測到的變異位點進行質(zhì)量過濾。SNP的質(zhì)量過濾標準：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)中的混合線性模型，對儋州雞體尺性狀與全基因組上SNP進行GWAS分析，混合線性模型為：

y=μ+Xb+u+e

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

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

1.5.3 多重檢驗 利用Bonferroni方法對SNPs的顯著性進行判定。若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位點前后各100 kb區(qū)域的上、下游基因，將所得到的基因與GO、KEGG等數(shù)據(jù)庫進行比對，根據(jù)基因注釋情況分析可能的候選基因。

2 結(jié) 果

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.4 候選基因功能注釋

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

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

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

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

3 討 論

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

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

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

4 結(jié) 論

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