基于基因芯片的TNF-α诱导骨关节炎细胞模型生物信息学分析

袁发浒，胡 松，黄丽霞，黄念芳，宋文剑，孙宾莲

(江汉大学，武汉 430056)

骨关节炎(osteoarthritis，OA)是一种常发于老年人的慢性退行性疾病，日益快速地成为世界人群主要致残因素[1-2]。据最新数据分析显示，目前在全球范围内，六十岁以上的人群中有18%的女性及9.6%的男性患有症状性的OA，其中四分之一的人更是无法进行正常的日常活动[3]。到2050年，预计将有1.3亿人患有OA，这将构成巨大的全球性社会负担，随着我国老龄化进程的加快，OA患者更是日益增加[4]。

OA的病理学因素是多方面的，具体发病机制尚不清晰，主要包含软骨下骨的重塑、滑膜炎症和关节软骨的丧失[5]。软骨细胞是软骨组织的唯一细胞类型，对维持和重塑软骨结构、保证软骨细胞外基质功能等方面起着至关重要的作用。近年来，由于软骨细胞凋亡在软骨退化中的作用，软骨细胞凋亡相关信号通路在OA方面的研究日益受到关注，抑制软骨细胞死亡目前已成为治疗软骨退化的一个潜在治疗靶点。此外炎性细胞因子也被认为是OA发生发展中的重要因素，特别是白介素1β(IL-1β)和肿瘤坏死因子α(TNF-α)，他们驱动分解代谢途径，并使OA进展持续下去[6-8]。

高通量基因表达分析目前已成为研究复杂疾病机制的重要方法，此外，高通量基因表达谱分析也有助于发现疾病新的生物标志物，监测疾病进程及评价药物疗效。目前，已有应用基因芯片技术分析OA患者相关基因表达的研究，这些研究成果有助于我们理解OA进展过程中所涉及的病理机制。然而，由于病人的个体差异，如患者的年龄、关节的分层等，这些因素往往导致研究结果产生一些偏倚。本研究通过培养原代小鼠软骨细胞，以TNF-α刺激建立OA细胞模型，应用小鼠全基因组表达谱芯片分析，从而识别差异基因表达谱和信号通路，为全面深入探寻OA致病机制，发现新的OA诊断标志物提供更多帮助。

1 材料和方法

1.1 实验动物

SPF级KM小鼠，体重 20 ～22 g，6 周龄，性别不限，来源于华中科技大学实验动物中心[SCXK(鄂)2016-0009]。实验于华中科技大学实验动物中心进行[SYXK(鄂)2016-0057]，本研究涉及的实验动物的使用及管理经华中科技大学实验动物伦理审查委员会核准，审批号：201803170028，研究按实验动物使用的3R原则给予实验动物人道的关怀。

1.2 主要试剂与仪器

重组TNF-α(西格玛奥德里奇(上海)贸易有限公司，货号：T7539)；DMEM培养基(赛默飞世尔科技(中国)有限公司，货号：11965084)；胎牛血清(赛默飞世尔科技(中国)有限公司，货号：12483020)；胰蛋白酶(Sigma)；Ⅱ型胶原酶(Gibco)；TRIzol试剂(Invitrogen)；PrimeScriptTMRT Reagent试剂盒(TaKaRa)；SYBR Green qPCR Master Mix试剂盒(TaKaRa)；小鼠基因表达微阵列芯片(Agilent，型号：G4846 A)。二氧化碳培养箱(力康生物医疗科技控股集团，型号：HF100)；体视显微镜(德国莱卡，型号：M205)；生物超净台(北京东联哈尔仪器制造有限公司，型号：BCN-1360)；Stepone Plus荧光定量PCR仪(ABI)。

1.3 实验方法

1.3.1 小鼠原代软骨细胞分离、培养及处理

小鼠原代软骨细胞于KM小鼠膝关节分离培养获得。选取6只同窝生SPF级KM小鼠，参照金旻等[9]报道的方法分离培养小鼠软骨细胞，无菌环境下分别剥离新生小鼠双侧膝关节软骨，剪碎后以2.5 g/L胰蛋白酶与0.1% Ⅱ型胶原酶联合消化后离心收集。原代软骨细胞以1×107cell/cm2密度接种于DMEM培养基(含10%胎牛血清)置于二氧化碳培养箱中培养，培养条件为37℃、5% CO2。原代软骨细胞稳定培养3 d后，将其分设为两组，一组为正常对照组(C组)，另一组为TNF-α诱导的OA模型组(T组)，每组三瓶细胞(分别来源于一只小鼠)，OA模型组以50 ng/mL的TNF-α孵育24 h，对照组以等积极PBS做相同处理。

1.3.2 RNA提取

两组细胞经24 h处理后，离心除去培养基，经PBS漂洗沉淀后立即加入总RNA抽提试剂TRIzol，依照产品说明书操作提取总RNA，提取所得的RNA用超微量分光光度计测定其浓度后，分装于-80℃冰箱备用。

1.3.3 微阵列芯片杂交与数据分析

微阵列芯片杂交实验委托于康成生物科技有限公司完成，实验流程参考Li等[10]所述。所用芯片为Cy3标记的小鼠全基因组表达谱芯片(Agilent Mouse 4x44 K Gene Expression Microarrays v2)。原始荧光数据利用Agilent公司的Gene Spring 11.0 软件进行均一化处理。以参照组为对照，以表达变化倍数大于2，P值小于0.001为条件筛选OA组差异表达基因(differentially expressed genes，DEGs)。获得的DEGs提交至DAVID平台(DAVID Functional Annotation Bioinformatics Microarray Analysis，https://david.ncifcrf.gov/home.jsp)进行基因功能聚类及通路分析。

1.3.4 荧光定量PCR验证

为评估芯片表达谱数据可靠性，从差异表达基因列表中挑选20个基因(上调、下调基因各十个)进行荧光定量PCR实验验证。软骨细胞总RNA用PrimeScriptTMRT Reagent试剂盒进行提取，采用随机六核苷酸引物合成第一链cDNA，采用SYBR Green qPCR Master Mix试剂盒配置反应体系，于荧光定量PCR仪扩增检测基因表达情况。将Gapdh设置为内参基因，基因相对表达量采用2-ΔΔCt方法进行分析。

1.4 统计学方法

2 结果

2.1 样本表达量聚类分析

对样本之间进行系统聚类，利用最广泛、最简便的层次聚类法(Hierarchical Clustering)，检查不同样品之间同一基因表达的相关性，如果大部分基因的表达量相似，相关系数高，说明数据的均一化程度高。在此，两组样本间(各三个样品)的差异表达基因进行了层次聚类，可见C组的三个样品和T组的三个样品分别聚类在一起，而两个组别也可区分开来(图1)，说明两组样本的基因表达有一定的差异性。

注：C组：正常对照组；T组：TNF-α诱导的OA模型组。下同。图1 层次聚类分析热图Note. Group C, Normal group. Group T, TNF-α-induced OA model group. The same below.Figure 1 Heat map and supervised hierarchical clustering

2.2 差异表达基因分析

火山图(Volcano Plots)是可视化表达两种不同模式条件的实用工具，它能统一体现组间变化倍数值和P值，可直观了解变化幅度和统计显著性之间的关系(考虑了变化幅度和变异性)。同时，基于这些值的基因子集也被分离表现出来。将C、T两组的基因表达量做火山图分析，如图2所示，纵向的左右两根绿线示意为表达差异倍数2.0，横向的绿线示意为P<0.05。图中的红点表示DEG(FC>2.0，P<0.05)，其中有8096个表达上调DEG和6413个表达下调DEG。

图2 T组与C组差异表达基因火山图Figure 2 Volcano plots of differentially expressed genes between the groups C and T

图3 T组与C组前35条差异表达基因分类注释条目Figure 3 The top 35 GO terms of differentially expressed genes between the groups C and T

为了获知这些DEGs的功能及所在的调控通路，分别将这些上调和下调的DEGs提交至DAVID在线平台进行功能聚类分析， DAVID的功能聚类数据库整合了Gene ontology， Interpro，KEGG等基因功能数据库，获得的基因注释与功能通路细致、可信度高。相较于正常对照组，TNF-α诱导处理的OA模型组的DEGs注释在细胞核、转录调节、细胞分化、ATP结合、细胞周期、CXC趋化因子等方面，图3详细列出了主要的前35个GO条目，图4列出了包含神经营养因子信号通路、赖氨酸降解、T细胞受体信号通路、ErbB信号通路、mTOR信号通路等前25条主要通路。其中下调表达差异最大前10个基因分别是P2ry1、Fjx1、Wnk1、Tln1、Lasp1、Pdia4、Prm2、Olfr214、Tnfrsf14、Spata2(表1)，而上调表达差异最大前10个基因分别是Dnahc14、Slc4a7、Cxcl3、Rasgrp1、Spata31d1c、Zfp866、Il1b、Nf1、Slc7a11、Cxcl2(表2)。

对骨关节炎进展相关的基因表达情况进行了整理(表3)，发现免疫调节和促炎类细胞因子相关基因(Il6、Il9、Il17a)显著上调。细胞外基质降解和重塑相关基因(Mmp3、Mmp9)显著下调。线粒体功能相关基因(Cyc1、Cox3、Uqcr10)表达显著上调。

为评估芯片表达谱数据可靠性，从差异表达基因列表中挑选20个基因(上调、下调基因各十个)进行荧光定量PCR实验验证。如图5所示，可见验证的绝大多数的基因表达情况与芯片数据吻合。

表1T组与C组前10个表达下调的差异表达基因

Table1The top 10 significantly down-regulated genes between the groups C and T

表2 T组与C组前10个表达上调的差异表达基因

表3 TNF-α诱导的OA相关基因表达情况

注：↑表示T组与C组相比，基因表达上调；↓表示T组与C组相比，基因表达下调。

Note. Compared with the group C, ↑ and ↓ indicate up-regulation and down-regulation of gene expression in the group T, respectively.

注：以Gapdh表达量做内参校正。Mmp9的基因芯片检测表达量与qPCR实验测试表达量相比， *P<0.05(t检验)。图5 基因芯片表达数据的荧光定量PCR验证Note. Normalized using Gapdh as an internal control. Microarray data of Mmp9 compared with qPCR data, *P<0.05 (t-test).Figure 5 Verification of microarray data through qPCR

3 讨论

以往对于OA病理机制的研究大多以某几个基因或某条信号通路进行，例如申成凯等[11]探讨了白介素-1β基因与晚期OA的关联，张荣凯[8]等则研究了脂蛋白基因在早期OA的表达情况，马春辉等[12]研究了凋亡基因通路在OA的作用，杨占东[13]等总结了SREBPs在调节骨关节炎关节软骨退变中的作用机制，这些研究对OA中某些显著的功能基因进行了深入细致的探讨，但总体而言，这些研究存在一定的片面性，极少有从整体上去挖掘OA相关的病理机制的研究， Qing等学者[14]曾以瘦素诱导的大鼠OA模型进行了表达谱芯片的分析，筛选出了如Bcl2l11等未被研究过的可能在OA起重要作用的基因，极大地丰富了OA相关领域的信息。本研究中，首次全面比较了TNF-α诱导的OA原代软骨细胞模型全基因组表达谱，用荧光定量PCR方法对一部分基因的芯片结果进行了验证，一些新的候选基因首次被发现与OA的发病机制有关。

最新的一些研究显示白介素IL-1和基质金属蛋白酶(MMP)的表达异常在OA的发病机理中有重要作用[15-17]。通过对GO和KEGG通路的分析，证实了这些基因在OA发病机制中的重要作用。IL-1是OA进展过程中软骨退变最重要的因素之一，也是基质稳态病理生理失衡的主要诱发因素。IL-1还可抑制软骨细胞外基质合成代谢，导致软骨组织[15]的破坏。本研究中，TNF-α处理后的软骨细胞IL-1α、IL-1β(表2)和IL-6、IL-9和 IL-10显著升高(表3)，提示TNF-α诱导OA的分子机制可能涉及炎症因子的调节。其次，基质金属蛋白酶系编码基因的表达(如MMP-3、9、12)同样由TNF-α显著调节(表3)，这进一步提示了TNF-α参与OA的发病机理。

氧化应激是自由基在机体产生的一种负面作用，是导致衰老和多种疾病的一个重要因素。越来越多的证据表明，线粒体功能障碍和线粒体DNA损伤在OA的发生发展中起重要作用[18]。Grishko等[19]研究发现OA患者软骨细胞中发现了线粒体DNA的损伤，活性氧自由基暴露后，线粒体DNA修复受阻，软骨细胞凋亡增加。本研究发现，经TNF-α诱导后，软骨细胞中涉及线粒体呼吸链相关的基因表达显著上调(表3)，印证了线粒体功能障碍在OA中的重要作用。

本研究发现小鼠原代软骨细胞经TNF-α处理后，大量细胞色素C氧化酶相关基因的表达显著上调(表3)，Wang等[20]发现可通过过氧化物酶体增殖激活受体γ(PPARγ)共激活剂1α可逆转OA中线粒体功能的衰减，此外，嘌呤受体基因P2ry1同样异常表达(表1)，这些都提示线粒体功能损伤在OA中可能有重要作用。此外，本研究发现一些神经功能相关的基因在本模型中的异常表达，如含嗅质蛋白(olfactomedin)结构域的基因Olfml1、Olfml3及嗅觉受体基因家族(olfactory receptor)，神经纤维瘤基因(neurofibromatosis 1，Nf1)都显著上调，这些全新的发现对OA中神经相关的病理机制提供了全新的线索。

综上所述，本研究利用小鼠基因表达芯片分析比较了TNF-α诱导的OA细胞模型基因表达的情况，并鉴定了一些DEGs和表达通路。此外，本研究还挖掘了一些新颖的OA相关的理论分子机制，后续这些病理假设还需要进一步探讨验证。