基于WGCNA分析挖掘调控莱芜猪肉色相关基因

Identification of genes regulating meat color in Laiwu pigs based on WGCNA

摘要

为筛选影响莱芜猪肉色的关键候选基因，本研究利用莱芜猪转录组数据，通过加权基因共表达网络分析（Weighted gene co-expression network analysis， WGCNA）对33头莱芜猪测定了45 min及24 h的肉色（包括L、a、b），对其背最长肌样本开展了转录组测序，挖掘影响猪肉色的关键基因。结合转录组数据和表型数据利用R语言WGCNA软件包开展了WGCNA分析，在剔除2个离群样本后利用31个样本数据构建基因共表达网络，筛选目标模块及模块中的关键基因。结果表明：1）通过WGCNA分析，共获得20个共表达模块，其中4个模块与2个及以上肉色性状显著相关（|相关性系数|≥0.35且显著性P<0.01），分别为Floral white、Dark green、Pink及Orangered 4模块。2）功能富集分析发现，目标模块内基因显著富集在线粒体结构和功能、蛋白修饰、脂肪沉积和代谢、肌肉结构和功能等相关的GO项和通路中。3）以模块基因表达量与模块特征值的相关性系数绝对值>0.85且基因表达量与至少2个肉色性状的相关性系数绝对值>0.3为标准筛选目标模块基因，并将筛选到的基因进行蛋白互作分析，最终从4个模块中筛选到8个影响肉色的候选基因MYO9A、NDUFA13、COX5B、MRPL12、RPL8、ROBO4、EGFL7和CIDEC。综上，本研究通过WGCNA分析鉴别到影响肉色的4个目标模块及8个候选基因，为猪肉色分子调控机制研究与肉色育种改良奠定了理论基础。

Abstract

To screen for key candidate genes that affect the color of pork in Laiwu pigs， this study utilized transcriptome data from Laiwu pigs and identified key genes affecting pork color using weighted gene co -expression network analysis （WGCNA）. The color of the meat， including L， a， and b， was measured at 45 m and 24 h postmortem for a total of 33 Laiwu pigs. RNA sequencing was conducted on the longissimus dorsi muscles of these pigs. Based on the transcriptome data and phenotype data， WGCNA analysis was conducted using the WGCNA package of R. After removing two outlier samples， 31 samples were used to construct a gene co-expression network for the identification of target modules and key genes. The results showed that： 1）Through WGCNA analysis， a total of 20 co-expression modules were identified. Among these modules， Floral white， Dark green， Pink， and Orangered 4 demonstrated significant correlations with two or more meat color traits， with a |correlation coefficient| ≥0.35 and a significance P value <0.01. 2）Functional enrichment analysis revealed that genes within the target module were significantly enriched in GO terms and pathways related to mitochondrial structure and function， protein modification， fat deposition and metabolism， muscle structure and function， etc. 3）Identification of key genes in the target module involved applying an absolute correlation coefficient greater than 0.85 between gene expression and module eigengene， as well as an absolute correlation coefficient greater than 0.3 between gene expression and meat color traits for at least two traits. Subsequently， these genes were carried out protein-protein interaction analysis. A total of eight candidate genes affecting meat color， including MYO9A， NDUFA13， COX5B， MRPL12， RPL8， ROBO4， EGFL7， and CIDEC， were obtained from the four target modules. In summary， this study identified four target modules and eight candidate genes regulating meat color in Laiwu pigs through WGCNA， laying the foundation for exploring the molecular regulatory mechanisms of pig meat color trait and making advancements in meat color improvement through breeding.

Graphical abstract

关键词

莱芜猪 / 肉色 / 转录组测序 / WGCNA

Key words

Laiwu pigs / meat color / RNA sequencing / WGCNA

中国是全球猪肉生产和消费大国，据报道，中国的猪肉消费量在总体肉类消费比例中超过50% ^［1］。猪肉在肉类食品市场中占据重要地位，而猪肉品质是养殖户最为关注的经济性状之一。肉色作为影响猪肉质性状的基本表型之一，是判定肉品质最直观的指标，能够在很大程度上决定消费者的购买欲望。因此，挖掘出影响肉色的关键基因对猪肉品质的改良具有重要意义。肉色性状作为一种微效多基因控制的数量性状，遗传力在0.28左右 ^［2］，属于中等遗传力。肌红蛋白作为肉色主要的显色物质，通常以氧肌红蛋白、脱氧肌红蛋白和肌红蛋白3种形式存在，这3种物质的比例以及氧化程度决定了肉质的颜色。李国辉等 ^［3］利用RNA-seq技术对14周龄的金茅黑鸡进行研究，挖掘出PLEKHS1、MUC等影响肉色的候选基因。Wang等 ^［4］对苏淮猪miR-22的多态性进行研究发现，该基因序列上游的G>A突变位点与猪肉颜色的红度显著相关；由于猪肉色与Ca²⁺的平衡密切相关，研究结果进一步证明miR-22对骨骼肌细胞内Ca²⁺的影响部分原因是由于抑制ELOVL6的表达引起的。廖印长等 ^［5］通过全基因组关联研究对宁乡猪的pH值以及肉色进行分析，发现15个与肉色显著相关的SNP位点，以及2个影响肉色的候选基因FGF10和GRM4，同时该研究分析发现候选基因对肉色和pH的调控主要通过肌肉发育以及脂肪沉积等路径。

加权基因共表达网络分析（Weighted gene co-expression network analysis， WGCNA）是一种系统的生物学算法，基于相同性状的基因在表达模式中具有共表达特征的假设，从构建的表达调控网络中找到并鉴定目标性状的相关基因，用于筛选候选生物标记或治疗的靶点 ^［6］。近年来，通过WGCNA分析探究畜禽重要经济性状调控分子机制的研究越来越多。荣丽云等 ^［7］运用WGCNA技术对大蒲莲猪和长白猪的肉质性状进行分析，共识别出两个关键模块和21个关键基因，其中COL5A3、COL6A2等7个基因在猪肉质性状形成中可能有重要作用。Zhao等 ^［8］通过WGCNA方法对28头杜洛克猪中选择的8个具有极高和极低滴水损失的个体进行分析，结果共筛选出AASS、BCKDHB、ALDH6A1、MUT和MCCC15个影响滴水损失的候选基因。Guo等 ^［9］基于WGCNA对2个黄羽鸡品种的肉色进行分析，鉴定出ATP5L、UQCR10、COX7C、CAV3、RBP4A、APOH等基因可能在肉色调控中起关键作用。

莱芜猪是我国华北型地方猪典型代表，具有繁殖力高，抗逆性强的优良特点，其肉色鲜红、肉质细嫩，与其他商品猪肉形成鲜明对比，视觉上直接体现出莱芜猪肉的高端品质。因此，在莱芜猪群体中开展挖掘影响肉色候选基因研究，解析猪肉色的分子机制，对商品猪等的肉色改良和分子育种具有重要意义。本研究以莱芜猪为研究对象，在RNA-seq的基础上利用WGCNA分析，筛选影响肉色的候选基因，解析猪肉色的分子机制，为猪肉色分子改良和分子育种奠定理论基础并提供科学依据。

1　材料与方法

1.1　试验动物

试验选用的33头莱芜猪由济南市莱芜猪种猪繁育有限公司提供。猪只采用相同饲粮和饲养管理措施进行饲养，当猪只体重达到约100 kg时，于2020年9—12月送至济南市东晟肉食品加工厂进行屠宰。屠宰前空腹24 h，采用电击麻醉、心脏放血的方式进行屠宰。取胸腰椎结合处背最长肌200 g用于肉色测定，另取0.2 g背最长肌样本，存放于液氮中用于mRNA转录组测序。

1.2　表型测定

参照国家标准NY/T 821-2019《猪肉品质测定技术规程》，采用Minolta CR-300色差仪测定肉色，测定指标包括红度（a）、黄度（b）和亮度（L）。屠宰后45~60 min内测定的肉色用L₁，a₁，b₁表示，将样本置于0~4 ℃保存至屠宰后（24 h ± 15 min）测定的肉色用L₂₄，a₂₄，b₂₄表示。每个样品在不同位点测量3次，最终测定结果为3次的平均值。

1.3　总RNA提取及mRNA转录组测序

取黄豆大小的背最长肌肉样，使用TRIZOL^® Reagent（Invitrogen，SanDiego，CA，USA）试剂盒提取组织中的总RNA。通过Nanodrop 2000测定RNA浓度和质量，保证样本OD₂₆₀/OD₂₈₀数值在1.8~2.1，OD₂₆₀/OD₂₃₀>2.0，并通过琼脂糖电泳检查RNA的完整性。转录组文库构建由诺禾致源科技股份有限公司完成，于Illumina HiSeq平台进行150 bp双端测序。

在分析之前首先对原始数据进行过滤（去除adaptor、含N和低质量的reads），并检查GC含量，得到Clean reads。通过HISAT2建立猪参考基因组（Sus scrofa 11.1）索引，并进行比对得到Reads在参考基因组的位置信息。采用subread软件中的Feature Counts工具对基因进行定量，并计算基因的FPKM（Fragments per kilobase of exon model per million mapped reads）值 ^［10］。在33个莱芜猪个体中，将FPKM 值在超过16个个体中均<1的基因剔除以提高数据的可靠性，利用剩余基因进行加权基因共表达网络构建。

1.4　加权基因共表达网络分析

利用R 4.1.2软件的WGCNA包，基于筛选后的FPKM值和肉色表型值进行加权基因共表达网络构建。首先对33个样本通过hclust函数进行层次聚类分析，以15 000的阈值剔除表达异常样本。基于筛选后的样本通过pick-Soft-Threshold（sft）函数计算两两基因的皮尔森相关系数以构建邻接矩阵，确定无尺度分布拓扑网络指数R²达到0.85作为最佳软阈值（β）。通过adjacency函数将邻接矩阵转换为拓扑重叠矩阵（Topology overlap measure，TOM），从而基于TOM矩阵构建加权基因共表达网络，将具有相似表达模式的基因聚类到相同的模块中，并为其分配颜色。

模块构建后，计算模块的特征值（Module eigengene，ME），并分别计算ME与肉色间的皮尔森相关性系数r和显著性P值，选择在2个及以上性状中|r|≥0.35且P<0.01的模块作为目标模块 ^［11］。随后对目标模块中的基因利用DAVID ^［12］数据库（https：//david.ncifcrf.gov/）进行GO和KEGG功能富集分析，以P<0.05为阈值筛选模块内基因显著富集的GO 条目（GO term）和通路（Pathway）。

1.5　蛋白互作网络构建及关键基因筛选

为了进一步从目标模块中筛选调控肉色的枢纽基因，计算模块中基因表达量和模块ME的相关系数MM值（Module membership）以及基因表达量与表型的皮尔森相关系数GS值（Gene significance），基于|MM|>0.85且|GS|>0.3的标准对目标模块的基因进行筛选，随后将筛选出的基因导入STRING数据库中（https：//cn.string-db.org），构建蛋白互作（Protein-protein interaction， PPI）网络，通过Cytoscape 3.9.1软件的cytoNCA插件计算PPI网络节点的Degree值并对top30的nodes进行可视化（不足30的全部展示），基于Degree值及基因功能筛选出各模块最终的关键基因。

2　结果与分析

2.1　表型数据统计

33头莱芜猪（胴体均值为（72.61±6.53）kg）的6个肉色性状测定结果的描述性统计见表1。结果显示，33头莱芜猪的L₁、a₁、b₁、L₂₄、a₂₄、b₂₄的均值依次为44.17、15.72、8.94、46.44、17.07和9.46，变异系数分别为7.26%、15.33%、12.35%、11.02%、13.05%和12.92%。相对于45~60 min的肉色指标，24 h的肉色指标均有提高。除L₁之外，其余5个肉色性状的变异系数均在10%以上，具有较大的遗传改良潜力。

2.2　加权基因共表达网络构建

33头莱芜猪背最长肌转录组测序共产生65.63 Gb的数据，平均值为6.56 Gb。将过滤后的 Clean reads 比对到猪参考基因组（Sus scrofa 11.1）中，得到的总比对率为92.98%~96.05%，唯一比对率为89.98%~92.84%，比对率较高，证明此次测序质量较好，可用于后续分析。删除了在超过半数以上个体中低表达量（FPKM<1）的基因后，利用剩下的10 235个基因对33个样本进行聚类（图1（a）），剔除2个离群样本LW115和LW122，获得31个样本用于后续分析（图1（b））。在无尺度分布拓扑网络指数R²为0.85时筛选出最佳软阈值12，而后基于动态剪切树模型对基因进行聚类，模块最小基因数设定为35，共得到了20个共表达模块（图2），其中Dark green模块基因数目最多，为3 974个；Dark orange 2模块基因数最少，为51个；其余模块的基因数量为66~1 814。

2.3　目标模块筛选及模块内基因的功能富集分析

为筛选目标模块，计算了模块特征值ME与肉色表型的相关性（r）及P值如图3所示。基于ME与至少2个肉色性状的|r| ≥0.35且P<0.01的标准，共筛选到4个目标模块：Floral white、Dark green、Pink及Orangered 4。其中Flora lwhite模块中包含66个基因，模块特征值ME与L₁（r=0.49，P=5×10^-3）呈显著正相关，与a₁（r=-0.57，P=8.1×10^-4）和a₂₄（r=-0.53，P=2.1×10^-3）均呈显著性负相关；Dark green模块包含1 814个基因，ME与L₁（r=-0.51，P=3.1×10^-3）呈显著性负相关，与a₁（r=0.48，P=6.8×10^-3）呈显著正相关；Pink模块含有388个基因，ME与a₁（r=0.73，P=2.6×10^-6）及a₂₄（r=0.59，P=5×10^-4）均呈显著性正相关，与L₁（r=-0.51，P=3.2×10^-3）呈显著性负相关；Orangered 4模块含有72个基因，ME与L₁（r=0.52，P=2.5×10^-3）、b₁（r=0.54，P=1.9×10^-3）、L₂₄（r=0.53，P=2.0×10^-3）及b₂₄（r=0.46，P=9.5×10^-3）均呈显著性正相关。

为进一步阐明目标模块与肉色性状的关系，对4个目标模块内的基因开展了功能富集分析。肉色受到线粒体结构和功能、蛋白修饰、肌内脂肪沉积和代谢、肌肉结构和功能等因素的影响，富集分析发现4个目标模块内的基因均有富集在与此相关的GO term中（图4）。Floral white模块的基因显著富集在肌动蛋白丝组织（GO：0007015，P=8.18×10^-3）、肌动蛋白丝结合（GO：0051015，P=2.27×10^-2）、肌动蛋白结合（GO：0003779，P=3.70×10^-2）等与肌肉结构和功能相关的GO项中。Dark green模块的基因显著富集在线粒体呼吸链复合物I组装（GO：0032981，P=8.87×10^-6）、线粒体NADH到泛醌的电子传递（GO：0006120，P=1.44×10^-2）、线粒体（GO：0005739，P=1.22×10^-19）、线粒体内膜（GO：0005743，P=3.50×10^-12）、氧化还原酶活性（GO：0016684，P=3.72×10^-2）等与线粒体结构及其功能相关的GO项；也富集在血红素分解代谢过程（GO：0042167，P=4.76×10^-2）等与肌肉色素代谢相关的GO项；以及肌肉结构成分（GO：0008307，P=8.03×10^-3）、Z盘（GO：0030018，P=1.20×10^-3）等与肌肉结构和功能相关的GO项。Pink模块的基因显著富集在肌动蛋白细胞骨架（GO：0030036，P=2.93×10^-3）、肌动蛋白细胞骨架重组的调控（GO：20000249，P=4.02×10^-3）、细胞骨架（GO：0005856，P=3.47×10^-3）等与肌肉细胞结构和功能相关的GO项中。Orangered 4的基因显著富集在脂肪细胞分化的正调控（GO：0045600，P=5.20×10^-3），棕色脂肪细胞分化（GO：0050873，P=4.42×10^-3）等脂肪沉积相关的GO项中。

对4个目标模块内的基因的KEGG富集分析发现（图5），Floral white、Dark green和Pink模块中的基因显著富集在氧化磷酸化（ssc00190，P=1.53×10^-10），代谢途径（ssc01100，P=2.94×10^-5），AMPK信号通路（ssc04152，P=2.05×10^-2），细胞粘附分子（ssc04514，P=1.05×10^-2）等与肉色相关的通路。而Orangered 4模块富集到了脂肪细胞脂解的调控（ssc04923，P=4.24×10^-7）、PPAR信号通路（ssc03320，P=3.07×10^-6）等与脂肪代谢相关的通路，该结果进一步表明肉色可能受到脂质沉积的影响。以上功能富集结果表明筛选到的4个目标模块与肉色性状密切相关。

2.4　PPI网络构建及关键基因鉴定

基于|MM|>0.85且|GS|>0.3的标准初步筛选目标模块的关键基因，对模块中关键基因进行PPI网络构建（图6），筛选基因编码蛋白与模块中其他蛋白存在高互作，且基因功能与肉色密切相关的基因作为模块的关键基因（表2）。Floral white模块筛选出26个基因，PPI网络分析中包含12个蛋白（节点），筛选到1个与肉色相关的关键基因MYO9A，其编码蛋白与2个蛋白存在互作；Dark green模块中筛选出229个基因，PPI网络分析中包含162个蛋白，共筛选到的4个关键基因NDUFA13、COX5B、MRPL12和RPL8所编码的蛋白与其他16~21个蛋白存在互作；Pink模块筛选出69个基因，所构建的PPI网络共包含29个基因，筛选出2个关键基因ROBO4、EGFL7编码的蛋白分别与其他6和3个蛋白存在互作；Orangered 4模块筛选出17个基因，所构建的PPI网络节点数最少，仅包含7个蛋白，筛选出1个关键基因CIDEC，该基因所编码的蛋白与其他5个蛋白存在互作。

3　讨 论

在猪肉的生产工业中，肉色是评价肉品质的重要指标之一，也是决定消费者购买行为的直接因素，因此肉色的改良对猪肉生产者及市场营销具有重大意义 ^［13］。目前，有关猪肉质性状的研究已经有了较大的进展，根据猪QTL数据库中的数据（https：//www.animalgenome.org/cgi-bin/QTLdb/SS/index），现已发现892个与肉色相关的QTL。但由于肉色性状的复杂性，且肉色的评价中利用了多个存在密切关系的指标（L₁、a₁、b₁、L₂₄、a₂₄和b₂₄）参与评价，致使目前尚未发现影响猪肉色主效基因。WGCNA分析可用于多性状的综合分析，通过系统性地整合大量基因表达数据，揭示基因模块之间的关联性，进而识别与多个性状紧密相关的关键基因网络。因此，本研究以莱芜猪的6个肉色指标为表型，通过WGCNA分析以挖掘影响肉色性状的候选基因，为探究调控猪肉色的关键基因和功能通路提供参考，以期为肉色调控机制的深入研究奠定基础。

3.1　WGCNA在单一和多个性状基因挖掘研究中的应用探索

WGCNA 是一种系统的生物学算法，基于相同性状的基因在表达模式中具有共表达特征的假设，从构建的表达调控网络中找到并鉴定目标性状的相关基因。肉色、IMF、pH 值、滴水损失等性状均是影响肉质的关键因素，综合分析各性状、挖掘莱芜猪肉质特色基因对解析莱芜猪肉质性状的遗传机理具有非常重要的意义。本研究也尝试利用多性状进行WGCNA分析，但由于性状过多，在分析过程中发现存在一些问题：一是性状较多的情况下，存在较多性状与某一模块相关性极高，但该模块的基因在进行分析时，发现模块内的基因一般为核酸合成、代谢等基础生命活动；二是性状与模块内的基因普遍存在相关性不高，基于|MM|>0.85 且|GS|>0.3 的标准筛选基因进行下一步 PPI 分析，筛选出的基因数量极少或者筛选不到基因，且由于基因数太少在 PPI分析时其编码的蛋白不存在互作。另外，在目前有关猪肉质性状的 WGCNA 研究文献报道中，一般为较少性状的分析，例如：Chen等 ^［14］基于WGCNA对宁乡猪的背膘组织发育进行研究，共挖掘出13个与背膘发育相关的枢纽基因；Wang等 ^［15］对147头杜长大商品猪的IMF及肉色a值进行分析，发现MYC、SOX9、CEBPB和PPAGRC1A可能对IMF及肉色的a值有共同作用；其原因可能与 WGCNA 分析本身的算法相关，较少的性状更有利于挖掘出与该性状高度相关的重要调控基因，而基于较多性状的共性基因挖掘方面该算法较为欠缺。

3.2　筛选到的目标模块与肉色密切相关

通过WGCNA分析将基因划分为20个模块，从中筛选出4个与肉色相关的关键模块：Floral white、Dark green、Pink和Orangered 4。Dark green模块基因显著富集到线粒体呼吸链复合物I，IV组装、线粒体电子转运NADH对泛醌的正调控、呼吸链等多条与线粒体相关的条目，说明该模块中包含了大量与线粒体相关的基因。Guo等 ^［9］对黄羽肉鸡肉色的WGCNA关键模块富集分析同样发现了大量与线粒体相关的条目。研究发现，线粒体与肌红蛋白的相互作用与肉色存在密切联系。肌红蛋白（Myoglobin，Mb）是肉色呈色的主要物质。Mb存在的形式主要有3种：脱氧肌红蛋白（Deoxymyoglobin，DeoMb）、氧合肌红蛋白（Oxymyoglobin，OxyMb）、高铁肌红蛋白（Metmyoglobin，MetMb），不同程度的肉色取决于Mb的含量。而线粒体的氧化还原作用与肌红蛋白存在很大的相关性。肌红蛋白可以促进线粒体吸收氧的能力，是由于线粒体基于NADH通过电子传递复合物III和IV传递的电子促进MetMb还原生成OxyMb ^［16］。另外，线粒体与Mb的相互作用也影响肉色变化的稳定性 ^［17］。Wu等 ^［18］在基于线粒体蛋白组学的肉色稳定性研究中发现，有21种线粒体蛋白在肉色较差的肌肉中过表达，这些蛋白限制了线粒体的功能并降低了Mb的稳定性。同时，不同类型的肌纤维由于肌纤维中的线粒体及Mb含量不同也与肉色密切相关，如氧化型的Ⅰ型和ⅡA型肌纤维富含线粒体且Mb含量相对较多，外观上呈现红色；酵解型的IIB型肌纤维Mb及线粒体含量相对较少，颜色更浅；IIX型肌纤维的肉色则介于两者之间 ^［19］。Ryu等 ^［20］报道称，猪背最长肌肉色的L值与Ⅰ型和ⅡA型肌纤维数量呈负相关。以上结果均表明，线粒体与肌红蛋白的含量及互相作用在决定肉色品质方面具有重要地位。

另外，Dark green模块基因除富集到与线粒体相关的条目，也富集到很多与蛋白质修饰的条目，如：蛋白质磷酸化的正调控、组蛋白去乙酰化的正调控及蛋白酶体泛素依赖性蛋白分解代谢过程的正调控等条目，表明Orangered 4模块中的基因也参与蛋白质修饰的过程。蛋白质修饰为共价键作用，通过对蛋白质或氨基酸残基添加修饰基团改变其本质。目前已知的蛋白质修饰有300多种，最常见的蛋白质修饰方式为：磷酸化，乙酰化，泛素化等 ^［21-22］。Ono等 ^［23］研究发现PSE肉中的磷酸化酶α活性较高，可能导致对应激敏感的猪糖原分解速度提升。Li等 ^［24］对绵羊胸肌和腰肌进行肉色稳定性研究，结果显示蛋白质磷酸化显著影响羊肉色泽和嫩度，肉色稳定性与肌浆蛋白的磷酸化呈负相关；这说明蛋白质磷酸化可能通过调节糖酵解和肌红蛋白的氧化还原稳定性参与肉色发育。Li等 ^［25］、Zhou等 ^［26］等研究了蛋白质乙酰化对肉质性状的调控，并提出AMPK（5'单磷酸腺苷激活蛋白激酶）可以靶向调控屠宰后的肌肉代谢，减少家禽出现PSE肉的概率。以上结果均表明，筛选到的Dark green模块是影响肉色的关键模块。

Orangered 4模块显著富集到多条与脂肪细胞分化及脂质氧化代谢相关的条目，及PPAR等与脂肪调控相关的通路，此外该模块基因还富集到与PSE形成相关的AMPK通路，表明IMF含量及脂质代谢与肉色有密切的关联。肌内脂肪含量也会显著影响肉色，主要是由于脂肪在外观上呈现白色，肌内脂肪含量越高，肉色越亮 ^［27］。本课题组前期的研究结果发现莱芜猪IMF含量与屠宰后45 min及24 h的肉色均呈极显著正相关（P<0.01） ^［28］。在GWAS研究中也发现了IMF含量和肉色可定位到相同的最强关联位点，如熊信威 ^［29］对莱芜猪的肉质性状进行GWAS分析，发现位于12号染色体上的MYH4基因存在一个突变位点rs222，不仅能够显著提升IMF的含量，而且对a等肉色也有正向调控作用；Cho等 ^［30］对西方商品猪和韩国本地猪杂交猪种进行GWAS分析，也发现了同时影响IMF含量和肉色指标的最强关联位点；Li等 ^［31］对杜长大进行肉色及IMF的GWAS分析，发现IMF与肉色表型呈现极显著相关。

Floral white模块基因、Pink模块的基因和Dark green模块基因均有在肌动蛋白结合、肌动蛋白细胞骨架、钙粘蛋白介导的细胞-细胞粘附负调控、肌肉结构成分、Z盘等与肌肉结构和功能相关的GO项中显著富集。研究已发现，肉色不仅受到生化因素的影响，还取决于肌肉的物理结构和光的散射，包含肌浆蛋白的组成、肌球蛋白变性、细胞骨架的完整性以及细胞粘附蛋白的变化等 ^［32］。完整的细胞骨架蛋白能够为肌原纤维收缩提供框架并带来更多的光散射，从而影响肉色 ^［33-34］。此外，Z盘降解和肌浆蛋白的活性可能会影响光的散射进而影响肉色。

3.3　筛选的候选基因与肉色形成密切相关

本研究基于|MM|>0.85且与至少2个性状|GS|>0.3 的选择标准筛选关键模块的基因随后进行PPI网络构建，最终4个模块共筛选出8个关键基因MYO9A、NDUFA13、COX5B、MRPL12、RPL8、ROBO4、EGFL7和CIDEC。Dark green模块筛选到的关键基因NDUFA13是线粒体复合物I的一个辅助亚基，对于线粒体膜电位的维持至关重要 ^［35-36］；COX5B是编码线粒体复合物IV（细胞色素C氧化酶）的一个结构亚基，线粒体可能通过细胞色素c与Mb相互作用。Guo等 ^［9］基于WGCNA分析发现COX7C是影响黄羽肉鸡肉色的重要基因，COX7C及COX5B均为细胞色素c氧化酶的家族成员，可能是物种的原因致使研究结果存在差别。MRPL12编码线粒体核糖体蛋白12，目前已有研究表明，敲除COX5B或MRPL12均会引起线粒体功能障碍 ^［37-40］，这些基因可能通过影响线粒体的功能影响肉色。Dark green模块筛选的关键基因RPL8编码的蛋白是铁死亡标志蛋白之一，铁死亡是一种铁依赖性脂质过氧化物积累引起的细胞变化过程 ^［41］，表明肉色的调控与脂质氧化过程存在一定的相互作用。另外Orangered 4筛选到的关键基因CIDEC也与脂质代谢密切相关，CIDEC为CIDE家族蛋白的一个成员，在脂肪组织分化和脂肪细胞脂质代谢的调节中起着重要作用 ^［42］；也有学者发现CIDEC与人体脂肪组织中的线粒体及脂质代谢相关基因PLIN1存在相互影响 ^［43］。Pink模块的关键基因EGFL7属于表皮生长因子（EGF）样家族成员，Yang等 ^［44］报道称，EGFL7可以增强细胞粘附、增强细胞骨架和减少细胞间聚集。MYO9A和EGFL7基因可能是通过影响肌肉结构影响肉色。Pink的另一关键基因ROBO4在维持血管稳定，抑制血管生成和血管通透性方面起到关键作用 ^［45］。目前为止，虽未见血管生成等在肉色调控中作用的相关报道，但血管自身的颜色可能对肉色造成影响。本研究基于WGCNA初步筛选出与肉色相关的基因，但由于样本有限，这些关键基因调节肉色形成的功能还需进行细胞甚至个体层面的功能验证，从而进一步指导肉质改良工作。

4　结 论

本研究利用莱芜猪背最长肌的RNA-seq数据，通过WGCNA构建了20个共表达网络模块，筛选出4个与肉色相关的关键模块并对模块基因进行富集分析，结合PPI网络等分析方法筛选出8个与肉色调控相关的潜在关键基因MYO9A、NDUFA13、COX5B、MRPL12、RPL8、ROBO4、EGFL7和CIDEC。本研究基于WGCNA初步筛选出与肉色相关的基因，为猪肉色分子调控机制研究与肉色育种改良奠定了理论基础。

参考文献

[1] 潘浩汉. 猪肉食品中微生物的控制途径［J］. 食品安全导刊， 2019， （18）： 37

[2] Pan H H. Control of microorganisms in pork food［J］. China Food Safety Magazine， 2019， （18）： 37 （in Chinese）

[3] 刘航， 侯黎明， 王彬彬， 刘晨曦， 陶伟， 张总平， 牛培培， 李强， 李平华， 黄瑞华. 显性效应对苏淮猪肉色性状遗传评估和基因组选择的影响［J］. 畜牧与兽医， 2021， 53（6）： 1-6

[4] Liu H， Hou L M， Wang B B， Liu C X， Tao W， Zhang Z P， Niu P P， Li Q， Li P H， Huang R H. The influence of dominance effects on the estimation of meat color genetic parameters and genomic selection in Suhuai pigs［J］. Animal Husbandry & Veterinary Medicine， 2021， 53（6）： 1-6 （in Chinese）

[5] 李国辉， 谢恺舟， 戴国俊， 张跟喜， 张涛， 杨国明， 陈兰， 韩威， 苏一军， 王金玉. 金茅黑鸡肉品质性状的关联转录组分析［J］. 中国畜牧杂志， 2019， 55（12）： 85-90

[6] Li G H， Xie K Z， Dai G J， Zhang G X， Zhang T， Yang G M， Chen L， Han W， Su Y J， Wang J Y. Correlation analysis on the meat quality traits and transcriptome in Jinmao Black chicken［J］. Chinese Journal of Animal Science， 2019， 55（12）： 85-90 （in Chinese）

[7] Wang H， Shen Z H， Huang R H， Zhao A Y， Jiang J N， Li P H， Zhou X L， Yang S B， Hou L M. A polymorphism in porcine miR-22 is associated with pork color［J］. Frontiers in Veterinary Science， 2022， 9： 939440

[8] 廖印长， 高虎， 张跃博， 殷诗舒， 徐康， 何俊. 宁乡猪宰后pH和肉色性状的全基因组关联分析［J］. 中国畜牧杂志， 2021， 57（S1）： 174-181

[9] Liao Y C， Gao H， Zhang Y B， Yin S S， Xu K， He J. Genome wide association analysis of pH and meat color traits after slaughter in Ningxiang pigs［J］. Chinese Journal of Animal Science， 2021， 57（S1）： 174-181 （in Chinese）

[10] Langfelder P， Horvath S. WGCNA： An R package for weighted correlation network analysis［J］. BMC Bioinformatics， 2008， 9： 559

[11] 荣丽云， 黄宁， 王孝义， 李明丽， 陈强， 鲁绍雄. 基于WGCNA方法识别影响猪肉质性状的关键模块和基因［J］. 云南农业大学学报：自然科学版， 2022， 37（5）： 769-779

[12] Rong L Y， Huang N， Wang X Y， Li M L， Chen Q， Lu S X. Identification of key gene co-expression modules and genes associated with pork quality based on WGCNA Method［J］. Journal of Yunnan Agricultural University：Natural Science， 2022， 37（5）： 769-779 （in Chinese）

[13] Zhao X， Wang C， Wang Y， Zhou L， Hu H， Bai L， Wang J. Weighted gene co-expression network analysis reveals potential candidate genes affecting drip loss in pork［J］. Animal Genetics， 2020， 51（6）： 855-865

[14] Guo X， Zhang H， Wang H， He X X， Wang J X， Wei W， Liu M， Xu J M， Liu Y N， Jiang R S. Identification of key modules and hub genes involved in regulating the color of chicken breast meat using WGCNA［J］. Animals： an Open Access Journal from MDPI， 2023， 13（14）： 2356

[15] Trapnell C， Pachter L， Salzberg S L. TopHat： Discovering splice junctions with RNA-Seq ［J］. Bioinformatics， 2009， 25（9）： 1105-1111

[16] 王文翔， 杜丽丽， 胡俊伟， 张琰翔， 马铭昊， 段蕊， 钱聪， 王欣玥， 李三禄， 张长庆， 张路培， 高雪， 徐凌洋， 李俊雅， 高会江. 基于转录组分析挖掘平凉红牛背最长肌肉质性状候选基因［J］. 畜牧兽医学报， 2023， 54（11）： 4589-4604

[17] Wang W X， Du L L， Hu J W， Zhang Y X， Ma M H， Duan R， Qian C， Wang X Y， Li S L， Zhang C Q， Zhang L P， Gao X， Xu L Y， Li J Y， Gao H J. Mining candidate genes related to meat quality traits of longissimus lumborum in Pingliang red steer based on transcriptome［J］. Acta Veterinaria et Zootechnica Sinica， 2023， 54（11）： 4589-4604 （in Chinese）

[18] Huang D W， Sherman B T， Lempicki R A. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources［J］. Nature Protocols， 2009， 4（1）： 44-57

[19] Altmann B A， Trinks A， Mörlein D. Consumer preferences for the color of unprocessed animal foods［J］. Journal of Food Science， 2023， 88（3）： 909-925

[20] Chen C， Ren H B， Li H L， Deng Y， Cui Q M， Zhu J， Zhang S Y， Yu J E， Wang H M， Yu X D， Yang S L， Hu X G， Peng Y L. Identification of crucial modules and genes associated with backfat tissue development by WGCNA in Ningxiang pigs［J］. Frontiers in Genetics， 2023， 14： 1234757

[21] Wang B B， Hou L M， Yang W， Men X M， Qi K K， Xu Z W， Wu W J. Construction of a co-expression network affecting intramuscular fat content and meat color redness based on transcriptome analysis［J］. Frontiers in Genetics， 2024， 15： 1351429

[22] Tang J L， Faustman C， Mancini R A， Seyfert M， Hunt M C. Mitochondrial reduction of metmyoglobin： Dependence on the electron transport chain［J］. Journal of Agricultural and Food Chemistry， 2005， 53（13）： 5449-5455

[23] Yamada T， Takakura H， Jue T， Hashimoto T， Ishizawa R， Furuichi Y， Kato Y， Iwanaka N， Masuda K. Myoglobin and the regulation of mitochondrial respiratory chain complex IV［J］. The Journal of Physiology， 2016， 594（2）： 483-495

[24] Wu S， Luo X， Yang X Y， Hopkins D L， Mao Y W， Zhang Y M. Understanding the development of color and color stability of dark cutting beef based on mitochondrial proteomics［J］. Meat Science， 2020， 163： 108046

[25] Picard B， Gagaoua M. Muscle fiber properties in cattle and their relationships with meat qualities： An overview［J］. Journal of Agricultural and Food Chemistry， 2020， 68（22）： 6021-6039

[26] Ryu Y C， Kim B C. The relationship between muscle fiber characteristics， postmortem metabolic rate， and meat quality of pig longissimus dorsi muscle ［J］. Meat Science， 2005， 71（2）： 351-357

[27] Khoury G A， Baliban R C， Floudas C A. Proteome-wide post-translational modification statistics： Frequency analysis and curation of the swiss-prot database ［J］. Scientific Reports， 2011， 1

[28] Jensen O N. Modification-specific proteomics： Characterization of post-translational modifications by mass spectrometry ［J］. Current Opinion in Chemical Biology， 2004， 8（1）： 33-41

[29] Ono K， Topel D G， Christian L L， Altheni T G. Relationship of cyclic-amp and phosphorylase a in stress-susceptible and control pigs ［J］. Journal of Food Science， 1977， 42（1）： 108-110

[30] Li M， Li X， Xin J Z， Li Z， Li G X， Zhang Y， Du M T， Shen Q W， Zhang D Q. Effects of protein phosphorylation on color stability of ground meat［J］. Food Chemistry， 2017， 219： 304-310

[31] Li Z W， Li X， Wang Z Y， Shen Q W， Zhang D Q. Antemortem stress regulates protein acetylation and glycolysis in postmortem muscle［J］. Food Chemistry， 2016， 202： 94-98

[32] Zhou B， Shen Z L， Liu Y S， Wang C T， Shen Q W. Proteomic analysis reveals that lysine acetylation mediates the effect of antemortem stress on postmortem meat quality development［J］. Food Chemistry， 2019， 293： 396-407

[33] Zhang X Y， Liu C Y， Kong Y Y， Li F D， Yue X P. Effects of intramuscular fat on meat quality and its regulation mechanism in Tan sheep［J］. Frontiers in Nutrition， 2022， 9： 908355

[34] 任一帆， 赵雪艳， 王彦平， 李菁璇， 武英， 张传生， 王继英. 莱芜猪与杜长大肉质性状比较与分析［J］. 养猪， 2021， （5）： 45-48

[35] Ren Y F， Zhao X Y， Wang Y P， Li J X， Wu Y， Zhang C S， Wang J Y. Comparison of carcass traits and meat quality between Laiwu pig and Duroc × Landrace × Yorkshire pigs ［J］. Swine Production， 2021， （5）： 45-48 （in Chinese）

[36] 熊信威. 猪12号染色体上增加肌内脂肪含量的因果基因及其因果突变的鉴别与作用机理研究［D］. 南昌：江西农业大学， 2017

[37] Xiong X W. Identification of the causative gene and mutation for increasing IMF content at SSC12 and its mechanism study［D］. Nanchang：Jiangxi Agricultural University， 2017 （in Chinese）

[38] Cho I C， Park H B， Ahn J S， Han S H， Lee J B， Lim H T， Yoo C K， Jung E J， Kim D H， Sun W S， Ramayo-Caldas Y， Kim S G， Kang Y J， Kim Y K， Shin H S， Seong P N， Hwang I S， Park B Y， Hwang S， Lee S S， Ryu Y C， Lee J H， Ko M S， Lee K， Andersson G， Ps affEnciso M， Lee J W. A functional regulatory variant of MYH3 influences muscle fiber-type composition and intramuscular fat content in pigs［J］. PLoS Genetics， 2019， 15（10）： e1008279

[39] Li H， Xu C N， Meng F M， Yao Z K， Fan Z F， Yang Y S， Meng X L， Zhan Y X， Sun Y， Ma F C， Yang J F， Yang M， Yang J， Wu Z F， Cai G， Zheng E Q Y. Genome-wide association studies for flesh color and intramuscular fat in （Duroc ciation studies for flesh crossbred commercial pigs ［J］. Genes， 2022， 13（11）： 2131

[40] Hughes J M， Clarke F M， Purslow P P， Warner R D. Meat color is determined not only by chromatic heme pigments but also by the physical structure and achromatic light scattering properties of the muscle［J］. Comprehensive Reviews in Food Science and Food Safety， 2020， 19（1）： 44-63

[41] Bee G， Anderson A L， Lonergan S M， Huff-Lonergan E. Rate and extent of pH decline affect proteolysis of cytoskeletal proteins and water-holding capacity in pork［J］. Meat Science， 2007， 76（2）： 359-365

[42] Lomiwes D， Farouk M M， Wu G， Young O A. The development of meat tenderness is likely to be compartmentalised by ultimate pH［J］. Meat Science， 2014， 96（1）： 646-651

[43] Fearnley I M， Carroll J， Shannon R J， Runswick M J， Walker J E， Hirst J. GRIM-19， a cell death regulatory gene product， is a subunit of bovine mitochondrial NADH： Ubiquinone oxidoreductase （complex I）［J］. The Journal of Biological Chemistry， 2001， 276（42）： 38345-38348

[44] Lu H， Cao X M. GRIM-19 is essential for maintenance of mitochondrial membrane potential［J］. Molecular Biology of the Cell， 2008， 19（5）： 1893-1902

[45] Gao S P， Sun H F， Jiang H L， Li L D， Hu X， Xu X E， Jin W. Loss of COX5B inhibits proliferation and promotes senescence via mitochondrial dysfunction in breast cancer［J］. Oncotarget， 2015， 6（41）： 43363-43374

[46] Tim 3741 ómez A， Bartley-Dier E L， Fontanesi F， Barrientos A. HIGD-driven regulation of cytochrome c oxidase biogenesis and function ［J］. Cells， 2020， 9（12） ： 2620

[47] Chen Y， Cairns R， Papandreou I， Koong A， Denko N C. Oxygen consumption can regulate the growth of tumors， a new perspective on the Warburg effect［J］. PLoS One， 2009， 4（9）： e7033

[48] Postnikova G B， Shekhovtsova E A. The effect of mitochondrial and artificial bilayer phospholipid membranes on conformation of myoglobin and its affinity for oxygen［J］.American Journal of Bilological chemistry， 2015，3（1）：16-32

[49] Dixon S J， Lemberg K M， Lamprecht M R， Skouta R， Zaitsev E M， Gleason C E， Patel D N， Bauer A J， Cantley A M， Yang W S， Morrison B， 3rd， Stockwell B R. Ferroptosis： An iron-dependent form of nonapoptotic cell death ［J］. Cell， 2012， 149（5）： 1060-1072

[50] Li F F， Gu Y， Dong W P， Li H， Zhang L Y， Li N L， Li W Z， Zhang L J， Song Y， Jiang L N， Ye J， Li Q. Cell death-inducing DFF45-like effector， a lipid droplet-associated protein， might be involved in the differentiation of human adipocytes ［J］. The FEBS Journal， 2010， 277（20）： 4173-4183

[51] Moreno-Navarrete J M， Ortega F， Serrano M， Rodriguez-Hermosa J I， Ricart W， Mingrone G， Fernández-Real J M. CIDEC/FSP27 and PLIN1 gene expression run in parallel to mitochondrial genes in human adipose tissue， both increasing after weight loss［J］. International Journal of Obesity， 2014， 38（6）： 865-872

[52] Yang R H， Qi S H， Ruan S B， Lin Z P， Lin Y， Zhang F G， Chen X D， Xie J L. EGFL7-overexpressing epidermal stem cells promotes fibroblast proliferation and migration via mediating cell adhesion and strengthening cytoskeleton［J］. Molecular and Cellular Biochemistry， 2016， 423（1/2）： 1-8

[53] Zhang F， Prahst C， Mathivet T， Pibouin-Fragner L， Zhang J S， Genet G， Tong R， Dubrac A， Eichmann A. The Robo4 cytoplasmic domain is dispensable for vascular permeability and neovascularization［J］. Nature Communications， 2016， 7： 13517

基金资助

山东省重点研发计划(2022LZGC003)

山东省重点研发计划(2022LZGCQY007)

山东省生猪产业技术体系(SDAIT-08-03)

山东省农业科学院农业科技创新工程(CXGC2023A23)