金茅黑鸡肉品质性状的关联转录组分析

李国辉，谢恺舟，戴国俊，张跟喜，张 涛，杨国明，陈 兰，韩 威，苏一军，王金玉*

（1.扬州大学动物科学与技术学院，教育部农业与农产品安全国际合作联合实验室，江苏扬州 225009；2.江苏省家禽科学研究所，江苏扬州 225125；3.江苏三得利牧业发展有限公司，江苏常州 213200）

肉品质性状是复杂性状，具有遗传力低、易受环境影响等特点。一般影响复杂性状的基因都是由多个基因控制，而每个基因对性状的贡献率相对较小，并且基因与基因之间、环境与基因之间还可能存在各种互作关系。研究人员通过不同途径来解析性状形成的分子机制。首先通过候选基因关联分析的方法鉴定了一些影响鸡肉品质的候选基因，如鸡肌内脂肪（Intramuscular Fat，IMF）含量的候选基因H-FABP和A-FABP[1]。随着测序技术的发展，科研人员应用GWAS 技术鉴定了影响鸡腹脂、体重和生长、免疫和繁殖性状的基因和染色体区域，初步揭示这些性状形成的分子机制[2-3]。Kong 等[4]利用全基因重测序技术对2 个品系的肉鸡肌肉进行研究，发现APOB、ATM、MIA3、MICAL2等基因表达可能影响肉色的形成，是影响肉鸡肌肉肉色的重要候选基因。然而，GWAS 发现的位点大多位于非编码区，如何从功能上解释这些位点的调控作用一直是人们关注的焦点。以芯片为基础的基因表达谱研究能够分析某一状态下的细胞或组织的mRNA 表达丰度，可以高通量地研究基因在转录组水平的表达量。通过全基因组表达谱研究发现了存在大量的差异表达基因（Differentially Expressed Genes，DEGs），但由于其覆盖度较低，难以检测低丰度的基因表达。

转录组测序（RNA-seq）可全面检测特定组织的整体转录水平，研究基因转录表达与性状之间的关系，在一定程度上弥补了表达谱芯片研究的不足，能够更大范围地准确检测基因表达水平。在此基础上，进行生物信息学分析，最后对大样本群体的表型和基因表达量进行关联分析，可以大规模挖掘基因组的特异分子标记。本研究对金茅黑鸡的胸肌进行转录组测序并和肉品质性状进行关联分析，筛选一系列影响肉品质性状的重要候选基因并分析这些基因的分子功能，为金茅黑鸡的品种选育提供分子理论依据。

1 材料与方法

1.1 实验材料 选取来自24 个家系同一批次的金茅黑鸡A 系母鸡72 只（由江苏三得利牧业发展有限公司金茅黑鸡资源场提供），保证鸡群在相同的饲养管理条件下自由采食，饲养至14 周龄，宰前禁食后称重，放血后迅速拨开皮肤，未用水烫毛，屠宰后取其胸肌肌肉组织，装入无酶冻存管后立即放入液氮中，随后保存于-80℃超低温冰箱中，直到RNA 的抽提，用于后续的转录组测序。在保鲜库中再分别取胸肌进行肉品质指标的测定。

1.2 胸肌肉品质测定 在屠宰现场测定胸肌肉色、系水力、剪切力、pH。肉色和pH：采用肉色测定仪（德国，OPTO-STAR）、MATTHAUS 酮体肉质PH 直测仪（德国，PH-STAR）直接插入胸肌肉表面0.5～1.0 cm 下测定，每个肉样在不同位置测定3 次，取其平均值。系水力：取新鲜肉样称重(W1)，将肉样上、下各垫16 层滤纸，滤纸外层各放一块硬质塑料板，置于铜环允许膨胀仪平台上，加压 35 kg，持续 5 min，撤除压力后称重（W2）。系水力=（W1－W2）/W1×100%。剪切力：取新鲜肉样块，修成宽1 cm、厚0.5 cm 长条肉样，用C—LM 3 型肌肉嫩度仪测定，每个肉样剪切3 次，取其平均值。

1.3 转录组测序与数据分析 使用TRI Reagent（Invitrogen，美国）提取肌肉样品的总RNA；Nanodrop ND-1000 分光光度计（ThermoFisher Scientific，Waltham，美国）和Agilent 2100 Bioanalyzer（Agilent Technologie，上海）测定RNA 浓度和质量。PCR 进一步扩增以构建cDNA 文库。使用百迈客生物科技有限公司（中国，北京）的HiSeq 2000 测序系统（Ilumina）进行每个文库的测序。在测序质量控制后，获得总Clean Data，用于后续的基因表达量分析；采用DESeq 软件包分析差异表达基因，在差异表达基因检测过程中，将Fold Change≥2 且FDR＜0.05（False Discovery Rate）作为筛选标准。

1.3.1 基因表达量的相关分析 72 个样本采用Cufflinks软件的Cuffdiff 组件，通过Mapped Reads 在基因上的位置信息，对转录本和基因的表达水平进行定量。为了让片段数目能真实反映转录本表达水平，需要对样品中Mapped Reads 的数目和转录本长度进行归一化。采用FPKM（Fragments Per Kilobase of transcript per Million fragments mapped）作为衡量转录本或基因表达水平的指标，FPKM 计算公式：

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

1.3.2 基因表达量和性状（GEM）的关联分析 利用R（version 3.1.1）中stat 软件包转录本表达量与肉品质性状做线性回归分析，并计算关联显著性。模型为 y=Xα+Qβ+e，其中X 为基因表达量，y 为表型值，Q 为广义线性模型的中的第一和第二主成分矩阵，e 为随机残差效应向量α、β为关联矩阵；最终将候选目标性状关联到某个转录本(基因)上。

1.3.3 生物信息学分析 群体的进化树分析采用ape 软件和neighbor-joining 算法；利用DAVID 在线分析软件（https://david.ncifcrf.gov/summary.jsp）对候选基因进行GO（Gene Ontology）（http://www.geneontology.org）和KEGG 注释以及Pathway 富集分析（Kyoto Encyclopedia Of Genes And Genomes）（http://www.genome.jp/kegg），以P＜0.05作为 term 显著富集的标准。

2 结果

2.1 RNA-seq 数据产出统计 本研究中，72 个胸肌样本的转录组经过测序和质量控制，共获得363.09 Gb Clean Data，各样品Clean Data 均达到4.22 Gb 以上，Q30 碱基百分比在85.53%及以上，满足后续基因表达量的分析要求。

2.2 基因表达量分析 根据单个样品基因表达水平分布的离散程度，可以直观比较不同样品的整体基因表达水平。如图1 所示，同一条件的不同生物学重复样品的箱子的离散程度比较接近，而不同条件的样品的箱子高度存在明显差别，说明该生物学重复实验比较成功，满足后续关联分析的要求。

图1 各样品FPKM 箱线图

2.3 遗传进化分析 基于72 个样品RNA-seq 测序，对所有的SNP 进行分析，通过ape 软件和neighborjoining 算法，计算得到群体的进化树（图2）。根据测序所得的进化树分析与实际系谱记录资料高度吻合，证明测序数据的准确性。

图2 群体进化树分析

2.4 基因表达量和性状的关联分析 关联分析发现（表1），与金茅黑鸡14 周龄肉色显著相关的基因5 个，系水力17 个，pH57 个，剪切力33 个。

表1 各性状显著关联基因数目统计

2.4.1 肉色、系水力候选基因分析 将肉色、系水力与所测全基因表达量进行关联分析，共得到22 个差异显著基因（表2）。对肉色、系水力显著相关的表达基因GO 富集发现，这些基因如TXK、LHX9、NOBOX、UNCX、SPDEF等多富集在发育过程（GO:0032502）、组织发育（GO:0009888）、神经发育（GO:0048666）、细胞芳香族化合物的代谢过程（GO:0006725）、细胞含氮化物的代谢过程（GO:0034641）等细胞组分和生物学过程。

2.4.2 pH 候选基因分析 对影响胸肌pH 显著表达的57个基因进行GO 和KEGG 富集分析发现，这些基因显著富集在生物学过程中的腺体发育（GO:0048732）、类固醇激素的应答（GO:004854）、脂类的应答（GO:0033993）、核酸的转录活性（GO:0001071）、含氮化合物的代谢调节（GO:0051171）等过程中。KEGG 富集发现，NPY、POMC、FGF20、MGAT4C和GUK1等基因显著富集于脂肪因子信号通路（gga04920）、神经受体-配体相互作用（gga04080）等通路（图3）。

表2 肉色、系水力显著相关的基因统计

2.4.3 剪切力候选基因分析 对影响肌肉剪切力显著表达的33 个基因进行GO 和KEGG 富集分析发现（图4），这些基因显著富集在内部组织构造的形成（GO:0048646）、细胞内信号传导调节（GO:1902531）、有机物的应答（GO:0010033）、细胞骨架（GO:0005856）、应激反应（GO:0006950）等生物学过程或细胞组分中。KEGG 富集发现PLA2G4E、EDAR、HS6ST3、HSPB7和HSPB2等基因富集在代谢通路（hsa01100）、细胞因子受体相互作用（hsa04060）、脂肪的消化和吸收（hsa04975）等通路。

图3 pHGO 和KEGG pathway 富集分析

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

3 讨 论

肉色是鸡重要的肉品质性状之一，肉色不但影响消费者的心理，还与鸡肉的pH、系水力和保存期等有关[5-6]。在肉色的候选基因中，Plekhsl（Plcckstrin Homology Domain Containing Family S Member l）是血小板-白细胞C 激酶底物同源区域的S 家族的l 号成员，虽然目前未能详细了解其全部功能，但已有文献表明Plekhsl 主要与信号传递、膜运输、细脆骨架的构成有关，许多参与信号传递的分子中都含有pH 结构域，它们能够直接介导其宿主蛋白靶向细胞膜，并与膜上的磷脂酰肌醇类相结合[7]。MUC 黏蛋白是大分子糖蛋白，富含载有大量O-低聚糖链的丝氨酸、苏氨酸、脯氨酸残的可变联重复序列，主要由黏膜杯状细胞分泌，具有润滑和保护消化道上皮作用[8]。MUC 还与消化道上皮细胞的更新、分化、完整性有密切关系[9]。徐春瑛[10]采用SNP 芯片技术发现MUC4基因与猪的系水力相关。此外，新发现的与肉色相关基因有MUM1、LYVE1，这些基因的功能作用未见在家禽肉品质的研究中报道，其功能有待进一步研究。

动物屠宰后，肌肉组织形态和化学构成发生改变，影响肌肉束缚水分的能力[11]。肌肉水分含量的维持对肉品质和风味有着重要影响。在系水力的候选基因中，LHX9是许多器官（性腺、四肢、心脏和神经系统）发育所必需的LIM-同源域基因家族的成员。LHX9被认为是肢体发育所需的成纤维细胞生长因子（FGF）和Sonic hedgehog 信号的整合子，与肌纤维的形成和发育密切相关[12-13]。SPDEF 是ETS 转录因子家族的新成员，SPDEF 转录因子受到抑制可导致与TGF-β相关的信号通路的基因的转录受到影响，而TGF-β及其超家族成员肌肉生长抑制素（Myostatin，MSTN）对骨骼肌的生长发育、损伤再生、运动生理及病理生理都有重要调节作用[14-15]，GWAS 关联分析表明SPDEF基因是猪的背膘厚的候选基因[16]。GRXCR1（Glutaredoxin Cysteine-rich 1 Protein）编码进化上保守的富含半胱氨酸的蛋白质，其具有与谷氧还蛋白家族蛋白质的序列相似性。GRXCR1 和GRXCR2 的结构、功能分析表明，这些蛋白质在肌动蛋白富含丝状结构和C 端富含半胱氨酸结构域介导多聚体形成中起保守作用，并调节纤毛细胞发育[17]。以往的报道表明，肌肉骨架蛋白降解对系水力也有重要影响[18]。

肌肉pH 是反映动物屠宰后肌糖原酵解速度和强度的指标，也是肉品质测定的重要指标之一。动物宰杀放血后由于呼吸作用和血液循环停止，糖原、葡萄糖和6-磷酸葡萄糖通过无氧糖酵解转化为乳酸，乳酸的积累以及三磷酸腺苷的水解释放氢离子，从而使pH 下降，影响肉的持水性、嫩度和色泽。本研究在影响pH 的候选基因中发现NPY、POMC基因被富集在脂肪细胞因子信号通路（Adipocytokine Signaling Pathway），与脂肪的形成和代谢密切相关。研究表明，该通路上NPY、POMC等一系列基因的表达变化趋向于刺激牛采食、抑制能量分解代谢、提高能量摄入和脂肪沉积的方向变化，进而影响肉牛胴体品质以及能量平衡与脂肪代谢之间的关系[17]。FGF20基因被富集在MAPK 信号通路，在此通路中，FGF20是成纤维细胞生长因子（FGFs）家族的成员[19]，主要用于调节骨骼肌的生长和发育。MGAT4C和GUK1基因被显著富集在代谢通路，MGAT4C 属于单酰甘油酰基转移酶MGAT 家族，MGAT 在调节机体的脂类代谢、糖代谢以及系统能量平衡中发挥重要作用[20]。GUK1（Guanylate Kinase 1）作为cGMP 循环的一部分，GUK1 催化GMP 向GTP的转化，是cGMP 磷酸二酯酶水解后再生所必需[21]，在生物体内对控制环磷酸核苷浓度，细胞内外信息传导以及调节生命活动等起着举足轻重的作用[22]。因此推测，NPY、POMC、MGAT4C、GUK1等基因在调节脂类代谢、糖代谢过程和磷酸二酯酶水解中进一步影响肌肉pH 的变化。

肌肉嫩度是肉品质的首要指标，常用剪切切力值大小来表示。肌纤维直径和肌原纤维大小以及肌间脂肪含量是影响剪切力的主要因素。外异蛋白受体基因（Ectodysplasin A receptor，EDAR）作为影响剪切力的候选基因。以前研究表明，Wnt/β-catenin 和/activin/BMP 信号通路调控EDA和EDAR的表达，Wnt10b 是EDA 通路上的靶基因。TGF-β信号通路上的FGF20 调控皮肤成纤维细胞的聚集，影响毛囊和毛发纤维的形成[23-24]。关联分析表明，EDAR基因在肌肉中表达和剪切力存在显著关系，推测可能也影响肌肉纤维的形成；HS6ST3基因（硫酸乙酰肝素6-O-磺基转移酶3）是磺基转移酶（Sulfotransferase）家族的一员，参与调节肌肉中的脂质代谢[25]。Wilsie 等[26]研究发现，硫酸乙酰肝素蛋白聚糖促进跨血浆脂肪酸跨膜转运脂肪细胞，从而有助于细胞内的脂质积累进而影响肌肉中脂肪的沉积。将EDAR基因和HS6ST3基因作为影响金茅黑鸡肌肉剪切力的主要候选基因，其功能有待进一步研究。