MAPK8、AR、RARA和FGF13是骨关节炎软骨退变过程中的潜在生物标志物

杨明义 苏亚妮 许珂 彭侃 刘林 鲁超 侯卫坤 井文森 冯磊 许鹏

(1西安交通大学附属红会医院关节外科,陕西 西安 710054;2延安大学)

骨关节炎(OA)是一种常见的与年龄相关的致残疾病,会引起疼痛、肿胀、僵硬和活动受限,严重影响患者的生活质量,其病理通常包括关节软骨、滑膜、软骨下骨、韧带、关节囊和关节周围肌肉的结构改变〔1〕。OA患者会经历关节软骨进行性退变、软骨骨赘形成、软骨下骨重塑(硬化)、滑膜炎和关节囊内炎性细胞浸润〔2～5〕。软骨退行性变是OA的主要问题,导致关节疼痛和功能障碍。随着年龄的增大,软骨磨损日益加重并进行性退化,引起并促进了OA的进展。对于OA的治疗,目前并没有有效的药物〔6,7〕,这主要是因为对OA起始阶段驱动病理过程的机制了解有限〔8〕。目前OA的发病机制尚不明确,软骨退化作为OA的始发因素,在OA的发病机制及病变进展中发挥着重要的作用。

微阵列在基因表达分析中发挥着重要作用,是医学研究的重要工具,具有重要的临床应用价值〔9〕。近年来,利用微阵列技术和生物信息学分析进行了大量的基因表达谱的研究,为研究疾病的分子机制提供了一定的依据。本研究通过微阵列基因表达谱和生物信息学分析来识别OA软骨退变过程中的生物标志物和相关通路。从生物学功能和通路的结果来看,这些参与OA软骨退变的潜在生物标志物和通路在OA软骨进行性退变的过程中发挥着重要的作用,对OA的发病机制和治疗研究具有一定的生物学指导意义。

1 材料与方法

1.1微阵列数据 从基因表达综合数据库(GEO,http://www.ncbi.nlm.nih.gov/geo)下载GSE16464和GSE10575基因表达数据集。GSE16464数据集中包含有3例OA软骨组织样本和3例正常软骨组织样本,其平台文件为GPL570。GSE10575数据集中包含有6例OA软骨组织样本和2例正常软骨组织样本,其平台文件为GPL570。

1.2鉴定差异表达基因(DEGs) GEO2R(http://www.ncbi.nlm.nih.gov/geo/geo2r)是一个交互式的网络工具,它允许用户比较GEO系列中的两组或两组以上的样本,以便识别在不同实验条件下表达的不同基因。GEO2R筛选OA软骨组织样本和正常软骨组织样本之间的DEGs,筛选标准:P<0.05、|log差异倍数(FC)|≥1。

1.3DEGs的基因本体论(GO)和京都基因与基因组百科全书(KEGG)富集分析 用于标注、可视化和集成发现的数据库(DAVID,http://david.ncifcrf.gov,version6.8)是一个整合了生物学数据和分析工具的生物信息数据库〔10〕。Cytoscape(version3.8.0)是一个专注于开源网络可视化和分析的软件,其插件ClueGO支持多个数据集注释,可以多集合富集结果比较。基因集富集分析(GSEA)根据对基因的定位、功能和生物学意义构建一个分子标签数据库,数据库中包含了多个功能基因集,通过分析基因表达数据,得到表达状况是否在某种功能上显著富集。GO是一个主要的生物信息学工具,用于注释基因和分析这些基因的生物学过程〔11〕。KEGG是一个整合了基因组、化学和系统功能信息的数据库。分别用DAVID和Cytoscape对DEGs进行GO和KEGG富集分析,GSEA对两个芯片中的全部基因信息做富集分析。DAVID选取标准为P<0.05。Cytoscape选取标准为P<0.05。GSEA选取标准为错误发现率(FDR)≤0.25,P≤0.05。

1.4蛋白互作网络(PPI)构建 检索相互作用基因/蛋白质的搜索工具(STRING,http://string-db.org,version11.0)数据库可以对互作基因网络进行预测,为了解疾病产生或发展的机制提供更加深入的科学依据。利用STRING构建DEGs的PPI网络,综合评分>0.4的交互作用被认为具有统计学意义。

1.5Hub基因的选择和分析 Cytoscape的插件Cytohubba(Degree算法)筛选在PPI中网络连接度最高的前10个基因为Hub基因,并对其进行可视化。对10个Hub基因利用Cytoscape的BINGO插件进行GO富集分析,FDR<0.000 1被认为有统计学意义。

1.6Hub基因-miRNAs相互作用网络的构建和关键miRNAs的预测 TargetScan、miRDB和miRWalk三个数据库预测Hub基因的靶向miRNAs,预测结果取交集。Cytoscape构建Hub基因-miRNAs的相互作用网络,针对两个以上Hub基因的miRNAs被认为是关键miRNAs〔12〕。

2 结 果

2.1DEGs的鉴定 GEO2R在线分析后,GSE16464鉴定出1 091个上调DEGs和1 539个下调DEGs,GSE10575鉴定出604个上调DEGs和849个下调DEGs。对两个数据集的上调DEGs和下调DEGs分别取交集,GSE16464和GSE10575以P<0.05,logFC≥1为选取标准,得到53个上调DEGs;GSE16464和GSE10575以P<0.05,logFC≤-1为选取标准,得到119个下调DEGs,最终得到172个DEGs。

2.2DEGs的GO和KEGG富集分析 DAVID对DEGs进行GO和KEGG富集分析,生物途径(BP)在信号转导、免疫反应、RNA聚合酶Ⅱ启动子转录的正调控和凝血显著富集;细胞定位(CC)在细胞膜、细胞质膜、细胞外泌体和胞外区显著富集;分子功能(MF)在蛋白质结合、蛋白质异二聚体活性、钙黏着蛋白结合参与细胞与细胞的黏附和核小体DNA结合显著富集(表1)。KEGG显著富集于癌症中的转录失调、系统性红斑狼疮、酗酒、cAMP信号通路和细胞因子与细胞因子受体的相互作用(表2)。Cytoscape对DEGs进行GO和KEGG富集分析,GO主要富集于类固醇激素的生物合成过程、吞噬作用、器官生长、核苷二磷酸生物合成过程、蛋白质解聚的负调控、上皮细胞增殖的负调控、白细胞介素7介导的信号通路和免疫受体活性(图1A);KEGG主要富集于胆固醇代谢、FcεRI信号通路、幽门螺杆菌感染中的上皮细胞信号传导、癌症中的转录失调、酗酒和系统性红斑狼疮(图1B)。GSEA对两个芯片中的全部基因做富集分析,发现也显著富集于蛋白质的生物合成、RNA聚合酶、类固醇激素的生物合成和核苷酸的生物合成。

表1 DEGs 的GO富集分析

表2 DEGs的KEGG富集分析

2.3PPI构建 利用STRING对DEGs进行PPI分析,并对其可视化红色为上调基因,绿色为下调基因。见图2。

图2 STRING构建DEGs的PPI

2.4Hub基因的选择和分析 在PPI中,用Cyto-scape的插件Cytohubba(Degree算法)筛选前10个网络连接度最高的基因为Hub基因,10个Hub基因在PPI中的网络连接度由高到低依次为丝裂原活化蛋白激酶(MAPK)8、雄激素受体(AR)、SRY-box转录因子(SOX)9、骨形态发生蛋白(BMP)2、视黄酸受体α(RARA)、RUNX家族转录因子(RUNX)1、细胞色素P450家族19亚家族A成员(CYP19A1)、成纤维细胞生长因子(FGF)13、维生素D受体(VDR)和超氧化物歧化酶(SOD)2,见图3A。10个Hub基因相互作用网络图见图3B。利用Cytoscape的BINGO插件对10个Hub基因进一步GO富集分析,显著富集于RNA聚合酶Ⅱ启动子转录的调控、氮化合物代谢过程的正调控、多细胞生物发展、骨骼发育、类固醇激素受体活性、配体依赖性核受体活性、调节细胞增殖、细胞生物合成过程的正调控和转录激活子活性。见图4。

深色为Hub基因,颜色越深,排名越靠前图3 从PPI中筛选的Hub基因

每一个圆圈代表一个GO Term,颜色越深,显著性越高,圆圈越大图4 Cytoscape的BINGO插件对Hub基因进行GO富集分析

2.5Hub基因-miRNAs相互作用网络分析的构建和关键miRNAs的预测 TargetScan、miRDB和miRWalk三个数据库预测Hub基因的靶向miRNAs,Cytoscape构建Hub基因和miRNAs的相互作用网络见图5。同时针对两个以上Hub基因的miRNAs被定义为关键miRNAs,共得到23个关键miRNAs:miR-1185-1-3p、miR-1185-2-3p(靶基因均为CYP19A1和FGF13),miR-125a-5p、miR-125b-5p(靶基因均为VDR和SOD2),miR-138-5p(靶基因为SOX9和RARA),miR-139-5p(靶基因为MAPK8和RUNX1),miR-27a-3p、miR-27b-3p(靶基因均为RARA和RUNX1),miR-302a-3p、miR-302d-3p、miR-302e(靶基因均为RUNX1和VDR),miR-3153、miR-3159(靶基因均为CYP19A1和SOD2),miR-3169(靶基因为FGF13和SOD2),miR-3675-5p(靶基因为CYP19A1和FGF13),miR-373-3p(靶基因为RUNX1和VDR),miR-4319(靶基因为VDR,SOD2),miR-432-3p、miR-512-5p(靶基因均为CYP19A1和SOD2),miR-520c-3p、miR-520d-3p(靶基因均为RUNX1和VDR),miR-6733-5p、miR-6744-5p(靶基因均为CYP19A1和SOD2)。

3 讨 论

本研究结果发现,GSE16464和GSE10575两个芯片中的OA软骨组织样本大多富集于RNA聚合酶Ⅱ启动子转录的调控、类固醇激素的生物合成和细胞的增殖。对这10个Hub基因进行了更深层次的miRNAs挖掘,共得到23个关键miRNAs,通过对Hub基因-关键miRNAs相互作用网络的分析我们发现这23个关键miRNAs大多靶向于SOD2、CYP19A1、RUNX1和VDR这4个Hub基因。可见SOD2、CYP19A1、RUNX1和VDR这4个Hub基因是OA软骨退变的过程中的重要生物标志物。

SOD2是铁/锰超氧化物歧化酶家族的成员。该基因的突变与特发性心肌病(IDC)、过早衰老、偶发性运动神经元疾病和癌症有关。人OA软骨中SOD1、SOD2和SOD3显著降低〔13〕。三种超氧化物歧化酶均在人软骨中大量表达,而在终末期OA软骨中,尤其是SOD2表达明显下调,SOD2启动子在OA软骨组织中有明显的DNA甲基化改变,软骨细胞中SOD2的缺失不但增加活性氧(ROS)和降低胶原酶的表达〔14〕,而且导致氧化损伤和线粒体功能障碍〔15〕,提示SOD2下调可能是OA发病机制的一个潜在因素。研究发现牛磺酸A通过增加SOD2含量并抑制丙二醛(MDA)和ROS在OA的治疗中具有很强的抗氧化活性〔16〕。本研究发现,SOD2在OA软骨细胞中下调,与上述研究结论相一致,可见SOD2是OA软骨退变过程中的重要生物标志物。

CYP19A1编码酶是细胞色素P450超家族的成员。细胞色素P450蛋白是单加氧酶,可催化涉及药物代谢和胆固醇、类固醇和其他脂质合成的许多反应,该蛋白定位于内质网并催化雌激素生物合成的最后步骤。性激素特别是雌激素,与关节软骨代谢和绝经后OA的发病机制有关,可调节软骨细胞的活性及涉及软骨基质合成代谢和分解代谢的多种因素的合成,包括金属蛋白酶一氧化氮和活性氧〔17～20〕。在绝经后的女性和男性中,周围组织中的雄激素前体是雌激素的主要来源,该反应由位于15号染色体上的CYP19A1基因产物芳香化酶催化〔21〕。通过CYP19A1将雄烯二酮转化为雌酮(E1),将睾酮转化为17β雌二醇(E2),在组织内源性合成雌激素中起关键作用〔22〕。可见CYP19A1与OA软骨退变密切相关,这与本研究结论相一致。

RUNX1编码的蛋白质代表核心结合因子(CBF)的α亚基,被认为与正常造血过程有关。CBF是一种异二聚体转录因子,可与许多增强子和启动子的核心元件结合。Runx1是决定造血功能的转录因子,存在于软骨膜和软骨细胞中。Runx1能够像润滑蛋白一样对机械负荷和软骨损伤作出反应,Runx1在浅层细胞中与其他因子一起具有通过刺激增殖或促进表面蛋白表达来维持组织表面完整性的潜在功能〔23〕。有研究表明针对RUNX1的靶向治疗为开发针对OA的疾病改良药物提供了希望〔24〕。本研究发现,RUNX1与OA软骨退变密切相关,其在OA软骨退变过程中的机制值得进一步研究。

VDR编码维生素D3受体,它是配体诱导型转录因子的核激素受体超家族的成员。维生素D3受体调节多种代谢途径,如参与免疫应答和癌症反应,但其下游靶标主要参与矿物质代谢。VDR被认为参与软骨和骨代谢,VDR基因多态性已被鉴定和分析,其与多种疾病相关,包括OA、癌症、糖尿病、心血管疾病、病毒感染、结核病、牙周炎、尿石和自身免疫性疾病〔25〕。有报道表明,VDR多态性与OA易感性相关〔26〕。这与本研究结论VDR是OA软骨退变过程中的潜在生物标志物相一致。

MAPK8在PPI中网络连接度最高,在10个Hub基因中排名第一。MAPK8编码的蛋白质是丝裂原活化蛋白(MAP)激酶家族的成员。MAP激酶充当多种生化信号的整合点,并参与多种细胞过程,例如增殖,分化,转录调节和发育。该激酶被各种细胞刺激激活,并靶向特定的转录因子,从而响应细胞刺激介导早期的基因表达。Rabl3的高表达与非小细胞肺癌患者MAPK8/9/10介导的自噬抑制的低生存率相关〔27〕。miR-190通过靶向MAPK8和调节MAPK8/ERK信号通路保护心肌细胞免受H2O2诱导的凋亡〔28〕。MAPK8通过MAPK信号通路介导对替莫唑胺的抗性和胶质母细胞瘤的细胞凋亡〔29〕。研究发现他莫昔芬通过MAPK8/叉头转录因子(FoxO)途径诱导乳腺癌的脂肪肝疾病〔30〕。然而目前尚无报道MAPK8与OA相关,本研究发现MAPK8在OA软骨退变过程中上调,是OA软骨退变的潜在生物标志物,具有较大的研究价值。

SOX9编码的蛋白质与HMG-box类DNA结合蛋白的其他成员一起识别CCTTGAG序列。 它在软骨细胞分化过程中起作用,缺乏可导致骨骼畸形综合征。据报道,OA软骨中SOX9表达降低〔31〕,SOX9过表达可减轻人类OA的进展〔32〕,microRNA-384-5p通过靶向SOX9经核转录因子(NF)-κB信号通路抑制软骨细胞的凋亡而减轻OA〔33〕。本研究发现,SOX9在OA软骨退变中下调。

BMP2编码蛋白质的(TGF)-β超家族的分泌配体。 编码的前蛋白经过蛋白水解生成二硫键连接的同型二聚体的每个亚基,后者在骨骼和软骨发育中发挥作用。幼鼠中诱导BMP2过表达可导致实验性OA中骨赘形成加重,但不改变软骨损伤〔34〕。OA病变周围的软骨细胞中BMP2水平升高〔34〕,本研究发现,BMP2在OA软骨退变中下调,可见BMP2与OA软骨具有密切的关系。

AR编码具有3个主要功能域的蛋白质:N末端域、DNA结合域和雄激素结合域。该蛋白质起类固醇激素激活的转录因子的作用。AR与膀胱癌发生发展的多种途径相互作用。RARA以配体依赖性方式调节转录。该基因座和其他几个基因座之间的易位与急性早幼粒细胞白血病有关。FGF13编码的蛋白质是FGF家族的成员。FGF家族成员具有广泛的促有丝分裂和细胞存活活性,并参与多种生物学过程,包括胚胎发育、细胞生长、形态发生、组织修复、肿瘤生长和侵袭。研究发现,FGF13在局灶性皮质异常病变中的表达增加〔35〕;FGF13是一种新的NF-κB调节因子,可增强病理性心肌肥厚〔36〕。AR、RARA和FGF13与OA的相关性目前尚无报道,本研究发现AR、RARA和FGF13在OA软骨退变中下调,是OA软骨退变中的潜在作用因子。