鸡胸肌条纹/木质表型相关基因的筛选

樊庆灿，郑路程，王倩倩，刘 犇,2，杨 雪，胡 威,2，郑文亚

（1.宜春学院生命科学与资源环境学院，江西宜春 336000；2.江西绿科农牧科技有限公司，江西宜春 336000）

过去几年中，快大型肉鸡的胸肌中白条纹肉（White Striping，WS）和木质肉（Woody Breast，WB）异常表型不断上升。WS 存在数量不等的与肌纤维平行的白色条纹[1]。WB 呈木质化，发病严重的肌肉表面可同时出现白色条纹和弥漫分布的脊状突起[2]，严重影响鸡肉的屠宰外观和肉品质。WS 和WB 表型有所不同，但肌肉组织均有所损伤。组织学评估表明，异常样品均呈现肌肉表层的退化，包括纤维坏死和溶解、异常纤维化和脂质增多[3]。为了阐明WS/WB 表型发展的分子机制，Zambonelli 等[4]比较了的Ross 708 肉鸡WS/WB胸肌和正常胸肌基因表达模式，鉴定出的差异表达基因（DEG）与肌肉发育、炎症、多糖代谢和钙信号传导有关。Malila 等[5]则筛选出钙信号传导和Ras 相关蛋白1信号传导等8 个信号通路与WS/WB 的产生有关，并发现细胞内Ca2+等离子、氧化应激和细胞程序性死亡引起WS/WB 表型。本研究利用GEO2R 软件的相同标准分别分析2 组表达谱数据，统计快大型肉鸡正常肌肉表型和WS/WB 肌肉表型的差异表达基因，以进行共同基因分析，并对筛选到的共同基因进行功能注释和功能富集分析，探讨WS/WB 肌肉表型的发生机理，为鸡屠宰和肉质性状的研究提供参考。

1 材料和方法

1.1 数据及样本 本研究数据来自于美国国家生物信息数据库中的GEO 数据库，编号分别为GSE79276 和GSE107362，2 个数据系列的实验动物分别为罗斯708和罗斯308 的公鸡，屠宰日龄均为49 d，样品均为胸大肌，样品表达量检测分别采用GPL21385 和GPL24307平台。具体信息见表1。

1.2 软件与方法 打开GEO 数据库中的分析工具GEO2R（https://www.ncbi.nlm.nih.gov/geo/geo2r/），输入GSE 79276，按照表1 设定正常表型组（NM）和WS/WB 表型组（WS/WB），选择Benjamini &Hochberg（False discovery rate）矫正，数据进行log2转换，进行差异表达分析。下载差异分析结果、数据表达量图和差异分析火山图，将Padj＜0.05 且|log2FC|＞1 的转录本设为差异表达转录本（Differential Expression Transcriptions，DET）。根据GPL21385 平台信息，查找DET 对应的基因，设为差异表达基因（Differential Expression Gene，DEG）。删除无基因信息的DETs。在DAVID V6.8（https://david.ncifcrf.gov/tools.jsp）中进行官方基因名转换。GSE107362 数据分析步骤如上。将GSE79276 和GSE107362 中的DEGs 输入到Bioinformatics &Evolut ionary Genomics（http://bioinformatics.psb.ugent.be/webtools/Venn/）中绘制韦恩图，下载共同差异表达基因（Identical Differential Expression Gene，IDEG）。利用David v6.8 进行IDEGs 的信号通路分析。将IDEGs 导入到String（https://string-db.org/cgi/input.pl?sessionId=NPdugdyy2MoP&input_page_show_search=on）中进行互作分析，并下载蛋白互作文件，导入Cytoscape v3.7.2 中绘制互作图，cytoHubba 插件MHC 算法计算核心基因，BINGO 插件进行生物过程的富集分析。

表1 数据及样本信息统计表

2 结果

2.1 2 个数据系筛选到的DEGs 2 个数据系的DETs 结果见图1。Log2FC＜0 表示WS/WB 胸肌中基因表达量下调。表2 中为2 个数据系的差异DEGs。GSE79276数据系共筛选出104 个DEGs（Padj＜0.05），其中17个上调，87 个下调。下调基因白介素3 调节核因子基因（Nuclear Factor,Interleukin 3 Regulated，NFIL3）是P值最小的基因，下调基因半胱氨酸和甘氨酸富集蛋白质3 基因（Cysteine and Glycine-Rich Protein 3，CSRP3）是差异倍数最大的基因。GSE107362 数据系共 筛 选 出953 个DEGs（Padj＜0.05），其中179 个 上调，774 个下调。上调基因细胞因子信号传导抑制因子1（Suppressor Of Cytokine Signaling 1，SOCS1）是P值最小的基因，上调基因B 细胞表面抗原CD40（CD40 Molecule，CD40）基因是差异倍数最大的基因。

表2 2 个数据系差异基因统计表

图1 2 个数据系差异表达转录本火山图

2.2 2 个数据系筛选到的IDEGs 2 个数据系IDEG 分析见图2。2 组DEGs 共筛选出23 个IDEGs，且所有基因均下调，基因相关信息见表3。NFIL3和CSRP3分别为GSE79276 数据系差异分析的P值最小和差异倍数最大的基因。这些基因中，纤连蛋白1（Fibronectin，FN1）、VI 型胶原蛋白α3 链编码基因（Collagen Type VI Alpha 3 Chain，COL6A3）和血小板反应蛋白1 基因（Thrombospondin 1，THBS1）为筛选到的核心基因，3 个基因间存在互作关系（图3）。

表3 2 个数据系的IDEGs 统计信息

图3 部分有互作关系的IDEGs

2.3 IDEGs 的聚类和信号通路分析 23 个IDEGs 的聚类分析结果见图4a。GO 分析共筛选出6 条“P＜0.05”的条目，其中生物黏附（Biological Adhesion）和细胞黏附（Cell Adhesion）是P值最小的条目，参与2个条目的基因均为FN1、COL6A3、THBS1和血小板反应蛋白2（Thrombospondin 2，THBS2）。剩余的4 个条目中，3 个条目为信号传导相关条目，最后一个为生物调节条目。23 个IDEGs 的信号通路分析结果见图4b。信号通路分析共筛选出3 个“P＜0.05”的信号通路，其中细胞外基质受体互作通路（ECM-receptor interaction）是P值最小的信号通路，其次为黏附斑信号通路（Focal Adhesion），参与2 个信号通路的基因均为软骨寡聚基质蛋白（Cartilage Oligomeric Matrix Protein，COMP），I 型胶原蛋白α3 链编码基因（Collagen Type VI Alpha 2 Chain，COL1A2）、FN1、COL6A3、THBS2和THBS1。最后一个通路为吞噬体信号（Phagosome）通路，参与该通路的基因为COMP、THBS2、THBS1基因。

图4 IDEGs 的GO 分析和信号通路分析结果

3 讨论

3.1 差异表达分析方法 本文使用GEO2R 软件进行分析，与Zambonelli 等[4]和Malila 等[5]分析方法有所差异。生物信息学尽管有多种方法，但并没有公认的最准确的方法，只是从不同角度尽量挖掘转录组数据中的更多信息，以便进行后续的功能富集分析。利用GEO2R 软件同样的标准和方法对2 组数据进行差异表达分析，一方面利用GEO2R 挖掘转录组数据中的更多信息，另一方面尽量减少分析方法等对结果的影响，使筛选到的共同基因更为准确可靠。

3.2 差异表达基因分析结果 本研究对2 个数据系进行转录本差异表达分析，发现WS/WB 表型与正常表型胸肌肉相比，下调的DEGs 的数目均远大于上调基因的数目，说明在WS/WB 表型鸡中，胸肌部分正常基因的表达降低。韦恩图分析筛选到的基因均为下调，也进一步证明了上述内容，即营养、氧化应激和缺氧等因素导致快大型鸡胸肌中正常的基因表达受到抑制[6]，从而导致鸡胸肌不能维持正常的生理状态和功能，出现肌纤维坏死和脂质增多等表型。

3.3 IDEGs 的聚类分析 IDEGs 的GO 聚类分析结果中，筛选出细胞黏附和信号传导相关的5 个条目。信号通路分析中发现基质受体互作和黏附斑信号通路，二者结果一致。细胞外基质（ECM）由结构和功能大分子的复杂混合物组成，细胞和ECM 之间的特异性相互作用是由跨膜分子整合素介导的[7]。黏附斑信号通路是一种依赖于黏附斑激酶（Focal Adhesion Kinase 1，FAK）的信号通路，该通路利用整合素介导ECM 和细胞之间的黏附作用，并刺激几种细胞内信号传导途径的活性[8]。研究表明，黏附斑信号通路与成肌细胞分化、肌肉纤维形成和肌肉大小的发育状态相关[9]，并可以通过整合素翻译跨细胞质膜传递的细胞骨架应力和应变信号，引起肌动蛋白细胞骨架的重组[10]。基质受体互作和黏附斑信号通路可能介导了胞外异常信号进入胞内，引入肌细胞的异常代谢。

3.4 IDEGs 功能分析 23 个IDEGs 中，COMP、FN1、COL1A2、COL6A3、THBS2、THBS1基因与细胞黏附和信号传导有关。研究发现这些基因还可能参与了肌肉发育和调控。COMP 蛋白是一种在肌肉骨骼和心血管系统中都存在的细胞外基质非胶原糖蛋白，研究发现，COMP

可维持线粒体稳态并控制平滑肌细胞的表型特征[11]。FN1纤连蛋白可以结合肌细胞和肌动蛋白表面，参与细胞运动和维持细胞形态[12]。COL1A2和COL6A3均编码胶原蛋白多肽链。胶原蛋白家族成员是在肌肉发育中具有重要调控作用的ECM 蛋白。XIII 型胶原蛋白能够调节神经肌肉接头正确组装[13]。XXV 胶原蛋白对于成肌细胞融合到肌纤维中是必需的[14]。血小板反应蛋白（THBS）是参与多种生物学过程的基质细胞蛋白。THBS 家族每个成员受组织生长、愈合和重塑等诱导。研究发现，肌纤维中THBS 的表达可稳定肌膜[15]。这些信号蛋白的下调可能会影响肌细胞的正常形态和功能。

IDEGs 中，膜联蛋白A1（Annexin A1，ANXA1）、CSRP3、肌球蛋白2（Myozenin 2，MYOZ2）、ATPase胞质/ 内质网Ca2+转运2（ATPase Sarcoplasmic/Endop lasmic Reticulum Ca2+，ATP2A2）直接参与肌肉发育与功能的调节。ANXA1 蛋白参与和调节糖皮质激素介导的炎症反应，促进肌动蛋白细胞骨架的重排[16]，并介导吞噬体和肌动蛋白细胞骨架之间的相互作用[17]。CSRP3 蛋白是肌发生的正调节剂，用作肌源性bHLH（Basic Helix-Loop-Helix）转录因子，直接结合肌动蛋白丝，将肌动蛋白丝交联成束，防止丝切蛋白介导的解聚反应[18]。MYOZ2 是一种细胞内结合蛋白，在肌原纤维形成中起重要作用，参与连接α-肌动蛋白，并将钙调神经磷酸酶信号转导至肌小节，参与肌肉活动的调节[19]。ATP2A2基因编码心肌肌网Ca2+-ATP 酶，它是位于骨骼肌的肌质或内质网中的细胞内泵。该酶催化ATP 的水解，并促进钙从胞质向肌浆网腔的转运，参与收缩/松弛周期的调节[20]。这4 个基因的下调，可能直接影响肌动蛋白的形态和正常功能，导致肌肉表型异常。

IDEGs 中，肌 醇5-磷 酸 酶2（Synaptojanin 2，SYNJ2）参与不同的膜运输和信号转导途径，对网格蛋白介导的内吞具有抑制作用[21-22]，可能会导致携带胆固醇的低密度脂蛋白（LDL）大量进入细胞。载脂蛋白A1（Apolipoprotein A1，APOA1）是血浆中高密度脂蛋白（HDL）的主要蛋白成分，HDL 具有抗氧化和保护LDL 的作用，并促进胞内胆固醇和甘油三酯等脂类外排[23]。ANXA1介导胆固醇从内质网到多囊泡体（MVB）的转运，并刺激MVB 分泌具有保护和运输脂质功能的外泌体[24]。推断SYNJ2、APOA1、ANXA1表达量的下降可能引起胆固醇等脂质在肌肉细胞的积累和破坏，诱导炎性反应[25]，产生大量吞噬LDL 和氧化型低密度脂蛋白（oxLDL）的巨噬细胞，最终形成泡沫细胞并破裂，在这个过程中将脂蛋白降解为各种脂质[26]，这可能是发生炎性细胞浸润和脂质沉积的因素之一。微管相关蛋白1B 基因（Microtubule Associated Protein 1B，MAP1B）编码微管相关蛋白。该家族的蛋白质被认为与微管组装有关，微管与其他纤维一起构成细胞骨架，微管双螺旋结构支撑着细胞生理形态[27]。该基因表达量的下调可能引起肌肉细胞形态的变化。

综上所述，WS/WB 表型的产生可能是由于各种因素导致肌肉组织接受到异常信号，经整合素介导后，通过黏附斑信号通路引起胞内相关基因表达的异常，最终导致胸肌细胞功能异常，肌肉纤维和肌动蛋白的结构和功能受到破坏，引起异常表型。

4 结论

本研究利用GEO 数据的2 个数据系筛选WS/WB表型的差异表达基因，筛选到23 个共同表达基因，功能富集分析发现部分基因富集在信号传导和细胞黏附通路，参与ECM 受体互作通路和黏附斑信号通路。此外，23 个IDEGs 中，ANXA1、CSRP3、MYOZ2和ATP2A2基因直接参与肌动蛋白结构组成和功能调节，COMP和FN1等6 个基因通过调节肌肉组织信号传导和细胞黏附活动，维持肌肉正常生理功能。