基于GEO数据库分析糖尿病心肌病的差异表达基因

2022-09-05陈嘉敏吴会会苏国海

中国医学科学院学报 2022年4期

陈嘉敏，李莹，，吴会会，刘鹏，郑燕，，苏国海

1山东大学齐鲁医学院济南市中心医院专科转化研究中心，济南 250013

2山东第一医科大学附属中心医院专科转化研究中心，济南 250013

糖尿病心肌病是指在没有冠状动脉疾病、高血压和心脏瓣膜病的情况下发生心力衰竭的病理生理疾病[1]。欧洲及美国的多项回顾性研究显示糖尿病增加了心力衰竭患者的住院率和死亡风险[2]，在我国，糖尿病导致的心力衰竭发病率在1993至2007 年从12.3%增加到22.1%，糖尿病是除冠心病、高血压外导致心力衰竭最主要的原因[3]。2017至2018年中国心力衰竭住院患者的横断面研究显示心力衰竭住院患者合并糖尿病的比例高达29.2%[4]，在包括中国人在内的多种族的研究显示，糖尿病增加了心力衰竭患者的发病率和死亡率[5]，近30年来我国糖尿病的患病人数逐年上涨[6]，使糖尿病心肌病的防治工作面临更加严峻的挑战。研究表明单纯将糖尿病患者的血糖恢复到基线水平并不能降低心力衰竭的发病率和死亡率[7]，但是新型降糖药钠葡萄糖协同转运蛋白2抑制剂如达格列净可以使心力衰竭合并糖尿病患者的心力衰竭恶化风险及因心血管事件死亡的风险降低[8]，心血管药物沙库巴曲-缬沙坦也可使糖尿病合并心力衰竭的患者获益[9]。多种机制参与了糖尿病心肌病的进展：(1)心肌细胞胰岛素信号通路受损：正常情况下，磷酸肌醇3激酶-蛋白激酶B信号通路激活，刺激葡萄糖转运体4(glucose transporter type 4，GLUT4)向质膜募集，促进葡萄糖进入心脏细胞。胰岛素信号通路受损时，GLUT4表达下降，心脏摄取葡萄糖的能力下降。(2)高血糖和糖毒性：高血糖可导致糖基化终末产物增加，糖基化终末产物是长寿命蛋白被糖化的产物，糖基化终末产物沉积使结缔组织交联、纤维化，导致心脏顺应性降低、舒张功能障碍。(3)脂肪酸利用增加和脂毒性：糖尿病患者葡萄糖利用下降，参与脂肪酸β氧化的信号转导分子过氧化物酶体增殖物激活受体(peroxisome proliferator-activated receptor，PPAR)α表达增加促进脂肪酸的摄取和利用，脂质和脂质代谢产物的增加导致胰岛素抵抗和心脏纤维化。(4)氧化应激：脂肪酸的β氧化促进活性氧的产生和心脏的氧化应激。(5)Ca2+处理受损：高血糖状态下，内质网应激使Ca2+处理异常，导致心肌细胞兴奋性收缩偶联障碍，心脏收缩能力下降[1，10- 12]。本研究旨在通过对现有数据库数据的再分析，寻找导致糖尿病心肌病的关键基因，对糖尿病心肌病的发病机制有更深一步的研究。

材料和方法

GEO数据库选取基因集进行差异基因分析从GEO数据库中选取与糖尿病心肌病研究相关的基因芯片。纳入标准：(1)种属：Rattus norvegicus；(2)同时具备正常心脏组织和糖尿病心肌病心脏组织。经过检索GSE4745[13](https：//www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE4745)和GSE5606[14](https：//www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE5606)两个基因集符合条件。GSE4745基因集的芯片来自GPL85(Affymetrix Rat Genome U34 Array)平台，数据集中包括正常组雄性大鼠心室组织和注射链脲佐菌素后3、28、42 d后糖尿病组雄性大鼠心室组织共24 份，注射链脲佐菌素42 d后的糖尿病组大鼠心室组织以及其对照组样本共8 份用于差异表达基因的分析。GSE5606基因集的芯片来自GPL1355(Affymetrix Rat Genome 230 2.0 Array)平台，数据集中包括注射链脲佐菌素16周的雄性大鼠正常组和糖尿病组心室组织样本共14 份，全部用于差异基因的分析。GEO数据库的在线分析工具GEO2R用于差异基因的分析，以P≤0.05，|Log|差异倍数(fold change，FC)|≥1为标准进行筛选，在R[15](3.6.3版本)中利用ggplot2[16]包(3.3.2版本)绘制差异基因的火山图，在线韦恩图绘制工具(http：//bioinformatics.psb.ugent.be/webtools/Venn/)绘制韦恩图，寻找两个数据集中共同表达的上调和下调差异基因。

富集分析将差异表达基因利用R中的AnnotationHub包[17]、org.Rn.eg.db包[18]进行基因ID的转换，clusterProfile包[19]进行基因本体论(gene ontology，GO)和京都基因和基因组百科全书(Kyoto encyclopedia of genes and genome，KEGG)富集分析，ggplot2绘图，矫正后P≤0.05为有统计学意义。基因集富集分析(gene set enrichment analysis，GSEA)是一种计算方法，用于确定一组预先确定的基因集在两种生物状态之间是否显示出差异有统计学意义，GSEA软件[20- 21](版本4.1.0，https：//www.gsea-msigdb.org/gsea/datasets.jsp)用于基因集富集分析，依据GSEA官网提供的文件格式[Data formats-GeneSetEnrichmentAnalysisWiki (broadinstitute.org)]构建文件进行GSEA，基因表达数据被分为正常和糖尿心肌病(diabetic cardiomyopathy，DCM)两种生物状态，PPAR信号通路基因集和脂肪酸代谢基因集[GSEA|MSigDB(gsea-msigdb.org)]作为参考基因集用于富集分析，名义P<0.05、错误发现率即q<0.25、标准化富集评分绝对值>1被认为有统计学意义。

蛋白质互作网络构建利用STRING(https：//string-db.org)进行蛋白质互作网络的构建，将差异表达基因导入STRING数据库，选择Rattus norvegicus进行蛋白互作网络的构建，将结果以TSV格式导出。将互作网络中节点出现的次数进行统计并排序，将STRING导出的TSV文件以“Network”方式导入Cytoscape，将排序文件以“Table”方式导入，将节点大小及颜色以“Mapping Type”方式依据节点出现次数多少进行颜色填充和节点大小的调整。利用插件Cytohubba的最大集团中心算法计算出排名前10的关键基因，以得分排序并调整节点大小绘图。

原代心肌细胞的提取取出生3 d内的Wistar大鼠乳鼠，75%的酒精将乳鼠短暂浸泡消毒，眼科剪剪开乳鼠前胸，镊子取出心脏置于含有1%青链霉素的Hank’s平衡盐溶液(货号：CC106，中科迈晨科技有限公司)中，剪碎心脏后在含有200 U/ml 胶原酶(货号：LS004176，美国Worthington)的Hank’s平衡盐溶液中进行消化，预消化5 min后弃去胶原酶，重新添加胶原酶后消化1 h，取上清加入等量的心肌细胞培养基终止消化并离心取沉淀，余下组织再次进行消化，重复以上步骤直至无组织团块，将离心取得的沉淀用培养基吹起后置于培养瓶中差速贴壁80 min后除去成纤维细胞，心肌细胞种于0.3%明胶(货号：9000708，默克公司)包被后的6孔板中。

心肌细胞培养原代心肌细胞贴壁18 h后换为含有100 μmol/L的5-溴脱氧尿嘧啶核苷(货号：S7918，Selleck)的心肌细胞培养基中以抑制成纤维细胞的生长，心肌细胞培养基为8%的马血清(货号：16050122，Gibco)、5%的胎牛血清(货号：10270106，Gibco)、1%的青链霉素和64%的低糖5.5 mmol/L DMEM的低糖培养基(货号：CM1004，中科迈晨科技有限公司)构成。原代心肌细胞在含有Brdu的低糖培养基培养72 h后正常组更换为不含Brdu的低糖培养基，实验组更换为不含Brdu的25 mmol/L高糖培养基(货号：C11995500BT，Gibco)，72 h后提取RNA。

实时荧光定量PCR从6孔板收集贴壁心肌细胞，用Trizol(货号：AG21102，艾科瑞生物工程有限公司)提取总RNA，使用反转录试剂盒(货号：AG11711，艾科瑞生物工程有限公司)反转录后利用SYBERGeeen qPCR试剂盒(货号：AG11718，艾科瑞生物工程有限公司)进行实时荧光定量PCR，每次实验设置3 个复孔，实验重复3 次。引物利用美国国家生物技术信息中心的Primer-BLAST工具设计，FASTA序列从美国国家生物技术信息中心的基因数据库获取，产物长度限制在100～300 bp，溶解温度设定为57～63 ℃，以β-actin作为内参基因，引物序列见表1。

统计学处理PCR以2-△△Ct法计算目的基因表达水平，利用GraphPad(8.0.1)进行统计分析，采用独立样本t检验，P≤0.05为差异有统计学意义。

表1 关键基因的引物序列

结果

GEO2R与韦恩图分析结果在GEO数据库中选择数据集GSE4745和GSE5606作为研究对象，利用GEO2R进行分析后，以|LogFC|≥1，P≤0.05为阈值筛选差异基因并利用ggplot2绘制火山图，其中LogFC≥1，P≤0.05为上调基因，LogFC≤-1，P≤0.05为下调基因。在GSE4745中以GSM107428、GSM107429、GSM107430、GSM107431作为对照组，GSM107432、GSM107433、GSM107434、GSM107435作为糖尿病心肌病组，经GEO2R分析、ggplot2绘图后显示上调基因164个、下调基因164个(图1A)。在GSE5606中以GSM130860、GSM130861、GSM130862、GSM130863、GSM130864、GSM130865、GSM130866作为对照组，以GSM130867、GSM130868、GSM130869、GSM130870、GSM130871、GSM130872、GSM130873作为糖尿病心肌病组，经GEO2R分析、ggplot2绘图显示上调基因286个，下调基因135个(图1B)。在线韦恩图绘制工具分析后显示共同表达的上调差异基因共35个(图1C)，共同表达的下调差异基因共15个(图1D)。GSE5606和GSE4745基因集中共同表达的差异基因的LogFC值和P值显示，GSE5606和GSE4745基因集中共同表达的上调基因共35个，共同表达的下调基因共15个(表2)。

基因富集分析结果在R中利用clusterProfile包，将50个差异基因进行GO和KEGG富集，矫正后P≤0.05为有统计学意义。这50个差异基因在GO富集中的分子功能、生物学过程有明显富集，而在细胞组分无富集。在生物学过程中，这50个差异基因在脂肪酸代谢过程、小分子分解过程、嘌呤核苷酸代谢过程、脂质修饰等生物过程中明显富集(图2A)。在分子功能中，差异基因主要在棕榈酰辅酶A水解酶活性、酰基辅酶A水解酶活性、辅酶A水解酶活性、脂肪酸衍生物的结合、硫酯水解酶活性这些条目中明显富集(图2B)。在KEGG，富集主要集中在PPAR信号通路、脂肪酸的延伸、不饱和脂肪酸的生物合成(图2C)。对脂肪酸代谢基因集和PPAR信号通路基因集进行GSEA富集显示，糖尿病组在脂肪酸代谢和PPAR信号通路的基因集中明显富集(图3)。

FC：差异倍数

STRING在线数据库蛋白互作网络分析结果利用STRING在线数据库(https：//www.string-db.org/)将GSE4745和GSE5606数据集共同表达的差异基因构建蛋白互作网络，利用Cytoscape绘制出蛋白质的互作关系网络(图4A)。利用插件Cytohubba的最大集团中心算法计算出排名前10的关键基因并绘制关键基因的互作网络图(图4B)。等级评分见表3。

关键基因在原代大鼠心肌细胞中的验证结果原代心肌细胞对照组在5.5 mmol/L的葡萄糖低糖环境、实验组在25 mmol/L的葡萄糖高糖环境中培养72 h后，提取RNA进行实时荧光定量PCR。t检验结果显示，Pdk4(t=23.520，P<0.001)、Hmgcs2(t=10.700，P=0.004)、Ucp3(t=5.184，P=0.035)、Acsl6(t=6.213，P=0.003)、Slc2a4(t=19.300，P=0.005)具有统计学意义且与表2的分析结果相符，其中Pdk4、Hmgcs2、Ucp3在高糖刺激后表达增加，Acsl6、Slc2a4在高糖刺激后表达降低(图5)。而Decr1(t=17.180，P<0.001)、Acot1(t=7.186，P=0.002)、Acot2(t=34.990，P<0.001)、Slc27a1(t=5.508，P=0.005)、Cpt1a(t=30.870，P<0.001)虽具有统计学意义但与表2基因芯片分析结果并不相符(图5)。

讨论

通过对GEO数据库中的两个糖尿病心肌病数据集进行分析后，共筛选出50个共同表达差异基因，对差异基因进行GO富集后显示，这些差异基因在脂肪酸代谢和脂肪酸修饰等生物功能中明显富集，KEGG富集分析提示差异基因在脂肪酸延伸和不饱和脂肪酸生物合成通路中明显富集。底物代谢紊乱是糖尿病心肌病的发病机制之一，在生理情况下，心脏可使用脂肪酸和葡萄糖作为主要的代谢底物，胰岛素抵抗会导致葡萄糖摄取减少而游离脂肪酸的摄取和氧化增加[22]，基因富集的分析结果提示这些差异基因可能参与了糖尿病心肌病的底物代谢紊乱。在差异基因的KEGG富集分析中，PPAR信号通路在富集第1位，PPAR在心脏糖脂代谢以及能量稳态中发挥关键作用，PPAR有多种亚型，其中PPARα在心脏中高表达，影响脂肪酸的摄取和线粒体脂肪酸的氧化[23]。PPARα不仅参与脂肪酸和葡萄糖的代谢，还参与支链氨基酸、酮体在糖尿病心肌病中的代谢[24]。差异基因在PPAR信号通路中的明显富集表明这些差异基因可能通过PPAR信号通路参与脂肪酸的代谢，进而参与心脏的底物代谢紊乱。GSEA富集分析提示糖尿病组在脂肪酸代谢和PPAR信号通路明显富集。

表2 基因集中差异基因的表达情况

GO：基因本体论；KEGG：京都基因和基因组百科全书；BP：生物学过程；MF：分子功能；PPAR：过氧化物酶体增殖物激活受体

表3 MCC算法计算出的排名前10的关键基因的分值

实时荧光定量PCR的验证结果显示，在高糖刺激下，Pdk4、Ucp3、Hmgcs2基因表达增加，Slc2a4、Acsl6基因表达下降。Pdk4是细胞能量代谢的关键激酶[25]，Pdk4特异性过表达的转基因小鼠表现出胰岛素抵抗和心肌葡萄糖氧化降低[26]，Ucp3是线粒体解偶联蛋白家族成员，Ucp3与脂肪酸的β氧化直接相关[27]。在一项利用胰岛β细胞遗传缺陷的糖尿病心肌病的模型研究中，高糖低胰岛素状态会导致Pdk4和Ucp3的表达增加[28]，这与本研究结果一致，Pdk4和Ucp3在糖尿病心肌病细胞模型中表达升高，但是目前的研究尚无针对Pdk4和Ucp3在糖尿病心肌病中的功能研究。Hmgcs2是酮生成过程中的限速酶[29]，与本研究一致，在大鼠糖尿病模型中，心脏Hmgcs2表达明显增加而肝脏中Hmgcs2无变化[30]。Slc2a4又被称作Glut4，Glut4是葡萄糖转运体最主要的亚型，占葡萄糖转运体的70%，在糖尿病心肌病中Glut4表达下降[31]。但是Glut4特异性过表达的糖尿病心肌病小鼠，非但不能挽救糖尿病心肌病的心脏代谢障碍，反而可能导致糖尿病心肌病的进一步发展[32]，Glut4在糖尿病心肌病中的作用有赖于进一步的研究。

Acsl6在糖尿病心肌病中尚无研究，在人类和大鼠骨骼肌中Acsl6调节脂质合成和线粒体的氧化能力[33]，Acsl6与小鼠大脑ω- 3脂肪酸DHA的富集相关[34]，也是小鼠神经保护所必需的[35]。在Acsl的各种亚型中，Acsl3、Acsl4敲减的细胞系中，葡萄糖刺激胰岛素释放降低了50%，但是这项研究中Acsl6对胰岛素的释放并未产生影响[36]。在糖尿病心肌病的细胞模型中Acsl6基因表达降低，在本研究中是首次报道，Acsl6基因在糖尿病心肌病模型中的功能可作为糖尿病心肌病的一个研究方向。

ES：富集分数；CTR：对照组；DCM：糖尿病心肌病组；NES：标准化富集分数；FDR：错误发现率

与现有的基于基因集GSE4745和GSE5606进行差异基因分析的文章相比，Dai等[37]通过对基因集GSE5606进行加权基因共表达网络分析显示，Angptl4、Acot1、Decr1、Hmgcs2、Pdk4在糖尿病心肌病的小鼠心脏中表达增加，这些基因同样是本研究中的关键差异表达基因，但是与Dai等[37]的研究结果不同，Decr1和Acot1基因在本研究中表达下降，考虑为实验选择动物种属不同。另一项关于GSE4745基因集的研究仅对GSE4745基因集的差异表达基因进行了分析，未在动物或细胞水平检测差异基因的表达水平[38]。与之前的研究相比，本研究同时使用了两个基因集的数据，关键差异表达基因与Dai等[37]的研究虽有重合但仍有新的发现，Acsl6基因不仅在现有的基因芯片的分析中报道较少，在糖尿病心肌病的研究中报道也较少，且与Acsl家族中的Acsl3和Acsl4胰岛素的释放相关，这为以后的研究提供了很好的理论依据。另外，在现有的与GSE4745和GSE5606基因集研究中仅对差异基因进行GO和KEGG富集不同，本研究同时进行了基因芯片整体基因代谢通路的研究，是对GO和KEGG仅对差异基因进行富集分析的补充。