小鼠急性心肌梗死后不同时间点的基因表达谱

李昊，贾孝，白雅琴，吴鹏，郭华林，贠克明，高彩荣，郭相杰

1.山西医科大学法医学院，山西 晋中 030600；2.证据科学教育部重点实验室（中国政法大学），北京 100088；3.内蒙古医科大学，内蒙古 呼和浩特010110

心肌梗死（myocardial infarction，MI）是指持续性缺血缺氧导致的心肌细胞坏死。绝大多数心肌梗死是由于冠状动脉粥样硬化引起的管腔狭窄或血栓栓塞导致，少数是由于冠状动脉痉挛引起血管壁严重收缩使管腔闭塞而致。急性心肌梗死（acute myocardial infarction，AMI）起病急骤，每年可导致900 万人死亡，是全球死亡事件发生的主要原因之一[1-3]。在法医学实践中，推断AMI 的发生时间存在困难，本研究通过分析美国国立生物技术信息中心（National Center for Biotechnology Information，NCBI）基因表达综合数据库（gene expression omnibus，GEO）公共数据中的GSE4648 数据集，分析数据集内小鼠AMI 后0.25、1、4、12、24、48 h 的差异表达基因，并对其进行功能通路富集分析、构建蛋白质-蛋白质相互作用关系（protein-protein interaction，PPI）网络、拓扑分析筛选关键基因并进行时序性聚类分析，以期为AMI 相关研究提供借鉴或参考，并为AMI 猝死的法医学鉴定提供潜在的生物标志物。

1 材料与方法

1.1 数据样本

从NCBI 的GEO 数据库（https://www.ncbi.nlm.nih.gov/geo/）下载GSE4648 数据集，该数据集是将C57BL/6J 小鼠冠状动脉左前降支进行手术结扎建立AMI 实验组，只穿线不结扎建立假手术对照组，分别于术后0.25、1、4、12、24、48 h 共6 个时间点收集左心室心肌样本，同时设立不做任何处理的空白对照组模型并收集左心室心肌样本，对上述3 组心肌组织进行基因芯片检测，得到转录本数据集[4]。

本研究从GSE4648 数据集中选取AMI 实验组和假手术对照组两组术后0.25、1、4、12、24、48 h 的左心室心肌组织样本数据，每组每个时间点6 个样本，空白对照组左心室心肌组织样本数据6 个，共78 个样本转录组数据进行分析。

1.2 方法

1.2.1 不同时间点的差异表达基因鉴定

基于芯片平台（GPL81、GPL82、GPL83）的基因注释信息使用perl语言进行转录本注释，将探针ID 转化为基因符号或者Entrez 编号（Entrez ID）。对于同一基因的多个探针，取差异表达最明显的数值进一步分析。删除测序结果中的非mRNA 探针，并对处理后的结果进行归一化处理。在R.4.0.4 环境下，使用R 软件包impute（R/Bioconductor package impute v3.15）对芯片原始数据进行背景校正和分位数标准化后，用R软件包limma（R/Bioconductor package limma v3.32.2）进行差异表达分析，将阈值设定为差异变化倍数（fold change，FC）上调或下调2 倍（|log2FC|＞1）且P＜0.05。使用R 软件包ggplot2（R/Bioconductor package ggplot2 v3.3.5）绘制6 个时间点差异表达基因的火山图。

1.2.2 不同时间点差异表达基因的功能通路富集分析

基因本体论（gene ontology，GO）富集从生物过程、分子功能和细胞组分对基因进行分类。京都基因与基因组百科全书（Kyoto encyclopedia of genes and enomes，KEGG）主要包含三类内容：基因信息数据、系统信息数据和化学信息数据。为了探寻不同时间点差异表达基因的功能和信号通路，利用R 软件包clusterProfiler、org.Hs.eg.db、org.Mm.eg.db、enrichplot、ggplot2 对不同时间点的差异表达基因进行GO 和KEGG富集分析。

1.2.3 筛选关键基因

STRING（https://cn.string-db.org/）是分析各种来源蛋白相互作用的数据库，包括已知和未知的预测关系。将各个时间点的差异表达基因导入STRING，以最小互作评分≥0.4 的蛋白构建网络图。随后利用Cytoscape v3.7.1（National Resource for Network Biology，https://cytoscape.org/）中的CytoHubba 插件分析PPI 网络的拓扑学参数，通过Degree 拓扑分析算法筛选各时间点排在前10位的关键基因。

1.2.4 基因聚类分析

对前述不同时间点的差异表达基因取交集，得到在不同时间点中均发生差异表达的交集基因，12～48 h基因数较多，为进一步分析这些基因的表达特点，使用R软件包Mfuzz（R/Bioconductor package Mfuzz v2.36.0）对该数据进行噪声抗变性软聚类分析。根据基因随时间的变化趋势将结果分为不同的模块（如持续上升、持续下降、先升后降），将结果可视化为基因聚类趋势图。

2 结果

2.1 不同时间点差异表达基因

通过R/Bioconductor package limma 差异分析共筛选出1 320个差异表达基因。与空白对照组相比，冠状动脉左前降支结扎术后各时间点差异表达基因的数目及表达情况（图1、表1）显示，随着心肌缺血时间的延长，DEG数量呈上升趋势，且表达以上调为主。

表1 结扎冠状动脉左前降支后不同时间点的差异表达基因数量Tab.1 The number of differentially expressed gene at different time points after the ligations of left anterior de‐cending coronary artery （个）

2.2 GO 功能富集分析

对不同时间点缺血心肌组织中的差异表达基因进行GO 富集，AMI 后0.25 h 仅2 个基因，未能产生富集结果。AMI 后1 h 差异表达基因富集到的生物功能包括细胞分解代谢过程的正向调节、心肌收缩的调节、跨膜转运的正向调控等；分子功能涉及转录辅阻遏物活性、阳离子通道活性及蛋白丝氨酸/苏氨酸激酶活性等；细胞组分有肌节、肌原纤维及收缩纤维等（图2A）。AMI 后4 h 差异表达基因富集到的生物功能涉及肌肉系统过程、横纹肌组织发育、心肌肥大等；分子功能涉及转录辅阻遏物活性、转录辅助调节剂活性、蛋白质自缔合等；细胞组分有肌原纤维、收缩纤维、肌动蛋白细胞骨架等（图2B）。AMI 后12 h 差异表达基因富集到的生物功能涉及磷酸化负调控、血管发育的调节、血管生成的调节等；分子功能包括肌动蛋白结合、蛋白质丝氨酸/苏氨酸激酶活性、转录辅助调节剂活性等；细胞组分有肌动蛋白细胞骨架、肌原纤维、收缩纤维等（图2C）。AMI 后24 h 差异表达基因富集到的生物功能涉及细胞因子产生正向调节、白细胞迁移、炎症反应的调节等；分子功能包括肌动蛋白结合、受体配体活性、细胞黏附分子结合等；细胞组分有肌动蛋白细胞骨架、膜区、膜微区等（图2D）。AMI 后48 h 差异表达基因富集到的生物功能包括创伤反应、白细胞迁移、炎症反应的调节等；分子功能包括肌动蛋白结合、酶抑制剂活性、受体配体活性等；细胞组分有含胶原的细胞外基质、膜区、肌动蛋白细胞骨架等（图2E）。

图2 不同时间点差异表达基因的GO 富集结果Fig.2 The GO enrichment analysis of differentially expressed genes at different time points

2.3 KEGG 通路富集分析

KEGG信号通路分析结果显示，AMI后0.25h仅2个差异基因未得到富集结果。从图3 可以看到，差异表达基因在AMI 1 h 后主要参与了Apelin 信号通路、丝裂原活化蛋白激酶（mitogen-activated protein kinase，MAPK）信号通路、心肌收缩、心肌细胞中的肾上腺素能信号传导等；AMI 后4 h 主要富集于癌症中的蛋白聚糖、钙信号通路、ErbB 信号通路等；AMI 后12 h 主要富集于MAPK 信号通路、PI3K-Akt 信号通路、流体剪切应力与动脉粥样硬化等；AMI 后24 h 主要参与细胞因子-细胞因子受体相互作用、MAPK 信号通路、趋化因子信号通路等；AMI 后48 h 主要参与沙门氏菌感染、吞噬体、肌动蛋白细胞骨架的调节等（图3E）。

图3 不同时间点差异表达基因的KEGG 富集结果Fig.3 The KEGG enrichment analysis of differentially expressed genes at different time points

2.4 PPI 网络和拓扑算法鉴定关键基因

通过Degree 拓扑分析算法筛选到各时间点排在前10 位关键基因的共表达及互作网络（图4），AMI 后0.25 h 基因数低于10 个，未产生结果。基因之间的连线代表其相互作用关系，基因的Degree 值越大，颜色越深，说明在网络中发挥的作用越重要。

图4 关键基因的共表达及互作网络Fig.4 The co-expression and interaction network of key genes

2.5 差异表达基因的时序性聚类分析

对AMI实验组各时间点差异表达基因取交集获得多时间点交集差异表达基因，将1～48 h 的基因交集进行可视化（图5），并分析其中存在3个时间点及以上交集基因的时序性变化特点，发现AMI后6个时间点中没有基因持续发生差异表达。根据基因发生差异时间点的早晚和数量，在AMI后0.25、1、4、12 h肌球蛋白轻链7（myosin light chain 7，MYL7）基因差异表达持续下调；AMI后1、4、12、24、48 h TSC22域家族成员2（TSC22 domain family member 2，TSC22D2）基因差异表达持续上调；AMI 后1、4、12、24 h 诱导型热激蛋白1A（heat shock protein 1A，HSPA1A）基因差异表达持续上调；AMI后1、4、12 h B 细胞异位基因2（B-cell translocation gene 2，BTG2）、孤儿核受体亚家族4A1（nuclear receptor subfamily 4 group A member 1，NR4A1）基因差异表达持续上调。兰尼碱受体2（ryanodine receptor 2，RYR2）仅在1、4 h表达上调，但该基因在2个时间点的关键基因共表达及互作网络中均发挥了重要作用（表2）。

表2 TSC22D2、NR4A1、MYL7、HSPA1A、BTG2 在AMI后的时序性表达Tab.2 The post-AMI sequential expression of TSC22D2，NR4A1，MYL7，HSPA1A and BTG2

图5 不同时间点差异表达基因交集韦恩图Fig.5 The Venn diagram of the intersection of differentially expressed genes at different time points

此外，在AMI 后12、24、48 h 这3 个时间点，共有82 个基因存在差异表达，对其进行时序性基因聚类分析以探索这些基因的表达变化，结果表明，82 个基因可分为3 个聚类，聚类结果随时间的改变呈现不同的变化趋势（图6）。聚类1 有33 个基因表达值随时间变化呈先上升后下降的趋势，为先上调后下调基因；聚类2 有26 个基因表达值随时间变化整体呈上升趋势，为上调基因；聚类3 有23 个基因表达值随时间变化整体呈下降趋势，为下调基因。

图6 AMI后12、24、48 h 差异表达基因的时序性表达Fig.6 The sequential expression of post-AMI of different expressed genes at 12 h，24 h and 48 h

3 讨论

在法医学鉴定实践中，关于心脏性猝死（sudden cardiac death，SCD），法医病理学家多参考美国心脏学会和美国心脏病学会的定义[5-6]“出现症状1 h 内发生的死亡，或指先前无任何症状但可能是由于心脏原因导致24 h 内发生的死亡”。冠状动脉粥样硬化性心脏病（coronary heart disease，CAD）引起的猝死约占SCD病因的80%[7]，CAD 引发的猝死死亡过程较短，缺乏特异的急性心肌梗死的病理学及形态学改变，尽管MI 是CAD 中的重要改变，但MI 一般在心肌受损6 h之后才能肉眼辨认。本研究对AMI 后6 个时间点心肌组织的差异表达基因进行生物信息学分析，锁定AMI 发生发展过程中的关键基因，以期为AMI 猝死的法医学诊断提供分子标志物，为法医学研究提供借鉴和参考。

本研究通过对AMI 的1 320 个差异表达基因进行GO 和KEGG 富集分析，发现随着AMI 后时间延长，差异表达基因数量整体呈上升趋势，且表达以上调为主。GO 富集结果显示，随着AMI 后的时间延长，不同时间点差异表达基因沿时序所参与的心肌收缩的调节、心肌肥大、血管发育的调节、炎症反应的调节等生物功能，符合AMI 后的病理生理过程。KEGG 结果显示，在AMI 后1～48 h 的5 个时间点均富集到MAPK 信号通路。MAPK 级联反应在MI 后可以促进心肌细胞生长，抑制心肌细胞间的缝隙连接，促进心肌纤维化进而推动病理性心肌纤维化的发展[8-9]，这提示MAPK信号通路在AMI 的发生发展过程中发挥关键作用。

本研究分析的多个AMI 关键基因中，MYL7、TSC22D2、HSPA1A、BTG2、NR4A1、RYR2的时序性差异表达有较明显的特异性，可能在AMI 的发生发展中发挥重要作用。SURESH 等[10]将发生AMI 的31 名患者与正常人的外周血进行比较，发现MYL7、TSC22D2、BTG2、NR4A1、RYR2均存在异常表达。MYL7参与肌动蛋白细胞骨架合成，调节心脏发育和心肌收缩[11]，已有学者通过监测AMI 患者MYL7的同系物MYL1评估心肌梗死组织面积，被认为是极具潜力的AMI 诊断标志物[12-13]。如能在AMI 猝死者中证实MYL7的特异性，则可以通过单克隆固相酶免疫测定和放射免疫测定将检测MYL7应用于AMI 猝死的法医学实践。TSC22D2隶属于TSC22家族，是转化生长因子-β（transforming growth factor-β，TGF-β）的主要效应因子。研究[14-15]结果表明，TSC22通过促进TGF-β 信号传导引起多种纤维化因子（如α-SMA、血浆纤溶酶原激活物抑制剂-1、胶原蛋白Ⅰ和纤连蛋白）的过表达，促进MI 后心肌纤维化。HSPA1A 是一种自我保护蛋白，过表达可以增强损伤心肌的耐受性，抑制心肌纤维化、改善心肌细胞坏死性改变[16]。有研究[17]结果表明，在AMI 后早期（第2 天）和晚期（第7 天）的炎症反应阶段，异常表达的HSPA1A 可能参与炎症反应，本研究结果与其相同。BTG2调节细胞周期、参与应激反应[18]，有研究[19]结果显示，BTG2在MI 后表达上调，本研究结果与其一致。BTG2可能是MI 关键调节因子。NR4A1可以有效抑制心肌肥大、血管狭窄，在MI 后早期可以参与单核细胞的炎症反应，晚期可能介导巨噬细胞的损伤修复反应[20-21]。RYR2的过度活化已被证实与恶性心律不齐、心肌收缩功能障碍和心脏中的不良结构重塑密切相关[22]。

综上所述，与AMI 关系密切的关键基因有MYL7、TSC22D2、HSPA1A、BTG2、NR4A1、RYR2，在AMI 发生后呈现一定的时序性变化，可为AMI 猝死的法医学鉴定及相关研究提供借鉴和参考。