轉錄組分析鑒定腰椎椎間盤退行性變中潛在免疫相關生物標志物

李志超，王 磊，薛景才，王文波，陳 煜，齊軍強，許 鵬

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 發展中的作用機制，并為治療靶點的開發提供思路。