牛快慢肌差异表达基因的筛选和分析

Screening and analysis of differentially expressed genes in fast and slow twitch muscles of cattle

摘要

为挖掘调控肌肉生长及骨骼肌纤维转化的关键基因，本研究利用差异分析、基因共表达网络分析和外显子选择信号分析的研究方法，通过对肌纤维类型差异明显的腓肠肌（G， n=18）和比目鱼肌（S， n=21）进行转录组测序分析。结果表明：1）共有1 576个差异表达基因，其中腓肠肌与比目鱼肌相比有170个基因上调，1 406个基因下调，富集分析发现主要富集在肌肉收缩、骨骼肌收缩、糖酵解过程、肌原纤维、肌钙蛋白复合体、AMPK信号通路和脂肪细胞脂解调节等通路中。2）基因共表达网络（WGCNA）分析显示，共被聚类为4个模块其中品红色、黑色、蓝色模块为显著模块。3）外显子选择信号分析（EWSS）发现26 467个高质量的SNP，最后联合多种方法共筛选到6个候选基因CACNA2D2、PFKFB4、ATP2B2、PRUNE2、PRKG1和FAM13A。4）对候选基因进行qPCR检测，结果与RNA-Seq测序结果一致，将候选基因输入到STRING数据库以及OmicsNet中进行蛋白质及代谢物互作网络分析，发现与候选基因互作的蛋白质之间共同作用促成了一个共享功能。综上，本研究最终筛选到6个相关差异基因，其可能是调控肌肉生长脂质代谢以及骨骼肌纤维转化的关键基因，为骨骼肌纤维转化的调控机制研究提供了基础，并为探索改善牛肉品质的新靶点和潜在机制提供了新的研究视角。

Abstract

The study was to explore the key genes regulating muscle growth and skeletal muscle fiber transformation.By differential expression analysis， gene co-expression network analysis （WGCNA）， and exon-based selection signal analysis， transcriptome sequencing was performed on the gastrocnemius （G， n=18） and soleus （S， n=21） muscles with significant differences in muscle fiber types. The results showed that： 1） There were 1 576 differentially expressed genes， with 170 genes upregulated and 1 406 genes downregulated in gastrocnemius compared to soleus muscle. Enrichment analysis revealed that these genes were mainly enriched in pathways such as muscle contraction， skeletal muscle contraction， glycolysis process， myofibrils， troponin complex， AMPK signaling pathway， and Regulation of lipolysis in adipocytes. 2） Weighted gene co-expression network analysis showed that the genes were clustered into four modules， with the blue， black， and blue modules being significant modules. 3） Exon wide selection signature （EWSS） identified 26 467 high-quality SNPs， and ultimately combined multiple methods to screen for 6 candidate genes， including CACNA2D2， PFKFB4， ATP2B2， PRUNE2， PRKG1， and FAM13A. 4） The candidate genes were validated by qPCR， and the results were consistent with the RNA-Seq data. Protein-protein and protein-metabolite interaction networks were constructed using the STRING database and OmicsNet， which revealed that the proteins that interact with the candidate genes work together contribute to shared biological processes. In conclusion， this study ultimately identified six key differentially expressed genes （DEGs） that are likely to regulate muscle growth， lipid metabolism， and skeletal muscle fiber-type transformation. These findings provide new insights into the molecular mechanisms underlying muscle fiber-type switching and offer potential targets for improving beef quality.

Graphical abstract

关键词

转录组 / 腓肠肌 / 比目鱼肌 / 肉质 / 肌纤维类型

Key words

transcriptome / gastrocnemius / soleus muscle / meat quality / muscle fiber types

牛肉由水（70%~76%）、蛋白质（18%~23%）、肌肉脂肪（0.7%~10%）、碳水化合物（0.5%~2%）、矿物质（0.5%~2%）、维生素等组成^［1-2］，是人类饮食中蛋白质的重要来源，也是具有抗氧化特性的氨基酸和人体健康肽的良好来源^［3］，因此成为世界上消费最广泛的肉类之一。我国牛肉消费持续增长，同时消费者对牛肉品质的需求越来越高。评价牛肉品质的指标包括大理石花纹（大理石花纹程度和大理石花纹细密度）、理化（剪切力、肉色、脂肪色、质地和成熟度）和感官（嫩度、风味、多汁性、味道、气味和外观）^［4］。研究发现，肌纤维类型与许多肉质性状密切相关，如pH值、肌内脂肪、肉色、持水力、嫩度等^［5］。

肌纤维类型分为I型和II型肌纤维^［6］，根据肌球蛋白重链（MyHC）的多样性，肌纤维可分为I型、IIa型、IIx型和IIb型^［7］，它们在收缩特性、氧化能力、代谢方式、线粒体含量、形态和肌球蛋白ATP酶活性方面存在差异。哺乳动物体内的许多不同骨骼肌含有从慢肌到快肌纤维的谱系，这些肌纤维具有不同水平的蛋白质异构体，这些异构体决定了它们独特的收缩性、代谢性以及其他特性^［8］。I型肌纤维富含线粒体和肌红蛋白，收缩速度相对较慢，但持续时间较长；II型肌纤维的线粒体和肌红蛋白含量相对较低，收缩速度较快，但持续时间较短^［9］。此外，不同类型的肌纤维会在病理或生理变化下相互转化，导致不同肌肉纤维类型比例变化，纤维类型之间的转换遵循I⇔IIa⇔IIx⇔IIb的原则^［10-12］。研究表明，IIb型纤维比例高的肌肉表面颜色较浅，滴水损失较高、屠宰后糖酵解快、肉的嫩度较高^［13-14］，而I型纤维高的肌肉风味性和多汁性更好^［15］。除此之外，Umek等^［16］研究发现在快肌腓肠肌、中肌跖肌和慢肌比目鱼肌中，与控制饮食的瘦小鼠相比，肥胖小鼠的快肌腓肠肌和中肌跖肌的IIa型和IIx/d型纤维的脂质积累最为显著，而在慢肌比目鱼肌中，肥胖小鼠没有显著的脂质积累。此外，肥胖小鼠的慢肌比目鱼肌在肌纤维直径没有显著变化的情况下，表现出最显著的向快型肌球蛋白重链表达的转变，相比之下，腓肠肌和跖肌的肌纤维直径增加，对肌球蛋白重链表达变化的抵抗性更强。证明了肥胖小鼠肌内脂质积累与肌肉和纤维类型相关。因此，阐明肌纤维类型形成和转化的机制对于改善牲畜肉质具有重要意义。

比目鱼肌主要由I型肌纤维组成，腓肠肌主要由Ⅱ型肌纤维组成^［17-19］。本研究通过差异分析、WGCNA和EWSS分析的联合比较，探索2种不同肌纤维类型的差异基因表达谱以揭示控制肌纤维形成的关键候选基因和途径。进一步利用相关候选基因进行蛋白质、代谢物互作网络预测，为日后肉牛肌纤维发育和转化的研究提供依据。

1　材料与方法

1.1　试验材料

本研究共采集39个安格斯肉牛的骨骼肌样本，均为24月龄。屠宰前禁食24 h，按照标准屠宰流程共收集腓肠肌（G）样本18个，比目鱼肌（S）样本21个；牛群均来自内蒙古牧牛海畜牧发展有限公司。

1.2　腓肠肌和比目鱼肌转录组数据比对

根据质控要求，利用FastQC和Fastp软件对39个原始转录组数据进行质量评估和过滤，最终获得高质量的Clean reads（数据过滤主要参数：-Q--thread=3--length_required=50--n_base_limit=6）；使用Hisat2软件与牛参考基因组（USAD ARS-UCD2.0）建立索引；其次使用Hisat2、samtools和stringtie软件将Clean reads与参考基因组进行比对以及基因的定量（比对主要参数：hisat2 -p 1--dta-x '$index'-1 PREFIX_filter_1.fq.gz-2 PREFIX_filter_2.fq.gz-S'$OUTPUT'/align/$fastq.align.sam 2>'$OUTPUT'/align/$fastq.align.log； samtools sort-@ 18-o'$OUTPUT'/align/$fastq.align.bam'$OUTPUT'/align/$fastq.align.sam； stringtie-e-B-p 18-G'$index'_genomic.gff）。最后，使用R中的Ballgown包计算样本的基因表达量，使用FPKM值来衡量基因表达水平。

1.3　腓肠肌和比目鱼肌基因表达分析

对转录组处理后得到的基因表达量文件，利用Web平台（https：//www.omicstudio.cn/tool）对腓肠肌（G）和比目鱼肌（S）的基因表达量以及样本聚类情况进行判断。基因表达箱线图纵坐标为lgFPKM，其对数转换旨在标准化数据分布以降低高表达基因对图形尺度的影响，已删除FPKM＝0的基因。

1.4　差异表达基因分析及富集分析

本文利用R软件中的DESeq2包计算样本之间的差异表达基因（DEGs），以log₂ FC>2，P<0.05作为筛选标准，使用DAVID数据库（https：//david.ncifcrf.gov/）对差异表达基因进行功能富集分析。

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

加权基因共表达网络分析（Weighted gene co-expression network analysis， WGCNA），是一种从高通量的表达数据中挖掘模块（module）信息的算法，在该方法中module被定义为一组具有类似表达谱的基因，它可以探索module内基因网络与研究者关注的表型之间的关联关系。在快肌和慢肌的39个样本中，利用Web平台（https：//cloud.oebiotech.com/）对FPKM基因表达矩阵文件进行WGCNA分析，过滤表达量低的基因（标准偏差≤0.5）、网络构建参数power值为30，基于power值进行模块识别；最后进行模块性状关联分析，相关系数绝对值≥0.3且P<0.05，筛选各个性状相关的模块。对每个性状相关的模块，分别计算模块基因表达量与对应性状的相关性（Gene significance， GS），同时计算模块基因表达量与模块特征基因的相关性。

1.6　外显子选择信号分析

外显子选择信号（Exon wide selection signature，EWSS）是一种基于转录水平上以群体分化指数（F_ST）作为衡量群体分化程度的指标，进行SNP检测和筛选的高效突变检测方法^［20］。本研究主要使用samtools和bcftools软件统计牛参考序列（GCF_002263795.3_ARS-UCD2.0_genomic.fna）中每个碱基位点的比对情况，并对腓肠肌和比目鱼肌之间均表达的SNP位点进行筛选并过滤掉DP>10的SNP位点。以F_ST作为筛选标准，利用vcftools软件计算，保留F_ST为正数的位点。F_ST值为0~0.05表示分化程度小，0.05~0.15表示在种群间中度分化，0.15~0.25为高度分化，>0.25为极高度分化^［21］，F_ST值越大表明2个群体间分化程度越高，亲缘关系越远；F_ST值越小表明群体间分化程度越低，亲缘关系越近。本研究采用滑动窗口法定位候选基因（2 kb滑动步长），使用R软件对SNP进行基因注释以及可视化，最后联合差异表达基因、WGCNA和EWSS分析确定与肌肉生长发育相关的候选基因。

1.7　相关基因表达的qPCR检测

用RNAiso Plus试剂盒（TaKaRa，大连）提取腓肠肌和比目鱼肌的总RNA，使用反转录试剂盒（TaKaRa，大连）将RNA反转录成cDNA。使用ABI7500Real-time PCR系统（Applied Biosystems，美国） 扩增细胞周期相关基因。PCR引物通过Primer Premier 5软件设计见表1。PCR反应体系（20 μL）：2×SYBR GreenⅠMix 10 μL，上下游引物（10 μmol/L）各0.8 μL，模板1 μL，ROX 0.4 μL，RNase-free 水 7 μL；PCR 反应程序：95 ℃ 30 s；95 ℃ 5 s，60 ℃ 34 s，循环40 次。以 GAPDH作为内参基因，以2^-ΔΔCt方法计算基因相对表达量^［22］。

1.8　蛋白质以及代谢物互作网络分析

本研究主要通过使用STRING（https：//cn.string-db.org/）数据库对候选基因进行蛋白质网络预测，同时利用OmicsNet（https：//www.omicsnet.ca/）平台中的STRING和KEGG数据库预测候选基因与蛋白质和代谢物之间的互作网络关系，以实现基因、蛋白质和代谢物之间的可视化。

1.9　统计分析

利用qPCR验证基因的表达，通过t检验分析对筛选到的关键基因进行差异显著性分析（P<0.05显著，P>0.05不显著）。

2　结果与分析

2.1　测序质量控制与参考基因组比对分析

利用牛的腓肠肌和比目鱼肌组织成功构建cDNA测序文库。共产生304.14 G原始数据，过滤后得到282.26 G有效数据，平均有效率为92.81%，Q30碱基百分比平均为97.78%，各样本Reads与参考基因组的比对效率为95.95%~97.45%。序列质量良好，适用于后续分析。

2.2　腓肠肌和比目鱼肌基因表达分析

对39个样本进行转录组分析，共鉴定到33 671个肌肉比较表达基因，FPKM表达量分布如图1（a）。对基因表达谱进行主成分分析（PCA），通过数据降维来观察腓肠肌和比目鱼肌样品之间的基因表达差异，发现比目鱼肌与腓肠肌样品间明显分离（图1（b）），表明可以通过这2种肌肉得到可靠的差异基因。

2.3　差异表达基因分析

为更好地研究腓肠肌和比目鱼肌这2种肌肉之间基因表达模式的差异，使用DESeq2来识别2种肌肉样本之间的差异基因。根据log₂ FC>2，P<0.05的标准，发现两组间有1 576个显著差异基因，腓肠肌与比目鱼肌相比有170个基因上调和1 406个基因下调（图2（a））。

对差异表达基因进行功能富集分析（P<0.05）发现，在下调基因中，主要富集在GO微管结合（Microtubule binding）、蛋白质磷酸化正向调控（Positive regulation of protein phosphorylation）和钙离子结合（Calcium ion binding）以及平行纤维-浦肯野细胞突触（Parallel fiber to purkinje cell synapse）；在KEGG中主要通路包括脂肪细胞脂解调节（Regulation of lipolysis in adipocytes）、丝裂原活化蛋白激酶信号通路（MAPK signaling pathway）、AMP活化蛋白激酶信号通路（AMPK signaling pathway）和脂肪细胞因子信号通路（Adipocytokine signaling pathway）（图2（b）~（c））；在上调基因中GO主要富集在糖酵解过程（Glycolytic process）、横纹肌收缩（Striated muscle contraction）、肌动蛋白丝结合（Actin filament binding）、肌球蛋白Ⅱ复合物（Myosin II complex）和肌肉收缩（Muscle contraction）等；KEGG中上调基因主要富集在糖酵解/糖异生（Glycolysis/Gluconeogenesis）、肌肉细胞骨架系统（Cyboskeleton in muscle cells）等通路里（图2（d）~（e））。

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

为了提高网络构建的准确性，过滤后最终共获得4 337个基因。根据power值构建共表达网络识别出5个模块，同一种颜色表示同一个模块，由于某些模块之间存在一定相关性，所以将对应的模块合并为同一个模块，即下方的Merged Dynamic为最终得到的模块，这些模块用于后续的分析，其中灰色模块（Grey）不能归属到任意模块的基因集合，没有参考意义（图3（a））。结果发现，聚集在蓝色（Blue）模块的基因数量最多 （2 667个），品红色（Magenta）模块包含的基因数量最少（140个），红色（Red）模块基因数量（571），黑色（Black）模块基因数量（198）。按照相关系数绝对值≥0.3且P<0.05，筛选各个性状相关的模块。品红色、黑色、蓝色模块为显著模块，与快慢肌的相关系数分别为0.93、0.86和0.80（图3（b））。

2.5　外显子选择信号分析

通过使用Samtoools和Bcftools软件对牛快慢肌进行SNP位点筛选（DP>10），共得到180 370个SNP位点，删除负值（F_ST）后得到26 467个SNP，其中F_ST在0.05~0.15内有1 655个SNP，F_ST在0.15~0.25有1 218个SNP，F_ST >0.25有11个SNP。曼哈顿图展示了全部的外显子选择信号位点，红色线以上代表F_ST>0.15的SNP（图4（a））。通过对DEGs、WGCNA以及EWSS的整合取交集，发现CACNA2D2、PFKFB4、ATP2B2、PRUNE2、PRKG1和FAM13A基因可能是影响肌肉生长发育的关键基因（图4（b）、表2）。筛选得到的关键基因F_ST值在0.03~0.05，SNP分化程度较低，原因是样品来自同一个群体，亲缘关系较近。

2.6　相关基因表达的qPCR检测

对候选基因进行qPCR检测。结果显示，下调基因CACNA2D2、PFKFB4和ATP2B2的表达量显著高于S组，而在上调基因PRUNE2、PRKG1和FAM13A中S组的表达量高于G组。结果与RNA-Seq测序结果一致（图5）。

2.7　候选基因的蛋白质以及代谢物互作网络分析

FAM13A受PIGY、HERC3和GPR61等蛋白的相互调控，从而可能促成FAM13A在脂肪组织中的高表达（图6（a））。ATP2B2编码的蛋白质属于P型初级离子转运ATP酶家族，在该蛋白的互作网络中发现，ATP2A1、ATP2A3、ATP2A2和MYO6等互相作用，MYO6是基于肌动蛋白的运动分子，该蛋白向肌动蛋白丝的负端移动，ATP2A1是肌质/内质网钙ATP酶1是影响横纹肌功能的关键因子（图6（b））。PRKG1受MAPK14和GUCY1A2等蛋白调控的同时与蛋白质代谢物相关（图6（c）、表3）。PFKFB4基因是糖代谢调控中重要的一个因子，与β-D-呋喃果糖-6-磷酸（Beta-D-Fructofuranose 6-phosphate）、果糖-6-磷酸（Neuberg ester）和D-果糖2，6-二磷酸（D-Fructose 2，6-bisphosphate）代谢物作用（图6（d）、表3）。总的来说，在蛋白质及代谢物互作网络中，发现与关键基因互作的蛋白质它们之间共同作用促成了一个共享功能，这为促进肌纤维转换和肌肉生长以及改善肉质提供一定见解。

3　讨 论

肉质对肉牛产业的经济价值至关重要，随着高通量测序技术的迅速发展，大量研究者应用该技术挖掘调控牛肉质的基因表达模式^［23］。研究样本多集中在背最长肌，Yan等^［24］以新疆棕牛和哈萨克牛为研究对象，对背最长肌的转录组和蛋白质组分析，鉴定影响肉质的功能基因和蛋白质。华西牛是我国自主培育的肉牛新品种，Deng等^［25］通过转录组学分析其背最长肌，探究补偿性生长下调控肉质的分子机制。持水量对肉质有很大影响，通过RNA测序鉴定中国西门塔尔肉牛背最长肌影响持水量的候选基因^［26］。其他评估肉质的因素，如嫩度、多汁性、大理石纹、剪切力、肉色、pH值、蒸煮损失等，研究者们以转录组为主要手段分析不同品种牛的背最长肌来探索各肉质性状的潜在分子机制与相关的候选基因和生物学指标^［27-33］。也有研究者分析除背最长肌外的其他部位肌肉来探究影响肉质差异的分子机制，包括腰肌、肋骨肌、胸肌、小腿、颈部、臀部、股二头肌、腹外斜肌等^［34-36］。本研究以肉牛腓肠肌和比目鱼肌为研究对象，构建2种肌肉的转录图谱，为肉牛全骨骼肌转录组与牛肉质研究补充数据。

研究发现，不同的肌纤维类型会影响肉质^［37］。根据肌球蛋白重链（Myosin heavy chain， MHC）的类型，肌纤维可分为I型纤维和II型纤维。I型纤维，又称“慢氧化纤维”，含有慢肌纤维（Myhc slow）的肌球蛋白，收缩速度慢，抗疲劳性强。II型纤维被称为“快速糖酵解纤维”，它们含有Myhc，主要从厌氧糖酵解中获取能量^［38］。迄今为止，一些基因已被阐明在肌纤维型转化中发挥重要作用，包括肌细胞增强因子2（MEF2）、肌生长抑制素（MSTN）、肌纤维调节因子（MRFs）、胰岛素样生长因子（IGFs）、过氧化物酶体增殖物激活受体γ（PPARγ）、forkhead转录因子1（FOXO1）、丝裂原活化蛋白激酶（MAPK）、MYH6、钙调神经磷酸酶等^［39-40］。钙调神经磷酸酶它作为一种钙介导的丝氨酸-苏氨酸磷酸酶，在许多类型的组织中表达。钙调神经磷酸酶参与许多生物学功能，如T淋巴细胞激活、血管和心脏发育以及骨骼肌发育^［41］。骨骼肌纤维型重塑也涉及多个关键的信号通路，包括TGF-β/Smad、MAPK/ERK、PI3K/AKT、CaN/NFAT等^［42-43］，这与本研究结果一致，在差异表达基因中也发现了调控慢肌的基因MYH7B又名MYH14以及在WGCNA品红色（Magenta）模块中发现快肌基因MYH2。

在本研究中筛选到的关键基因CACNA2D2、PFKFB4、ATP2B2、PRUNE2、PRKG1和FAM13A，发现FAM13A在全基因组关联分析显示，它与慢性阻塞性肺疾病（Chronic obstructive pulmonary disease，COPD）有关^［44］，之前一项研究表明，在喂食高脂肪细胞饮食的小鼠和对照小鼠中进行DNA微阵列分析^［45］，确定了表达下调的前20个基因，其中包括FAM13A基因。此外还分析了FAM13A在成年牛不同组织中的表达模式，构建FAM13A氨基酸序列系统发育树，发现除肺组织外，皮下脂肪组织在所有组织中的表达最高。证明FAM13A基因在牛脂肪细胞中的提供了功能基础，以改善肉牛肌内脂肪沉积^［46］。此前研究证明在人类脂肪细胞模型中，前脂肪细胞中FAM13A的敲除加速了脂肪细胞的分化。在小鼠中，与野生型相比，敲除（KO）FAM13A基因的小鼠在高脂肪饮食饲喂后，其内脏脂肪与皮下脂肪（VAT/SAT）的比例较低，KO小鼠的皮下脂肪细胞，发现大小分布更多向小型脂肪细胞转移的趋势，同时伴随着脂肪生成能力的提高^［47］。CACNA2D2基因在GO注释中包括电压门控钙通道活性（Voltage-gated calcium channel activity）和KEGG中MAPK信号通路（MAPK signaling pathway）；在肌肉和骨骼中均有表型影响动物生长，先前研究发现，CACNA2D2 突变对癫痫性脑病（EE）的贡献尚不清楚，在对一名患有癫痫、运动障碍、小脑萎缩、精神运动迟缓和畸形特征的患者进行了全外显子组测序（WES），发现2个新的纯合候选突变：CACNA2D2中的c.1295delA（p.Asn432fs）和CELSR3中的c.G6407A（p.Gly2136Asp），基因优先排序表明CACNA2D2是最突出的候选基因^［48］。靶向破坏CACNA2D2纯和小鼠表现出共济失调，生长迟缓，对癫痫发作的易感性增加，寿命缩短^［49］；由于小脑功能障碍引起的共济失调导致各种基本神经功能缺损，从而引起行走运动迟缓、肌阵挛和肌张无力等症状。在骨骼肌中，糖酵解和氧化磷酸化的相对表达量和平衡对于肌肉和骨骼肌代谢至关重要，在一项研究环状RNA对快肌和慢肌在肌肉纤维代谢的作用中发现，circMYLK4与CACNA2D2相互作用导致肌浆网中Ca²⁺释放被抑制从而抑制糖酵解，但是脂肪酸β-氧化增强了三羧酸循环和线粒体氧化磷酸化^［50］，而在上述的富集分析中也发现基因显著富集在糖酵解和脂肪酸代谢过程中这与本研究结果一致。除此之外，PFKFB4基因在AMP活化蛋白激酶信号通路（AMPK signaling pathway）显著富集，研究表明，AMPK在调节脂肪和胆固醇代谢方面发挥着重要作用，尤其是在肌肉、脂肪组织和肝脏等特定组织中，通过促进脂肪酸转运蛋白（FAT）/CD36向细胞膜转移来增加脂肪酸的摄取^［51］。PRKG1是cGMP依赖性蛋白激酶1，PRKG1对RHOA的磷酸化阻断了该蛋白在无数过程中的作用，比如调节RHOA易位；减少收缩；控制囊泡运输，减少肌球蛋白轻链磷酸化等；PRKG1蛋白除了放松平滑肌张力、防止血小板聚集和调节细胞生长外，还在调节心血管和神经元功能方面发挥着核心作用。在所有类型的平滑肌、血小板、小脑浦肯野细胞、海马神经元和外侧杏仁核中表达最强烈在平滑肌细胞中^［52］，cGMP和PRKG1活性的增加会影响平滑肌特异性收缩蛋白的表达。有人通过利用GWAS证实PRKG1基因在脂肪细胞中通过三酰甘油的水解来释放脂肪酸和甘油，在脂肪分解过程中发挥关键作用，同时发现PRKG1对奶牛牛奶脂肪酸有显著影响^［53］。上述富集分析也发现PRKG1基因在脂肪细胞脂解调节（Regulation of lipolysis in adipocytes）通路中富集，通过水解释放脂肪酸（FA）和甘油，供其他器官用作能量底物，这是白色脂肪组织的独特功能^［54］。ATP2B2属于P型初级离子转运ATP酶家族，其在细胞内钙稳态中发挥关键作用，富集在钙离子结合和钙调蛋白结合中，精确的钙调节对神经发育和小脑功能的重要性，研究发现ATP2B2变体在引起异质性神经发育和运动障碍综合征中作用^［55］。cAMP作为一种第二信使，在脂解过程中具有重要作用，发现ATP2B2参与cAMP通路且在脂解途径中被下调，cAMP通过影响脂肪组织中的脂质代谢，参与了对猪饲料效率的调节^［56］。除此之外，在互作网络关系中，发现ATP2B2与ATP2A1相互作用，ATP2A1是位于肌肉细胞肌浆或内质网中的细胞内泵，横纹肌性能的关键调节因子，有助于参与肌肉兴奋/收缩的钙封存^［57］。总之，这些候选基因及其参与的代谢调控通路在骨骼肌生长、发育和脂质代谢中发挥着不可或缺的作用。

4　结 论

本研究通过对腓肠肌和比目鱼肌的转录组数据进行差异分析、WGCNA以及EWSS联合分析，最终确定了CACNA2D2、PFKFB4、ATP2B2、PRUNE2、PRKG1和FAM13A候选基因，以及富集到的调节脂肪和胆固醇代谢的AMPK信号通路，这些基因和信号通路可能是调控肌肉生长脂质代谢以及骨骼肌纤维转化的关键基因。研究结果将为探索肌纤维类型转化的分子机制提供有效参考。

参考文献

[1] Park S J， Beak S H， Jung D J， Kim S Y， Jeong I H， Piao M Y， Kang H J， Fassah D M， Na S W， Yoo S P， Baik M. Genetic， management， and nutritional factors affecting intramuscular fat deposition in beef cattle： A review ［J］ Asian-Australasian.Journal of Animal Sciences， 2018，31（7）：1043-1061

[2] Sakowski T， Grodkowski G， Gołebiewski M， Slósarz J， Kostusiak P， Solarczyk P， Puppel K. Genetic and environmental determinants of beef quality-a review［J］. Frontiers in Veterinary Science， 2022，24（9）：819605

[3] Wu G， Cross H R， Gehring K B， Savell J W， Arnold A N， McNeill S H. Composition of free and peptide-bound amino acids in beef chuck， loin， and round cuts ［J］. Journal of Animal Science， 2016，94（6）：2603-2613

[4] Baik M， Lee J， Kim S Y， Ranaweera K K T N. Invited review： Factors affecting beef quality and nutrigenomics of intramuscular adipose tissue depositions ［J］. Animal Bioscience， 2023，36（2）：350-363

[5] 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

[6] Schiaffino S， Reggiani C. Fiber types in mammalian skeletal muscles［J］. Physiological Reviews， 2011，91（4）：1447-1531

[7] Song S M， Ahn C H， Kim G D. Muscle fiber typing in Bovine and porcine skeletal muscles using immunofluorescence with monoclonal antibodies specific to myosin heavy chain isoforms［J］. Food Science of Animal Resources， 2020，40（1）：132-144

[8] Lloyd E M， Pinniger G J， Murphy R M， Grounds M D. Slow or fast： Implications of myofibre type and associated differences for manifestation of neuromuscular disorders［J］. Acta Physiologica， 2023，238（4）：e14012

[9] Picard B， Lefaucheur L， Berri C， Duclos M J. Muscle fibre ontogenesis in farm animal species［J］. Reproduction，Nutrition Development， 2002，42（5）：415-431

[10] Campos G E， Luecke T J， Wendeln H K， Toma K， Hagerman F C， Murray T F， Ragg K E， Ratamess N A， Kraemer W J， Staron R S. Muscular adaptations in response to three different resistance-training regimens： Specificity of repetition maximum training zones ［J］. European Journal of Applied Physiology， 2002，88（1）：50-60

[11] Schiaffino S， Reggiani C. Fiber types in mammalian skeletal muscles［J］. Physiological Reviews， 2011，91（4）：1447-1531

[12] Ying F， Zhang L， Bu G W， Xiong Y Z， Zuo B. Muscle fiber-type conversion in the transgenic pigs with overexpression of PGC1α gene in muscle［J］. Biochemical and Biophysical Research Communications， 2016，480（4）：669-674

[13] Choe J H， Choi Y M， Lee S H， Shin H G， Ryu Y C， Hong K C， Kim B C. The relation between glycogen lactate content and muscle fiber type composition， and their influence on postmortem glycolytic rate and pork quality［J］. Meat Science， 2008，80（2）：355-362

[14] Lee S H， Choe J H， Choi Y M， Jung K C， Rhee M S， Hong K C， Lee S K， Ryu Y C， Kim B C. The influence of pork quality traits and muscle fiber characteristics on the eating quality of pork from various breeds［J］. Meat Science， 2012，90（2）：284-291

[15] Sibut V， Hennequet-Antier C， Le Bihan-Duval E， Marthey S， Duclos M J， Berri C. Identification of differentially expressed genes in chickens differing in muscle glycogen content and meat quality［J］. BMC Genomics， 2011，12：112

[16] Umek N， Horvat S， Cvetko E.. Skeletal muscle and fiber type-specific intramyocellular lipid accumulation in obese mice［J］. Bosnian Journal of Basic Medical Sciences， 2021，21（6）：730-738

[17] Arany Z， Lebrasseur N， Morris C， Smith E， Yang W L， Ma Y H， Chin S， Spiegelman B M. The transcriptional coactivator PGC-1beta drives the formation of oxidative type IIX fibers in skeletal muscle［J］. Cell Metabolism， 2007，5（1）：35-46

[18] Johnson M A， Polgar J， Weightman D， Appleton D. Data on the distribution of fibre types in thirty-six human muscles： An autopsy study ［J］. Journal of the Neurological Sciences， 1973，18（1）：111-129

[19] Pette D， Spamer C. Metabolic properties of muscle fibers［J］. Federation proceedings， 1986，45（13）：2910-2914

[20] 刘在霞，刘永斌，石彩霞，谷明娟，满达，王宇，付绍印，何小龙，戴伶俐，白音巴图，鲍艳春，张文广，娜日苏.不同繁殖力绵羊卵巢转录组比较分析［J］.中国农业大学学报， 2023，28（12）：132-143

[21] Liu Z X， Liu Y B， Shi C X， Gu M J， Man D， Wang Y， Fu S Y， He X L， Dai L L， Bai Y， Bao Y C， Zhang W G， Na R S. Comparative analysis of ovary transcriptome of sheep with different fecundity［J］. Journal of China Agricultural University， 2023，28（12）：132-143 （in Chinese）

[22] Hall S J G. Genetic differentiation among livestock breed-values for F_st ［J］. Animals， 2022，12（9）：1115

[23] 高丽，谷明娟，杨苗苗，刘云鹏，卜彩虹，程鹏，白春玲. MSTN基因突变通过SMAD2/SMAD3-CDKN1C调控牛肌肉卫星细胞增殖分化［J］.农业生物技术学报， 2022，30（6）：1128-1139

[24] Gao L， Gu M J， Yang M M， Liu Y P， Bu C H， Cheng P， Bai C L. MSTN mutation regulates proliferation and differentiation of bovine （ Bos taurus ） muscle-derived satellite cells through SMAD2/SMAD3 CDKN1C ［J］. Journal of Agricultural Biotechnology， 2022，30（6）：1128-1139 （in Chinese）

[25] Zhang F W， Hanif Q， Luo X Y， Jin X D， Zhang J C， He Z X， Lei C Z， Liu J Y， Huang B Z， Qu K X. Muscle transcriptome analysis reveal candidate genes and pathways related to fat and lipid metabolism in Yunling cattle［J］. Animal Biotechnology， 2023，34（4）：1022-1029

[26] Yan X M， Wang J， Li H B， Gao L， Geng J， Ma Z， Liu J M， Zhang J S， Xie P G， Chen L. Combined transcriptome and proteome analyses reveal differences in the longissimus dorsi muscle between Kazakh cattle and Xinjiang brown cattle［J］. Animal Bioscience， 2021，34（9）：1439-1450

[27] Deng T Y， Liang M， Du L L， Li K， Li J N， Qian L， Xue Q Q， Qiu S Y， Xu L Y， Zhang L P， Gao X， Li J Y， Lan X Y， Gao H J. Transcriptome analysis of compensatory growth and meat quality alteration after varied restricted feeding conditions in beef cattle［J］. International Journal of Molecular Sciences， 2024，25（5）：2704-2725

[28] Du L L， Chang T P， An B X， Liang M， Duan X H， Cai W T， Zhu B， Gao X， Chen Y， Xu L Y， Zhang L P， Li J Y， Gao H J. Transcriptome profiling analysis of muscle tissue reveals potential candidate genes affecting water holding capacity in Chinese Simmental beef cattle［J］. Scientific Reports， 2021，11（1）：11897

[29] Leal-Gutiérrez J D， Elzo M A， Carr C， Mateescu R G. RNA-seq analysis identifies cytoskeletal structural genes and pathways for meat quality in beef［J］. PLoS One， 2020，15（11）：e0240895

[30] Li G H， Yang R J， Lu X， Liu Y， He W， Li Y， Yu H B， Qin L H， Cao Y， Zhao Z H， Fang X B. RNA-seq analysis identifies differentially expressed genes in the longissimus dorsi of Wagyu and Chinese Red Steppe Cattle［J］. International Journal of Molecular Sciences，2022，24（1）：387

[31] Li Y X， Yang M S， Lou A G， Yun J Y， Ren C Y， Li X C， Xia G J， Nam K C， Yoon D， Jin H G， Seo K S， Jin X. Integrated analysis of expression profiles with meat quality traits in cattle［J］. Scientific Reports， 2022，12（1）：5926

[32] Meng X R， Gao Z W， Liang Y S， Zhang C L， Chen Z， Mao Y J， Huang B Z， Kui K X， Yang Z P. Longissimus dorsi muscle transcriptomic analysis of Simmental and Chinese native cattle differing in meat quality［J］. Frontiers in Veterinary Science， 2020，12（7）：601064

[33] Silva D B S， Fonseca L F S， Pinheiro D G， Magalhães A F B， Muniz M M M， Ferro J A， Baldi F， Chardulo L A L， Schnabel R D， Taylor J F， Albuquerque L G. Spliced genes in muscle from Nelore Cattle and their association with carcass and meat quality［J］. Scientific Reports， 2020，10（1）：14701

[34] Wu G S， Qiu X J， Jiao Z Z， Yang W J， Pan H J， Li H， Bian Z Y， Geng Q， Wu H， Jiang J M， Chen Y Y， Cheng Y W， Chen Q L， Chen S， Man C， Du L， Li L B， Wang F Y. Integrated analysis of transcriptome and metabolome profiles in the longissimus dorsi muscle of buffalo and cattle［J］. Current Issues in Molecular Biology， 2023，45（12）：9723-9736

[35] Yu H W， Yu S C， Guo J T， Wang J F， Mei C G， Abbas Raza S H， Cheng G， Zan L S. Comprehensive analysis of transcriptome and metabolome reveals regulatory mechanism of intramuscular fat content in beef cattle［J］. Agriculture Food Chemistry， 2024，2（6）：2911-2924

[36] Wu X Y， Zhou X L， Chu M， Guo X， Pei J， Xiong L， Ma X M， Bao P J， Liang C N， Yan P. Whole transcriptome analyses and comparison reveal the metabolic differences between oxidative and glycolytic skeletal muscles of yak［J］. Meat Science， 2022，194：108948

[37] Yu H W， Wang J F， Zhang K， Cheng G， Mei C G， Zan L S. Integrated multi-omics analysis reveals variation in intramuscular fat among muscle locations of Qinchuan cattle［J］. BMC Genomics， 2023，24（1）：367

[38] Zhang T L， Wang T Z， Niu Q H， Zheng X， Li H P， Gao X， Chen Y， Gao H J， Zhang L P， Liu G E， Li J Y， Xu L Y. Comparative transcriptomic analysis reveals region-specific expression patterns in different beef cuts［J］. BMC Genomics， 2022，23（1）：387

[39] Pette D， Staron R S. Transitions of muscle fiber phenotypic profiles［J］. Histochemistry and Cell Biology， 2001，115（5）：359-372

[40] Shi S， Xia Y H， Chen S C， Li M， Chen D H， Liu F， Zheng H L. The relationship between structural/MHC changes in upper airway palatopharyngeal muscle morphology and obstructive sleep apnea/hypopnea syndrome［J］. European Archives of Oto-Rhino-Laryngology， 2014，271（1）：109-116

[41] Park S K， Sheffler T L， Spurlock M E， Grant A L， Gerrard D E. Chronic activation of 5'-AMP-activated protein kinase changes myosin heavy chain expression in growing pigs ［J］. Journal of Animal Science， 2009，87（10）： 3124-3133

[42] Zhu J Y， Shi X E， Lu H Z， Xia B， Li Y F， Li X， Zhang Q L， Yang G S. RNA-seq transcriptome analysis of extensor digitorum longus and soleus muscles in large white pigs［J］. Molecular Genetics Genomics， 2016，291（2）：687-701

[43] Liu Y F， Zhang M， Shan Y J， Ji G G， Ju X J， Tu Y J， Sheng Z W， Xie J F， Zou J M， Shu J T. miRNA-mRNA network regulation in the skeletal muscle fiber phenotype of chickens revealed by integrated analysis of miRNAome and transcriptome ［J］. Scientific Reports， 2020，10（1）：10619

[44] Chin E R， Olson E N， Richardson J A， Yang Q， Humphries C， Shelton J M， Wu H， Zhu W， Bassel-Duby R， Williams RS. A calcineurin-dependent transcriptional pathway controls skeletal muscle fiber type［J］. Gene & Development， 1998，12（16）：2499-2509

[45] Vainshtein A， Sandri M. Signaling pathways that control muscle mass［J］. International Journal of Molecular Sciences， 2020，21（13）：4759

[46] Zhang Y N， Qiu J， Zhang P， Zhang J， Jiang M， Ma Z B. Genetic variants in FAM13A and IREB2 are associated with the susceptibility to COPD in a Chinese rural population： A case-control study ［J］. International Journal of Chronic Obstructive Pulmonary Disease. 2018，13：1735-1745

[47] Wardhana D A， Ikeda K， Barinda A J， Nugroho D B， Qurania K R， Yagi K， Miyata K， Oike Y， Hirata K I， Emoto N. Family with sequence similarity 13， member A modulates adipocyte insulin signaling and preserves systemic metabolic homeostasis［J］. Proceedings of the National Academy of Sciences of The United States of America. 2018，115（7）：1529-1534

[48] Liang C C， Li A N， Raza S H A， Khan R， Wang X Y， Wang S H， Wang G H， Zhang Y， Zan L S. The molecular characteristics of the FAM13A gene and the role of transcription factors ACSL1 and ASCL2 in its core promoter region［J］. Genes， 2019，10（12）：981

[49] Fathzadeh M， Li J H， Rao A， Cook N， Chennamsetty I， Seldin M， Zhou X， Sangwung P， Gloudemans M J， Keller M， Attie A， Yang J， Wabitsch M， Carcamo-Orive I， Tada Y， Lusis A J， Shin M K， Molony C M， McLaughlin T， Reaven G， Montgomery S B， Reilly D， Quertermous T， Ingelsson E， Knowles J W. FAM13A affects body fat distribution and adipocyte function［J］. Nature Communications， 2020，11（1）：1465

[50] Pippucci T， Parmeggiani A， Palombo F， Maresca A， Angius A， Crisponi L， Cucca F， Liguori R， Valentino M L， Seri M， Carelli V. A novel null homozygous mutation confirms CACNA2D2 as a gene mutated in epileptic encephalopathy［J］. Plos One， 2013，8（12）：e82154

[51] Ivanov S V， Ward J M， Tessarollo L， McAreavey D， Sachdev V， Fananapazir L， Banks M K， Morris N， Djurickovic D， Devor-Henneman D E， Wei M H， Alvord G W， Gao B N， Richardson J A， Minna J D， Rogawski M A， Lerman M I. Cerebellar ataxia， seizures， premature death， and cardiac abnormalities in mice with targeted disruption of the CACNA2D2 gene ［J］. American Journal of Pathology， 2004，165（3）：1007-1018

[52] Cao H G， Li C C， Sun X H， Yang J J， Li X， Yang G S， Jin J J， Shi X E. Circular RNA circMYLK4 shifts energy metabolism from glycolysis to OXPHOS by binding to the calcium channel auxiliary subunit CACNA2D2［J］. The Journal of Biological Chemistry， 2024，300（7）：107426

[53] Sanchez A M J， Candau R B， Csibi A， Pagano A F， Raibon A， Bernardi H. The role of AMP-activated protein kinase in the coordination of skeletal muscle turnover and energy homeostasis［J］. American Journal of Physiology Cell Physiology， 2012，303（5）：C475-485

[54] Sellak H， Lincoln T M， Choi C S. Stabilization of cGMP-dependent protein kinase G （PKG） expression in vascular smooth muscle cells： Contribution of 30UTR of its mRNA ［J］. Journal of Biochemistry， 2011，149（4）：433-441

[55] Shi L J， Lv X Q， Liu L， Yang Y Z， Ma Z， Han B， Sun D X. A post-GWAS confirming effects of PRKG1 gene on milk fatty acids in a Chinese Holstein dairy population［J］. BMC Genetics， 2019，20（1）：53

[56] Cho C H， Patel S， Rajbhandari P. Adipose tissue lipid metabolism： Lipolysis［J］. Current Opinion in Genetics & Development， 2023，83（12）：102114

[57] Poggio E， Barazzuol L， Salmaso A， Milani C， Deligiannopoulou A， Cazorla Á G， Jang S S， Juliá-Palacios N， Keren B， Kopajtich R， Lynch S A， Mignot C， Moorwood C， Neuhofer C， Nigro V， Oostra A， Prokisch H， Saillour V， Schuermans N， Torella A， Verloo P， Yazbeck E， Zollino M， Jech R， Winkelmann J， Necpal J， Calì T， Brini M， Zech M. ATP2B2 de novo variants as a cause of variable neurodevelopmental disorders that feature dystonia， Ataxia， intellectual disability， behavioral symptoms， and seizures ［J］. Genetics in Medicine， 2023，25（12）：100971

[58] Xu Y Y， Qi X L， Hu M Y， Lin R Y， Hou Y， Wang Z X， Zhou H H， Zhao Y X， Luan Y， Zhao S H， Li X Y. Transcriptome analysis of adipose tissue indicates that the cAMP signaling pathway affects the feed efficiency of pigs［J］. Genes， 2018，9（7）：336

[59] Odermatt A， Barton K， Khanna V K， Mathieu J， Escolar D， Kuntzer T， Karpati G， MacLennan D H. The mutation of Pro789 to Leu reduces the activity of the fast-twitch skeletal muscle sarco （endo） plasmic reticulum Ca²⁺ ATPase （SERCA1） and is associated with Brody disease ［J］. Human Genetics， 2000，106（5）：482-491

基金资助

内蒙古自治区揭榜挂帅项目(2022JBGS0025)

内蒙古自治区重大专项(2021ZD0009-2)

内蒙古自治区科技计划(2023YFHH0058)

内蒙古农业大学高水平成果培育专项(CG202432)