转录组分析鉴定腰椎椎间盘退行性变中潜在免疫相关生物标志物

李志超，王 磊，薛景才，王文波，陈 煜，齐军强，许 鹏

1.山东中医药大学第二附属医院骨科，济南 250002

2.山东中医药大学研究生院，济南 250001

3.海军军医大学研究生院，上海 200433

4.海军军医大学长征医院骨科，上海 200003

以椎间盘退行性变（IDD）为特征的椎间盘疾病，给患者带来了生理和心理上的痛苦，给家庭和社会造成巨大的经济负担［1］。IDD 的致病因素包括营养物质减少、基质降解、细胞凋亡、炎性因子增加、免疫应答、遗传因素和机械负荷等［2］，但其具体机制目前仍不清楚。IDD 与免疫之间有很强的相关性［3］，椎间盘的驻留细胞（纤维环和髓核）和浸润的免疫细胞分泌的炎性细胞因子是引发IDD 的关键因素［4］。发生退行性变的椎间盘使巨噬细胞倾向于促炎性极化，这似乎加剧了IDD［5］。

竞争性内源RNA（ceRNA）是一种新的基因表达调控模式。在调控过程中，一些长链非编码RNA（lncRNA）和微RNA（miRNA）结合，lncRNA 起到miRNA 海绵的作用，降低miRNA 对其靶基因的抑制作用，使靶基因在细胞中的表达量增加［6］。有研究［7-8］表明，IDD 的病理变化与lncRNA、miRNA 及其靶点的失衡以及许多不同细胞生物过程有关。失衡的ceRNA 相互作用轴可能是治疗IDD 的关键目标［9］。本研究通过加权基因共表达分析（WGCNA）构 建IDD 中lncRNA-miRNA-mRNA 的ceRNA 网 络，筛选与IDD 免疫相关的关键基因，为寻找IDD 的新靶点和治疗策略提供参考。

1 材料和方法

1.1 数据来源

从基因表达汇编（GEO）数据库（https：//www.ncbi.nlm.nih.gov/geo）GSE67567 数 据 集 下 载IDD 人类退行性髓核组织和对照的髓核组织微阵列芯片数据，GSE67567 数据集包含GSE56081（lncRNA 和mRNA 数据）和GSE63492（miRNA）2 个子集，各包括来自5名健康志愿者和5例IDD患者的髓核样本。

1.2 数据预处理和差异表达分析

使用GEO 的在线分析工具GEO2R 获得标准化的差异表达基因。根据GPL15314、GPL19449文件和GEO中的注释信息，将探针转换为相应的基因符号。对于差异表达的miRNA，将阈值设置为|log2FC| ＞ 2、调整P＜ 0.05；对于差异表达的lncRNA，将阈值设置为|log2FC| ＞ 1、调整P＜ 0.05。

1.3 共表达模块建立和基因集富集分析（GSEA）

使用Image GP工具（https：//www.ehbio.com/Cloud Platform/front/#/）分 析GEO 中的mRNA 表达矩阵。首先构建差异表达基因的共表达网络，将其转换为拓扑重叠矩阵，通过分层聚类和动态树切割算法，使用模块特征基因和模块关系来确定与表型相关的共表达模块，最小模块大小为30［10］。使用Enrichr数据库（https：//maayanlab.cloud/Enrichr/）分析每个模块中的mRNA，进行京都基因与基因组百科全书（KEGG）通路和基因本体（GO）生物过程的GSEA，筛选出与免疫调节关系最为密切的模块。

1.4 ceRNA 网络建立和GSEA

通过Targetscan 数据库（https：//www.targetscan.org）和miRTarbase 数据库（http：//mirtarbase.mbc.nctu.edu.tw/php/index.php）筛 选miRNA 靶 向 的mRNA，将这些mRNA 与上述所选模块中的差异表达mRNA 取交集获得免疫相关mRNA。同理，通过StarBase 数据库（http：//starbase.sysu.edu.cn/index.php）预测差异表达miRNA 调控的lncRNA。最后，使用Cytoscape 3.7.2建立lncRNA-miRNA-mRNA ceRNA网 络。使 用Enrichr 对ceRNA 网 络 中 的mRNA 的KEGG 通路和GO 生物过程进行GSEA。

1.5 免疫浸润分析

使用ImmuCellAI（http：//bioinfo.life.hust.edu.cn/web/ImmuCellAI/）［11］比较GEO 数据库中IDD 组和对照组的免疫浸润特征。ImmuCellAI 是一种基于转录组测序或芯片工具估计免疫细胞浸润丰度数据的工具，可以判断18 种T 细胞和6 种其他类型免疫细胞的比例。

1.6 与CD4 幼稚T 细胞（CD4 naive）相关的ceRNA子网络的构建

使用Pearson 相关分析确定ceRNA 网络中CD4 naive 浸润与mRNA 和lncRNA 表达的相关性，以相关系数绝对值＞ 0.8、P＜ 0.01 为条件筛选mRNA和lncRNA［12］。同法分析mRNA 和lncRNA 之间的相关性，筛选出关键lncRNA。利用筛选出的mRNA 和lncRNA，在之前的ceRNA网络基础上构建与CD4 naive相关的ceRNA 子网络。

1.7 确定关键靶点的通路

使用关键lncRNA 表达量的中位数对所有表达量进行分组，并使用进行GO 生物过程和KEGG 通路的GSEA。通过获得的标准化富集评分（NES）及其相应的P值（NOMP-val）筛选关键的生物过程和途径。

1.8 验证

使用GEO 数据库中的另一个IDD 数据集GSE124272（包括8 例IDD 患者和8 名健康志愿者），通过IDD 和对照组基因表达的比较，分析差异表达基因，对获得的关键基因进行验证。

2 结 果

2.1 IDD 中mRNA、lncRNA 和miRNA的差异表达

通过GEO2R 分析GSE56081 和GSE63492 数据集中的数据，经过校正，得到15 220 个mRNA、370个miRNA 和623 个与IDD 相关的lncRNA。进一步根据|log2FC| ＞ 2、调整P＜ 0.05 筛选出35 个差异表达miRNA，均为低表达；根据|log2FC| ＞ 1、调整P＜ 0.05筛选出87个差异表达lncRNA，其中45个为低表达，42 个为高表达（图1）。

图1 差异表达基因Fig.1 Differentially expressed genes

2.2 共表达模块构建及GSEA

通过WGCNA 完成IDD 和对照组的共表达模块构建，通过聚类分析和动态切树算法，得到5 个不同颜色的模块，并识别出与临床特征相关的共表达模块（图2）。

图2 WGCNA 结果Fig.2 Results of WGCNA

对每个模块的KEGG 通路和GO 生物过程进行GSEA，结果表明，与免疫调节最相关的是蓝绿色（turquoise）模块，总共包含1 808 个mRNA。蓝绿色模块的KEGG 通路主要包括剪接体、溶酶体、RNA降解、亨廷顿病、人T 细胞白血病病毒1 型感染、氧化磷酸化等，GO 生物学过程主要包括参与免疫反应的中性粒细胞活化、纤溶酶原激活的调控、巨自噬调节、TORC1 信号的调节、中性粒细胞介导的免疫、中性粒细胞脱颗粒等（图3）。

图3 蓝绿色模块的KEGG通路和GO生物过程Fig.3 KEGG pathways and GO biological processes of turquoise module

2.3 ceRNA 网络构建和GSEA

使用Targetscan 和miRTarbase 数据库预测受差异表达miRNA 调节的mRNA，Targetscan 数据库中获得2 617 个、miRTarbase 数据库中获得3 940 个mRNA，取交集后获得160 个mRNA，再与免疫相关模块中的mRNA 取交集得到25 个mRNA。使用StarBase 数据库预测出差异表达miRNA 调节的797个lncRNA，将其与差异表达lncRNA 取交集后获得8 个lncRNA。使 用Cytoscape 3.7.2 构 建lncRNAmiRNA-mRNA ceRNA 网络（图4a），获得ceRNA 网络中的mRNA，并对其KEGG 通路和GO 生物过程进行GSEA。结果表明，GO 生物过程主要包括T 细胞耐受诱导调控、细胞对粒细胞-巨噬细胞集落刺激因子刺激的反应、CD4+α-β T细胞活化的正调节、细胞对转化的反应生长因子刺激的反应、α-β T细胞分化的正调节等（图4b），KEGG通路主要包括Fox O信号通路、细胞衰老、SNARE 在囊泡转运中的相互作用、MAPK 信号通路（图4c）。

图4 免疫相关ceRNA 网络及GSEAFig.4 Immune-related ceRNA network and GSEA

2.4 与CD4 幼稚T 细胞相关的ceRNA 子网络构建及GSEA

使用ImmuCellAI 进行免疫细胞浸润分析，结果显示，CD4 naive、CD8 幼稚T 细胞（CD8 naive）、衰竭T 细胞（Tex）、自然调节性T 细胞（nTreg）、诱导调节性T 细胞（iTreg）、辅助性T 细胞1（Th1）、辅助性T 细胞17（Th17）、滤泡辅助性T 细胞（Tfh）、中央记忆细胞（Tcm）、自然杀伤细胞（NKT）、黏膜相关不变T 细胞（MAIT）、树突状细胞（DC）、单核细胞（monocyte）、γ-σ T 细胞（Tgd）、CD4 T 细胞在IDD组与对照组间差异有统计学意义（P＜ 0.05，图5a），其中以CD4 naive T 细胞差异最为显著，IDD 组丰度显著低于对照组（P＜ 0.05）。

图5 免疫浸润分析及lncRNA 与mRNA 的Pearson 相关分析Fig.5 ImmuCellAI immune infiltration analysis and Pearson correlation analysis of lncRNAs and mRNAs

利 用ceRNA 网 络 中mRNA 和lncRNA 的 表 达量筛选与CD4 naive T 细胞免疫浸润相关的关键基因，经Pearson 相关分析，发现8 个lncRNA 和25 个mRNA 与CD4 naive T细胞浸润具有强相关性，其中1 个lncRNA（LINC00641）与10 个mRNA 呈强 正相关性（图5b）。

在ceRNA网络基础上构建一个与CD4 naive T细胞相关的ceRNA子网络（图6a），关键mRNA可能是UHMK1、ZFP36L2、ZCCHC3、ZBTB20。在GSE56081数据集的GSEA 中，使用LINC00641 的中位表达值进行分组，明确关键因子的调控途径。结果显示，许多机制与免疫密切相关，这些相关的KEGG 通路包括嗅觉转导、转化生长因子信号通路、溶酶体、胞质DNA 感应通路、剪接体等（图6b）；GO 生物过程主要包括负调控细胞基质黏附、白细胞迁移、细胞对细胞因子刺激的反应、炎性反应、白细胞迁移的调节等（图6c）。

图6 与CD4 naive相关的ceRNA 子网络及GSEAFig.6 ceRNA sub-network associated with CD4 naive and GSEA

2.6 验证结果

通过GEO2R 分析GSE124272 数据集的差异表达基因，结果显示，LINC00641、UHMK1、ZFP36L2、ZCCHC3和ZBTB20这5个关键基因在IDD组和对照组中均存在差异表达（表1）。

表1 GSE124272数据集中关键基因的差异表达Tab.1 Differential expression of key genes in GSE124272 dataset

3 讨 论

IDD 是临床常见病，发生机制复杂［13］。一项荟萃分析［14］表明，突出的椎间盘会自发重吸收，即IDD 患者未经手术治疗或侵入性治疗，脱出的髓核自发消失或明显缩小，但现有文献报道多为病例报告或回顾性研究，其机制可能与炎性反应和免疫调节有关［15］。因此，寻找新的免疫靶点和新的治疗方法具有重要意义。近年来，基因表达谱的研究为越来越多的疾病筛选新型生物标志物带来了新思路与新方法［16］，本研究旨在探索IDD 与免疫细胞的相关性并筛选相关的特征基因。

本研究通过基因表达谱鉴定差异表达基因，使用WGCNA 和GSEA 筛选免疫相关模块，并构建了 由8 个lncRNA、21 个miRNA 和25 个mRNA 组成的ceRNA 网络。mRNA 的GO 生物过程与KEGG通路分析表明其与免疫关系密切，如T 细胞耐受诱导的调节、CD4+的正向调节、T 细胞活化与分化的正向调节等生物过程，表明IDD 与T 细胞的调节有关，T 细胞参与IDD 的病理过程。而结果中的其他相关通路，如Fox O 信号通路、MAPK 信号通路与炎性反应及氧化应激关系密切［17-18］，鉴于炎性反应是IDD 或椎间盘再生机制的中心，其平衡是维持组织稳态的关键［19］，抑制炎性反应可能有助于治疗IDD［20］。

本研究使用免疫细胞浸润分析的方法探索与IDD 具有强相关性的免疫细胞，结果发现，IDD患者与对照组的CD4 naive 表达差异性最为显著。Naive T 细胞也称非致敏T 细胞，在暴露于抗原刺激之前处于相对静止的状态。CD4 T 细胞是人体免疫系统中重要的免疫细胞［21］。最近的研究［22］表明，CD4 naive 的功能比以前认为的要多，并且它们在表型、分化阶段、持久性、功能和解剖位置方面是多种多样的。动物实验［19］表明，椎间盘内促炎治疗可使血液和淋巴结中的T 细胞亚群增加。另一项研究［23］表明，在髓核突出的椎间盘中存在CD4 和CD8 T 细胞。在本研究中，CD4 naive 和CD4 T 细胞在免疫细胞浸润分析后均存在显著差异。Pearson 相关分析发现，UHMK1、ZFP36L2、ZCCHC3、ZBTB20和LINC00641 可能是IDD 免疫相关的关键靶点，即IDD 的特征基因。本研究还筛选了1 个与CD4 naive相关的ceRNA 子网络，并将LINC00641 识别为关键因素，可将其作为IDD的潜在治疗靶点。GSEA结果显示，转化生长因子信号通路可能是IDD 潜在的信号通路。上述关键基因的表达在GSE124272 数据集中得到了验证。

目前大多数与LINC00641相关的研究都与恶性肿瘤有关［24］，也与糖尿病足的一些生物学过程有关［25］。在IDD 的研究中，Wang 等［26］发现LINC00641 可以直接与miR-153-3p 结合，作为miR-153-3p 的天然海绵，通过miRNA介导的沉默途径靶向自噬相关基因5，并调节髓核细胞的自噬和死亡。本研究结果表明，LINC00641 可能通过免疫相关竞争性调节基因机制影响IDD 的发生与发展。

本研究发现的与CD4 naive 相关的mRNA 目前研究尚少，这些mRNA与肿瘤和免疫密切相关［27-29］。UHMK1 可能参与RNA 加工，也可能参与白细胞的形成［30］。ZFP36L2 在CD4 naive 中高表达，可抑制调节性T 细胞的功能［31］。本研究的免疫浸润分析结果显示，CD4 和CD4 naive 与IDD 密切相关，间接验证了ZFP36L2 是IDD 的重要靶点之一。同时，ZBTB20 可抑制CD8+T 细胞免疫代谢［32］，并能促进巨噬细胞炎性反应，启动先天免疫反应［33］。这些结果为IDD 的诊断和靶向治疗提供了新思路。

本研究通过GSEA 来识别LINC00641 的关键调控途径，发现与免疫相关的关键途径是转化生长因子信号通路、胞质DNA 感应通路、溶酶体通路等。胞质DNA 感应通路与先天免疫反应密切相关［34］，溶酶体是真核细胞的主要分解代谢区［35］，许多炎性疾病的根本原因可能是自噬-溶酶体途径的功能障碍［36］。最近有研究［37］表明，溶酶体与椎间盘细胞的死亡有关。近年来，有研究［38］发现转化生长因子在调节细胞生长分化和免疫功能方面具有重要作用。转化生长因子信号通路对椎间盘的发育和生长至关重要［39］，但过度激活转化生长因子信号对椎间盘是有害的，抑制异常的转化生长因子信号可以延缓椎间盘的退行性变［40］。因此，通过lncRNA 抑制转化生长因子信号通路来减少炎性反应可能是治疗IDD 的有效方法。

本研究尚有一定的局限性，受现有数据库中数据的限制，本研究选取的样本量有限，以后将增加样本量进一步研究探索；此外，本研究组未开展细胞实验或动物实验，未来将通过基础实验研究进一步验证关键lncRNA 和mRNA 在IDD 中的作用。

总之，本研究在WGCNA 的基础上，将ceRNA、IDD 和免疫调节联系起来，构建了一个与免疫调节相关的IDD ceRNA 网络。在这个网络中，lncRNA LINC00641 是一个关键靶点，它干预了mRNA（UHMK1、ZFP36L2、ZCCHC3、ZBTB20）的表达。同时，CD4 naive 表达减少，这可能与IDD 的发生机制密切相关。lncRNA LINC00641 可能通过转化生长因子信号通路调节IDD。这些结果将有助于从ceRNA 的角度更好地探索免疫调节在IDD 发展中的作用机制，并为治疗靶点的开发提供思路。