神经源性异位骨化相关差异基因的筛选及生物信息学分析

陈梓杰，陈柏行，许隆，张严，冯俊铭，李世杰，高怡加，曾展鹏

神经源性异位骨化（neurogenic heterotopic ossification，NHO）是指中枢神经系统损伤后板层骨在肌肉、肌腱、韧带等非骨骼组织中的异位生长，常见于创伤性颅脑损伤和脊髓损伤。异常的骨组织常引起疼痛，关节活动受限等并发症，严重影响患者的生活质量[1]。非甾体类抗炎药和双膦酸盐可在一定程度上预防NHO，然而一旦异位骨形成后，药物干预的收效甚微，手术切除是唯一的治疗方式，但术后高复发率成为困扰医生的难题[2]。目前，NHO确切的病理生理机制仍未明确，了解其发病机制对预防和治疗NHO具有重要的意义。

Chalmers[3]最早指出发生异位骨化需要满足3个条件：存在骨诱导因子、成骨前体细胞和允许成骨的环境。NHO的形成可能通过损伤的神经组织分泌诱导因子与软组织中成骨前体细胞相互作用，使细胞趋化、增殖，并成骨分化，最终激活异位成骨。实验发现脑损伤后的血清促进骨髓间充质干细胞（bone marrow mesenchymal stem cells，BMSCs）的增殖，并增强细胞内血清碱性磷酸酶的表达[4]，提示BMSCs可能是损伤中枢神经释放骨诱导因子的作用部位，但目前对于BMSCs参与NHO形成的分子机制研究仍较少。基因微阵列技术已广泛应用于人类疾病研究，可用于寻找疾病相关基因、治疗靶点和生物标志物[5]。本文试通过分析基因表达数据库（Gene Expression Omnibus，GEO）中NHO基因芯片，寻找关键基因及信号通道，为深入了解发病机制和发现新的治疗靶点提供参考。

1 材料与方法

1.1 芯片数据的获取与预处理

下载GEO中基因表达谱芯片系列GSE94683。该系列包含16个样本，分成NHO组（n=7）及健康对照组（n=9）。所有样本提取的来源均为BMSCs。我们通过R软 件（v3.6.1）的preprocess Core软 件 包 中 的normalize.quantiles（）函数对原始数据集进行分位数归一化，得到的标准化表达矩阵用Log2转换计算表达值。对原始数据进行注释，当存在多个探针对应同一个基因名时对其合并取均值。

1.2 差异基因的筛选

使用R软件的Limma包筛选差异表达基因（Differentially expressed genes，DEG）。采用Benjamini-Hochberg法进行多重比较和控制假阳性率，以调整后的P＜0.05和差异倍数|Log 2(fold-change)|＞2作为阈值，并绘制火山图及热图。

1.3 基因本体论与KEGG通路富集分析

基因本体论（Gene Ontology，GO）分析和京都基因和基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）分析是进行基因功能注释和通路富集分析的常用方法。本研究将筛选得到的DEG提交到Metascape在线工具（http://metascape.org/gp/index.html#/main/step1）进行分析，确定这些DEG的关键生物学过程（biological process，BP）、分子功能（molecular function，MF），细胞成分（cellular component，CC）和KEGG通路后，运用R软件生成散点图。

1.4 蛋白质-蛋白质相互作用网络和hub基因识别

通过Metascape数据库中的蛋白质-蛋白质相互作用信息（protein-protein interaction，PPI）分析功能，可识别NHO组和对照组之间的核心基因。用来构建PPI网络的数据主要来源于BioGrid[6]，InWeb_IM[7]，OmniPath[8]数据库。将Metascape导出的数据利用Cytoscape（V3.8.0）软件构建PPI网络，并以节点连接度为标准排序PPI网络中的hub基因。利用分子复杂检测（MCODE）算法对Cytoscape中的PPI网络的中的典型模块进行筛选。

2 结果

2.1 差异表达分析

将芯片数据标准化后，共识别出276个DEG，包括150个上调基因和126个下调基因，见图1。对前25个上调基因和25个下调基因进行可视化，蓝色表示基因表达上调，白色表示表达无差异，粉色表示表达下调，结果显示差异基因表达水平在2组间具有较好的组间一致性，见图2。

2.2 差异基因的富集分析

Metascape工具用于DEG的功能和途径富集分析，图中点的大小代表该术语/通路富集的差异基因数；点的颜色代表调整后的P值，表示通路富集的显著性程度；横坐标GeneRatio为对应术语/通路的差异基因占该术语/通路包含所有基因的比例，见图3。GO富集分析包含BP、MF与CC 3个方面，BP主要在化学突触传递的调节（modulation of chemical synaptic transmission）、跨突触信号的调控（regulation of trans-synaptic signaling）和骨化（ossification）等过程富集，见图3A。MF主要富集在细胞外基质结构成分（extracellular matrix structural constituent）、肝素结合（heparin binding）和糖胺聚糖结合（glycosaminoglycan binding），见图3 B。DEG的CC主要富集在胶原三聚体（collagen trimer）、含胶原的细胞外基质（collagen-containing extracellular matrix）和细胞外基质（extracellular matrix），见图3 C。KEGG通路分析表明，DEG主要在补体与凝血级联（complement and coagulation cascades）、IL-17信号通路（IL-17 signal pathway）、Wnt信号通路（Wnt signaling pathway）、蛋白质消化吸收（protein digestion and absorption）等通路富集，见图3 D。

图1表达谱差异分析的火山图

图2差异度大小排名前25的上、下调基因热图

为了进一步探索术语之间的关系，我们选择了富集的术语作为节点，构建术语富集网络：图4（A）展示了不同术语之间的相关性，越接近则说明相似性越高，颜色表示不同的聚类；图4（B）根据P值大小对节点进行着色，P值越小，表示该术语的富集程度越高。

图5展示了各通路与相关差异基因的网络图，黄色箭头代表富集的通路，节点代表在该通路的差异基因。

2.3 PPI网络构建与关键基因筛选

通过Metascape在线工具对DEG进行分析寻找蛋白质之间的相互作用，并使用Cytoscape软件对PPI网络进行可视化，得到1个包含76个节点，96个连接的网络，见图6。通过Cytohubba工具筛选出连接度最高的前10个hub基因，分别为EGF、GPR18、ADRA2C、CXCL1、C3、CXCL6、ADRB2、PENK、CXCL2、CDH1，见图7。Metascape网站内置MCODE算法，根据网络拓扑结构分析对基因网络进行聚类分析，共筛选出4个模块，见图6。根据Metascape节点中的邻接度及密度，每个模块包含1个种子节点作为该模块的关键基因[9]，分别为：CXCL2、PTH2、EphB1、HTR2B。从hub基因和模块基因比对中可以发现，连接度前10的hub基因中包含了所有模块1的基因，说明模块1可能在参与疾病调节起到关键作用。

3 讨论

神经源性异位骨化多发生在中枢神经损伤的患者，早期发病隐匿，药物及物理疗法的疗效不理想，手术干预后复发率高，且目前对NHO确切发病机制知之甚少。

本研究通过对基因芯片GSE94683进行生物信息学分析，筛选出150个上调基因和126个下调基因，GO富集分析表明NHO参与的分子机制主要富集在化学突触传递的调节、跨突触信号的调控、骨化、细胞外基质结构成分、肝素结合、糖胺聚糖结合、胶原三聚体、含胶原的细胞外基质、细胞外基质。富集的KEGG通路有：补体与凝血级联、IL-17信号通路、Wnt信号通路、蛋白质消化吸收。构建PPI网络筛选10个hub基因：EGF、GPR18、ADRA2C、CXCL1、C3、CXCL6、ADRB2、PENK、CXCL2、CDH1。同时筛选出4个具有生物学相关性的模块，对应种子基因分别为CXCL2、PTH2、EPHB1、HTR2B。

表皮生长因子（Epidermal growth factor，EGF）被鉴定为PPI网络中连接度最高的hub基因。EGF是一种生长因子，也是一种有丝分裂原，与表皮生长因子受体（Epidermal growth factor receptor，EGFR）结合，激活受体酪氨酸激酶活性，促进信号级联，参与细胞增殖、有丝分裂。Laflamme等[10]证明了EGF参与了促成骨细胞分化的早期阶段，然而Liu等[11]研究显示EGF本身不能诱导未分化MSCs的成骨分化，但可增强BMP9介导的MSCs的异位成骨能力。同时，Lee等[12]也发现用BMP-2和EGF处理的BMSCs，成骨标志物ALP、BSP、OPN和OST的表达比单纯BMP-2处理组明显增加，提示EGF在与BMP协同作用下可诱导骨矿化。而BMP信号通路被公认为异位骨化的共同信号通路，在一种脊髓损伤的NHO动物模型中，BMP2、BMP4、BMP7和BMP9的表达显著增强[13]。这些证据指向了EGF在NHO中可能协同BMP调节促进MSCs的成骨分化与矿化。

图3富集分析散点图

图5 KEGG富集通路-差异基因网络

图6 PPI蛋白互作网络

图7连接度最高的前10个hub基因

CXC基序趋化因子配体2（CXC Motif Chemokine Ligand 2，CXCL2）为模块1的种子基因，CXCL2是一种细胞因子，主要通过趋化中性粒细胞来调节免疫功能[14]。研究显示CXCL2在破骨细胞生成过程中也发挥了作用，Ha等[15]的研究证明CXCL2可增强破骨前体细胞的迁移，促进RANKL介导破骨前体细胞向破骨细胞分化。骨重塑依赖于骨形成与骨吸收结合的骨稳态调控平衡，抑制破骨细胞分化可减少骨吸收，间接增强骨形成，这与本研究中CXCL2表达下调的结果相一致。我们认为CXCL2在NHO形成中参与的机制可能是通过抑制MSCs旁分泌从而减弱RANKL介导的破骨细胞迁移和分化。

作为模块2的种子基因，甲状旁腺激素2（Parathyroid Hormone 2，PTH2）是与甲状旁腺激素具有编码序列相似性的肽激素前体，通过与甲状旁腺激素2型受体（Parathyroid hormone type 2 receptor，PTH2R）结合，参与G蛋白偶联受体相关途径。Oki等[16]发现PTH2干预可明显增强兔股骨干骺端内间隙的骨形成。模块2中又同时存在PTH1R，然而PTH1R可被PTHrP激活却不能被PTH2激活。Zhang等[17]的实验证明了通过激活PTH1R增强SOX9的表达，促进MSCs的软骨分化。而本研究中PTH2与PTH1R均上调，推测这2个基因可能通过协同促进MSCs分化为软骨表型，并以软骨内成骨方式产生异位骨。

模块3的种子基因是促红细胞生成素产生肝细胞受体B1（erythropoieti-producing hepatocellular receptor B1，EphB1），是受体酪氨酸激酶的最大亚族Eph受体家族的一员，与跨膜蛋白Ephrins配体以细胞-细胞接触介导的方式传递双向信号，Eph受体下游的信号传导途径称为正向信号传导，而Ephrins配体下游的信号传导途径为反向信号传导[18]。体外实验观察到颅骨成骨前体细胞中EphB2和EphB1的表达，并且EphB1基因敲除的小鼠颅骨骨量减少[19]。与此同时，刘军等[20]发现EphrinB2调节BMSCs的OPG/RANKL比率上调，抑制了对骨微环境中破骨细胞的分化，减少了骨吸收。这些发现提示NHO发生的调控机制可能为EpHB1激活的双向信号介导成骨-破骨耦联，促进了MSCs的成骨分化，同时抑制破骨细胞分化。

模块4的种子基因5-羟色胺受体2B（5-Hydroxytryptamine receptor 2B，HTR2B）为5-HT的受体，5-HT是一种神经递质，除了在大脑中发挥重要的作用外，5HT-HTR2B通路也参与了骨形成过程。据报道，HTR2B激活可以调节胚胎小鼠后肢MSCs的增殖[21]；同时，Locker等[22]发现通过HTR2B可激活下游的NO信号通路和磷脂酶A2信号通路，释放NO及花生四烯酸，促进成骨前体细胞粘附、成骨分化及矿化，这提示HTR2B可能在NHO形成中发挥了作用。

综上所述，NHO的发生存在多靶点、多途径调控机制，这些基因和作用通路可能成为未来治疗NHO的方向，值得进一步的研究验证。