高糖诱导阴茎海绵体内皮细胞损伤关键基因的筛选和生物信息学分析*

马 帅 陈 颖 齐 涛 李世雄 杨文涛 陈 俊**

1.广西壮族自治区南宁市广西中医药大学研究生学院（广西南宁 530200）

2.广东省广州市中山大学附属第三医院不育与性医学科（广东广州 510631）

3.广西壮族自治区南宁市广西中医药大学附属瑞康医院男性科（广西南宁 530012）

糖尿病性勃起功能障碍（diabetic erectile dysfunction,DED）是糖尿病最主要的并发症之一，是难治性ED，严重影响着患者及其配偶的生活质量和心理健康[1]。尽管目前ED的治疗方法有很多，如口服PDE-5抑制剂、阴茎海绵体药物注射以及假体植入等，但治疗DED的效果均欠佳，究其原因在于DED的发病机制不清[2,3]。因此，研究DED的发病机制，探索新的、有效的DED治疗方法迫在眉睫。

随着生物信息学技术的不断发展和应用，对基因芯片数据的挖掘和利用已成为生物医学的研究热点，为DED的基础研究提供了新的思路[4,5]。研究表明，DED的发生和发展与阴茎海绵体内皮细胞受损密切相关[6]。本研究拟对模拟血糖紊乱病理状况下的小鼠阴茎海绵体内皮细胞的转录谱RNA-seq测序数据进行分析，构建差异基因的蛋白质相互作用网络，进行GO功能富集分析，筛选核心差异基因，以期为DED治疗开发新的药物靶点提供思路和理论基础。

材料方法

一、芯片数据的选择与分析

在GEO数据库中检索编号为GSE146078的芯片表达谱为研究对象。该芯片是模拟持续高糖环境导致阴茎内皮细胞损伤，正常（5 mmol/LD-葡萄糖）对照组和高糖（30 mmol/LD-葡萄糖）处理72小时组的小鼠阴茎海绵体内皮细胞的转录谱RNAseq测序数据[7]。采用检测平台GPL17021，使用Illumina HiSeq 2500(Mus musculus)进行测序，利用R语言工具包edgeR和DEseq2对表达矩阵进行归一化分析，并进行差异分析，|log FC|＞0.80和P＜0.05。

二、蛋白质相互作用网络构建

通过STRING（https://string-db.org/）将筛选到的差异基因进行蛋白质与蛋白质相互作用网络验证。限制物种为小鼠，所使用的参数阈值为中度可信数据：0.400，使用的相互作用类型包括7种（Textmining，Experiments，Databases，Co-expression，Neighborhood，Gene Fusion，Co-occurrence）。

三、差异基因的GO和KEGG富集分析

使用R软件的clusterProfiler包和reactome（https://reactome.org/）对筛选的差异基因进行GO富集分析和通路富集分析；使用DAVID数据库（https://david.ncifcrf.gov）对差异基因进行GO和KEGG通路富集分析。设置P<0.05，分析结果以气泡图展示。

四、核心差异基因分析

利用Cytoscape的Hub插件分析核心基因，限定物种为小鼠，获取蛋白相互作用关系，并导入Cytoscape软件绘制相互作用网络。利用cytoHubba插件分析节点的大小和颜色表示Degree值的大小，节点越大，颜色越红对应的Degree值越大。

结 果

一、阴茎海绵体内皮细胞高糖处理损伤芯片分析

GSE146078芯片中初始数据共检测覆盖了23282个基因，其中20296个基因表达矩阵满足初筛要求（总计数counts>1，即至少在其中一个样本有表达）。经过R语言工具包edgeR和DEseq2对表达矩阵进行归一化分析，并进行差异分析，最终得到了17139个基因的表达数据。其中总的差异基因为367个，包含有150个上调以及217个下调基因（图1）。

图1 小鼠阴茎海绵体内皮细胞高糖处理损伤芯片分析的火山图（A）和热图（B）

二、高糖处理后内皮细胞差异基因通路富集分析

通过STRING在线分析工具对差异基因进行蛋白质与蛋白质相互作用网络验证，最终得到的相互作用网中包含结点数据332个、边界数505条、平均结点度为3.04、平均聚类系数0.415、期待的边界数据为243（图2）。

图2 差异基因的蛋白质相互作用网络分析

在此基础上对差异基因进行GO富集分析，结果显示差异基因的分子功能主要集中在细胞因子活性、趋化因子活性、受体信号转导、受体调控活性等过程（图3A）。细胞组分富集分析发现其主要在细胞外胞质空间、质膜蛋白复合物（图3B）。生物学功能通路富集可以发现差异表达基因在炎症反应通路极显著的富集，如急性时相反应、急性炎症应答、白细胞迁移、白细胞趋化以及促炎症因子白介素-1的应答等。同时也有一些代谢调控通路发生变化，如活性氧的代谢过程、氮氧化物的合成过程等（图3C）。进一步进行Reactome相互作用组的富集分析，则发现这些差异表达基因主要富集在DNA损伤过程，如DNA碎片化激活因子、凋亡诱导的DNA碎片化等（图3D）。

图3 差异基因的GO富集分析

KEGG通路富集分析结果如图4所示，上调基因主要富集在细胞色素P450对异生物质（xenobiotics）的代谢、类固醇激素的合成、谷胱甘肽代谢、军团菌病、药物代谢-细胞色素P450、血清素能突触、沙门氏菌感染、类风湿关节炎、药物代谢、IL-17信号通路以及化学致癌作用（图4）。下调基因富集结果显示其与Hippo通路相关（图5）。

图4 上调表达差异基因的KEGG通路富集分析

图5 下调表达差异基因的KEGG通路富集分析

三、高糖处理后内皮细胞核心差异基因分析

利用Cytoscape的Hub插件分析发现，高糖处理小鼠阴茎海绵体内皮细胞差异基因中有两个关键基因群，图6左侧为炎症应答为代表的免疫调控网络（代表的细胞因子有上调表达基因Ccl19、Cxcl1、IL-6和下调表达基因Cxcl12；代表的转录因子有下调表达基因Sox9等），右侧基因群为结合DNA的组蛋白家族可能在DNA损伤凋亡以及片段化时期发生的变化（代表有上调表达基因Hist1h4d）（图6）。我们在差异基因中发现多个基因家族特征性表达，如锌指蛋白Zscan4（Zinc finger and SCAN domain containing 4）家族的Zscan4b、Zscan4c、Zscan4d；谷 胱 甘 肽S-转 移 酶GSTA（Glutathione S-Transferase Alpha）家族Gsta1、Gsta3以及Gsta4；细胞色素P450家族成员CYP2D包括Cyp2d10、Cyp2d11、Cyp2d12、Cyp2d13、Cyp2d34、Cyp2d37-ps、Cyp2d40以及Cyp2d9；血清淀粉样蛋白A（Serum amyloid A,Saa）家 族Saa1、Saa2、Saa3，且Saa1、Saa2、Saa3均显著上调表达（图7）。核心差异基因中，Cxcl12、IL-6等已被证明参与ED进展[8,9]，Zscan4家族、GSTA家族及P450家族成员在内皮细胞中的功能作用均未见报道，而Saa家族成员被发现与内皮损伤、炎症反应等相关疾病过程密切相关[10,11]，但其在DED中的作用目前尚未见报道。因此，本研究将对Saa家族成员与DED的关联进行进一步分析。

图6 高糖处理后内皮细胞通路核心差异基因分析

图7 高糖处理后内皮细胞中

四、Saa家族相关蛋白质分析

Saa1-3均定位于细胞外，属于分泌型蛋白质[12]。我们经过STRING（https://string-db.org/）搜索得到Saa家族相关蛋白质相互作用网络，发现其与多个载脂蛋白或者脂相关蛋白有关联。其中包括载脂蛋白Lcn2（Lipocalin-2）、血清类粘蛋白Orm2（Orosomucoid 2）、激素原转化酶1抑制剂Pcsk1n（Proprotein convertase subtilisin/kexin type 1 inhibitor）、氧磷脂酶1Pon1（Paraoxonase 1）及载脂蛋白（Apolipoprotein）家族多个成员Apoa1，Apoa2，Apoc3，Apoe（图8）。

图8 STRING分析Saa家族相关蛋白质相互作用网络

五、Saa家族功能富集分析

GO分析结果显示Saa家族基因的分子功能主要集中在包括脂酶活性，热休克蛋白结合，脂蛋白颗粒结合，磷脂结合等过程（9A）。其生物学功能主要富集于甘油三脂的代谢，胆固醇代谢，胆固醇内流，脂解反应，脂蛋白代谢等过程（9B）。其细胞组分主要集中在极低密度脂蛋白，低密度脂蛋白，高密度脂蛋白，细胞外胞质等（9C）。KEGG分析发现其主要参与PPAR信号通路以及胆固醇代谢（9D）。

讨 论

ED在糖尿病患者中的发病率是普通人群的3~4倍，在诊断为糖尿病后的前10年，超过一半的男性确诊为DED，为患者和其家庭带来巨大痛苦[13]。然而，由于DED发病机制尚未阐明，临床上缺乏针对该病的特效药物以及有效的治疗靶点。近年来，随着微阵列和测序技术的快速发展，基因芯片数据分析已被广泛应用于研究多种疾病的发生发展以及靶基因的筛选等。如，Chen等[14]分析3个包含子宫内膜异位症组织和正常子宫内膜组织的微阵列数据集，发现上皮-间充质转化在子宫内膜异位症中具有重要作用；Xue等[15]生物信息学分析发现CDC45、GINS2、MCM2和PCNA可能与宫颈癌患者预后相关，可作为宫颈癌患者的潜在治疗靶点。因此，筛选与DED相关的靶标基因，对于开发治疗DED的靶标药物具有重要意义。

图9 Saa家族功能富集分析

研究表明，DED的发生和发展与阴茎海绵体内皮细胞受损密切相关[6]。因此，本研究通过模拟血糖紊乱病理状况致使阴茎内皮细胞损伤的芯片数据，利用GPL17021检测平台筛选出正常组小鼠和高糖条件处理组小鼠之间的差异基因分析，共发现367个差异表达基因，对这些差异基因进行GO功能富集分析，结果显示这些差异表达基因主要涉及急性炎症答应、活性氧化的代谢、DNA损伤等生物学过程，与已报道的炎症因子及活性氧参与DED的结果一致。如炎性细胞因子IL-6和IL-1β含量降低可以有效改善糖尿病大鼠勃起功能[16]；血浆IL-8（sIL-8）水平与男性不育和ED有密切联系[17]；抑制大鼠海绵体中活性氧（ROS）的产生可以缓解ED进展[18]。KEGG分析结果显示，上调基因富集通路在细胞色素P450对异生物质（xenobiotics）的代谢、类固醇激素的合成、类风湿关节炎、药物代谢和IL-17信号通路等。其中IL-17信号通路与内皮细胞病理功能相关，如IL-17影响肺动脉内皮细胞的增殖，血管新生以及细胞黏附[19]；其它如类固醇激素对内皮细胞功能也有维持作用[20]。下调基因与Hippo通路相关，Hippo通路的关键因子如YAP能够调控内皮细胞激活并参与血管炎症反应[21]。

在高糖诱导阴茎海绵体内皮细胞损伤相关基因分析中发现，其核心差异网络中包括多个细胞因子如Ccl19、Cxcl1、Cxcl12、IL-6等；转录因子Sox9等；及Zscan4家族、GSTA家族、Saa家族成员等。据报道，基因工程改造后表达Cxcl12的间充质干细胞改善链脲佐菌素诱导的糖尿病大鼠ED症状[8]；抑制IL-6可减轻神经保留性根治性前列腺切除术大鼠模型的勃起功能障碍[9]。本研究分析发现，Saa家族成员Saa1-3在高糖诱导的小鼠海绵体内皮细胞中均显著上调表达，提示Saa家族可能在DED过程中具有重要作用。

Saa即血清淀粉样蛋白A（serum amyloid A），是炎症性载脂蛋白家族，定位于细胞外，是一种分泌型蛋白质，属于急性期反应过程中血清蛋白的降解与调控基因群，可引起全身系统性炎症[12,22]。大量研究证明，Saa家族在内皮损伤、炎症反应等相关疾病过程具有重要作用，如，Saa能够增强人颈动脉内皮细胞的迁移、增殖，并增加内皮细胞管腔形成[10]；IL-1β可以诱导人类冠状动脉内皮细胞中Saa的表达，最终加速冠状动脉粥样硬化进程[11]等。然而，其与DED的关联尚未见报道。本研究进一步分析发现，发现其与Lcn2、Orm2、Pcsk1n、Pon1及Apo家族成员等多个载脂或者脂相关蛋白质有关联。GO和KEGG功能富集分析结果也显示其主要集中于低密度脂蛋白、高密度脂蛋白、细胞外胞质等，参与脂蛋白颗粒结合、脂解反应，胆固醇代谢、脂蛋白代谢等过程。据报道，脂质代谢的异常与ED存在密切关联，脂质代谢混乱和氧化应激都是DED的独立危险因素[23,24]。部分研究表明，Pon1活性的降低可能与ED的致病性有关，并且由于抗动脉粥样硬化酶Pon1的活性降低，患者的动脉粥样硬化发展可能更快[25]；Apoe缺失会引起小鼠海绵体钙化，最终导致小鼠勃起功能障碍的发生[26]。而本研究分析结果显示Saa家族与脂相关蛋白存在相互作用，提示Saa家族可能通过载脂或者脂相关蛋白质参与DED的发生和发展过程。

综上所述，本研究通过生物信息学筛选发现正常组小鼠和高糖条件处理组小鼠之间的差异基因基因主要涉及急性炎症答应、活性氧化的代谢过程、DNA损伤过程等生物学过程。核心差异基因中Saa1-3在高糖处理后小鼠阴茎海绵体内皮细胞中显著上调表达，可能是DED过程中的重要基因，Saa1-3与多个载脂蛋白或者脂相关蛋白有关联，并涉及脂解反应，胆固醇代谢、脂蛋白代谢等过程。这对我们揭示DED发病机制以及药物靶点研究提供了重要的理论参考。然而，Saa家族对DED的发生和发展过程是否具有重要的作用，有待进一步的功能验证。