骨关节炎患者滑膜组织 miRNA-mRNA 相互调控网络及分子作用机制的研究

徐扬 吴驷东 钱臣 朱斌杰

骨关节炎 ( osteoarthritis，OA ) 作为临床工作中最常见的关节退行性病变，其主要的临床特征为关节疼痛、活动受限以及功能障碍。流行病学调查报告显示，随着年龄的增加，OA 的患病率也逐渐增加，75 岁以上年龄的患病率可高达 80%，严重影响了老年人的生活质量。如何有效预防 OA 的发生发展，从 OA 的发生机制中及早进行干预是目前临床研究的热点与难点。

基础研究结果显示，OA 发病的主要病理学变化为关节软骨基质的减少、骨赘形成以及滑膜无菌性的炎症。而滑膜作为关节中内部结构与肌骨骼相互连接的桥梁，在维持关节的正常生理功能中发挥着重要的作用。研究显示，滑膜组织可以通过分泌透明质酸从而有效地保持关节面的润滑程度，同时也可以有效地较少关节面的相互摩擦，为关节面的软骨提供较好的保护作用。而随着滑膜中巨噬细胞的形态功能改变，会进一步诱发甚至加重 OA的疾病进程。既往研究显示，当滑膜组织中的巨噬细胞可以产生多种炎性因子并使其参与到 OA 的病程发展过程中，同时还与患者的疼痛、关节僵硬等临床特征密切相关。但是，在 OA 患者中滑膜组织的具体分子调控机制尚未完全明晰，仍须进一步的深入研究。

近年来随着基因芯片以及二代测序技术的普及，越来越多的基因测序被应用到 OA 患者中，从而对于明确 OA 发生发展的具体分子调控机制提供了更多的理论工具。Shen 通过整合基因表达数据库( gene expression omnibus，GEO ) 数据库中关于 OA患者膝关节滑膜组织的基因芯片并通过生物信息学分析的方法对数据进行统计分析，结果显示 OA 主要的分子调控机制是通过丝裂原活化蛋白激酶、低氧诱导因子等信号通路的开关从而影响 OA 的疾病发展。但是，上述研究仅处于对 mRNA 的研究，具体 mRNA 在 OA 疾病中是如何被上级信号所调控的仍未有相关研究报道。

随着研究的进一步深入，微小 RNA ( microRNA，miRNA ) 的概念也被相应提出与证实。miRNA 作为细胞内一类小的非编码 RNA，可以通过调控其对应靶基因的表达以及信号通路的开关，从而影响疾病的发生发展。既往研究显示，miRNA 在 OA 的发生发展中起到了重要的调控作用。如 miR-146 可以通过降低 OA 滑膜成纤维细胞内血管黏附分子-1 的表达，从而调控 OA 的发生；miR-602 可以通过调控滑膜成纤维细胞的血管生成活动，从而有效阻止OA 的发生发展。诸多研究结果都显示，miRNA在 OA 的发生发展过程中都起到了重要的调控作用，但是既往研究往往仅针对某一单一 miRNA，缺乏整体性。而本研究拟通过生物信息学分析的方法有效筛选出 OA 患者滑膜组织内差异表达的 mRNA和 miRNA 并以此构建相应的调控网络，明确引起OA 发生的具体分子调控机制。

材料与方法

一、芯片数据的获取

本研究中数据来源为从 GEO 获取。mRNA 数据来源于 GEO 数据库中的基因芯片数据 GSE554457( 芯片平台 GPL96 )，GSE554457 包含 33 例患者的滑膜组织数据，其中 10 例为健康滑膜组织，10 例为OA 患者膝关节滑膜组织；miRNA 数据来源于 GEO数据库中的 RNA 测序数据 GSE126677 ( 芯片平台GPL167 )，一共包含 5 例患者的滑膜组织数据，其中 3 例为健康滑膜组织，2 例为 OA 患者膝关节滑膜组织。

二、获取差异表达 miRNA / mRNA

差异表达 miRNA 和差异表达 mRNA 的获取是基于 R 4.0.3 进行分析。使用 limma 加载包对获取的mRNA 和 miRNA 数据进行标准化处理后通过设标准为＜ 0.05，|| ＞ 1 从而分别获得差异表达的 miRNA 和 mRNA。通过 R 软件中的 pheatmap 加载包绘制热图，ggplot2 加载包绘制火山图可视化差异基因的表达情况。将差异表达的 mRNA 通过 string( string-db.org ) 建立蛋白互作网络，并将获得的网络导入 Cytoscape 3.8.0 中进行可视化调整，最终呈现出差异表达 mRNA，通过 MCODE 插件筛选排名前三的关键模块。

三、差异基因的功能富集分析

本研究中针对差异表达 mRNA 的功能富集分析主要是通过 Fun Rich 3.1.3 和 R 4.0.3 进行分析并绘图。差异基因的功能富集分析主要是从三个方面对差异基因进行注释：生物学过程 ( biological process，BP )、细胞组分 ( cell component，CC ) 和分子功能( molecular function，MF )，差异基因的富集通路分析主要是通过 DAVID ( ncifcrf.gov ) 进行。

四、差异 miRNA 的靶基因预测

通过 Fun Rich 3.1.3 对差异表达的 miRNA 进行靶基因的预测。

五、miRNA-mRNA 调控网络的可视化构建

由于 miRNA 对于 mRNA 具有负向调控作用，本研究中首先通过差异表达的 miRNA 预测靶基因。将上调的 miRNA 预测靶基因与下调的差异表达mRNA，下调的 miRNA 与上调的差异表达 mRNA 分别取交集绘制韦恩图 ( Veen ) 获得在最终靶基因。根据靶基因和对于 miRNA 的相互作用关系在 Cytoscape中进行可视化分析。

六、统计学处理

本研究统计学分析通过 R 4.0.3 进行，差异基因的筛选采用 R 4.0.3 的 Wilcox 检验。＜ 0.05 为差异有统计学意义。设置差异基因的筛选标准为＜0.05，|| ＞ 1。

结 果

一、差异表达 mRNA 和 miRNA 的筛选

通过差异基因表达分析结果显示，GSE126677包含 71 个差异表达的 miRNA，其中包含 9 个上调差异表达 miRNA 和 62 个下调差异表达 miRNA。GSE55457 中包含 874 个差异表达的 mRNA，其中包含 113 个下调差异表达 mRNA 以及 761 个上调差异表达 mRNA ( 图 1 )。通过 Cytoscape 构建差异表达mRNA 的蛋白可视化网络同时通过 MCOD 筛选出三个关联程度最高的模块 ( 图 2 )。

图1 mRNA 与 miRNA 的差异表达基因分析图 a：GSE55457 展示 OA 患者滑膜细胞差异表达 mRNA 热图分析；b：GSE126677 展示 OA患者滑膜细胞差异表达 miRNA 热图；c：GSE55457 展示 OA 患者滑膜细胞差异表达 mRNA 火山图；d：GSE126677 展示 OA 患者滑膜细胞差异表达 miRNA 火山图 ( 注：a、b 中的红色表示高表达，蓝色表示低表达；c、d 中的红色表示上调基因，蓝色表示下调基因，灰色表示无统计学差异的基因 )Fig.1 Analysis of differentially expressed genes between mRNA and miRNA a: mRNA thermomap of differential expression of synoviocytes of patients with osteoarthritis in GSE55457; b: miRNA thermomap of differentially expressed synoviocytes of patients with osteoarthritis in GSE126677; c: Volcanic map of differentially expressed mRNA of synoviocytes of patients with osteoarthritis in GSE55457; d: Volcanic map of differentially expressed miRNA of synoviocytes of patients with osteoarthritis in GSE126677 ( Notice: Red in figure a and b indicated high expression, while blue indicated low expression; red in figure c and d indicated up-regulated genes, blue indicated down-regulated genes, and gray indicated genes with no statistical differences )

图2 差异表达的 mRNA 之间的蛋白互作网络图 ( PPI ) [ 注：a：差异表达基因的蛋白互作网络，b、c、d：Cytoscape 3.8.0 中MCODE 插件筛选出来的关联强度最高的 3 个模块 ( 红色表示上调，绿色表示下调 ) ]Fig.2 Protein interaction network map between differentially expressed mRNA ( PPI ) [ Notice: a: Protein interaction network of differentially expressed genes. b, c, d: The three modules with the highest correlation intensity selected by MCODE plug-ins in Cytoscape 3.8.0 ( red for up-regulation and green for down-regulation ) ]

二、差异表达 mRNA 的富集分析结果

通过 FunRich 对差异基因进行 GO 富集分析。结果显示差异表达 mRNA 在细胞组分中主要存在于细胞外基质蛋白质、血小板颗粒、质膜、细胞核、细胞膜、细胞外基质、细胞膜小窝；在分子功能中主要参与转膜受体活动、转录调控活动、转录因子活动、酪氨酸蛋白、金属离子结合、生长因子活性、细胞外基质形成、细胞骨架蛋白结合、细胞骨架蛋白锚定活动、细胞粘连分子活动等过程中；在生物学进程中参与信号转导、碱基的调控、基因表达的调控、胰酶表达的调控、细胞迁移、细胞增殖、细胞交流以及细胞的锚定等进程中 ( 图 3 )。而 KEGG富集通路分析结果显示差异基因主要富集于 b1-整合素细胞表面相互作用、整合素家族细胞表面相互作用、血管内皮生长因子 ( VEGF ) 和 VEGF 受体调控网络以及 CDC42 活动调控等信号通路中 ( 图 4 )。

图3 差异表达的 mRNA GO 功能富集分析 a：GO 富集分析展示差异基因富集于细胞组分的功能变化；b：GO 分析展示差异基因富集于分子功能的功能变化；c：GO 分析展示差异基因基于生物学过程的功能变化；d：气泡图描述细胞组分的变化情况；e：气泡图描述分子功能的功能变化；f：气泡图描述生物学过程的功能变化情况 ( 注：颜色表示差异基因的统计学意义，红色表示统计学意义明显，绿色表示统计学意义较低；大小表示富集于相应功能的基因数目，越大表示富集的基因越多 )Fig.3 Functional enrichment analysis of mRNA a: GO functional enrichment analysis of differential expression showed that the functional changes of differential genes were enriched in cellular components; b: GO analysis showed functional changes of differential genes enriched in molecular function; c: GO analysis showed functional changes of differential genes based on biological processes; d: Bubble diagram described the changes of cell components; e: The bubble diagram described the functional changes of molecular functions; f: The bubble chart described the functional changes of biological processes ( Notice: Color indicated the statistical significance of differential genes. Red indicated statistically significant and green indicated statistically significant. Size indicates the number of genes enriched in the corresponding function. The larger the number of genes enriched, the more genes were enriched )

图4 差异表达 mRNA 的富集通路分析Fig.4 Enrichment pathway analysis of differentially expressed mRNA

三、差异表达 miRNA 靶基因预测

利用 FunRich 3.1.3 中 miRNA Enrichment 的 Find targets 功能对获取的差异表达的 miRNA 进行靶基因的预测，最终显示 9 个 miRNA 可以预测出相应的靶基因，而预测的靶基因总数为 2110 个。

四、构建 miRNA-mRNA 调控网络

通过统计学分析发现最终获取的 9 个 miRNA 都为上调的 RNA，则其预测的靶基因应与下调的差异表达 mRNA 取交集，通过 Venn 图绘制交集后共获得19 个下调的差异表达 mRNA ( 图 5 )，相关的 miRNA和 mRNA 表达情况详见表 1。基于 miRNA 与 mRNA的调控关系将各基因导入 Cytoscape 3.8.0 中构建miRNA 差异 mRNA 调控网络，具体网络详见图 6。

图6 miRNA-RNA 相互调控作用网络 ( 红色表示上调，蓝色表示下调，菱形为 miRNA，圆形为 mRNA，大小表示 logFC，logFC )Fig.6 miRNA-RNA interaction network ( Notice: Red indicated up-regulation and blue indicated down-regulation. The diamond was miRNA and the circle was mRNA. The size indicated that the bigger the logFC, logFC is, the bigger the shape is )

表1 具有潜在靶向作用关系的 miRNA 与 mRNA 表达情况Tab.1 miRNA and mRNA expression with potential targeting relationship

图5 miRNA 预测靶基因与 GSE55457 下调基因的 Venn 图Fig.5 Venn map of target genes predicted by miRNA and downregulated genes of GSE55457

讨 论

既往研究结果显示，miRNA 与 mRNA 在 OA 的发病过程中都起到了至关重要的调控作用。而本项研究通过分别分析两组高通量测序数据并进行整合分析后，构建了 OA 发病过程中的 miRNA-mRNA 调控网络，对于 OA 的发病机制研究提供了更加可靠的理论依据。

一、mRNA 与 OA 的发生

GSE554457 为基因芯片技术获得的 OA 患者与健康人膝关节滑膜组织的 mRNA 差异表达情况。GO 富集分析结果显示，差异基因主要分布于细胞外基质、细胞核的组分，转录因子、生长因子的活化以及细胞外基质的形成上。细胞外基质作为细胞外主要的屏障，可有效地保护细胞的活性，同时抵御机械压力对于细胞的刺激、损伤，而细胞外基质的合成与降解的异常会导致细胞的活性、功能明显受损，最终引起OA 的发生发展，加重 OA 的病情。

KEGG 富集通路分析结果显示差异基因主要富集于 b1-整合素细胞表面相互作用、整合素家族细胞表面相互作用、VEGF 和 VEGF 受体调控网络以及CDC42 活动调控等信号通路中。整合素作为黏附分子家族中的主要组成部分，其可介导细胞外的机械信号传导至细胞内，从而有效地抑制细胞外基质的降解，促进细胞外基质的合成，最终达到修复软骨损伤的作用。而整合素的降低会明显影响机械信号刺激的传导，促进细胞外基质的降解，最终诱导 OA 的发生。而 VEGF 可以通过其特异性 VEGF受体作用于血管内皮细胞，促进血管内皮细胞的增生以及新生血管的形成，既往研究表明，VEGF 表达升高可以促进 OA 炎症的发生发展以及关节软骨重构，加重 OA 患者的病情。

二、miRNA 与 OA 的发生

GSE126677 为 RNA 测序技术获得的健康人膝关节滑膜组织以及 OA 患者膝关节滑膜组织的差异表达 miRNA，通过统计分析共获得 72 个差异表达miRNA。通过 FunRich 预测靶基因功能显示共 9 个上调的 miRNA 可以有效地预测出其靶基因 ( 表 1 )。将预测的靶基因与差异表达中下调的 mRNA 取交集后显示 hsa-miR-199a-5p、hsa-miR-23b-3p、hsa-miR-31-5p、hsa-miR-140-5p 的靶向作用较多，同时差异性表达较强。既往研究结果显示，has-miR-23b-3p在 OA 中表达明显上升，同时 miR-23b-3p 的表达升高会增加血清中炎症因子的表达、增加软骨细胞的凋亡、抑制二型胶原蛋白的表达。Costa Viviana通过对比健康患者和 OA 患者膝关节组织 miR-31-5p的表达发现，miR-31-5p 可以通过调控缝隙连接蛋白 ( Cx43 ) 的影响成骨细胞以及软骨细胞的表达，促进炎性因子的释放，最终促进 OA 病程的发展。miR-140-5p 在 OA 的发生发展过程中的影响调控作用已得到多方研究验证，有学者通过基础实验证实miR-140-5p 可以通过调控软骨细胞内 NF-κB 的开关从而调控细胞外基质的合成、降解，最终导致 OA的发生发展。hsa-miR-199a-5p 与 OA 发生的具体调控关系尚未被研究报道，需要进一步探索。

三、miRNA 靶向调控的 miRNA 与 OA 的发生

通过 miRNA-mRNA 调控网络可以发现，GAP43、HGF、SCRG1 与 miRNA 的关联程度较高，为miRNA-mRNA 调控网络中的关键节点。HGF 作为上游分子，其与 c-Met 结合后可以快速调控下游信号通路 ( PI3K / AKT、Ras / MAPK、JAK / STAT )的开关，从而影响细胞的功能与活性，最终引起OA 的发生发展。有学者通过大鼠实验研究发现，c-Met / HGF 的激活可以促进 OA 大鼠滑膜组织功能的恢复。SCRG1 作为关节软骨相关基因，在OA 的发生发展中也起到了主要的调控作用，Ochi Kensuke 通过基础实验显示 SCRG1 可以促进软骨细胞的形成，促进细胞外基质的合成，抑制其降解。GAP43 与 OA 发生的具体调控关系尚未被研究报道，须进一步探索。

在 OA 的疾病发生中存在 miRNA-mRNA 相互调控网络作用，本研究初步建立了 OA 发生过程中miRNA-mRNA 调控网络，对揭示 OA 发生发展的具体分子作用机制提供了强有力的理论参考依据。