基于生物信息学的椎间盘退变miRNA-mRNA 调控关系的分析

杜文佳 晏丽 李延宏

1.兰州大学第二医院/甘肃省骨关节病重点实验室，甘肃 兰州730030；2.兰州现代职业技术学院，甘肃 兰州730000

下腰痛（low back pain，LBP）是世界范围内重要的医学问题，给社会和患者造成了巨大的经济负担［1］。据相关研究表明，超过80%的成年人在生命中的某个阶段会遭受LBP［2］。其中，约10%的腰痛患者会发展为慢性残疾［3］。LBP 的病因是多方面的，包括遗传原因，生活方式和老化［4-7］。尽管LBP 的病因很复杂，但椎间盘退变（intervertebral disc deqeneration，IDD）被认为是主要原因［8］。正常椎间盘（IVD）是由髓核（NP）、纤维环（AF）和软骨终板（CEP）组成的无血管纤维软骨垫［9］。NP 组织通过产生细胞外基质（ECM）在维持IVD 完整性方面发挥重要作用。越来越多的研究表明，异常的NP 细胞功能在IDD 的发病机制中起关键作用［10，11］。现有证据表明，IDD 的发生与多种因素相关，包括身体损伤、遗传易感性、微环境变化和炎症［12，13］。目前对于IDD 的治疗主要侧重于药物治疗和手术，暂时缓解疼痛症状。然而，这些治疗并没有完全解决IDD［14］。因此，应从分子层面开展针对IDD 的分子靶向治疗，从根本上预防IDD或恢复IVD 功能［15］。既往研究表明，许多编码基因在IDD 中存在差异表达，其中一些已被证明在IDD中起重要作用［16］。

作为基因表达的关键调节因子，miRNA 是一种内源性基因编码的小型非编码RNA，通过与靶mRNA 分子的3′-非翻译区（3′-UTR）结合来调节基因转录后的表达［17］。通常，初级miRNA 在细胞核中产生，并通过Dicer 酶在细胞质中进一步加工为成熟的miRNA。在掺入RNA 诱导的沉默复合物（RISC）后，它们可以靶向抑制多个mRNA 的翻译，从而调节大约30%的人类蛋白质编码基因和多个细胞内过程，包括细胞增殖、凋亡和细胞因子释放［18］。此外，已知miRNA 与其他内源性RNA（如长链非编码RNA、环状RNA 和mRNA）相互作用，在体外诊断细胞中形成广泛的基因调控网络［19］。这激发了人们对miRNA 在IDD 中的作用以及作为新型生物标志物和治疗剂的潜力的兴趣。随着高通量测序技术的发展，如今基因芯片技术和生物信息学方法已被广泛用于鉴定疾病潜在的生物学标志物。本研究使用生物信息学方法分析IDD 相关的高通量基因芯片（GSE56081 和GSE116726），以鉴定差异表达基因（DEGs）和差异表达miRNA（DE-miRNAs），并分析了与椎间盘退变相关的通路。

1 材料和方法

1.1 芯片数据的获取 原始数据集GSE56081 和GSE116726 从GEO（Gene Expression Omnibus，http：//www.ncbi.nlm.nih.gov/geo）数据库下载。GSE56081 数据集基于Agilent 基因芯片GPL15314 平台（Agilent033010 Probe Name version），总共包括10 个样本，5 个来自IDD患者的髓核（NP）样本和5 个来自正常的NP 样本，GSE116726 数据集基于Agilent 基因芯片GPL20712平台（Agilent-070156 Human miRNA［miRNA version］），总共包括6 个样本，3 个来自IDD 患者的髓核（NP）样本和3 个来自正常的NP 样本。

1.2 数据预处理和鉴别DEGs 使用R 软件的limma包对芯片数据进行差异表达分析，limma 包是一个非常流行的统计计算和统计映射工具［20］。分别从GSE 56081 和GSE116726 芯片数据中提取了mRNA 数据和miRNA 数据，并对它们进行基因注释。然后使用Limma 数据包进行DEGs 和DE-miRNAs 的表达分析。差异表达的筛选标准为：P＜0.05，|log fold-change|＞2。

1.3 功能富集分析 为了识别DEGs 中的生物学功能并了解基因的重要性，通过GO 分析进行验证，其中包括以下类别：BP_Fat（生物学过程）；CC_Fat（细胞成分）；和MF_Fat（分子功能）。KEGG（The Kyoto Encyclopedia of Genes and Genomes；http：//www.genome.jp/kegg）数据库是由基因组，酶促途径和生物化学组成的在线数据库的集合［21］。WebGestalt（WEB-based Gene SeT Ana-Lysis Toolkit，http：//www.webgestalt.org）是一款基于网络的功能丰富的分析工具，分析结果支持12 种生物和来自各种数据库和技术平台的324 个基因标识符以及来自公共数据库的150937 个功能类别［22］。并采用WebGestalt 对上调DEGs 和下调DEGs 分别进行GO分析。KOBAS3.0（http：//kobas.cbi.pku.edu.cn/kobas3）是一款用于基因/蛋白质功能注释和基因集富集的网页分析工具，分析结果提供了5944 种物种的KEGG 途径信息以及71 种物种的基因GO 注释［23］。本研究利用KOBAS3.0 对上调DEGs 和下调DEGs 进行KEGG分析。P＜0.05 被认为筛选的GO 术语和KEGG 途径具有统计学差异。

1.4 构建miRNA-mRNA 调控关系网络 miRTarBase是一个经过实验验证的miRNA-靶标相互作用数据库［24］。本研究分别将上调和下调的DE-miRNAs 通过查询miRTarBase 数据库，获取置信度相对可靠的靶基因。随后将上调DE-miRNAs 的靶基因与下调的DEGs取交集，下调DE-miRNAs 的靶基因与上调的DEGs 取交集，分别构建出miRNA-mRNA 调控关系网络，并且由Cytoscape 软件进行可视化处理（版本3.2.0）［25］。

2 结果

2.1 筛选出DEGs 在GSE56081 的椎间盘样本中总共鉴定了523 个DEGs，其中存在3 上调DEGs74 个，下调DEGs149 个（见图1A-1B）。热图中，红色代表上调的基因，蓝色代表下调的基因；火山图中红色代表上调基因，蓝色代表下调基因。

2.2 GO 富集分析 为了了解变性和非变性椎间盘样本中DEGs 的生物学作用，对上调和下调的DEGs 分别进行了GO 富集分析。上调和下调的DEGs 涉及的生物进程均包括生物调节、代谢过程、对刺激的反应、多细胞生物过程、细胞通讯、定位等。上调和下调的DEGs 分子功能均参与蛋白质结合、离子结合、核酸结合、转移酶活性、水解酶活性、结构分子活性等（见图1C）。上调和下调的DEGs 细胞成分（CC）均涉及细胞膜、细胞核、内膜系统、含蛋白质复合物、细胞外基质、细胞投射、细胞骨架等（见图1D）。

图1 GSE56081 椎间盘样本中鉴定DEGs 和GO 富集分析上调和下调的DEGs

2.3 KEGG 富集分析 通过KOBAS3.0 软件系统对KEGG 通路分析后，分析结果发现上调DEGs 参与的信号通路包括人T 细胞白血病病毒1 型感染、松弛素信号通路、疟疾、代谢途径、TGF-β 信号通路、人乳头瘤病毒感染、泛素介导的蛋白质水解、雌激素信号通路等通路；下调DEGs 的KEGG 途径：包括蛋白质消化吸收、阿米巴病、幽门螺杆菌感染中的上皮细胞信号转导、类风湿性关节炎、IL-17 信号通路、AGE-RAGE 信号通路、NOD 样受体信号通路、肿瘤坏死因子信号通路、细胞因子-细胞因子受体相互作用等通路。

2.4 筛选出DE-miRNAs 对GSE116726 的椎间盘髓核样本中的DE-miRNAs 进行鉴定，从GSE116726 的椎间盘髓核样本中共鉴定出858 个DE-miRNAs，其中上调的DE-miRNAs 有362 个，下调的DE-miRNAs 有496 个。

2.5 构建miRNA-mRNA 调控关系 DE-miRNAs 通过miRTarBase 数据库预测了靶基因，通过韦恩图将靶基因与DEGs 取交集。与上调miRNA 的靶基因有交集的DEGs 共11 个，与下调miRNA 的靶基因有交集的DEGs 共有12 个。通过Cytoscape 软件可视化后，本研究分别得到了上调DE-miRNAs 与下调DEGs 和下调DE-miRNAs 与上调DEGs 组成的miRNA-mRNA 调控关系（见图2）。

图2 miRTarBase 数据库预测了靶基因与DEGs 取交集

3 讨论

椎间盘位于椎骨之间，由髓核、纤维环和软骨终板组成。因其独特的解剖结构、力学特性和有限的血液供应，在生活习惯、职业和外伤的影响下，椎间盘易发生变性［26-28］。随着椎间盘的退变，许多脊柱退行性疾病（颈椎病、腰椎间盘突出症、腰椎滑脱、腰椎管狭窄）亦随之而来［29-31］。目前，机械应力，年龄，创伤和遗传因素被普遍认为是导致IDD 的重要原因［32-33］。

尽管目前对椎间盘退变的发病机制和病因尚不十分清楚，但随着科学技术不断进步和基因组学研究技术的成熟，生物信息学方法的出现将加速了人类疾病机理的研究进展，并极大地促进疾病发病机制研究的突破［34-35］。本研究通过对变性椎间盘样本和对照样本之间的基因芯片进行比较分析，分别确定了374 个上调DEGs 和149 个下调DEGs；362 个上调DE-miRNAs和496 个下调DE-miRNAs，确定了4 条信号通路，并且鉴定出32 对miRNA-mRNA 调控关系。

根据GO 分析结果，上调和下调的DEGs 主要集中在生物调节、代谢过程、刺激反应、多细胞生物过程等生物学过程中。而KEGG 富集显示，上调DEGs 与IDD可能相关的通路包括TNF 信号通路、TGF-β 信号通路、PI3K-Akt 信号通路、AGE-RAGE 信号通路、雌激素信号通路等。而下调的DEGs 的KEGG 富集结果表明，下调的DEGS 与IDD 可能相关的途径包括IL-17信号通路、AGE-RAGE 信号通路、TNF 信号通路、NFkappa B 信号通路等。众所周知，细胞外基质（ECM）降解是椎间盘退变的重要特征之一。TNF 信号通路已在许多情况下显示出可以引起ECM 变性的改变。有研究表明，TNF-α 可抑制多种胶原蛋白类型（聚集蛋白聚糖，纤维调节蛋白）从而改变了椎间盘的生物力学，最终导致IDD 的发生［36-37］。TGF-β 信号通路在IDD 过程中起着重要作用。抑制ECM 降解和增加ECM 合成可以激活TGF-β 信号通路，进而促进细胞增殖和抑制细胞死亡以及减轻炎症反应，对椎间盘具有保护作用。但是，TGF-β 信号转导的过度激活亦可能会导致IDD 的发生［38］。雌激素信号通路在椎间盘退变过程中有着重要意义。雌激素信号传导途径遍布人体组织，并参与许多病理过程。相关研究表明，17-β-雌二醇可以通过激活雌激素β 受体来促进纤维环细胞的增殖［39］。随着年龄的增长，雌激素的分泌将逐渐减少，这会导致雌激素信号通路的激活，随之炎症因子的表达和基质的降解增加，从而加速IDD 的进程［40］。PI3K-Akt 信号通路可能对人髓核具有保护作用［19］。PI3K-Akt 途径的激活，可通过增加ECM 含量减少细胞凋亡和自噬来预防IDD。已有的研究表明，激活PI3K-Akt 途径可以抑制IL-1β 介导的NP 细胞凋亡从而抑制IDD 的发生［41-42］。

本研究一共鉴别出了32 对miRNA-mRNA 调控关系。在下调miRNA 和上调mRNA 的调控关系中，HMBOX1 受到hsa-miR-93-5p 和hsa-miR-373-3p 调控，而hsa-miR-93-5p 同时调控ANKN 和LZIC，hsamiR-373-3p 同时调控KREMEN1 和CCND2。在上调miRNA 和下调mRNA 的调控关系中，LOCR 同时受hsa-miR-107，hsa-miR-196b-5p，hsa-miR-27b-3p，hsa-miR-30b-5p 和hsa-miR-30c-5p 五个miRNA 调控，值得进一步深入研究，LOCR 目前报道主要与劳氏眼脑肾综合征有关［43］；ATG9A 同时受hsa-miR-15a-5p 和hsa-miR-497-5p 调控；SBNO1 和SYNJ1 同时受hsa-miR-15b-5p 调控。构成miRNA-mRNA 调控关系的miRNA 和mRNA 多报道与癌症进程有关，鲜有报道与IDD 有关，这或许为研究IDD 提供了新的思路和研究方向。

综上所述，本研究通过生物信息学方法对IDD 相关通路和基因进行分析，鉴别了523 个DEGs 和858个人DE-miRNAs。对DEGs 进行GO 和KEGG 的富集分析，进一步发现细胞外基质，TNF 信号通路，TGF-β 信号通路，雌激素信号通路和PI3K-Akt 信号通路在IDD过程中起了重要作用。此外，还构建了32 对miRNAmRNA 调控关系，其潜在的机制值得我们进一步研究。总之，本研究鉴别了IDD 相关的信号通路和可信的miRNA-mRNA 调控关系，为后续的研究提供了新的理论基础和方向。