类风湿关节炎亚临床滑膜炎转录组学特点

孙晓莹，邓雪蓉，李光韬，谢文慧，樊勇，张慧娟，杨新蕾，张卓莉

随着达标治疗(treat to target, T2T)理念的提出，类风湿关节炎(rheumatoid arthritis, RA)患者的治疗目标为疾病缓解或低疾病活动度。在临床缓解的关节中，通过关节超声检查仍可以发现亚临床滑膜炎症[1]。亚临床滑膜炎的存在可预测RA复发及关节结构破坏进展[2]。

很少有研究关注RA亚临床滑膜炎相关的血清学生物标志物。先前研究发现，钙卫蛋白与人碱性成纤维细胞生长因子(basic fibroblast growth factor, bFGF)可作为在RA临床缓解患者中识别超声定义的活动性滑膜炎的生物标志物[3-4]。然而，上述标志物仍需进一步验证，目前尚未见关于RA亚临床滑膜炎转录组学变化的研究。在既往研究中，RA的转录组分析通常集中在外周血单个核细胞(peripheral blood mononuclear cell, PBMC)[5]。本研究旨在探索RA亚临床滑膜炎患者PBMC转录组学特点，研究亚临床滑膜炎发生机制。

1 材料与方法

1.1 研究对象

本研究前瞻性入选2019年5月至2019年7月就诊于北京大学第一医院风湿免疫科临床缓解的RA患者共10例，其中亚临床滑膜炎组及超声下缓解组患者各5例。定义所有关节超声下能量多普勒(power doppler, PD)分级均为0级的患者为超声下缓解组，任意一个关节PD≥1级的患者为亚临床滑膜炎组。

RA患者纳入标准：(1)符合1987年美国风湿病学会(American College of Rheumatology, ACR)或2010年ACR/欧洲抗风湿病联盟(European League Against Rheumatism, EULAR) RA分类标准；(2)年龄>18岁；(3)基于C反应蛋白的28个关节疾病活动度评分(disease activity score based on 28 joint count and CRP， DAS28-CRP)<2.6。排除标准：(1)合并其他关节炎性疾病；(2)存在关节明显畸形或残毁，或有缺如及外伤等。

本研究方案经本院伦理委员会审批，获得所有入组患者的口头或书面知情同意。

1.2 临床数据及全血标本收集

所有患者在同一天完成临床评估、实验室及超声检查，并用EDTA抗凝采血管采集所有患者外周静脉血10 mL。记录患者的临床资料，包括姓名、性别、年龄、28个关节(包括双侧肩、肘、腕、膝、掌指关节[metacarpophalangeal, MCP]1-5和近端指间关节[proximal interphalangeal, PIP]1-5关节)的肿胀关节计数(swollen joint counts, SJC)及压痛关节计数(tender joint counts, TJC)、患者整体评估(patient’s global assessment, PGA)、评估者整体评估(evaluator’s global assessment，EGA)以及实验室资料，包括红细胞沉降率(erythrocyte sedimentation rate, ESR)及C反应蛋白(C-reactive protein, CRP)。通过DAS28-ESR、DAS28-CRP、简化临床疾病活动指数(simplified disease activity index, SDAI)、临床疾病活动指数(clinical disease activity index, CDAI)等指标评估患者的疾病活动度。应用美国GE LOGIQ-E9型彩色超声波诊断仪和ML探头(6-15MHz)，分别从掌侧及背侧对28个关节及双踝、双跖趾关节(metatarsophalangeal, MTP)1-5关节共40个关节进行超声检查评估滑膜炎。灰阶(grey scale, GS)及PD超声半定量分级标准均按照Szkudlarek等[6]提出的0～3级半定量评分系统。

1.3 提取PBMC总RNA及检测RNA

从患者新鲜的外周静脉血中分离PBMC，用Trizol法进行细胞总RNA提取。应用Agilent 2100生物分析仪检测总RNA的浓度、RIN/RQN值、28S/18S及片断大小。

1.4 转录组学分析

使用BGISEQ平台测序，测序所得数据为原始测序数据，随后对其进行质控，以确定测序数据是否适合后续分析。质控后，经过滤得到过滤后的数据比对到参数序列。比对完，通过统计比对率、读长在参考序列上的分布情况，判断比对结果是否通过第二次质控。若通过，进行基因定量分析、基于基因表达水平的各项分析，筛选出两组间差异表达基因(differentially expressed genes,DEGs)，进行GO功能显著性富集分析、KEGG通路显著性富集分析、聚类、蛋白网络互作(protein-protein interaction,PPI)分析及转录因子深入挖掘分析。本研究使用DEseq2算法进行差异基因检测，采集差异倍数在两倍以上且校正P值≤0.05筛选差异基因。应用DIAMOND将DEGs对比至STRING数据库，利用与已知蛋白的同源性获得DEG编码蛋白的互作关系。

1.5 统计学方法

计量资料先进行正态分布检验，符合正态分布计量资料用均值±标准差表示，非正态分布的计量资料用中位数(四分位间距)表示；计数资料以率表示。组间计量资料均数比较采用两个独立样本的t检验(正态分布)。非正态分布计量资料的比较使用Wilcoxon符号秩和检验，计数资料两组间比较用卡方检验或Fisher精确检验。用 SPSS 21.0 软件进行统计分析。以P<0.05为差异有统计学意义。

2 结果

2.1 临床缓解RA存在亚临床滑膜炎及超声下缓解患者人口学资料及临床特征

临床缓解RA超声下缓解组与亚临床滑膜炎组患者人口学、实验室指标及疾病活动度间无显著差异(表1)。40个关节GS总评分在两组患者之间也无显著差异。亚临床滑膜炎组患者40个关节总PD评分为2(1，4.5)，超声下缓解组患者40个关节总PD评分为0(0，0)，两组间差异有统计学意义(P<0.05)(表1)。

表1 临床缓解RA患者亚临床滑膜炎组及超声下缓解组人口学及临床特征Table 1 Demographic and clinical characteristics of RA patients with ultrasound-defined subclinical synovitis and patientswith ultrasound determined remission

2.2 DEGs分析结果

本项目使用BGISEQ平台一共测了10个样本，亚临床滑膜炎组及超声下缓解组各5个，样本比对基因组的平均比对率为 96.43%，比对基因集的平均比对率为67.29%；一共检测到18 397个基因。按照|log 2 FC|>1.0及P<0.05的标准，与超声缓解组患者对比，亚临床滑膜炎组患者DEGs共 304个，其中106个基因上调，198 个基因下调。

2.3 PPI分析结果

从筛选出的DEGs 中，利用 STRING 平台构建PPI网络图(图1)。其中筛选连接节点数排序前20个枢纽基因(表2)， 包括2′-5′-寡腺苷酸合成酶(2′-5′-oligoadenylate synthetase, OAS)1、OAS2、OAS3、OASL、干扰素诱导蛋白(interferon-induced protein, IFI)44、IFI44L、黏液病毒/流感病毒抗性1(myxovirus resistance 1, MX1)、MMP9、干扰素诱导蛋白四肽重复序列(interferon-induced protein with tetratricopeptide repeats, IFIT)1、IFIT2、IFIT3、ISG15泛素样修饰剂(ISG15 ubiquitin like modifier, ISG15)、XAF1、干扰素调节因子7(interferon regulatory factor 7, IRF7)、含2个自由基的S-腺苷甲硫氨酸域(radical S-adenosyl methionine domain containing 2, RSAD2)、泛素特异性肽酶18(Ubl carboxyl-terminal hydrolase 18, USP18)、受体转运蛋白4(receptor-transporting protein 4, RTP4)、鸟苷酸结合蛋白1(guanylate-binding protein 1, GBP1)、CXCL10、DDX60。在RA亚临床滑膜炎组患者的PBMC中，除MMP9为上调基因外，其余均为下调基因。

图 1 PPI 网络分析图

表2 PPI网络图中连接点数排序前20个差异表达基因Table 2 Top 20 differentially expressed genes (DEGs)in PPI network

2.4 GO和KEGG富集分析结果

基于生物过程(biological process, BP)、细胞组成(cellular component, CC)、分子功能(molecular function, MF)的GO富集分析如图2所示，与超声下缓解组相比，亚临床滑膜炎组PBMC显示出多种富集，如BP类别中的免疫功能、Ⅰ型干扰素(interferon,IFN)组信号通路，CC类别中的细胞外间隙，MF类别中的丝氨酸型内肽酶活性及OAS活性。KEGG富集结果显示这些靶基因在抗原处理和递呈、IL-17信号通路及自然杀伤细胞介导的细胞毒性等通路高度富集(图3)，且这些富集通路中的DEGs 主要影响运输和分解代谢过程以及细胞生长和死亡过程。

图 2 GO生物过程、细胞组成及分子功能分析

图 3 KEGG通路富集分析差异表达基因

3 讨论

在RA患者的日常临床实践中，主要通过患者的体格检查及急性期反应物等临床指标评估疾病活动性[7-10]。然而，部分处于临床缓解状态的患者仍会出现疾病复发及放射学进展[1,11]。近年来，随着超声技术的发展，处于临床缓解的RA患者使用敏感的超声检查可发现临床检查及常规实验室检查不容易检测到的残留炎症，即亚临床滑膜炎[12]。多项研究证实，超声下亚临床滑膜炎与疾病复发及影像学进展密切相关[2,12]。与大多数高脂血症或糖尿病等慢性疾病不同，RA没有金标准的血清生物标志物。因此，鉴定RA亚临床滑膜炎的敏感血清生物标志物具有很大的挑战性。

本研究是第一个关于RA亚临床滑膜炎患者PBMC转录组学特点的研究。在达临床缓解的RA患者中初步筛选出可能与亚临床滑膜炎发生相关的20个DEGs，其中MMP9基因在RA亚临床滑膜炎组患者的PBMC中上调，而下调的DEGs主要为Ⅰ型IFN信号通路相关的基因。筛选出的上述DEGs主要与RA发病机制和药物代谢有关，可能是阐明RA亚临床滑膜炎分子发病机制的关键基因标志。

MMP9属于MMPs家族的明胶醇类，在滑膜被覆细胞以及滑膜下组织中浸润的炎性细胞中表达，可以溶解明胶、Ⅰ型、Ⅴ型和Ⅹ型胶原以及粘合素，而这些成分也是关节软骨的重要组成成分。此外，MMP9能够激活和调节与血管新生相关的生长因子。Stojanovic等[13]发现MMP9水平在RA患者血浆及滑液中均高于骨关节炎患者。Hakim等[14]研究显示MMP9基因的表达水平随DMARDs药物治疗而下降。本研究显示相对超声下缓解组患者，MMP9基因在RA亚临床滑膜炎组患者的PBMC中上调，提示MMP9表达上调可能参与亚临床滑膜炎的发生，是一个潜在的鉴别亚临床滑膜炎的血清学标志物。MMP9参与亚临床滑膜炎的发生可能与其激活和调节血管新生相关的生长因子等机制有关，目前尚不清楚，需要在未来的研究中进一步阐明。

本研究发现在亚临床滑膜炎组患者中，下调的DEGs大部分为Ⅰ型IFN相关的基因，主要参与Ⅰ型IFN介导的信号通路。IFN是细胞因子超家族，分为Ⅰ型和Ⅱ型，Ⅰ型包括IFNα和IFNβ。IFN均通过Janus激酶(JAK)信号转导子和STAT通路触发细胞内信号转导级联反应，从而在生物学和病理学条件下发挥多种免疫功能。Ⅰ型IFN已被证明起着将先天免疫和适应性免疫联系在一起的信号作用，其对免疫系统具有双重作用，表现为免疫刺激和免疫抑制功能[15]。Ⅰ型IFN在RA发病机制中的作用目前尚不清楚。van der Pouw Kraan等[16]学者曾发现RA患者外周血细胞的基因表达分析显示了大量与Ⅰ型IFN活性相关的基因过表达，Ⅰ型IFN介导的免疫也明显上调。而He等[17]提出Ⅰ型IFN相关基因可能是RA的特征基因，可能通过影响免疫反应而与RA发生相关。Ⅰ型IFN相关基因还与雷公藤、TNF抑制剂、IL-6抑制剂及CD20单抗等药物治疗反应相关，其表达水平的变化对不同药物治疗有所区别[5,18-20]。另外，Ⅰ型IFN与TNF-α之间存在交叉调节。IFNβ在体外以及关节炎动物模型体内均可下调TNF-α生成[21-22]。相反，TNF-α在体外抑制浆细胞样树突状细胞产生的IFN-β，在系统性幼年关节炎患者及干燥综合征患者体内TNF-α的阻断可导致Ⅰ型IFN通路激活增加[23-24]。总之，先前研究的结果表明Ⅰ型IFN相关基因可能是RA特征基因，可能通过影响免疫反应而与RA的发展相关。根据其免疫双重作用，它们在RA的作用可能有弊或有利。既往研究显示在RA患者的滑膜组织中可以检测到IFNβ蛋白，其中FLS是导致RA滑膜中IFNβ水平升高的原因[25-26]。另外，有学者提出RA滑膜中炎性细胞可能是从预先激活的外周细胞中募集的，而不是或不仅是通过局部产生[27]。因此，外周血PBMC中Ⅰ型IFN表达水平可反应局部滑膜组织中的表达水平。本研究Ⅰ型IFN相关的基因在超声下缓解组PBMC中相对上调可能暗示Ⅰ型IFN在RA中主要起到抗炎作用，亦或者由于Ⅰ型IFN高水平患者滑膜增殖主要以炎症活动为主，使其对药物治疗反应好，滑膜炎可很快消退，这值得进一步探讨。

此外，亚临床滑膜炎组患者下调的DEGs还包括CXCL10。CXCL10由IFN-γ刺激后分泌产生，具有强大的招募中性粒细胞、促进细胞因子分泌的功能，可介导炎症反应，通常在RA患者外周血及滑膜组织表达增高[28]。然而，CXCL10在RA整个病程中是否持续发挥作用及其与RA疾病活动度及骨破坏等之间的关系目前尚未阐明。有研究发现，CXCL10在部分RA患者的滑液和滑膜组织中水平很高，但血清中几乎检测不到CXCL10[29]。Kraan等[30]研究显示，RA患者滑膜组织中CXCL10在临床受累及临床未受累关节中的表达无明显差异，表明CXCL10水平与RA临床活动性可能无明显相关。Strieter等[31]体外和体内研究均显示，CXCL10是一种血管生成抑制剂, 其作用机制可能与清除一些其他的生长因子, 如IL-8和bFGF等有关。提示在本研究中CXCL10基因水平下调与亚临床滑膜炎发生有关可能是由于其对血管形成相关因子的抑制作用下降有关，值得进一步深入研究。

本研究存在一些局限性：第一，本研究纳入的样本量小，两组人群性别、年龄等基本信息存在一定差异，可能导致偏倚。第二，两组人群均为RA达临床缓解的患者，目前疾病均处于相对稳定状态。部分超声下缓解组患者仍有轻度滑膜增生，使得两组患者之间差异较小。第三，本研究采用全外周血单个核细胞，未区分细胞群。白细胞亚群的动态变化也可能会混淆结果[32]。第四，本研究为初步探索性研究，仅研究了亚临床滑膜炎相关的转录组学，而mRNA水平未必能代表蛋白水平，本研究结果有待在基因及蛋白水平上进一步验证。另外，由于RA亚临床滑膜炎患者的滑膜组织及滑液不容易获取，本研究仅对患者PBMC进行转录组学分析。未来将收集患者滑膜组织及滑液进行转录组学检测。第五，本研究仅纳入临床及超声均缓解的RA患者作为对照，缺乏健康人群及活动期RA患者作为对照，未来研究将纳入这两组人群作为对照组进行验证，探索特异于RA亚临床滑膜炎患者的血清学标志物。

综上所述，在RA达临床缓解的患者，外周血MMP9升高以及Ⅰ型干扰素相关通路、CXCL10下调可能提示亚临床滑膜炎存在的潜在生物学标志物，或许可以为揭示亚临床滑膜炎的发生提供新的视角和潜在的新靶点。但其分子发病机制目前尚不清楚，仍需在基因及蛋白水平进一步验证，并深入探索。