基于双硫死亡的心肌缺血再灌注损伤相关风险基因筛选及分析

【摘 要】目的：筛选基于双硫死亡的心肌缺血再灌注损伤（myocardial ischemia reperfusion injury，MIRI）相关风险基因及探索其可能的作用途径。方法：用R语言limma分析并筛选MIRI相关风险基因24 h组与Control组之间的差异表达基因。基于GO和KEGG数据库利用 R语言cluterProfiler包对其进行功能富集分析。根据双硫死亡基因集以及GSE160516的表达数据，利用R语言GSVA包中的ssgsea法对其进行富集分析，并利用Venn图获得心肌缺血再灌注中双硫死亡相关基因的表达。用R语言psych包中的corr.test计算双硫死亡相关基因与凋亡因子、线粒体、铁死亡及炎症等相关基因的相关性，并绘制热图。利用R语言Mfuzz包来进行不同时间点心肌缺血再灌注转录组数据的表达模式聚类分析并利用Venn图获取时序性差异表达的双硫死亡基因。结果：本研究共筛选出17个差异表达双硫死亡相关基因，经过双硫死亡相关基因与炎症、凋亡、铁死亡、线粒体基因相关性分析，以及双硫死亡相关基因与心肌缺血再灌注时序性分析的差异表达基因筛选，最终得到双硫死亡相关基因Flna、Myl6、Tln1 在MIRI不同时间的特异性表达。结论：Flna、Myl6、Tln1 在MIRI中可能起着关键作用，其基于双硫死亡同时与线粒体相关基因密切相关，可作为MIRI患者风险因素的筛选指标。

【关键词】心肌缺血再灌注；双硫死亡；生信分析；风险因素筛选指标

【中图分类号】D919.1；DF795.1 【文献标志码】A 【收稿日期】2023-12-06

心肌梗死（myocardial infarct，MI）是心源性猝死的主要类型，心肌急性缺血缺氧是其重要的发生机制。急性心肌梗死得到及时治疗后，血管损伤反而呈现短暂加重的现象，被称为心肌缺血再灌注损伤（myocardial ischemia reperfusion injury，MIRI）。

目前，关于心肌缺血再灌注的病理生理学机制方面的研究已有不少。有研究表明，当心肌血管梗阻再通时，可以通过加重血管炎症导致这种现象的发生[1]。线粒体动力学被认为在MIRI扮演着中心环节的作用，线粒体功能障碍可能是导致梗阻血管再通后血管损伤加重的潜在机制[2]。线粒体在心肌梗死后促进炎症消退和心肌组织修复，当其出现功能障碍时，梗死组织炎症加重并抑制其修复[3]。此外，氧化应激也被认为参与心肌缺血再灌注的病理过程[4]。虽然已经明确了心肌缺血后血管再通会导致更严重的心肌及血管损伤，但是其具体的机制仍然有待研究。

低血糖是心血管疾病的危险因素之一，常常导致心律失常、心肌梗死甚至死亡[5]。最新的研究发现，在葡萄糖缺乏的条件下，SLC7A11的高表达会导致细胞的双硫死亡[6]。双硫死亡（disulfidptosis）是一种新的细胞死亡方式，它是由细胞葡萄糖缺乏诱导的SLC7A11高表达的细胞内二硫化物的异常积累导致的细胞死亡[6]。心肌梗死后，常常因为胰岛素分泌不足而导致低血糖。在葡萄糖缺乏的情况下引起细胞双硫死亡导致心肌细胞进一步损伤，这可能是MIRI的机制之一，目前很少相关研究。因此，筛选双硫死亡相关基因（disulfidptosis-relatedgenes，DRG）的表达对于诊断和治疗MIRI 非常重要。本文旨在筛选MIRI特征基因并分析其生物学功能，探索DRG与心肌梗死后不同时间的特征基因的关联性，提供潜在的MIRI风险因素筛选指标，为MIRI的临床防治提供参考依据。

1 资料与方法

1.1 数据获取及差异表达分析

从GEO 数据库（https：//www.ncbi.nlm.nih.gov/gds/）下载小鼠MIRI数据。数据集GSE160516。根据文献中最佳灌注时间，使用R包limma分析并筛选IR24h组与Control组之间的差异表达基因，筛选标准为：FCgt;1.5，Plt;0.05[7]。

1.2 功能分析

根据得到的差异基因，利用cluterProfiler包进行功能富集分析。基于Gene Ontology数据库，以及KEGG PATHWAYDATABASE生物化学通路的数据库，对候选基因进行了功能富集分析[8-9]。利用统计学算法（Fisher’s exact test）寻找出一组基因和哪些具体的功能条目联系最大，分析结果中每个条目都对应一个统计值P 来表示显著性。

1.3 ssGSEA分析。

根据双硫死亡基因集以及GSE160516的表达数据，利用GSVA包中的ssgsea法进行富集分析，最终得到各样本对应各基因集的ssgsea得分[10]。利用wilcox.test计算组间各基因集得分差异情况，多组间差异则利用kruskal.test计算各基因集得分差异情况。

1.4 双硫死亡相关基因的相关性分析

根据双硫死亡相关基因在GSE160516中的表达量，用psych中的corr.test计算双硫死亡相关基因与凋亡因子、线粒体、铁死亡及炎症等相关基因的相关性，并绘制热图[11]。凋亡因子、线粒体、铁死亡及炎症等相关基因筛选条件：|logFC|gt;1。

1.5 双硫死亡相关基因时序分析

对心肌缺血再灌注样本的表达矩阵以6、24、72 h的时间顺序研究表达差异。利用Mfuzz包来进行不同时间点转录组数据的表达模式聚类分析。参数选择α=0.5，基因筛选的标准偏差为0.3。

1.6 统计学方法

本研究数据处理和分析均使用R 语言进行。不同的分析用相应的R包进行。功能富集分析用Fisher’s exact test对对应条目进行注释。在ssGSEA 分析中，利用R 包中的wilcox.test和kruskal.test比较组间差异。检验水准α=0.05。

2 结果

2.1 差异分析

从GSE160516数据集中筛选Control与IR 24 h之间的差异表达基因，结果表明有2 291个基因表达上调，2 251个基因表达下调（图1A）。其中，Fn1、Chil3、ll1rn、serpina3n、Ccr2、Ccl12 等20个基因上调最明显；Ciart、Pkd212、Dusp18、Kcnip2、Lgals4、lgi1 等20 个基因表达水平下降幅度最大（图1B）。

2.2 功能富集分析

为了确定差异表达基因的生物学功能，对差异表达基因（不分上调与下调）进行了GO/KEGG 功能富集分析。在生物学过程层面主要包括Cellular macromolecule metabolic pro⁃cess、Multicellular organismdevelopment、Organelle organiza⁃tion、Positive regulation of cellular metab...、Protein metabolicprocess 5类功能变化，并且展示其功能变化对应的基因名称（图2A）。在分子功能的分析中，发现其主要富集于purinenucleotide binding、purine ribonucleotide binding、ribonucleo⁃tide binding、adenyl nucleotide binding、adenyl ribonucleotidebinding以及ATP binding处（图2B）。在细胞组成的分析中，分析出Intracellular membrane − bounded organelle、Nucleo⁃plasm、Nucleus等细胞定位功能变化（图2C）。表明MIRI后，Biological Process（BP）、Molecular Function（MF）、CellularComponent（CC）3类基因功能均发生变化。为了探索这些功能变化所属Gene Ontology数据库数据集，将所有表达上、下调的基因进行了富集分析，发现富集数量最多的是GO：0043231、GO：0005829基因注释条目，其代表的细胞组成分别是Intracellular membrane−bounded organelle、Cytosol。功能变化程度最明显的是GO：0005634（nucleus）、GO：0031981（nuclear_lumen）、GO：0046872（metal_ion_binding） 3 个条目（图2D）。

为了确定差异表达基因对应的生物学通路，利用cluter⁃Profiler包基于KEGG PATHWAY DATABASE数据库分析这些基因的生物学通路。分析发现，差异表达基因主要集中在Proteoglycans in cancer、Chemokine signaling pathway通路，而上调最显著的基因Spp1、Thc 对应Focal adhesion 通路（图3A），还分析MAPK signaling pathway、Metabolic pathways、Pathways in cancer、PI3K − Akt signaling pathway、Proteogly⁃cans in cancer 5个条目对应得的上、下调表达基因（图3B）。与上面的分析显示出类似的结果，MIRI后差异表达基因对应的生物学通路条目热图显示，上调程度最明显的Spp1、Thc基因都对应PI3K−Akt信号通路、Focal adhesion通路，下调程度最显著的Zbtb16 基因对应Transcriptional misregulation incancer、Pathways in cancer通路（图3C）。

2.3 ssGSEA分析及双硫死亡基因筛选

为了确定双硫死亡基因与心肌缺血再灌注之间的内在联系，再根据双硫死亡基因集以及GSE160516的表达数据，利用GSVA 包中的ssgsea 法对每个样本进行打分，确定在MIRI不同时间后双硫死亡相关基因的表达活性。随着血管恢复灌流的时间推移，双硫死亡相关基因的活性也逐步下降（图4A）。随后，将差异表达基因与双硫死亡基因集取交集，得到17个基因（图4B）。

2.4 双硫死亡基因表达相关性分析

为了探索双硫死亡导致损伤的可能机制，根据双硫死亡相关基因在GSE160516中的表达量，用psych中的corr.test计算双硫死亡相关基因与凋亡因子、线粒体、铁死亡及炎症等相关基因的相关性并绘制热图。筛选出的17个双硫死亡基因均与差异表达的凋亡因子存在不同程度的相关性，其中jam3 基因与凋亡因子相关性均较弱（图5A）。另外，除Myh10 基因与大部分铁死亡基因相关性较弱外，其余双硫死亡相关基因与铁死亡基因相关性较强（图5B）。在炎症相关基因中，双硫死亡基因与ccl 基因家族呈显著的正相关或负相关关系，少量双硫死亡基因与其他炎症相关基因相关，例如：Gys1、Myh9、Oxsm、Slc3a2（图5C）。随后分析了双硫死亡相关基因与线粒体相关基因的相关性，jam3、Slc7a11 与大部分线粒体相关基因无相关性，其他基因与线粒体相关基因相关性较强（图5D）。与凋亡因子、铁死亡及炎症等相关基因相比，线粒体相关基因与双硫死亡相关基因的相关性明显更强。

2.5 双硫死亡基因时序分析

为了确定在MIRI后不同时间段均有表达的双硫死亡基因，利用Mfuzz包对心肌缺血再灌注样本的表达矩阵以6 h、24 h、72 h的时间顺序研究基因表达差异（图6A）。将聚类得到的基因与差异表达双硫死亡基因取交集，筛选出3个基因（Flna、Myl6、Tln1）（图6B）。

3 讨论

心肌梗死发病率高，常引起严重后果。有研究表明，心肌梗死再通后血液再灌注会造成更严重的细胞损伤，被称为MIRI。研究MIRI的致病机制，可为患者的治疗及预后提供参考依据和治疗靶点。鉴于心肌缺血常因胰岛素分泌不足而导致低血糖、葡萄糖获取不足，葡萄糖代谢障碍，此种情况下引起的细胞双硫死亡可能进一步导致心肌细胞损伤。因此，本研究设计本生信分析实验，旨在筛选基于双硫死亡的MIRI的相关风险基因，为MIRI的预防和治疗提供参考。

本研究基于Gene Ontology 数据库的功能分析显示，心肌缺血再灌注后，细胞大分子代谢过程、多细胞生物的发育、细胞器组织、 细胞代谢的正向调节、蛋白质代谢过程等生物学过程变化明显，包括分子功能以及细胞定位均发生明显的变化，可能是导致血液再通后更严重损伤的原因之一。基于KEGG PATHWAY DATABASE 数据库分析差异表达基因的生物学通路，发现高度上调差异表达基因可以激活对应的Focal adhesion通路，抑制该通路可减轻心肌纤维化并保护心肌功能[12]。MAPK信号通路的激活也可能是双硫死亡导致MIRI的原因，有研究表明抑制该信号通路可以预防MIRI[13]。然而，也有研究表明MAPK家族成员MAPK3的上调会减轻心肌缺血再灌注导致的细胞凋亡，换句话说，MAPK3 在MIRI 后起保护心肌的作用[14]。因此，双硫死亡在MIRI中的作用及其机制仍有待深入研究。

趋化因子信号通路在癌症等多种疾病中发挥作用已得到验证，研究也发现，CXCL16的激活会导致更严重的缺血再灌注损伤[15]。同时，趋化因子通常通过加重机体的炎症而导致病情进一步加重。此外，富集在PI3K−Akt 信号通路上的差异表达基因也可以在MIRI 后起到保护心肌的作用，Maciel等[16]的研究也验证这点。在缺氧/复氧的细胞模型中，抑制PI3K−Akt信号通路降低中药对心肌细胞的保护作用[17]。除此之外，MIRI 还上调癌症发病途径、癌症蛋白聚糖信号通路的表达，这可能导致患者癌症的风险增高，但其在心肌损伤中的作用机制如何仍有待深入研究。

有研究表明，可通过PI3K−Akt 信号通路和MAPK信号通路抑制GSK3表达水平[18]，GSK3在多种细胞和组织中发挥作用，包括心脏。GSK3的抑制通过增加糖原合成和减少糖酵解在心肌组织中起保护作用[19]。GSK3是一种高度保守的细胞内丝氨酸/苏氨酸激酶，通过磷酸化抑制糖原合成酶并调节葡萄糖代谢[20]。另外，PI3K−Akt 信号通路过度激活导致糖酵解显著加速，阻止丙酮酸进入线粒体，抑制线粒体生物发生和功能，导致能量代谢紊乱[21]。因此在本研究的分析中，PI3K−Akt 信号通路和MAPK信号通路的富集可能是心肌缺血再灌注后导致葡萄糖代谢紊乱进一步引起心肌细胞双硫死亡的机制之一。

随后，本研究在差异表达基因中筛选出17个双硫死亡相关基因，表明MIRI过程中存在引起细胞双硫死亡基因激活，从而导致MIRI。为了探索MIRI基于双硫死亡可能的机制，本研究分析了双硫死亡相关基因与凋亡因子、铁死亡相关基因、炎症因子以及线粒体相关基因之间的相关性。有研究发现，在癌症中Jam3 基因的抑制可以促进肾癌细胞的迁移并抑制其凋亡[22]。但本研究表明，Jam3 基因与大部分凋亡因子及线粒体相关基因无明显相关性，提示Jam3 基因表达通过细胞凋亡及改变线粒体功能导致MIRI可能性小。文献表明，Myh10 基因与大部分铁死亡及炎症相关基因无明显相关性，相关研究寥寥无几。最新的研究表明，Myh10 基因参与受损线粒体自噬的诱导[23]；Myh10 基因在MIRI中的具体机制，仍有待深入研究阐明。

综上，差异表达双硫死亡相关基因与凋亡因子、铁死亡相关基因、炎症因子以及线粒体相关基因之间的相关性强，但其囊括范围较大，难以精确定位MIRI风险因素筛选指标。因此，分析MIRI后6 h、24 h、72 h 3个时段均显著表达的基因，将它们与差异表达双硫死亡基因取交集，得到Flna，Myl6，Tln13个基因。Tln是细胞骨架蛋白，可以将整合素和肌节联系起来[24]。整合素的活化在内皮细胞表型转化中起重要作用，整合素亲和力的调节在体内和体外导致动脉粥样硬化内皮活化[25]。此外，Tln在心脏成纤维细胞中的流失引起心肌压力超负荷导致心肌细胞肥大[26]。Tln1 还能调节局灶性黏附动力学、细胞迁移和侵袭，因此，Tln1在肿瘤中表达并促进肿瘤进展[27]。同样，Flna是一种肌动蛋白结合蛋白，参与细胞骨架的形成，将多种蛋白质锚定在细胞骨架中并调节细胞黏附和迁移[28]，因此在肿瘤和心血管疾病中发挥作用。Flna 基因的突变导致先天的心血管、脑和肺的疾病[29]，在MIRI中的作用有待研究阐明。

本研究仅对基于双硫死亡的MIRI相关风险基因在生信分析水平上进行探索，其具体的临床实用价值还需要动物实验以及临床试验予以验证和确认。

总之，本研究筛选出了3个在MIRI表达明显的基因，它们与细胞双硫死亡密切相关。时序性分析提示，双硫死亡可能贯穿MIRI的全程，其特异表达的基因可能成为MIRI 的风险因素筛选指标，为MIRI的临床防治提供新的参考和借鉴。

参考文献

[1] Algoet M，Janssens S，Himmelreich U，et al. Myocardial ischemiareperfusioninjury and the influence of inflammation[J]. Trends Cardio⁃vasc Med，2023，33（6）：357-366．

[2] Zhou ML，Yu YF，Luo XX，et al. Myocardial ischemia-reperfusioninjury：therapeutics from a mitochondria-centric perspective[J]. Cardiol⁃ogy，2021，146（6）：781-792．

[3] Cai SS，Zhao MY，Zhou B，et al. Mitochondrial dysfunction in mac⁃rophages promotes inflammation and suppresses repair after myocardialinfarction[J]. J Clin Invest，2023，133（4）：e159498．

[4] Xiang M，Lu YD，Xin LY，et al. Role of Oxidative Stress in Reper⁃fusion following Myocardial Ischemia and Its Treatments[J]. Oxid MedCell Longev，2021，2021：6614009．

[5] Freeman J. Management of hypoglycemia in older adults with type2 diabetes[J]. Postgrad Med，2019，131（4）：241-250．

[6] Liu XG，Nie LT，Zhang YL，et al. Actin cytoskeleton vulnerabilityto disulfide stress mediates disulfidptosis[J]. Nat Cell Biol，2023，25（3）：404-414．

[7] Ritchie ME，Phipson B，Wu D，et al. Limma powers differential ex⁃pression analyses for RNA-sequencing and microarray studies[J].Nucleic Acids Res，2015，43（7）：e47．

[8] Ashburner M，Ball CA，Blake JA，et al. Gene ontology：tool for theunification of biology. The Gene Ontology Consortium[J]. Nat Genet，2000，25（1）：25-29．

[9] Kanehisa M，Goto S. KEGG：Kyoto encyclopedia of genes and ge⁃nomes[J]. Nucleic Acids Res，2000，28（1）：27-30．

[10] Hänzelmann S，Castelo R，Guinney J. GSVA：gene set variationanalysis for microarray and RNA-seq data[J]. BMC Bioinformatics，2013，14：7．

[11] Lancaster PA，Davis ME，Tedeschi LO，et al. A mathematical nu⁃trition model adequately predicts beef and dairy cow intake and biologi⁃cal efficiency[J]. Transl Anim Sci，2022，6（1）：txab230．

[12] Valiente-Alandi I，Potter SJ，Salvador AM，et al. Inhibiting fibro⁃nectin attenuates fibrosis and improves cardiac function in a model ofheart failure[J]. Circulation，2018，138（12）：1236-1252．

[13] Pan J，Chen SY，Guo WZ，et al. Alpinetin protects against hepaticischemia/reperfusion injury in mice by inhibiting the NF-κB/MAPK sig⁃naling pathways[J]. Int Immunopharmacol，2021，95：107527．

[14] Wang H，Pang W，Xu X，et al. Cryptotanshinone Attenuates Isch⁃emia/Reperfusion- induced Apoptosis in Myocardium by UpregulatingMAPK3[J]. J Cardiovasc Pharmacol，2021，77（3）： 370-377.

[15] Shen B，Li J，Gao L，et al. Role of CC-chemokine receptor 5 onmyocardial ischemia-reperfusion injury in rats[J]. Mol Cell Biochem，2013，378（1/2）：137-144．

[16] Maciel L，de Oliveira DF，Mesquita F，et al. New cardiomyokinereduces myocardial ischemia/reperfusion injury by PI3K-AKT pathwayvia a putative KDEL-receptor binding[J]. J Am Heart Assoc，2021，10（1）：e019685．

[17] Zhang ZB，Qin XT，Wang ZH，et al. Oxymatrine pretreatment pro⁃tects H9c2 cardiomyocytes from hypoxia/reoxygenation injury by modu⁃lating the PI3K/Akt pathway[J]. Exp Ther Med，2021，21（6）：556．

[18] Omar MA，Wang LG，Clanachan AS. Cardioprotection by GSK-3inhibition：role of enhanced glycogen synthesis and attenuation of cal⁃cium overload[J]. Cardiovasc Res，2010，86（3）：478-486．

[19] Frame S，Cohen P. GSK3 takes centre stage more than 20 yearsafter its discovery[J]. Biochem J，2001，359（Pt 1）：1-16．

[20] Embi N，Rylatt DB，Cohen P. Glycogen synthase kinase-3 fromrabbit skeletal muscle. Separation from cyclic-AMP-dependent proteinkinase and phosphorylase kinase[J]. Eur J Biochem，1980，107（2）：519-527．

[21] Vaupel P，Multhoff G. Revisiting the Warburg effect：historicaldogma versus current understanding[J]. J Physiol，2021，599（6）：1745-1757．

[22] Li XD，Yin AP，Zhang WJ，et al. Jam3 promotes migration andsuppresses apoptosis of renal carcinoma cell lines[J]. Int J Mol Med，2018，42（5）：2923-2929．

[23] Jun YW，Lee S，Ban BK，et al. Non-muscle MYH10/myosin IIBrecruits ESCRT-Ⅲ to participate in autophagosome closure to maintainneuronal homeostasis[J]. Autophagy，2023，19（7）：2045-2061．

[24] Manso AM，Okada H，Sakamoto FM，et al. Loss of mouse cardio⁃myocyte talin-1 and talin-2 leads to β-1 integrin reduction，costamericinstability，and dilated cardiomyopathy[J]. Proc Natl Acad Sci U S A，2017，114（30）：E6250-E6259．

[25] Al-Yafeai Z，Pearson BH，Peretik JM，et al. Integrin affinitymodulation critically regulates atherogenic endothelial activation in vitroand in vivo[J]. Matrix Biol，2021，96：87-103．

[26] Azizi L，Cowell AR，Mykuliak VV，et al. Cancer associated talinpoint mutations disorganise cell adhesion and migration[J]. Sci Rep，2021，11（1）：347．

[27] Noll NA，Riley LA，Moore CS，et al. Loss of talin in cardiac fibro⁃blasts results in augmented ventricular cardiomyocyte hypertrophy in re⁃sponse to pressure overload[J]. Am J Physiol Heart Circ Physiol，2022，322（5）：H857-H866．

[28] Zhou J，Kang XM，An HX，et al. The function and pathogenicmechanism of filamin A[J]. Gene，2021，784：145575．

[29] Meliota G，Vairo U，Ficarella R，et al. Cardiovascular，brain，andlung involvement in a newborn with a novel FLNA mutation：a case re⁃port and literature review[J]. Adv Neonatal Care，2022，22（2）：125-131．

（责任编辑：周一青）

基金项目：国家自然科学基金资助项目（编号：82060340）；贵州省科技计划资助项目[编号：黔科合基础-ZK（2023）一般330]；贵州省普通高等学校青年科技人才成长资助项目[编号：黔教合KY字〔2021〕157）]；山东省证据鉴识重点实验室开放课题资助项目（编号：2018KFKT6）；国家级大学生创新创业训练资助项目（编号：202210660025、202110660004）。