软骨细胞单细胞转录组数据的再挖掘揭示人骨关节炎病程中转录调控因子和代谢的时序变化

施剑明,杜桂花

(1.景德镇市第一人民医院创伤骨科,江西 景德镇 333000;2.南昌大学公共卫生学院劳动卫生与毒理学教研室,江西省预防医学重点实验室,南昌 330006)

人骨关节炎(human osteoarthritis,hOA)是导致老年人运动能力丧失的主要原因之一,给社会和家庭造成巨大的经济损失[1-2]。诱导hOA发生的因素众多且关系复杂,有研究[3]表明,代谢、遗传、衰老和损伤等都与hOA的发生发展相关。目前,临床对hOA的治疗主要以非类固醇类抗炎治疗和物理治疗为主,但都无法逆转疾病进程,人工关节置换通常是患者的最终选择[4]。探究hOA的发病及进展机制,可为逆转病情进展治疗方案的开发提供方向和数据支持,是目前的研究热点。

临床上,根据国际软骨修复协会(International Cartilage Repair Society,ICRS)软骨损伤分级系统和国际骨关节炎研究协会(Osteoarthritis Research Society International,OARSI)评分标准,骨关节炎处软骨细胞受损程度可分为由轻微到严重的4个等级:正常无损伤或软骨表面完整处为0级(Stage 0,S0);表浅的钝性损伤或软骨表面不平整处为I级(Stage I,S1);软骨损伤小于软骨厚度一半或软骨表面不连续处为Ⅱ级(Stage Ⅱ,S2);软骨损伤大于软骨厚度一半但未达软骨下骨或软骨出现垂直裂变处为Ⅲ级(Stage Ⅲ,S3);软骨损伤可见软骨下骨暴露或出现软骨侵蚀处为Ⅳ级(Stage Ⅳ,S4)。目前,尽管大量研究[5-7]报道了hOA的病理表现和潜在分子机制,但并未对其病情进程中的分子动态变化进行探索。因此,对hOA进展过程中的调控网络进行研究可为临床阻止或逆转hOA进程提供理论基础。

单细胞转录组测序(single-cell RNA sequencing,scRNA-seq)利用有限数量的样本在单细胞水平上获得转录组信息,可为发育生物学的研究和生理病理机制的探索提供重要信息[8-9]。JI等[10]通过scRNA-seq分析了hOA患者软骨细胞的内部异质性及分子特征,但并未从hOA病程分级角度上分析疾病进展过程中基因表达及其生物学过程改变。而本研究通过对人膝关节软骨细胞的单细胞测序数据进行再挖掘,从单细胞水平上对hOA各发展阶段中基因及生物学过程变化进行连续性的研究,以期为临床开展hOA的阶段针对性治疗和新疗法的开发建立基础。

1 数据和方法

1.1 原始数据来源

人软骨细胞单细胞测序数据集(序列号:GSE104782)[10]下载于GEO数据库。该数据集中软骨细胞样本来源于10例膝关节置换手术患者的膝关节,依据每位患者关节损伤情况分别采集S0—S4分级的软骨细胞进行单细胞转录组建库及测序。

10例患者中,男3例、女7例,平均年龄(66.50±1.86)岁,体重(66.90±2.44)kg,美国特种外科医院膝关节评分(hospital for special surgery knee score,HSS)(60.20±2.56)分,具体见表1。

1.2 数据挖掘方法

基于R软件对单细测序数据集进行再挖掘。

1.2.1 质量控制

进行数据分析之前,利用单细胞检测mRNA分子总数(nCount_RNA)和基因数量(nFeature_RNA)剔除测序深度不足或因细胞粘连导致测序质量不佳的单细胞测序数据。据图用于后续分析的准入标准为:“2000

1.2.2 数据分析和可视化

基于R 3.6软件,通过软件包Seurat 3.0和Monocle 2进行数据分析和可视化的实现,其中在Seurat 3.0中进行测序数据的初始标准化和归一化,在Monocle 2中进行细胞拟时轨迹分析。通过软件包CellCycleSoring进行细胞周期检测,通过DAVID Bioinformatics Resources 6.8(https://david.ncifcrf.gov/)进行差异基因功能(gene ontology,GO)分析,通过PANTHER分类系统(http://pantherdb.org/)寻找差异基因中编码转录调控因子的基因。

使用GraphPad Prism 8.0、R 3.6和Adobe Illustrator CC 2018等软件进行图表制作。

1.3 统计学方法

2 结果

2.1 单细胞测序数据概况

对原始测序数据进行质控分析后,将数据过滤指标设定为“2000

对筛选后的细胞进行标准化和归一化分析后,进行均匀流形近似和投影(uniform manifold approximation projection,UMAP)降维分析(图1A)。在此基础上进一步验证软骨细胞已知标志基因(如ACAN、COL2A1、ITM2A和VEGFA)的表达(图1B),结果明确了筛选后的细胞均为软骨细胞。虽然来自于10例患者的软骨细胞并未各自成簇,但样本间仍有异质性(图1C)。为检测单细胞处理的批次差异效应,通过修正Pearson相关系数评定各样本间相似程度[11],结果显示Pearson相关系数范围在0.899～0.983,样本间虽存在异质性,但一致性程度较高,可合并用于后续分析(图1D)。

A:软骨细胞scRNA-seq文库UMAP;B:软骨细胞标志蛋白UMAP;C:10例hOA患者软骨细胞scRNA-seq文库UMAP;D:10例患者文库间Pearson相关系数热图。

2.2 hOA发展进程中软骨细胞的拟时轨迹分析

利用Monocle 2从单细胞水平上拟合hOA进展轨迹,并对随轨迹变化的相关基因进行分析,以从单细胞水平上分析hOA发展的调控网络。

各分级软骨细胞UMAP结果显示,虽然同一分级内的细胞聚集成簇且S0至S4排列具有一定方向性,但并非所有细胞都按该规律排布(图2A)。单细胞测序数据来自10例患者,虽然样本间相关性分析提示样本间的一致性较好,但不同分级软骨细胞的采集是依据关节软骨病损程度进行的。为明确各样本分级,有学者[12-13]检测了与骨关节炎进展呈正相关的细胞外基质降解相关基因,结果显示,ADAMTS12、ADAMTS4和ADAMTS16在hOA患者软骨细胞中的表达升高,且S4软骨细胞中以上基因表达水平均最高(图2B),这表明从hOA进展的角度看各样本一致性仍较好,可基于此进行拟时轨迹分析(pseudo-time trajectory analysis)。

A:不同分级(S0-S4)软骨细胞UMAP;B:细胞外基质降解相关标志基因在不同分级软骨细胞中的表达水平,其中点的大小表示基因表达的细胞百分比(Percent Expressed),灰度梯度表示基因均化表达量(Average Scaled Expression);C:各分级软骨细胞沿着拟时轨迹轴排列;D:沿拟时轨迹轴变化的差异基因及其表达热图;E:各组差异基因GO分析部分条目;F:各组差异基因KEGG通路分析;G:根据各组差异基因及其表达峰值总结hOA各发展阶段生物学过程。

利用各分级细胞前1500个差异基因(基于P值)将细胞进行排序,除个别S4的细胞排列在轨迹一端起点外,其余细胞根据其分级很好地依拟时轨迹轴排列,表明拟时轨迹分析模型构建较成功(图2C)。在此基础上,进一步将沿着拟时轨迹轴变化的差异基因(P<0.01)及其表达量拟合成热图(图2D),结果显示差异基因可分为4个组,组1的基因表达在拟时轨迹轴最开始处表达最高,组2的基因在拟时轨迹轴前中段表达最高,组3基因在拟时轨迹轴后中段表达最高,而组4基因在拟时轨迹轴后段表达最高。

对各组差异基因进行GO分析和KEGG分析,细胞对TNF的应答/TNF信号通路、对胰岛素的应答/MAPK信号通路、MAPK级联反应/PI3K-Akt信号通路分别出现在组1,组2和组3,而在组4中则出现细胞黏附和蛋白的消化、吸收过程(图2E—F)。差异基因的表达峰值呈现从组1至组4的前后排列,因此基于差异基因可归纳总结出hOA发展4个阶段的相关生物学过程(图2G):在hOA发展的早期阶段,其生物学过程主要涉及细胞因子的产生和一氧化氮的生物合成等;而hOA发展的晚期阶段,关节软骨组织中出现血管生成、胶原纤维形成和基质破坏等;炎症反应则贯穿hOA进展的整个过程。

2.3 调控hOA发展的转录调控因子

包括RUNX2、MAFF和ATF3等在内转录调控因子的在hOA的发生和发展过程中起重要作用[14-16]。为进一步从转录调控水平解析hOA进展中的关键因子,本研究利用PANTHER分类系统对随着hOA软骨细胞拟时轨迹变化的各组基因进行了分析(图3)。已知的对hOA存在关键调控作用的转录因子ATF3、MAFF和RUNX2分别出现在组1,组2和组4。而FOXO3、KLF10、KMT2C和ZMYM2在hOA的发展初期表达最高(图3A),转录调控因子HES1、ID1、ID3、JUN、KLF2和MEF2C表达则随hOA拟时轨迹先下调,随后又在拟时轨迹后期表达上调(图3B);转录调控因子ETS2、FOXA3、JUND、NFATC1和SOX8在hOA拟时轨迹的前中段表达最高(图3C),而CITED4、EGR2、ETV5、FHL2和SOX11则在拟时轨迹的中后段表达最高(图3D);转录调控因子AEBP1、CKLF和STAT1在hOA拟时轨迹后期表达上调,这可能与hOA晚期软骨组织自我修复可能性的彻底丧失有关(图3E)。

A—B:组1;C:组2;D:组3;E:组4。图3 hOA软骨细胞中各组关键转录调控因子沿拟时轨迹轴的表达模式

2.4 hOA发展过程中软骨细胞的能量代谢改变

对沿hOA进展拟时轨迹变化的软骨细胞差异基因进行KEGG分析,结果显示氧化磷酸化相关基因出现在组4中,且相关基因的表达随hOA进展升高(图4A)。关节软骨是一种无血管组织,这意味着软骨细胞对营养物质的摄取及其代谢产物的排出主要依赖于软骨表面的渗透,而软骨细胞主要通过糖酵解过程获得能量[17]。细胞能量代谢与细胞周期相关,增殖活跃的细胞(如肿瘤细胞)更倾向于利用糖酵解获取能量[18]。为排除hOA软骨细胞中能量代谢改变是由于细胞增殖变化导致的,本研究利用CellCycleSoring对软骨细胞的细胞周期进行了分析。结果表明,随着hOA进展,处于增殖状态(S/G2/M期)的软骨细胞百分比明显减少(图4B),排除了hOA进展后期软骨细胞氧化磷酸化水平增加是由细胞增殖活力改变而引起的可能性。考虑到个别S4细胞同大部分S0细胞一样排列于拟时轨迹轴一端起点,将拟时轨迹分析差异基因中与氧化磷酸化或糖酵解有关基因的均化表达量根据hOA进展分级制成热图,并对其均化表达量进行累计以直观地比较。结果显示,与软骨细胞中氧化磷酸化相关基因的表达水平随hOA进展升高(图4C—D),而糖酵解相关基因表达水平随着hOA进展而降低(图4E—F)。

D:各级细胞氧化磷酸化相关基因均化表达累计值柱状图;E—F:各级细胞糖酵解相关基因均化表达水平热图及均化表达累计值柱状图。

A:氧化磷酸化相关基因沿拟时轨迹轴的表达模式;B:各级软骨细胞处于S/G2/M期的百分比;C:各级细胞氧化磷酸化相关基因均化表达水平热图。

3 讨论

关节软骨是一种相对简单的组织,仅由软骨细胞组成。尽管如此,大量研究[10,19-20]揭示了软骨细胞内部异质性的存在。scRNA-seq通过高效降维和拟时轨迹分析,可精确地对细胞内部异质性进行识别,在细胞谱系的建立和疾病发生发展机制的探索中发挥重要作用。有相关研究[10,21]利用scRNA-seq对hOA软骨细胞内部异质性及其分子特征进行了探索,但对病情进展中基因表达水平及关键细胞生物学过程的变化进行分析。因此,本研究以骨关节炎连续性分级为基础,联合单细胞测序拟时轨迹分析,挖掘并总结归纳出hOA病情进展中4类不同的生物学过程(图2G),并对hOA病情进展中转录调控因子的变化趋势进行了解析。此外,本研究还发现,随着hOA病情进展软骨细胞代谢发生改变,表现为软骨细胞氧化磷酸化水平加强,而细胞糖酵解程度降低。

骨关节炎是一种以进行性软骨退化和关节炎症为特征的疾病[22-23]。除炎症反应外,一氧化氮的生物合成、软骨细胞分化、骨矿化、血管生成和入侵等生物学过程都与骨关节炎相关[22,24],本研究的分析结果也支持这一结论。与既往体内动物实验或体外细胞实验不同的是,本研究基于骨关节炎分级进行了scRNA-seq数据再挖掘,不仅揭示了hOA发生发展过程中的相关生物学过程,还描绘出不同生物学过程发生的先后顺序及交叉现象,为临床开展针对性治疗以阻止或减缓hOA病情进展提供了相对精准的时间轴。

转录因子通过调控基因表达控制着细胞周期、细胞分化、细胞衰老和代谢平衡等重要生物学过程[25]。在hOA的发生发展过程中,转录因子调控因子的作用同样至关重要。如RUNX2调控成骨细胞和软骨细胞分化过程,在hOA发展过程中发挥重要作用[26],而SOX11过表达与软骨损伤的恶化有关[27],本研究发现RUNX2或SOX11在关节软骨损伤严重的后期表达显著上调。IEZAKI等[16]发现转录因子ATF3在hOA组织中表达上调,且ATF3敲除可减缓OA进展,但本研究结果提示ATF3的表达随着OA进展而下调。产生这一矛盾的原因可能是组织与单细胞分析的异质性。YANG等[28]对不同时期hOA的差异基因进行分析,发现ATF3随疾病进展表达下调,与本研究结论一致的。此外,除了已知的转录调控因子(如RUNX2、ATF13、SOX11和STAT1)[14-16,27,29],本研究发现了大量此前未知的在hOA进展中发生表达水平变化的转录调控因子,它们可能对疾病发展有调控作用。相较于传统的基于对照组和hOA疾病组比较发现的差异转录因子,本研究利用患者软骨细胞scRNA-seq数据发现的hOA不同发展阶段的差异转录因子对hOA阶段性靶向治疗策略的研发有更加积极的意义。

除了生物学过程和转录因子,本研究还发现在hOA的发展过程中伴随着细胞代谢的变化。细胞代谢不仅仅是细胞能量的来源,它还与细胞信号通路和表观遗传调控密切相关。此外,细胞代谢的变化也满足了干细胞自我更新和分化等多方面的需求[30-31]。hOA的最终结果通常被认为是软骨祖细胞失去再生能力,逐渐朝向终端分化并失去自我修复能力[32]。这一现象是否与软骨细胞的代谢改变相关目前尚不清楚。有研究[33]表明,hOA患者软骨细胞中存在氮氧化物,这些氮氧化物会耗尽氧化磷酸化所需的物质。此外,炎症过程中产生的活性氧分子会影响呼吸链的运行,导致线粒体活性下降,从而影响软骨细胞的氧化磷酸化水平[34]。还有研究[35]表明,hOA中软骨细胞功能失调导致PFKFB3的表达下调,而PFKFB3的表达水平则与糖酵解代谢强度呈正相关。以上证据表明,随着hOA的发展,氧化磷酸化相关基因的表达水平上调,而糖酵解相关基因的表达下调。另一项研究[36]发现,线粒体DNA的单倍体G变异将导致细胞代谢从糖酵解转向氧化磷酸化,从而增加了hOA的风险。此外,hOA后期阶段通常发生软骨组织中血管生成,这会导致软骨细胞所处环境中氧含量的升高,促使软骨细胞从糖酵解代谢向氧化磷酸化代谢转变。而相对于氧化磷酸化代谢,糖酵解代谢更有利于维持干细胞的自我更新能力[37]。本研究结果表明,hOA的发展伴随着软骨细胞中氧化磷酸化代谢的增强和软骨祖细胞干性的逐渐丧失,这可能是hOA病程难以逆转的关键原因之一,虽然这一假说尚需证据支持,但它表明代谢干预可能是阻遏或逆转hOA进展的途径,可作为后续临床研究的方向。

综上所述,基于hOA分级,本研究对软骨细胞单细胞测序数据进行了再挖掘,探讨了随hOA病情进展的转录调控因子变化情况并归纳总结出hOA病情进展中的4种不同生物学过程演变,发现了疾病进展过程中伴随的软骨细胞代谢变化。本研究结果有望为临床分阶段的针对性治疗hOA提供了基础,同时也为新疗法的开发指明了方向。