李國(guó)輝，謝愷舟，戴國(guó)俊，張跟喜，張 濤，楊國(guó)明，陳 蘭，韓 威，蘇一軍，王金玉*

（1.揚(yáng)州大學(xué)動(dòng)物科學(xué)與技術(shù)學(xué)院，教育部農(nóng)業(yè)與農(nóng)產(chǎn)品安全國(guó)際合作聯(lián)合實(shí)驗(yàn)室，江蘇揚(yáng)州 225009；2.江蘇省家禽科學(xué)研究所，江蘇揚(yáng)州 225125；3.江蘇三得利牧業(yè)發(fā)展有限公司，江蘇常州 213200）

肉品質(zhì)性狀是復(fù)雜性狀，具有遺傳力低、易受環(huán)境影響等特點(diǎn)。一般影響復(fù)雜性狀的基因都是由多個(gè)基因控制，而每個(gè)基因?qū)π誀畹呢暙I(xiàn)率相對(duì)較小，并且基因與基因之間、環(huán)境與基因之間還可能存在各種互作關(guān)系。研究人員通過(guò)不同途徑來(lái)解析性狀形成的分子機(jī)制。首先通過(guò)候選基因關(guān)聯(lián)分析的方法鑒定了一些影響雞肉品質(zhì)的候選基因，如雞肌內(nèi)脂肪（Intramuscular Fat，IMF）含量的候選基因H-FABP和A-FABP[1]。隨著測(cè)序技術(shù)的發(fā)展，科研人員應(yīng)用GWAS 技術(shù)鑒定了影響雞腹脂、體重和生長(zhǎng)、免疫和繁殖性狀的基因和染色體區(qū)域，初步揭示這些性狀形成的分子機(jī)制[2-3]。Kong 等[4]利用全基因重測(cè)序技術(shù)對(duì)2 個(gè)品系的肉雞肌肉進(jìn)行研究，發(fā)現(xiàn)APOB、ATM、MIA3、MICAL2等基因表達(dá)可能影響肉色的形成，是影響肉雞肌肉肉色的重要候選基因。然而，GWAS 發(fā)現(xiàn)的位點(diǎn)大多位于非編碼區(qū)，如何從功能上解釋這些位點(diǎn)的調(diào)控作用一直是人們關(guān)注的焦點(diǎn)。以芯片為基礎(chǔ)的基因表達(dá)譜研究能夠分析某一狀態(tài)下的細(xì)胞或組織的mRNA 表達(dá)豐度，可以高通量地研究基因在轉(zhuǎn)錄組水平的表達(dá)量。通過(guò)全基因組表達(dá)譜研究發(fā)現(xiàn)了存在大量的差異表達(dá)基因（Differentially Expressed Genes，DEGs），但由于其覆蓋度較低，難以檢測(cè)低豐度的基因表達(dá)。

轉(zhuǎn)錄組測(cè)序（RNA-seq）可全面檢測(cè)特定組織的整體轉(zhuǎn)錄水平，研究基因轉(zhuǎn)錄表達(dá)與性狀之間的關(guān)系，在一定程度上彌補(bǔ)了表達(dá)譜芯片研究的不足，能夠更大范圍地準(zhǔn)確檢測(cè)基因表達(dá)水平。在此基礎(chǔ)上，進(jìn)行生物信息學(xué)分析，最后對(duì)大樣本群體的表型和基因表達(dá)量進(jìn)行關(guān)聯(lián)分析，可以大規(guī)模挖掘基因組的特異分子標(biāo)記。本研究對(duì)金茅黑雞的胸肌進(jìn)行轉(zhuǎn)錄組測(cè)序并和肉品質(zhì)性狀進(jìn)行關(guān)聯(lián)分析，篩選一系列影響肉品質(zhì)性狀的重要候選基因并分析這些基因的分子功能，為金茅黑雞的品種選育提供分子理論依據(jù)。

1 材料與方法

1.1 實(shí)驗(yàn)材料 選取來(lái)自24 個(gè)家系同一批次的金茅黑雞A 系母雞72 只（由江蘇三得利牧業(yè)發(fā)展有限公司金茅黑雞資源場(chǎng)提供），保證雞群在相同的飼養(yǎng)管理?xiàng)l件下自由采食，飼養(yǎng)至14 周齡，宰前禁食后稱(chēng)重，放血后迅速撥開(kāi)皮膚，未用水燙毛，屠宰后取其胸肌肌肉組織，裝入無(wú)酶凍存管后立即放入液氮中，隨后保存于-80℃超低溫冰箱中，直到RNA 的抽提，用于后續(xù)的轉(zhuǎn)錄組測(cè)序。在保鮮庫(kù)中再分別取胸肌進(jìn)行肉品質(zhì)指標(biāo)的測(cè)定。

1.2 胸肌肉品質(zhì)測(cè)定 在屠宰現(xiàn)場(chǎng)測(cè)定胸肌肉色、系水力、剪切力、pH。肉色和pH：采用肉色測(cè)定儀（德國(guó)，OPTO-STAR）、MATTHAUS 酮體肉質(zhì)PH 直測(cè)儀（德國(guó)，PH-STAR）直接插入胸肌肉表面0.5～1.0 cm 下測(cè)定，每個(gè)肉樣在不同位置測(cè)定3 次，取其平均值。系水力：取新鮮肉樣稱(chēng)重(W1)，將肉樣上、下各墊16 層濾紙，濾紙外層各放一塊硬質(zhì)塑料板，置于銅環(huán)允許膨脹儀平臺(tái)上，加壓 35 kg，持續(xù) 5 min，撤除壓力后稱(chēng)重（W2）。系水力=（W1－W2）/W1×100%。剪切力：取新鮮肉樣塊，修成寬1 cm、厚0.5 cm 長(zhǎng)條肉樣，用C—LM 3 型肌肉嫩度儀測(cè)定，每個(gè)肉樣剪切3 次，取其平均值。

1.3 轉(zhuǎn)錄組測(cè)序與數(shù)據(jù)分析 使用TRI Reagent（Invitrogen，美國(guó)）提取肌肉樣品的總RNA；Nanodrop ND-1000 分光光度計(jì)（ThermoFisher Scientific，Waltham，美國(guó)）和Agilent 2100 Bioanalyzer（Agilent Technologie，上海）測(cè)定RNA 濃度和質(zhì)量。PCR 進(jìn)一步擴(kuò)增以構(gòu)建cDNA 文庫(kù)。使用百邁客生物科技有限公司（中國(guó)，北京）的HiSeq 2000 測(cè)序系統(tǒng)（Ilumina）進(jìn)行每個(gè)文庫(kù)的測(cè)序。在測(cè)序質(zhì)量控制后，獲得總Clean Data，用于后續(xù)的基因表達(dá)量分析；采用DESeq 軟件包分析差異表達(dá)基因，在差異表達(dá)基因檢測(cè)過(guò)程中，將Fold Change≥2 且FDR＜0.05（False Discovery Rate）作為篩選標(biāo)準(zhǔn)。

1.3.1 基因表達(dá)量的相關(guān)分析 72 個(gè)樣本采用Cufflinks軟件的Cuffdiff 組件，通過(guò)Mapped Reads 在基因上的位置信息，對(duì)轉(zhuǎn)錄本和基因的表達(dá)水平進(jìn)行定量。為了讓片段數(shù)目能真實(shí)反映轉(zhuǎn)錄本表達(dá)水平，需要對(duì)樣品中Mapped Reads 的數(shù)目和轉(zhuǎn)錄本長(zhǎng)度進(jìn)行歸一化。采用FPKM（Fragments Per Kilobase of transcript per Million fragments mapped）作為衡量轉(zhuǎn)錄本或基因表達(dá)水平的指標(biāo)，F(xiàn)PKM 計(jì)算公式：

FPKM=cDNA Fragments/Mapped Fragments×Transcript Length（bp）

1.3.2 基因表達(dá)量和性狀（GEM）的關(guān)聯(lián)分析 利用R（version 3.1.1）中stat 軟件包轉(zhuǎn)錄本表達(dá)量與肉品質(zhì)性狀做線(xiàn)性回歸分析，并計(jì)算關(guān)聯(lián)顯著性。模型為 y=Xα+Qβ+e，其中X 為基因表達(dá)量，y 為表型值，Q 為廣義線(xiàn)性模型的中的第一和第二主成分矩陣，e 為隨機(jī)殘差效應(yīng)向量α、β為關(guān)聯(lián)矩陣；最終將候選目標(biāo)性狀關(guān)聯(lián)到某個(gè)轉(zhuǎn)錄本(基因)上。

1.3.3 生物信息學(xué)分析 群體的進(jìn)化樹(shù)分析采用ape 軟件和neighbor-joining 算法；利用DAVID 在線(xiàn)分析軟件（https://david.ncifcrf.gov/summary.jsp）對(duì)候選基因進(jìn)行GO（Gene Ontology）（http://www.geneontology.org）和KEGG 注釋以及Pathway 富集分析（Kyoto Encyclopedia Of Genes And Genomes）（http://www.genome.jp/kegg），以P＜0.05作為 term 顯著富集的標(biāo)準(zhǔn)。

2 結(jié)果

2.1 RNA-seq 數(shù)據(jù)產(chǎn)出統(tǒng)計(jì) 本研究中，72 個(gè)胸肌樣本的轉(zhuǎn)錄組經(jīng)過(guò)測(cè)序和質(zhì)量控制，共獲得363.09 Gb Clean Data，各樣品Clean Data 均達(dá)到4.22 Gb 以上，Q30 堿基百分比在85.53%及以上，滿(mǎn)足后續(xù)基因表達(dá)量的分析要求。

2.2 基因表達(dá)量分析 根據(jù)單個(gè)樣品基因表達(dá)水平分布的離散程度，可以直觀比較不同樣品的整體基因表達(dá)水平。如圖1 所示，同一條件的不同生物學(xué)重復(fù)樣品的箱子的離散程度比較接近，而不同條件的樣品的箱子高度存在明顯差別，說(shuō)明該生物學(xué)重復(fù)實(shí)驗(yàn)比較成功，滿(mǎn)足后續(xù)關(guān)聯(lián)分析的要求。

圖1 各樣品FPKM 箱線(xiàn)圖

2.3 遺傳進(jìn)化分析 基于72 個(gè)樣品RNA-seq 測(cè)序，對(duì)所有的SNP 進(jìn)行分析，通過(guò)ape 軟件和neighborjoining 算法，計(jì)算得到群體的進(jìn)化樹(shù)（圖2）。根據(jù)測(cè)序所得的進(jìn)化樹(shù)分析與實(shí)際系譜記錄資料高度吻合，證明測(cè)序數(shù)據(jù)的準(zhǔn)確性。

圖2 群體進(jìn)化樹(shù)分析

2.4 基因表達(dá)量和性狀的關(guān)聯(lián)分析 關(guān)聯(lián)分析發(fā)現(xiàn)（表1），與金茅黑雞14 周齡肉色顯著相關(guān)的基因5 個(gè)，系水力17 個(gè)，pH57 個(gè)，剪切力33 個(gè)。

表1 各性狀顯著關(guān)聯(lián)基因數(shù)目統(tǒng)計(jì)

2.4.1 肉色、系水力候選基因分析 將肉色、系水力與所測(cè)全基因表達(dá)量進(jìn)行關(guān)聯(lián)分析，共得到22 個(gè)差異顯著基因（表2）。對(duì)肉色、系水力顯著相關(guān)的表達(dá)基因GO 富集發(fā)現(xiàn)，這些基因如TXK、LHX9、NOBOX、UNCX、SPDEF等多富集在發(fā)育過(guò)程（GO:0032502）、組織發(fā)育（GO:0009888）、神經(jīng)發(fā)育（GO:0048666）、細(xì)胞芳香族化合物的代謝過(guò)程（GO:0006725）、細(xì)胞含氮化物的代謝過(guò)程（GO:0034641）等細(xì)胞組分和生物學(xué)過(guò)程。

2.4.2 pH 候選基因分析 對(duì)影響胸肌pH 顯著表達(dá)的57個(gè)基因進(jìn)行GO 和KEGG 富集分析發(fā)現(xiàn)，這些基因顯著富集在生物學(xué)過(guò)程中的腺體發(fā)育（GO:0048732）、類(lèi)固醇激素的應(yīng)答（GO:004854）、脂類(lèi)的應(yīng)答（GO:0033993）、核酸的轉(zhuǎn)錄活性（GO:0001071）、含氮化合物的代謝調(diào)節(jié)（GO:0051171）等過(guò)程中。KEGG 富集發(fā)現(xiàn)，NPY、POMC、FGF20、MGAT4C和GUK1等基因顯著富集于脂肪因子信號(hào)通路（gga04920）、神經(jīng)受體-配體相互作用（gga04080）等通路（圖3）。

表2 肉色、系水力顯著相關(guān)的基因統(tǒng)計(jì)

2.4.3 剪切力候選基因分析 對(duì)影響肌肉剪切力顯著表達(dá)的33 個(gè)基因進(jìn)行GO 和KEGG 富集分析發(fā)現(xiàn)（圖4），這些基因顯著富集在內(nèi)部組織構(gòu)造的形成（GO:0048646）、細(xì)胞內(nèi)信號(hào)傳導(dǎo)調(diào)節(jié)（GO:1902531）、有機(jī)物的應(yīng)答（GO:0010033）、細(xì)胞骨架（GO:0005856）、應(yīng)激反應(yīng)（GO:0006950）等生物學(xué)過(guò)程或細(xì)胞組分中。KEGG 富集發(fā)現(xiàn)PLA2G4E、EDAR、HS6ST3、HSPB7和HSPB2等基因富集在代謝通路（hsa01100）、細(xì)胞因子受體相互作用（hsa04060）、脂肪的消化和吸收（hsa04975）等通路。

圖3 pHGO 和KEGG pathway 富集分析

圖4 剪切力GO 和KEGG pathway 富集分析

3 討 論

肉色是雞重要的肉品質(zhì)性狀之一，肉色不但影響消費(fèi)者的心理，還與雞肉的pH、系水力和保存期等有關(guān)[5-6]。在肉色的候選基因中，Plekhsl（Plcckstrin Homology Domain Containing Family S Member l）是血小板-白細(xì)胞C 激酶底物同源區(qū)域的S 家族的l 號(hào)成員，雖然目前未能詳細(xì)了解其全部功能，但已有文獻(xiàn)表明Plekhsl 主要與信號(hào)傳遞、膜運(yùn)輸、細(xì)脆骨架的構(gòu)成有關(guān)，許多參與信號(hào)傳遞的分子中都含有pH 結(jié)構(gòu)域，它們能夠直接介導(dǎo)其宿主蛋白靶向細(xì)胞膜，并與膜上的磷脂酰肌醇類(lèi)相結(jié)合[7]。MUC 黏蛋白是大分子糖蛋白，富含載有大量O-低聚糖鏈的絲氨酸、蘇氨酸、脯氨酸殘的可變聯(lián)重復(fù)序列，主要由黏膜杯狀細(xì)胞分泌，具有潤(rùn)滑和保護(hù)消化道上皮作用[8]。MUC 還與消化道上皮細(xì)胞的更新、分化、完整性有密切關(guān)系[9]。徐春瑛[10]采用SNP 芯片技術(shù)發(fā)現(xiàn)MUC4基因與豬的系水力相關(guān)。此外，新發(fā)現(xiàn)的與肉色相關(guān)基因有MUM1、LYVE1，這些基因的功能作用未見(jiàn)在家禽肉品質(zhì)的研究中報(bào)道，其功能有待進(jìn)一步研究。

動(dòng)物屠宰后，肌肉組織形態(tài)和化學(xué)構(gòu)成發(fā)生改變，影響肌肉束縛水分的能力[11]。肌肉水分含量的維持對(duì)肉品質(zhì)和風(fēng)味有著重要影響。在系水力的候選基因中，LHX9是許多器官（性腺、四肢、心臟和神經(jīng)系統(tǒng)）發(fā)育所必需的LIM-同源域基因家族的成員。LHX9被認(rèn)為是肢體發(fā)育所需的成纖維細(xì)胞生長(zhǎng)因子（FGF）和Sonic hedgehog 信號(hào)的整合子，與肌纖維的形成和發(fā)育密切相關(guān)[12-13]。SPDEF 是ETS 轉(zhuǎn)錄因子家族的新成員，SPDEF 轉(zhuǎn)錄因子受到抑制可導(dǎo)致與TGF-β相關(guān)的信號(hào)通路的基因的轉(zhuǎn)錄受到影響，而TGF-β及其超家族成員肌肉生長(zhǎng)抑制素（Myostatin，MSTN）對(duì)骨骼肌的生長(zhǎng)發(fā)育、損傷再生、運(yùn)動(dòng)生理及病理生理都有重要調(diào)節(jié)作用[14-15]，GWAS 關(guān)聯(lián)分析表明SPDEF基因是豬的背膘厚的候選基因[16]。GRXCR1（Glutaredoxin Cysteine-rich 1 Protein）編碼進(jìn)化上保守的富含半胱氨酸的蛋白質(zhì)，其具有與谷氧還蛋白家族蛋白質(zhì)的序列相似性。GRXCR1 和GRXCR2 的結(jié)構(gòu)、功能分析表明，這些蛋白質(zhì)在肌動(dòng)蛋白富含絲狀結(jié)構(gòu)和C 端富含半胱氨酸結(jié)構(gòu)域介導(dǎo)多聚體形成中起保守作用，并調(diào)節(jié)纖毛細(xì)胞發(fā)育[17]。以往的報(bào)道表明，肌肉骨架蛋白降解對(duì)系水力也有重要影響[18]。

肌肉pH 是反映動(dòng)物屠宰后肌糖原酵解速度和強(qiáng)度的指標(biāo)，也是肉品質(zhì)測(cè)定的重要指標(biāo)之一。動(dòng)物宰殺放血后由于呼吸作用和血液循環(huán)停止，糖原、葡萄糖和6-磷酸葡萄糖通過(guò)無(wú)氧糖酵解轉(zhuǎn)化為乳酸，乳酸的積累以及三磷酸腺苷的水解釋放氫離子，從而使pH 下降，影響肉的持水性、嫩度和色澤。本研究在影響pH 的候選基因中發(fā)現(xiàn)NPY、POMC基因被富集在脂肪細(xì)胞因子信號(hào)通路（Adipocytokine Signaling Pathway），與脂肪的形成和代謝密切相關(guān)。研究表明，該通路上NPY、POMC等一系列基因的表達(dá)變化趨向于刺激牛采食、抑制能量分解代謝、提高能量攝入和脂肪沉積的方向變化，進(jìn)而影響肉牛胴體品質(zhì)以及能量平衡與脂肪代謝之間的關(guān)系[17]。FGF20基因被富集在MAPK 信號(hào)通路，在此通路中，F(xiàn)GF20是成纖維細(xì)胞生長(zhǎng)因子（FGFs）家族的成員[19]，主要用于調(diào)節(jié)骨骼肌的生長(zhǎng)和發(fā)育。MGAT4C和GUK1基因被顯著富集在代謝通路，MGAT4C 屬于單酰甘油?；D(zhuǎn)移酶MGAT 家族，MGAT 在調(diào)節(jié)機(jī)體的脂類(lèi)代謝、糖代謝以及系統(tǒng)能量平衡中發(fā)揮重要作用[20]。GUK1（Guanylate Kinase 1）作為cGMP 循環(huán)的一部分，GUK1 催化GMP 向GTP的轉(zhuǎn)化，是cGMP 磷酸二酯酶水解后再生所必需[21]，在生物體內(nèi)對(duì)控制環(huán)磷酸核苷濃度，細(xì)胞內(nèi)外信息傳導(dǎo)以及調(diào)節(jié)生命活動(dòng)等起著舉足輕重的作用[22]。因此推測(cè)，NPY、POMC、MGAT4C、GUK1等基因在調(diào)節(jié)脂類(lèi)代謝、糖代謝過(guò)程和磷酸二酯酶水解中進(jìn)一步影響肌肉pH 的變化。

肌肉嫩度是肉品質(zhì)的首要指標(biāo)，常用剪切切力值大小來(lái)表示。肌纖維直徑和肌原纖維大小以及肌間脂肪含量是影響剪切力的主要因素。外異蛋白受體基因（Ectodysplasin A receptor，EDAR）作為影響剪切力的候選基因。以前研究表明，Wnt/β-catenin 和/activin/BMP 信號(hào)通路調(diào)控EDA和EDAR的表達(dá)，Wnt10b 是EDA 通路上的靶基因。TGF-β信號(hào)通路上的FGF20 調(diào)控皮膚成纖維細(xì)胞的聚集，影響毛囊和毛發(fā)纖維的形成[23-24]。關(guān)聯(lián)分析表明，EDAR基因在肌肉中表達(dá)和剪切力存在顯著關(guān)系，推測(cè)可能也影響肌肉纖維的形成；HS6ST3基因（硫酸乙酰肝素6-O-磺基轉(zhuǎn)移酶3）是磺基轉(zhuǎn)移酶（Sulfotransferase）家族的一員，參與調(diào)節(jié)肌肉中的脂質(zhì)代謝[25]。Wilsie 等[26]研究發(fā)現(xiàn)，硫酸乙酰肝素蛋白聚糖促進(jìn)跨血漿脂肪酸跨膜轉(zhuǎn)運(yùn)脂肪細(xì)胞，從而有助于細(xì)胞內(nèi)的脂質(zhì)積累進(jìn)而影響肌肉中脂肪的沉積。將EDAR基因和HS6ST3基因作為影響金茅黑雞肌肉剪切力的主要候選基因，其功能有待進(jìn)一步研究。