儋州雞體尺性狀全基因組關聯分析

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

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

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

1 材料與方法

1.1 試驗動物與飼養管理

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

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。質量合格的樣品送杭州聯川生物技術股份有限公司進行10×深度的全基因組重測序。

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

1.4 測序數據的質控

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

1.5 GWAS分析

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

y=μ+Xb+u+e

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

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

1.5.3 多重檢驗 利用Bonferroni方法對SNPs的顯著性進行判定。若SNP的值<0.01/N的值(N表示質控后的SNP數目)，則SNP與性狀顯著關聯；若SNP的值<0.1/N的值，則表示SNP與性狀潛在顯著關聯。

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

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

2 結 果

2.1 儋州雞體尺性狀測定結果

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

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

2.2 儋州雞遺傳進化與群體結構分析

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

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

圖2 親緣關系熱圖

2.3 儋州雞體尺性狀GWAS結果

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

圖3 脛長的全基因組關聯分析曼哈頓圖

圖4 脛圍的全基因組關聯分析曼哈頓圖

圖5 體斜長的全基因組關聯分析曼哈頓圖

圖6 胸寬的全基因組關聯分析曼哈頓圖

圖7 髖骨寬的全基因組關聯分析曼哈頓圖

圖8 胸深的全基因組關聯分析曼哈頓圖

圖9 龍骨長的全基因組關聯分析曼哈頓圖

2.4 候選基因功能注釋

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

表2 脛長性狀顯著相關的SNPs位點

表3 脛圍性狀顯著相關的SNPs位點

表4 影響儋州雞脛長和脛圍的關鍵候選基因

3 討 論

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

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

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

4 結 論

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