金茅黑雞肉品質(zhì)性狀的關(guān)聯(lián)轉(zhuǎn)錄組分析

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

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

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

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

1 材料與方法

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

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

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

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

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ì)性狀做線性回歸分析，并計算關(guān)聯(lián)顯著性。模型為 y=Xα+Qβ+e，其中X 為基因表達(dá)量，y 為表型值，Q 為廣義線性模型的中的第一和第二主成分矩陣，e 為隨機殘差效應(yīng)向量α、β為關(guān)聯(lián)矩陣；最終將候選目標(biāo)性狀關(guān)聯(lián)到某個轉(zhuǎn)錄本(基因)上。

1.3.3 生物信息學(xué)分析 群體的進(jìn)化樹分析采用ape 軟件和neighbor-joining 算法；利用DAVID 在線分析軟件（https://david.ncifcrf.gov/summary.jsp）對候選基因進(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)計 本研究中，72 個胸肌樣本的轉(zhuǎn)錄組經(jīng)過測序和質(zhì)量控制，共獲得363.09 Gb Clean Data，各樣品Clean Data 均達(dá)到4.22 Gb 以上，Q30 堿基百分比在85.53%及以上，滿足后續(xù)基因表達(dá)量的分析要求。

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

圖1 各樣品FPKM 箱線圖

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

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

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

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

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

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

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

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

圖3 pHGO 和KEGG pathway 富集分析

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

3 討 論

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

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

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

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