结合WGCNA鉴定与猪肌纤维和肌内脂肪相关的中枢基因

刘娟，王舒，左周，路畅，杨阳，蔡春波，赵燕，郭晓红，曹果清，李步高，高鹏飞

（山西农业大学 动物科学学院，山西 太谷030801）

猪肉是人类蛋白质和脂肪的主要来源之一，约占世界肉类消费的30%。影响猪肉品质的因素有遗传、营养、应激等，评定肉质的感官指标有肉色、pH值、嫩度、系水力、风味等［1］。肌纤维组织学特性包括肌纤维直径、肌纤维密度、肌纤维类型、肌节长度、结缔组织和肌内脂肪含量，是评价肉品质的重要经济性状。肌纤维直径（Diameter of muscle fiber，DiaMF）与肌肉嫩度相关，直径越小，肉质越嫩［2，3］。肌纤维密度（Density of muscle fiber，DenMF）和肌内脂肪（Intramuscular fat，IMF）含量影响猪肉风味、嫩度和多汁性［4］。研究表明，肌内脂肪含量低于2.5%的肉类会影响适口性，较高的肌内脂肪含量通常被认为对肉质特征具有积极的影响［5，6］。马身猪是一个中国北方猪品种，先前研究分析了马身猪和大白猪之间与生长速度和肉质相关的表型差异。马身猪比大白猪具有较低的生长速度和更好的肉质，马身猪的肌内脂肪含量高于大白猪［7］。进行马身猪和大白猪的转录组分析，阐明2个猪种不同生长发育阶段的分子机制，发现马身猪和大白猪的主要区别在于氨基酸代谢，脂肪酸代谢和骨骼肌发育［8］。

高通量测序技术能够准确获得生物体或组织中特定时期的基因表达特征，近年来，国内外研究者利用RNA-Seq技术从转录组水平探析了不同品种肌肉发育［9］、肌内脂肪沉积［10］的基因表达模式，构建了相关基因调控网络，发现并鉴定了PPARGC1B、TRIM 63、CSRP3、ACTN2、MYL 1、MYH 6等一系列基因是影响肌肉生长和脂肪沉积的关键基因［11］，胰岛素、丝裂原活化蛋白激酶（MAPK）、Wnt等信号通路能够显著影响肌纤维类型和肌肉发育［12］。WGCNA是一种应用于基因组学、蛋白质组学和代谢组学中描述相应基因、蛋白或代谢物表达相关性的系统生物学分析方法［13］。Matthew等［14］首次将WGCNA应用于代谢组学数据，基于WGCNA分析番茄代谢组学数据，确定了与番茄成熟和遗传背景相关性状的3个主要代谢模块。Priscila等［15］利用WGCNA揭示了miRNA在内洛尔牛脂肪酸组成中的潜在调控作用；宋万勤等［16］使用从欧洲生物信息研究所获得的RNA-Seq数据鉴定与间充质干细胞成骨分化有关的长非编码RNA，发现32种新颖的LncRNA是 包 括miR-689，miR-640，miR-601和miR-544等在内的miRNA前体。

随机森林（Random forest，RF）［17］是一种机器学习领域的集成方法，它利用随机重采样技术和节点随机分裂技术构建多棵决策树，通过投票得到最终分类结果。程淑萍等［19］基于ncRNA和蛋白质的序列信息提取特征，训练包括随机森林算法在内的3个机器学习模型，预测ncRNA与蛋白质的相互作用，预测的准确率较高，证明机器学习方法能够预测ncRNA与蛋白质是否存在相互作用。

课题组前期以马身猪和大白猪为研究对象，对不同生长阶段肌肉组织进行转录组测序，构建了肌肉组织中19条长非编码RNA（Long noncoding RNA，LncRNA）539条mRNA组成的分子调控网络［20］。本研究采用WGCNA分析技术并结合随机森林算法，对2个品种肌肉转录组测序数据深入分析，以期获得影响猪肌肉发育及肌内脂肪沉积的关键基因。

1 材料和方法

1.1 实验动物与样品采集

所有动物手术都严格按照《世界医学协会道德守则》（http：//ec.europa.eu/environment/chemicals/lab_animals/legislation_en.htm）进行，实验方案经山西农业大学动物伦理委员会批准。本研究选用马身猪和大白猪为研究对象，饲养于山西省大同市种猪场，按照猪饲养标准（NY/T 65-2004）进行饲喂。分别在1、30、60、90、120、150、180 d随机选择3头体重相近的去势公猪屠宰，采集位于最后肋骨处的背最长肌样品，液氮保存备用。按照背最长肌纤维方向垂直切割，采集1 cm3肌肉组织放入4%多聚甲醛固定。

1.2 肌纤维的检测直径、密度和肌内脂肪含量的测定

肌肉组织在4%甲醛中固定，经乙醇脱水、二甲苯透明后用石蜡包埋，6μm切片后采用苏木精—伊红（H&E）染色。每个采样点平行样品不少于9张切片，每张切片取3～5张图像，采用马尾松三色染色测定肌纤维的直径和密度。在10×10倍显微镜下，对每个样品作20个视野大方格内肌纤维根数的测定（大方格面积约为0.25 mm2），换算后得到每平方毫米纤维根数，得出肌纤维密度。在10×40倍显微镜下，采用不规则多边形对每个样本进行100根肌纤维横断面积的描画，统计出所有的横截面积，取其平均值，得到肌纤维面积。根据肌纤维面积计算，假设肌纤维呈圆柱形，由圆面积公式求得肌纤维直径［21］。采用索氏脂肪抽提法分析肌内脂肪含量［22］。

1.3 测序数据来源与分析方法

数据来源于本课题组前期猪肌肉转录组数据（NCBI登录号为SRP068558（https：//trace.ncbi.nlm.nih.gov/Traces/sra_sub/sub.cgi？login=pda）），共18个样本，包含马身猪和大白猪2个品种3个日龄（1、90、180 d）的测序数据。使用R包limma［23］进行差异基因分析，DEG的阈值设置为|log2（fold-change）|＞0.8且P＜0.05。

1.4 WGCNA分析

将18个样本随机分为训练集和测试集，其中训练集含10个样本，测试集含有8个样本。训练集数据采用R-Studio 4.0.0软件的WGCNA包［13］进行分析。首先计算所有基因的表达相关系数，根据基因间相关性确定软阈值（power），构建基因聚类树。根据基因表达相关系数确定表达模块，每个表达模块最少基因数设定为30个。将每种基因表达模块作为整体对象，计算模块间的相关系数，当模块间相关系数高于0.85，则合并相关模块。

1.5 GO和KEGG富集分析

基因功能注释分析采用DAVID［24］（http：//david.abcc.ncifcrf.gov/）在线数据库进行。提取GO分析和KEGG富集分析的结果，分别基于阈值P＜0.05和P＜0.01进行分析。

1.6 模块核心基因的筛选及随机森林法验证

核心基因由模块内基因连通性决定。根据基因的表达趋势计算皮尔逊相关系数，取相关系数绝对值大于0.2的基因进行共表达分析。利用cytohubba插件［25］筛选每个模块的核心基因，使用Cytoscape（version 3.5.1）［26］对 模 块 构 建 互 作 网络。利用R语言Random Forest模块［17］进行随机森林分析，计算模块核心基因的平均基尼指数减少值（Mean Decrease in the Gini index，MDG），筛选hub基因。采用Graphpad 8.0软件可视化作图。

1.7 qPCR检测hub基因发育性表达变化

采用qPCR技术检测hub基因在不同月龄2个品种背最长肌的发育性表达变化。在NCBI数据库中获得筛选得到的hub基因序列，通过Primer3plus在线软件进行引物设计，选取18S作为检测内参基因，引物信息详见表1。采用Prime ScriptTM RT reagent Kit with gDNA Eraser试剂盒（TaKaRa，大连，中国）合成cDNA，反转录时总RNA浓度为500 ng·μL-1。使用SYBR Green RTPCR（Taκara，Dalian，China）进行qPCR检测。每个样本进行3次技术重复，分别以无模板和无逆转录酶的体系作为阴性对照。反应程序如下：预变性阶段，95℃30 s；循环阶段，95℃30 s，60℃30 s，设定循环数为45；溶解阶段，95℃30 s，60℃1 min，95℃30 s。

表1 关键基因的引物列表Table 1 The list of quantitative primers of hub genes

1.8 统计分析

qPCR结果采用2-△△Ct法计算hub基因的相对表达量。运用SPSS 22.0软件进行基因表达差异分析，显著水平选择0.05和0.01。采用Graphpad 8.0软件进行基因表达量作图。

2 结果与分析

2.1 WGCNA分析结果

利用R语言WGCNA包计算训练集样本中基因表达相关系数并聚类，表达趋势性相似的基因合并为一组，构建基因共表达网络（图1）。训练集样本聚为两大类，马身猪和大白猪1日龄样本聚为第一大类，90日龄和180日龄样本聚为第二大类，每一大类下各品种分别聚一簇。样本聚类树下方为猪背最长肌肌纤维直径、密度和肌内脂肪含量的检测热图，可以看出2个品种猪背最长肌肌纤维直径随着日龄均显著增加，肌纤维密度逐渐降低，肌内脂肪含量增加。同一日龄2个品种间相比，马身猪的肌纤维密度显著大于大白猪，肌内脂肪含量也显著高于大白猪。

图1 样本聚类图和性状热图Fig.1 Cluster diagram of samplesand heatdiagram of traits

采用平均链接层次聚类将表达模式相似的差异基因进行分组聚类，建立第一层模块聚类树状图（图2A Dynamic Tree Cut），计算模块间的相关系数，建立模块特征值聚类树状图（图2B），当模块相关系数大于0.85，软阈值大于5时，对模块进行二次合并，建立第二层模块聚类树状图（图2A Merged dynamic）。本研究中7262条差异基因首先获得19个表达模块，根据模块相关性优化后最终获得16个模块基因集，之后的分析按照合并后的模块基因集进行。16个模块平均含有450条基因，蓝色（blue）模块中含有的基因数最多（2960条），灰色60（grey60）和浅绿色（lightgreen）模块基因最少，仅含50条。

图2 差异表达基因基于相关系数聚类树状图Fig.2 Dendrogram of all differentially expressed genes clustered based on correlation coefficient

2.2 表型相关模块筛选

通过计算模块特征与表型间的相关系数，获得模块性状相关性热图（图3）。肌内脂肪与MEblack模块成极显著正相关（R=0.78，P=0.008），与MEsalmon模块（R=-0.73，P=0.02）呈显著负相关。肌纤维密度与MEblack模块（R=-0.78，P=0.007）和MElightgreen模块（R=-0.66，P=0.04）存在显著负相关。只有MEmagenta模块与肌纤维直径存在极显著正相关（R=0.9，P=0.0004）。

图3 模块特征基因与肌肉性状的相关性热图Fig.3 Heat map of correlation between module characteristic genes and muscle traits

2.3 GO和KEGG富集结果

对MEblack、MEsalmon、MElightgreen和MEmagenta模块基因集分别进行功能富集（图4）。GO分析结果表明，MEblack模块基因功能富集于生物合成、分子修饰及细胞内各种代谢进程，MElightgreen模块基因功能富集于核酸代谢、蛋白分子复合物组装及细胞内代谢进程，MEsalmon模块基因功能富集于细胞质、细胞器、细胞骨架、转录因子结合及一些细胞代谢进程，MEmagenta模块基因显著富集于细胞运动、细胞形态发生等进程。

图4 GO和KEGG富集图Fig.4 GO and KEGG enrichment plots

MEblack和MElightgreen模块与肌纤维密度存在极显著的负相关，这2个模块共同富集于代谢和大分子修饰等进程，主要包含总代谢进程（GO：0008152），细胞代谢（GO：0044237）、初级 代谢（GO：0044238）等进程。MEblack模块与肌内脂肪含量存在极显著的正相关，其中碳水化合物代谢（GO：0005975）、小分子代谢（GO：0044281）、氧化还原活性（GO：0016491）等生物进程可能是显著影响肌内脂肪沉积的关键分子进程。MEmagenta模块与肌纤维直径存在显著的正相关，主要集中在 细 胞 骨 架（GO：0007010）、细 胞 质（GO：0005829）和细胞器（GO：0006996）等生物进程。有趣的是MEsalmon模块没有与其它3个模块共有的生物进程，该模块与肌内脂肪含量存在极显著负相关，主要集中于激酶活性（GO：0016301）、离子结合（GO：0043167）、细胞及胞内组成结构形态发生（GO：0000902）等生物进程，说明这些进程对肌内脂肪沉积存在负调控作用。

对于KEGG通路富集分析发现，MEblack模块 显 著 富 集 于cGMP-PKG（ko04022）、MAPK（ko04010）、基 质 粘 附（ko04510）和 代 谢 途 径（ko01100）等信号通路。lightgreen模块基因显著富集于配体识别（ko04371）、C型凝集素受体（ko04625）、cAMP（ko04024）以 及 生 物 合 成（ko04722）等信号通路。MEmagenta模块基因极显 著 富 集 于cGMP-PKG（ko04022）、Wnt（ko04310）、Rap1（ko04015）、cAMP（ko04024）、配体识别（ko04371）、甲状腺激素（ko04919）、C型凝集素受体（ko04625）、液体剪切应力和动脉粥样硬化（ko05418）等信号通路。salmon模块基因显著富集于液体剪切应力和动脉粥样硬化（ko05418）、cAMP（ko04024）、Wnt（ko04310）等信号通路。

2.4 共表达子网络的构建及hub基因的筛选

根据基因聚类结果，MEblack、MElightgreen、MEsalmon和MEmagenta4个模块分别包含156条、50条、83条和113条差异基因，利用Cytoscape分别构建4个模块基因互作网络。图5A为MEblack模块共表达网络，156个基因共形成11800条连接，根据cytohubba程序获得各基因间的连通度，最 终 确 定CENPI、MAP3K 3、XRCC6、FAM 131B、RIC8B、PKP2和PLEKHA 5等基因为该模块的hub基因。图5B、图5C和图5D分别是MElightgreen、MEsalmon和MEmagenta模块基因形成的表达网络。与MEblack模块类似，基于基因连通度结果，最终确定SNRNP48和HK3为MElightgreen模块的hub基因，FAM 76B、TNXB、AGTRAP和CSDC2为MEsalmon模块hub基因，AKT 3、CDK 16、RAB21、CBX7、ZNF792和PARP6为MEmagenta模块hub基因。

图5 模块基因的共表达网络Fig.5 The coexpression network of module genes

2.5 随机森林对核心基因的检验结果

根据WGCNA分析结果，共获得19个核心基因。将这些核心基因作为样本特征属性，建立随机森林分类模型对训练集样本进行深度学习，通过测试集对模型进行检验，从而确定显著影响肌肉肌纤维直径、肌纤维密度和肌内脂肪沉积3个性状的关键基因。检验结果见表2，训练集中共有10个样本，只有1 d样本出现预测错误，90 d和180 d的样本全部预测准确，预测准确率为90%；测试集预测结果与训练集相似，1 d样品出现预测错误，其余样本均准确预测，预测准确率为87.5%。

表2 训练集和测试集的分类结果Table 2 Classification results of training set and test set

本研究中，随机森林的预测准确率达到87.5%。基因对样本分类预测的贡献率，可根据基因的MDG筛选（图6）。CBX7（MDG=1.58）、PLEKHA 5（MDG=0.83）、RIC8B（MDG=0.53）、FAM 131B（MDG=0.52）和PARP6（MDG=0.29）5个基因对样本分类贡献最大，最终确定这5个基因为影响肌纤维发育和肌内脂肪沉积的hub基因。

图6 基因MDG直方图Fig.6 Histogram plot of MDG of Gene

2.6 hub基因在猪肌肉组织中的时序性表达变化

为了验证hub基因在猪生长发育过程中的变化趋势，以马身猪和大白猪为对象，选择1、30、60、90、120、150、180 d共7个时间点进行采样，采用qPCR检测hub基表达（图7）。大白猪中，CBX7、PARP6、FAM 131B和RIC8B四个基因均随着日龄表达量逐渐升高。与大白猪相比，CBX7和PARP6在马身猪中基因表达趋势一致，随着日龄表达量逐渐升高；FAM 131B和RIC8B基因出现先升高后降低的趋势，在90 d时表达量最高。PLEKHA 5基因表达趋势较特殊，在2个品种中均为先降低后升高的变化，90 d时表达量降到最低。

图7 2个猪品种中5个关键基因的时序性表达Fig.7 Time-sequence expression line diagram of 5 hub genes in 2 pig breeds

3 讨论

肌肉中肌纤维的直径、密度和肌内脂肪含量是影响猪肉风味、嫩度、多汁性等肉质性状的关键因素［27］。我国地方品种和大白猪、长白猪等国外引进猪品种相比，肌纤维直径小，肌内脂肪含量高，系水力强，鲜嫩多汁，肉味香浓［27，28］。近年来随着高通量测序技术的快速发展，国内外学者在转录组水平对影响肌肉发育和脂肪沉积的关键基因和分子调控机制进行了大量研究，获得了一系列关键基因和信号通路。

Sollero等［29］对皮奥伊州猪和长大二元猪胎儿肌肉组织转录本进行分析，发现了能够显著影响胎儿期肌肉发育的1300条差异表达的mRNA。Hu等［30］对沙子岭猪和大白猪背最长肌进行转录组和蛋白组分析，获得了488条显著差异的基因转录本和23个差异蛋白，这些基因对猪肌肉脂肪酸代谢、糖酵解和肌肉发育等具有显著的影响。最近课题组利用基因表达模式分析揭示了影响大白猪发育的基因（GHSR，EZR，FOXS1，DRD2，SH 3PXD2B，CSF1和TSHR），并 利 用GO和KEGG寻找了与脂肪酸合成相关的差异表达基因（ACACA，ACSF3，OXSM，CBR4和HSD17B8），发现肌内脂肪在猪发育的90～180 d沉积［31］。本研究将WGCNA网络分析和随机森林深度学习相结合，从大数据多维角度探寻影响肌肉性状的关键基因，最终筛选到5条hub基因和一些关键调控通路。李鑫等［32］基于网络分析和随机森林方法对肝癌分期进行研究，挑选出BUB1B、CCNA 2、CCNB1、CCNB2、CDC20、MAD2L1、MCM 4、PC-NA、RFC4和TOP2A等10条对肝细胞癌有重要影响的关键基因。有研究者利用WGCNA挖掘牛皮癣的关键控制因素和驱动因素，使用候选基因构建基于随机森林的二元分类器，成功将验证组中的正常样本与牛皮癣样本区分开，准确度为0.95～1，筛选获得的hub基因作为潜在的牛皮癣病发进程的的治疗靶标［33］。

本研究发现CBX7和PARP6基因归类于MEmagenta模块中，主要作用细胞运动、细胞骨架形成、细胞形态发生等功能，对肌纤维直径有促进作用。研究发现CBX7通过基因重组蛋白（Dickkopf-1，DKK-1）介导的Wnt/β-catenin途径能够抑制乳腺肿瘤［34］，胶质瘤［35］等肿瘤细胞的增殖。结合网络分析，CBX7在MEmagenta模块中，MEmagenta模块基因和脂肪沉积性状高度相关，随机森林分析发现，CBX7对样本分类贡献最大。本研究中发现CBX7基因对肌内脂肪沉积没有产生显著影响，但有研究发现，在CBX7基因敲除小鼠的脂肪组织含量显著高于野生型小鼠［36］，说明了CBX7对脂肪沉积也具有重要作用，具体作用机制有待深入研究。因此可以推断，CBX7基因是影响脂肪沉积性状和肌肉发育的关键基因。研究表明PARP6是MAPK的信号靶点，影响肌肉发育和脂肪沉积［37］，对大白猪皮下和肌内脂肪组织进行测序发现MAPK信号通路中PARP6显著减少，说明该基因在脂肪细胞分化过程中具有重要的调节作用［22］。PARP6能够与胰岛素反应性氨基肽酶（insulin-responsive amino peptidase，IRAP）和GRB14结合，调节脂肪细胞中葡萄糖转运体4（glucose transporter 4，GLUT 4）的活性［37］，猜测PARP6可能是影响肌肉发育和肌内脂肪沉积的关键基因。

在本研究中，FAM 131B、RIC8B、PLEKHA 5均属于MEblack模块，能够促进肌内脂肪沉积并同时抑制肌纤维密度增加。FAM 131B基因参与MAPK信号通路，促进细胞增殖，机体免疫的形成。Huriye Cin等［38］发现FAM 131B与BRAF结合，激活MAPK信号通路，导致肿瘤发生。RIC8B基因编码GEF，属于G蛋白介导的信号传导途径，在细胞分裂过程中发挥重要作用［39］。研究发现在成骨细胞增殖过程中C/EBPβ调控RIC8B基因的高表达，而在成骨细胞分化过程中RIC8B基因表达量显著下降［40］。PLEKHA 5在在脂质和脂蛋白代谢通路中起到重要作用，对肿瘤细胞的增殖具有重要作用也有相关报道［41］。

4 结论

本研究通过WGCNA分析获得了与肌纤维发育和肌内脂肪沉积的基因集，构建了重要模块的共表达网络，利用随机森林算法对与性状相关的功能基因进行深度学习与验证，最终获得5条hub基因，包括CBX7、PARP6、PLEKHA 5、FAM 131B和RIC8B。这些基因对肌肉发育和脂肪代谢存在显著作用，但具体的作用机制还需进一步研究。