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

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

我國是禽肉生產(chǎn)大國，同時(shí)也是禽肉消費(fèi)大國。據(jù)統(tǒng)計(jì)，我國肉類總產(chǎn)量中禽肉產(chǎn)量所占比重近30%，僅次于豬肉的產(chǎn)量。隨著居民生活水平的提高，消費(fèi)者對(duì)肉品質(zhì)的需求也越來越高，肉色作為評(píng)價(jià)肉品質(zhì)最直觀的感官特征，已成為影響消費(fèi)者購買的主要因素之一。因此，在選育種及養(yǎng)殖生產(chǎn)過程中，篩選出與動(dòng)物肉色形成相關(guān)的候選基因或分子標(biāo)記，對(duì)提高畜禽肉品質(zhì)和經(jīng)濟(jì)價(jià)值有重要意義。全基因組關(guān)聯(lián)研究（Genome-Wide Association Study，GWAS）能夠檢測到與目標(biāo)性狀相關(guān)的單核苷酸多態(tài)性（Single Nucleotide Polymorphism，SNP）標(biāo)記、基因組區(qū)域和候選基因，在畜禽育種工作中已取得顯著效果。針對(duì)肉色性狀的全基因組關(guān)聯(lián)分析研究已經(jīng)在豬、牛身上廣泛開展。Wu 等對(duì)青峪豬群體進(jìn)行GWAS 分析，發(fā)現(xiàn)MYC 靶基因1（MYC Target 1，）和BCL2 相互作用蛋白3（BCL2 Interacting Protein 3，）基因同時(shí)影響著青峪豬的pH 和肉色性狀。Zhang 等利用全基因組關(guān)聯(lián)分析對(duì)雜交豬的肉色性狀進(jìn)行研究，發(fā)現(xiàn)6 個(gè)基因組區(qū)域與雜交豬肉色性狀相關(guān)。Luo 等研究發(fā)現(xiàn)泛素特異性肽酶43（Ubiquitin-Specific Peptidase 43，）是影響豬肉色性狀的基因。Sun 等利用雞60 K SNP 芯片對(duì)“北京油雞×科寶肉雞”F代資源群體的肉質(zhì)性狀進(jìn)行全基因組關(guān)聯(lián)分析，篩選出以I 型膠原蛋白A 2（Collagen Type I Alpha 2 Chain，）為代表的候選基因與肉色顯著相關(guān)。然而，利用全基因組關(guān)聯(lián)分析對(duì)地方特色雞種肉色性狀的研究報(bào)道較少。儋州雞作為海南省特有的海島型地方雞遺傳資源，具有肉質(zhì)好、營養(yǎng)高等諸多優(yōu)異種質(zhì)特性，市場前景廣闊，在儋州雞產(chǎn)業(yè)化開發(fā)過程中，應(yīng)用生物育種技術(shù)開展肉質(zhì)性狀的研究和相關(guān)分子標(biāo)記的開發(fā)具有重要意義。本研究利用10×重測序技術(shù)獲得的基因型數(shù)據(jù)對(duì)儋州雞群體進(jìn)行GWAS 研究，以期識(shí)別影響儋州雞肉色性狀的分子標(biāo)記，探討遺傳標(biāo)記的分布規(guī)律，尋找相關(guān)的候選基因，為儋州雞選種育種及開發(fā)利用奠定理論基礎(chǔ)。

1 材料與方法

1.1 實(shí)驗(yàn)動(dòng)物與飼養(yǎng)管理 隨機(jī)選取600 只同一批次孵化的1 日齡儋州雞（中國熱帶農(nóng)業(yè)科學(xué)院熱帶作物品種資源研究所提供），按照《NY/T33-2004 雞飼養(yǎng)標(biāo)準(zhǔn)》進(jìn)行飼養(yǎng)。實(shí)驗(yàn)期分為育雛期（0～35 日齡）和育成期（35～90 日齡），基礎(chǔ)日糧組成及營養(yǎng)成分見表1。

表1 實(shí)驗(yàn)日糧組成及營養(yǎng)水平（干物質(zhì)基礎(chǔ)）

1.2 樣品采集和分析

1.2.1 血樣采集及DNA 提取 60 日齡時(shí)，翅靜脈采血5 mL，EDTA 抗凝處理，其中1 mL 使用動(dòng)物血液/細(xì)胞/組織基因組DNA 提取試劑盒（天根生化科技（北京）有限公司）提取基因組DNA，紫外分光光度計(jì)（NanoDrop 2000，美國賽默飛世爾Thermo 公司）測量OD 值，用1%瓊脂糖凝膠電泳檢測完整性，Qubit 定量的濃度≥15 ng/μL，由杭州聯(lián)川生物技術(shù)股份有限公司進(jìn)行10×深度的全基因組重測序。

1.2.2 肉色性狀測定 90 日齡時(shí)，隨機(jī)挑選200 只雞進(jìn)行屠宰。胸肌顏色以CIELAB 體系進(jìn)行測定，使用M-600D 色度計(jì)（柯尼卡美能達(dá)辦公系統(tǒng)（中國）有限公司）檢測肌肉的L（亮度）、a（紅度）、b（黃度），測定3 個(gè)不同區(qū)域的色差值，最后求平均值。測得數(shù)據(jù)經(jīng)過Excel 2016 整理，進(jìn)行GWAS 分析。

1.3 測序數(shù)據(jù)的質(zhì)控 使用FastP（0.19.3）對(duì)原始數(shù)據(jù)（RawData）進(jìn)行預(yù)處理。主要步驟為：去除接頭（Adapter）；去除含有N（N 表示無法確定堿基信息）的比例大于5% 的reads；去除低質(zhì)量reads；利用BWA 將CleanData 的數(shù)據(jù)比對(duì)到參考基因組上，并使用GATK 軟件進(jìn)行SNP 的檢測，并對(duì)檢測到的變異位點(diǎn)進(jìn)行質(zhì)量過濾。SNP 的質(zhì)量過濾標(biāo)準(zhǔn)：QD>2.0、FS>60.0、MQ>40.0、MQRankSum>-12.5、ReadPosRankSum>-8.0、SOR>3.0。群體過濾標(biāo)準(zhǔn)：去除最小等位基因頻率小于5%（MAF<0.05）的位點(diǎn)和缺失率大于5%的位點(diǎn)（geno>0.05）。

1.4 全基因組關(guān)聯(lián)分析

1.4.1 關(guān)聯(lián)分析模型 使用EMMAX（version beta-07）中的混合線性模型，對(duì)儋州雞肉色性狀與全基因組上SNP 進(jìn)行GWAS 分析，混合線性模型如下：

=+++

其中，表示儋州雞肉色性狀表型值；表示性狀的平均值；表示固定效應(yīng)矩陣，為固定效應(yīng)向量；表示剩余多基因效應(yīng)；表示表型值的隨機(jī)殘差，和均服從正態(tài)分布。

1.4.2 群體結(jié)構(gòu)分析 為消除連鎖不平衡對(duì)個(gè)體間遺傳相關(guān)估計(jì)的偏差，選用獨(dú)立的SNPs，利用Plink（1.9）進(jìn)行SNPs 篩選，以50 個(gè)位點(diǎn)為一個(gè)窗口，以5 個(gè)位點(diǎn)為步長，設(shè)為0.2。利用主成分分析（Principal Component Analysis，PCA）方法對(duì)群體結(jié)構(gòu)進(jìn)行分析，利用R 作圖。

1.4.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.5 基因注釋及候選基因篩選 利用Ensembl（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ù)庫比對(duì)，根據(jù)基因注釋情況分析可能的候選基因。

2 結(jié)果與分析

2.1 儋州雞表型和基因分型數(shù)據(jù)質(zhì)控 儋州雞肉色性狀統(tǒng)計(jì)結(jié)果見表2。質(zhì)控前有200 個(gè)個(gè)體，20 897 371 個(gè)SNP 位點(diǎn)，按質(zhì)控標(biāo)準(zhǔn)進(jìn)行質(zhì)控后，最終有200 個(gè)個(gè)體、12 424 543 個(gè)SNP 位點(diǎn)保留下來可用于后續(xù)的GWAS關(guān)聯(lián)分析。SNP 分布如圖1 所示，SNP 主要分布在轉(zhuǎn)錄本和內(nèi)含子區(qū)。

圖1 SNP（A）和InDel（B）位置分布

表2 儋州雞肉色性狀描述性統(tǒng)計(jì)

2.2 群體結(jié)構(gòu)分析 由圖2 可知，實(shí)驗(yàn)群體不存在分層現(xiàn)象，群體未出現(xiàn)家系分化，適宜進(jìn)行后續(xù)的GWAS 分析。

圖2 儋州雞群體主成分分析

2.3 儋州雞體尺性狀全基因組關(guān)聯(lián)分析結(jié)果 本研究采用EMMAX 軟件對(duì)儋州雞肉色性狀進(jìn)行GWAS 分析，校正方法為Bonferroni，確定肉色性狀全基因組水平顯著性閾值為8.049×10（0.01/12 424 543），潛在閾值為8.049×10（0.1/12 424 543）。由圖3 可知，與肉色性狀相關(guān)的全基因組和潛在顯著水平的SNP 位點(diǎn)分別為6 個(gè)和13 個(gè)。其中，與肉色b關(guān)聯(lián)的SNP 有10 個(gè)，分布在1、2、4 號(hào)染色體上；與肉色a關(guān)聯(lián)的SNP 有4 個(gè)，分布在1、2、5 號(hào)染色體上；與肉色L關(guān)聯(lián)的SNP 有5 個(gè)，分布在1、2 號(hào)和Z 染色體上。

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

2.4 候選基因功能注釋 由表3 可知，利用NCBI 網(wǎng)站，對(duì)SNPs 進(jìn)行注釋。在各顯著和潛在顯著關(guān)聯(lián)的SNPs上下游區(qū)域內(nèi)，共發(fā)現(xiàn)67 個(gè)基因，其中20 個(gè)是未表征基因。使用GO 數(shù)據(jù)庫和KEGG 數(shù)據(jù)庫對(duì)其進(jìn)行生物信息學(xué)分析，發(fā)現(xiàn)基因可作為儋州雞肉色性狀候選基因（表4）。

表3 肉色性狀的全基因組和潛在顯著關(guān)聯(lián)SNP 位點(diǎn)

表4 重要候選基因

3 討論

肉色是衡量肉品質(zhì)的重要指標(biāo)之一，能夠反映肌肉的生理生化特點(diǎn)及微生物學(xué)變化。肌肉中色素含量、分布及其化學(xué)狀態(tài)均決定了肉色的深淺及分布的均勻度，同時(shí)肌肉色素對(duì)肉色深淺的影響還受到pH 的影響。肌肉色素包括肌紅蛋白、血紅蛋白和微量有色代謝物（維生素B、黃素蛋白、細(xì)胞色素等），其中肌紅蛋白（Mb）的含量最高，是肌肉色素的基本組成成分，肌紅蛋白有3 種存在形式：還原型肌紅蛋白（紫紅）、氧合肌紅蛋白（鮮紅）、高鐵肌紅蛋白（褐色），因此肌紅蛋白的化學(xué)結(jié)構(gòu)對(duì)肉色有較大的影響。另外，線粒體活性、肌肉結(jié)構(gòu)完整性、脂質(zhì)氧化、糖酵解和pH有關(guān)的肌肉代謝機(jī)制影響著肌紅蛋白的化學(xué)結(jié)構(gòu)。這些機(jī)制通過影響肌紅蛋白的氧化、氧化和還原反應(yīng)，從而影響肌肉維持顏色穩(wěn)定的能力。在本研究中，基于10×重測序?qū)僦蓦u肉色性狀進(jìn)行GWAS 分析，在全基因組及潛在顯著水平中發(fā)現(xiàn)19 個(gè)SNP 與肉色性狀相關(guān)。

肌肉中的肌紅蛋白作為肌肉儲(chǔ)存氧的結(jié)構(gòu)物質(zhì)，是由一條肽鏈和一個(gè)血紅素輔基組成的結(jié)合蛋白，肌紅蛋白的構(gòu)象及合成代謝與肌肉氧化還原密切相關(guān)。Joseph 等研究發(fā)現(xiàn)在宰后肌肉中添加蘋果酸、丙酮酸、琥珀酸鹽和乳酸鹽等糖酵解或三羧酸循環(huán)中間代謝產(chǎn)物，能夠刺激肌肉再生煙酰胺腺嘌呤二核苷酸（Nicotinamide Adenine Dinucleotide，NADH），進(jìn)而增加高鐵肌紅蛋白還原程度，影響著肉色。

本實(shí)驗(yàn)GWAS 預(yù)測到位于儋州雞1 號(hào)染色體上?；o酶A 合成酶短鏈家族成員3（Acyl-CoA Synthetase Short-Chain Family Member 3，）候選基因作為一種線粒體酶可催化丙酸等短鏈脂肪酸轉(zhuǎn)化為CoA 復(fù)合物。能夠激活丙酰輔酶A 合成酶將丙酸鹽轉(zhuǎn)化為丙酰輔酶A，丙酰輔酶A 羧化酶將丙酰輔酶A轉(zhuǎn)化為甲基丙二酰輔酶A，甲基丙二酰輔酶A 再通過甲基丙二酸單酰輔酶A 變位酶異構(gòu)化為琥珀酰輔酶A，琥珀酰輔酶A 進(jìn)入三羧酸循環(huán)并釋放能量。由此推測是與肉色相關(guān)的候選基因，但相關(guān)功能性仍需進(jìn)一步驗(yàn)證。

研究發(fā)現(xiàn)，肌內(nèi)脂肪含量也影響著肉色，主要是因?yàn)橹炯?xì)胞是白色的，肌內(nèi)脂肪含量越高，肌肉亮度越大，即L值越大；pH 值下降過程中，肉表面對(duì)光的反射增加，也會(huì)產(chǎn)生淺色的外觀；同時(shí)肌內(nèi)脂肪含量受到脂代謝調(diào)控。本實(shí)驗(yàn)GWAS 預(yù)測到位于儋州雞Z 號(hào)染色體上乙酰輔酶A ?；D(zhuǎn)移酶2（Acetyl-Coenzyme A Acyltransferase 2，）是脂肪酸-氧化途徑關(guān)鍵酶，可通過催化脂肪酸-氧化參與脂肪酸降解及延伸，在脂肪酸代謝中起著重要作用。脂肪酸在-氧化時(shí)產(chǎn)生乙酰輔酶A，在作用下縮合生成乙酰乙酰輔酶A，經(jīng)過一系列酶促反應(yīng)進(jìn)入三羧酸循環(huán)徹底氧化，并釋放能量，維持細(xì)胞活性。李文娟利用基因組芯片技術(shù)對(duì)肉品質(zhì)相關(guān)的脂肪性狀功能基因進(jìn)行篩選和分析，發(fā)現(xiàn)基因是影響脂類代謝的關(guān)鍵基因。De Boer 等采用基因組芯片對(duì)不同小鼠脂質(zhì)代謝組織進(jìn)行了表達(dá)譜分析，檢測到了等基因存在差異表達(dá)。Doi 等研究發(fā)現(xiàn)，脂蛋白脂肪酶（）受體基因可以使mRNA 表達(dá)水平增加，引起大鼠脂肪代謝異常。以上研究表明，三羧酸循環(huán)和脂肪代謝對(duì)肌肉色素沉積和代謝等生物過程起的重要作用進(jìn)而影響肌肉色澤的改變。綜上推測基因可列為影響儋州雞肉色性狀的重要候選基因。

4 結(jié)論

綜上所述，利用10× 重測序?qū)僦蓦u肉色性狀進(jìn)行全基因組關(guān)聯(lián)分析，發(fā)現(xiàn)與肉色性狀相關(guān)的全基因組和潛在顯著水平的SNP 位點(diǎn)分別為6 個(gè)和13 個(gè)。經(jīng)過KEGG 通路分析和GO 注釋分析后，推測和基因可能是與肉色相關(guān)的候選基因。這些結(jié)果將為儋州雞的育種工作提供候選的分子標(biāo)記，為地方品種雞標(biāo)記輔助選擇提供參考依據(jù)。