类风湿关节炎滑膜组织线粒体功能相关差异表达基因的分析

巴燕娜，袁 晴，贺 彤，穆 楠，吴振彪，袁志栋,3※

(1.第四军医大学西京医院临床免疫科，西安 710032； 2.第四军医大学基础医学院生理与病理生理学教研室，西安 710032；3.赣南医学院基础医学院，江西 赣州 341000)

类风湿关节炎(rheumatoid arthritis，RA)是一种以关节滑膜炎性病变为主的系统性自身免疫性疾病，主要表现为慢性、对称性和进行性的多关节炎[1]。该病影响了近1%的人群，女性发生率是男性的2～3倍[2]。滑膜组织中成纤维细胞样滑膜细胞在RA中显示出异常的行为，被认为与关节损害有关[3-4]，这种细胞对凋亡的抵抗性是RA的一个主要特征[5]。因此，系统分析RA滑膜组织细胞的基因表达水平，不但可加深对RA病理性过程及临床症状的理解，而且能更好地理解药物治疗的作用机制[6]。

线粒体是真核细胞中的一种重要细胞器，拥有独立的基因组DNA，共37个基因。这些基因和核基因组编码的定位于线粒体内蛋白的基因对线粒体的氧化磷酸化、能量转换等功能活动是必需的。RA患者高度失调的微血管结构造成了缺氧微环境，异常的滑膜细胞代谢和线粒体功能障碍随之而来，增加活性氧类，诱导炎症[7]。另外，在细胞应激或损伤时，线粒体成分可被释放到胞质或细胞外，可作为危险信号发挥作用[8]。因而，研究线粒体相关基因的异常表达对探讨RA的发病机制非常关键。

转录组芯片分析已被用来研究RA病理发生机制、治疗效果和生物标志物[9-13]。虽然线粒体与炎症、免疫密切相关，但线粒体在RA病理发生机制中的功能一直未受到研究者的重视。为探索这些问题，本研究通过系统分析RA与正常人滑膜组织的基因表达芯片数据，以期确定RA病理发生与线粒体之间的关系及鉴定关键的影响滑膜组织病变的线粒体相关基因。

1 材料与方法

1.1滑膜活检组织的公共芯片数据 RA滑膜活检组织的基因表达数据下载于美国国家生物技术信息中心的GEO数据库的芯片数据，ID为GSE77298。该数据集包含16例末期RA患者和7例无关节疾病的受试者的滑膜活检组织，提取总RNA后在Affymetrix芯片上检测，然后获得滑膜组织的基因表达数据[10]。但由于该芯片设计上未考虑线粒体基因组编码的基因，因此只能检测线粒体相关基因而不能检测线粒体编码基因的表达变化。

1.2分析方法

1.2.1基因表达芯片数据的分析 下载芯片数据后，用R软件安装simpleaffy、affyPLM和limma等程序包后对芯片基因表达数据依次处理，利用主成分分析、热图聚类分析评估样本间的总体相似度，剔除错误分组的样本，鉴定差异表达基因，绘制差异表达基因的火山图和热图。

1.2.2线粒体相关基因功能注释与蛋白网络分析 本研究定义核基因组编码但其蛋白产物定位于线粒体，进而与线粒体结构、功能相关的基因为线粒体功能相关基因。线粒体功能相关基因信息来源于MitoMiner(MitoMiner 4.0v2018JUN)数据库[14]，该数据库收录了1 626个基因，包括人的1 184个已知其蛋白产物定位于线粒体的基因和442个预测其蛋白产物定位于线粒体的基因。将其与鉴定的全转录组差异表达基因列表取交集，即提取显著性差异表达的线粒体功能相关基因。差异表达的线粒体相关基因的功能注释与蛋白质网络分析采用在线分析工具Metascape[15]，而Metascape应用了图论聚类算法——分子复合物检测MCODE(Molecular Complex Detection)[16]。

2 结 果

2.1滑膜组织差异表达的基因 首先对芯片数据进行标准化，并根据主成分分析和热图聚类分析结果评估样本间的总体相似度，发现RA1、RA2、RA9和RA15样本数据与总体的差异较大，进而剔除这4例样本，并分析基因差异表达，获得1 986个在RA滑膜组织中显著差异表达的基因(校正后的P<0.05且表达值变化2倍以上)，其中表达水平2倍以上变化的上调基因有1 490个，下调的基因有496个，见图1、图2。这些差异表达基因中的IGLC1、SPP1、MMP1、JCHAIN和CHI3L2为表达水平上调倍数变化6以上的基因；PLIN1和ADIPOQ(脂联素)是表达水平下调倍数变化6以上的基因。

Up：上调(红色)；Down：下调(绿色)；logFC：RA患者与正常人标准化后的基因表达值比值，即以2为底倍数变化的对数；-log10(adj.P.Val)：以10为底校正后P值的对数的相反数；NotSig：无显著变化

图1 基因表达的火山图

HC1-7：正常对照；RA：类风湿关节炎患者；蓝色：基因表达下调；红色：基因表达上调

图2 差异表达基因的层级聚类热图

2.2差异表达的线粒体功能相关基因 本研究结果显示MitoMiner 4.0数据库收录的1 355个线粒体功能相关基因在芯片表达数据中有表达，但由于芯片设计的局限性，线粒体基因组编码的基因未检测到表达。在前面鉴定的1 986个显著差异表达基因中发现有169个线粒体功能相关基因，占目前所知的蛋白产物定位于线粒体的基因总数量的10.39%，其中表达上调的基因数为131个，下调的基因数为38个。

2.3差异表达线粒体相关基因的功能注释 为了解169个差异表达线粒体相关基因的功能，对这些基因进行注释后获得了以条形图形式列出了P值排列前20的富集通路，离散色阶代表不同的统计学差异。P值最显著的通路是核苷酸代谢过程通路(GO:0009117)，其次是蛋白质在内质网中的蛋白质加工(KEGG通路：HSA04141)和线粒体结构组织通路(GO:0007005)。另外，还包括凋亡信号的调控(GO:2001233)、内质网应激反应(GO:0034976)、氧化应激反应(GO:0006979)和活性氧类代谢过程(GO:0072593)通路等(图3A)。差异表达线粒体相关基因的功能通路网络中颜色代表了聚类簇的ID，聚类簇ID与图3A富集的通路一致，每个节点代表了一个富集的术语，但功能通路网络聚类簇更体现了富集的术语间的功能相关性(图3B)。这些异常表达线粒体功能相关基因的功能富集通路体现了RA相关功能通路出现了扰动。

2.4差异表达的线粒体相关基因的蛋白相互作用网络 利用Metascape生成了差异表达的线粒体相关基因的蛋白质互作网络(图4A)。Metascape检测的蛋白-蛋白相互作用网络的密集连接区，它们代表了分子复合物。每一个MCODE成分分别进行了通路和过程富集分析，这些MCODE成分中根据P值打分最高的术语是KEGG通路HSA04141[内质网蛋白质加工，log10(P)=-17.1]、GO:0009117[核苷酸代谢过程，log10(P)=-14.1]和GO:0006753[核苷磷酸盐代谢过程，log10(P)=-14.0]。蛋白质互作网络中4种不同的颜色表示识别到的4个MCODE成分(图4A、图4B)，分别是MCODE_1：Reactome 通路R-HSA-1799339[靶向膜的SRP依赖的共翻译蛋白，log10(P)=-11.1]、Reactome通路R-HSA-72766[翻译，log10(P)=-8.6]和HSA04141[内质网蛋白质加工，log10(P)=-4.2]；MCODE_2：Reactome通路R-HSA-446203 [天冬酰胺N-连接糖基化，log10(P)=-9.5]、KEGG通路HSA04141[内质网蛋白质加工，log10(P)=-8.8]和KEGG通路HSA00510[N-聚糖生物合成，log10(P)=-6.2]；MCODE_3：Reactome通路R-HSA-382551[小分子转运，log10(P)=-3.1]和MCODE_4：GO:0050878[体液水平的调控，log10(P)=-3.8]。

4A：线粒体相关基因的蛋白质互作网络图； 4B：蛋白质网络中密集连接的蛋白质组

图4 差异表达的线粒体相关基因编码的蛋白质互作用网络

3 讨 论

线粒体在细胞内不仅参与物质的代谢、能量转化和钙离子存储等生命活动，还在病理状态下由损坏的线粒体导致先天性免疫系统异常激活而引起自身性免疫疾病[17]。然而，在RA中的线粒体功能一直未引起研究者的重视。为了探究线粒体功能相关的基因在RA中的表达变化和功能，本研究分析了GEO数据库基因表达芯片数据GSE77298，获得了1 986 个在RA滑膜组织中显著差异表达的基因，这些基因与DisGeNET数据库[18]收录的1 832个RA相关基因的交集为322个，与Liu等[19]报道的RA患者中显著差异表达上调的228个(实际是227个)基因和下调的36个基因相比，相一致的上调基因110个，相一致的下调基因10个；另一项研究报道在滑膜组织中RA组比正常组高表达的基因1 018个，下调的基因893个[20]，与本研究结果相同的分别为686个和257个。另外，在RA滑膜组织中显著差异表达、倍数变化6以上的上调基因SPP1、MMP1在以前的研究中[21-23]也得到了验证。而ADIPOQ是倍数变化6以上的表达下调基因，Harshan等[20]的研究结果显示脂联素在RA患者滑膜组织中表达水平较正常人显著降低，本研究结果与上述文献一致。在这1 986个显著差异表达基因中，本研究鉴定了169个与正常对照存在统计学差异表达的线粒体功能相关基因。而这些线粒体功能相关基因部分也得到了DisGeNET数据库[18]及RA相关研究报道[19-20,24-28]的数据支持与验证。169个线粒体相关基因的25个与DisGeNET数据库中的一致，其中FHL1、G0S2和SLIT3为表达下调的3个基因，其他的22个为上调基因。SLIT3能够抑制RA中Robo3诱导的滑膜成纤维细胞的入侵[24]，因此SLIT3的表达下调可能促进RA的病情恶化。另外，G0S2可作为生物标志物用来预测响应抗肿瘤坏死因子治疗RA患者[27]。而表达上调倍数变化最大的为TYMS，其遗传多态性与RA患者的甲氨蝶呤临床治疗的效果相关[25-26]。

本研究发现的部分差异表达基因不仅与已往研究报道[18-28]一致，而且通过功能通路注释和互作网络分析，结果显示GO:0009117(核苷酸代谢过程)和KEGG通路HSA04141(内质网蛋白质加工)不仅分别是差异表达的线粒体相关基因功能富集通路中P值最显著的通路之一，也分别是蛋白互作网络MCODE成分中根据P值得分最高的生物学过程术语之一。将蛋白质互作数据与通路和生物学过程富集数据相结合显示了差异表达的线粒体功能相关基因的蛋白质互作网络可能在RA中存在异常。而有报道发现青蒿素能通过抑制PDK1(GO:0009117共有34个差异表达的线粒体功能相关基因之一)诱导的蛋白激酶B和核糖体S6激酶家族成员2磷酸化激活来抑制RA成纤维样滑膜细胞的迁移和入侵[28]。

另外，富集的通路还包括线粒体结构组织通路(GO:0007005)、凋亡信号的调控(GO:2001233)、内质网应激反应(GO:0034976)、氧化应激反应(GO:0006979)和活性氧类代谢过程(GO:0072593)通路，以及蛋白质互作网络中识别到的MCODE成分R-HSA-1799339 (靶向膜的SRP依赖的共翻译蛋白)等。以上功能与通路的富集注释结果以及已有相关研究报道[18-28]的验证显示了这些线粒体功能相关基因的异常表达可能与RA滑膜组织的病理变化有关。

需要指出的是，本研究分析的是芯片数据，由于该芯片探针未覆盖线粒体基因组编码的基因，所以分析结果仅鉴定了核基因组编码的线粒体功能相关基因，并且这些基因均为编码蛋白基因。因此，后续更完整的研究应囊括1 626个定位线粒体内的相关编码基因、线粒体基因组编码的基因，以及非编码RNA基因，以期能够全面系统地揭示线粒体对RA的发生、发展及滑膜组织病变所起的作用，找到有效的RA疾病诊断分子标志物和治疗的药物靶标。