基于基因芯片的TNF-α誘導骨關節炎細胞模型生物信息學分析

袁發滸，胡 松，黃麗霞，黃念芳，宋文劍，孫賓蓮

(江漢大學，武漢 430056)

骨關節炎(osteoarthritis，OA)是一種常發于老年人的慢性退行性疾病，日益快速地成為世界人群主要致殘因素[1-2]。據最新數據分析顯示，目前在全球范圍內，六十歲以上的人群中有18%的女性及9.6%的男性患有癥狀性的OA，其中四分之一的人更是無法進行正常的日常活動[3]。到2050年，預計將有1.3億人患有OA，這將構成巨大的全球性社會負擔，隨著我國老齡化進程的加快，OA患者更是日益增加[4]。

OA的病理學因素是多方面的，具體發病機制尚不清晰，主要包含軟骨下骨的重塑、滑膜炎癥和關節軟骨的喪失[5]。軟骨細胞是軟骨組織的唯一細胞類型，對維持和重塑軟骨結構、保證軟骨細胞外基質功能等方面起著至關重要的作用。近年來，由于軟骨細胞凋亡在軟骨退化中的作用，軟骨細胞凋亡相關信號通路在OA方面的研究日益受到關注，抑制軟骨細胞死亡目前已成為治療軟骨退化的一個潛在治療靶點。此外炎性細胞因子也被認為是OA發生發展中的重要因素，特別是白介素1β(IL-1β)和腫瘤壞死因子α(TNF-α)，他們驅動分解代謝途徑，并使OA進展持續下去[6-8]。

高通量基因表達分析目前已成為研究復雜疾病機制的重要方法，此外，高通量基因表達譜分析也有助于發現疾病新的生物標志物，監測疾病進程及評價藥物療效。目前，已有應用基因芯片技術分析OA患者相關基因表達的研究，這些研究成果有助于我們理解OA進展過程中所涉及的病理機制。然而，由于病人的個體差異，如患者的年齡、關節的分層等，這些因素往往導致研究結果產生一些偏倚。本研究通過培養原代小鼠軟骨細胞，以TNF-α刺激建立OA細胞模型，應用小鼠全基因組表達譜芯片分析，從而識別差異基因表達譜和信號通路，為全面深入探尋OA致病機制，發現新的OA診斷標志物提供更多幫助。

1 材料和方法

1.1 實驗動物

SPF級KM小鼠，體重 20 ～22 g，6 周齡，性別不限，來源于華中科技大學實驗動物中心[SCXK(鄂)2016-0009]。實驗于華中科技大學實驗動物中心進行[SYXK(鄂)2016-0057]，本研究涉及的實驗動物的使用及管理經華中科技大學實驗動物倫理審查委員會核準，審批號：201803170028，研究按實驗動物使用的3R原則給予實驗動物人道的關懷。

1.2 主要試劑與儀器

重組TNF-α(西格瑪奧德里奇(上海)貿易有限公司，貨號：T7539)；DMEM培養基(賽默飛世爾科技(中國)有限公司，貨號：11965084)；胎牛血清(賽默飛世爾科技(中國)有限公司，貨號：12483020)；胰蛋白酶(Sigma)；Ⅱ型膠原酶(Gibco)；TRIzol試劑(Invitrogen)；PrimeScriptTMRT Reagent試劑盒(TaKaRa)；SYBR Green qPCR Master Mix試劑盒(TaKaRa)；小鼠基因表達微陣列芯片(Agilent，型號：G4846 A)。二氧化碳培養箱(力康生物醫療科技控股集團，型號：HF100)；體視顯微鏡(德國萊卡，型號：M205)；生物超凈臺(北京東聯哈爾儀器制造有限公司，型號：BCN-1360)；Stepone Plus熒光定量PCR儀(ABI)。

1.3 實驗方法

1.3.1 小鼠原代軟骨細胞分離、培養及處理

小鼠原代軟骨細胞于KM小鼠膝關節分離培養獲得。選取6只同窩生SPF級KM小鼠，參照金旻等[9]報道的方法分離培養小鼠軟骨細胞，無菌環境下分別剝離新生小鼠雙側膝關節軟骨，剪碎后以2.5 g/L胰蛋白酶與0.1% Ⅱ型膠原酶聯合消化后離心收集。原代軟骨細胞以1×107cell/cm2密度接種于DMEM培養基(含10%胎牛血清)置于二氧化碳培養箱中培養，培養條件為37℃、5% CO2。原代軟骨細胞穩定培養3 d后，將其分設為兩組，一組為正常對照組(C組)，另一組為TNF-α誘導的OA模型組(T組)，每組三瓶細胞(分別來源于一只小鼠)，OA模型組以50 ng/mL的TNF-α孵育24 h，對照組以等積極PBS做相同處理。

1.3.2 RNA提取

兩組細胞經24 h處理后，離心除去培養基，經PBS漂洗沉淀后立即加入總RNA抽提試劑TRIzol，依照產品說明書操作提取總RNA，提取所得的RNA用超微量分光光度計測定其濃度后，分裝于-80℃冰箱備用。

1.3.3 微陣列芯片雜交與數據分析

微陣列芯片雜交實驗委托于康成生物科技有限公司完成，實驗流程參考Li等[10]所述。所用芯片為Cy3標記的小鼠全基因組表達譜芯片(Agilent Mouse 4x44 K Gene Expression Microarrays v2)。原始熒光數據利用Agilent公司的Gene Spring 11.0 軟件進行均一化處理。以參照組為對照，以表達變化倍數大于2，P值小于0.001為條件篩選OA組差異表達基因(differentially expressed genes，DEGs)。獲得的DEGs提交至DAVID平臺(DAVID Functional Annotation Bioinformatics Microarray Analysis，https://david.ncifcrf.gov/home.jsp)進行基因功能聚類及通路分析。

1.3.4 熒光定量PCR驗證

為評估芯片表達譜數據可靠性，從差異表達基因列表中挑選20個基因(上調、下調基因各十個)進行熒光定量PCR實驗驗證。軟骨細胞總RNA用PrimeScriptTMRT Reagent試劑盒進行提取，采用隨機六核苷酸引物合成第一鏈cDNA，采用SYBR Green qPCR Master Mix試劑盒配置反應體系，于熒光定量PCR儀擴增檢測基因表達情況。將Gapdh設置為內參基因，基因相對表達量采用2-ΔΔCt方法進行分析。

1.4 統計學方法

2 結果

2.1 樣本表達量聚類分析

對樣本之間進行系統聚類，利用最廣泛、最簡便的層次聚類法(Hierarchical Clustering)，檢查不同樣品之間同一基因表達的相關性，如果大部分基因的表達量相似，相關系數高，說明數據的均一化程度高。在此，兩組樣本間(各三個樣品)的差異表達基因進行了層次聚類，可見C組的三個樣品和T組的三個樣品分別聚類在一起，而兩個組別也可區分開來(圖1)，說明兩組樣本的基因表達有一定的差異性。

注：C組：正常對照組；T組：TNF-α誘導的OA模型組。下同。圖1 層次聚類分析熱圖Note. Group C, Normal group. Group T, TNF-α-induced OA model group. The same below.Figure 1 Heat map and supervised hierarchical clustering

2.2 差異表達基因分析

火山圖(Volcano Plots)是可視化表達兩種不同模式條件的實用工具，它能統一體現組間變化倍數值和P值，可直觀了解變化幅度和統計顯著性之間的關系(考慮了變化幅度和變異性)。同時，基于這些值的基因子集也被分離表現出來。將C、T兩組的基因表達量做火山圖分析，如圖2所示，縱向的左右兩根綠線示意為表達差異倍數2.0，橫向的綠線示意為P<0.05。圖中的紅點表示DEG(FC>2.0，P<0.05)，其中有8096個表達上調DEG和6413個表達下調DEG。

圖2 T組與C組差異表達基因火山圖Figure 2 Volcano plots of differentially expressed genes between the groups C and T

圖3 T組與C組前35條差異表達基因分類注釋條目Figure 3 The top 35 GO terms of differentially expressed genes between the groups C and T

為了獲知這些DEGs的功能及所在的調控通路，分別將這些上調和下調的DEGs提交至DAVID在線平臺進行功能聚類分析， DAVID的功能聚類數據庫整合了Gene ontology， Interpro，KEGG等基因功能數據庫，獲得的基因注釋與功能通路細致、可信度高。相較于正常對照組，TNF-α誘導處理的OA模型組的DEGs注釋在細胞核、轉錄調節、細胞分化、ATP結合、細胞周期、CXC趨化因子等方面，圖3詳細列出了主要的前35個GO條目，圖4列出了包含神經營養因子信號通路、賴氨酸降解、T細胞受體信號通路、ErbB信號通路、mTOR信號通路等前25條主要通路。其中下調表達差異最大前10個基因分別是P2ry1、Fjx1、Wnk1、Tln1、Lasp1、Pdia4、Prm2、Olfr214、Tnfrsf14、Spata2(表1)，而上調表達差異最大前10個基因分別是Dnahc14、Slc4a7、Cxcl3、Rasgrp1、Spata31d1c、Zfp866、Il1b、Nf1、Slc7a11、Cxcl2(表2)。

對骨關節炎進展相關的基因表達情況進行了整理(表3)，發現免疫調節和促炎類細胞因子相關基因(Il6、Il9、Il17a)顯著上調。細胞外基質降解和重塑相關基因(Mmp3、Mmp9)顯著下調。線粒體功能相關基因(Cyc1、Cox3、Uqcr10)表達顯著上調。

為評估芯片表達譜數據可靠性，從差異表達基因列表中挑選20個基因(上調、下調基因各十個)進行熒光定量PCR實驗驗證。如圖5所示，可見驗證的絕大多數的基因表達情況與芯片數據吻合。

表1T組與C組前10個表達下調的差異表達基因

Table1The top 10 significantly down-regulated genes between the groups C and T

表2 T組與C組前10個表達上調的差異表達基因

表3 TNF-α誘導的OA相關基因表達情況

注：↑表示T組與C組相比，基因表達上調；↓表示T組與C組相比，基因表達下調。

Note. Compared with the group C, ↑ and ↓ indicate up-regulation and down-regulation of gene expression in the group T, respectively.

注：以Gapdh表達量做內參校正。Mmp9的基因芯片檢測表達量與qPCR實驗測試表達量相比， *P<0.05(t檢驗)。圖5 基因芯片表達數據的熒光定量PCR驗證Note. Normalized using Gapdh as an internal control. Microarray data of Mmp9 compared with qPCR data, *P<0.05 (t-test).Figure 5 Verification of microarray data through qPCR

3 討論

以往對于OA病理機制的研究大多以某幾個基因或某條信號通路進行，例如申成凱等[11]探討了白介素-1β基因與晚期OA的關聯，張榮凱[8]等則研究了脂蛋白基因在早期OA的表達情況，馬春輝等[12]研究了凋亡基因通路在OA的作用，楊占東[13]等總結了SREBPs在調節骨關節炎關節軟骨退變中的作用機制，這些研究對OA中某些顯著的功能基因進行了深入細致的探討，但總體而言，這些研究存在一定的片面性，極少有從整體上去挖掘OA相關的病理機制的研究， Qing等學者[14]曾以瘦素誘導的大鼠OA模型進行了表達譜芯片的分析，篩選出了如Bcl2l11等未被研究過的可能在OA起重要作用的基因，極大地豐富了OA相關領域的信息。本研究中，首次全面比較了TNF-α誘導的OA原代軟骨細胞模型全基因組表達譜，用熒光定量PCR方法對一部分基因的芯片結果進行了驗證，一些新的候選基因首次被發現與OA的發病機制有關。

最新的一些研究顯示白介素IL-1和基質金屬蛋白酶(MMP)的表達異常在OA的發病機理中有重要作用[15-17]。通過對GO和KEGG通路的分析，證實了這些基因在OA發病機制中的重要作用。IL-1是OA進展過程中軟骨退變最重要的因素之一，也是基質穩態病理生理失衡的主要誘發因素。IL-1還可抑制軟骨細胞外基質合成代謝，導致軟骨組織[15]的破壞。本研究中，TNF-α處理后的軟骨細胞IL-1α、IL-1β(表2)和IL-6、IL-9和 IL-10顯著升高(表3)，提示TNF-α誘導OA的分子機制可能涉及炎癥因子的調節。其次，基質金屬蛋白酶系編碼基因的表達(如MMP-3、9、12)同樣由TNF-α顯著調節(表3)，這進一步提示了TNF-α參與OA的發病機理。

氧化應激是自由基在機體產生的一種負面作用，是導致衰老和多種疾病的一個重要因素。越來越多的證據表明，線粒體功能障礙和線粒體DNA損傷在OA的發生發展中起重要作用[18]。Grishko等[19]研究發現OA患者軟骨細胞中發現了線粒體DNA的損傷，活性氧自由基暴露后，線粒體DNA修復受阻，軟骨細胞凋亡增加。本研究發現，經TNF-α誘導后，軟骨細胞中涉及線粒體呼吸鏈相關的基因表達顯著上調(表3)，印證了線粒體功能障礙在OA中的重要作用。

本研究發現小鼠原代軟骨細胞經TNF-α處理后，大量細胞色素C氧化酶相關基因的表達顯著上調(表3)，Wang等[20]發現可通過過氧化物酶體增殖激活受體γ(PPARγ)共激活劑1α可逆轉OA中線粒體功能的衰減，此外，嘌呤受體基因P2ry1同樣異常表達(表1)，這些都提示線粒體功能損傷在OA中可能有重要作用。此外，本研究發現一些神經功能相關的基因在本模型中的異常表達，如含嗅質蛋白(olfactomedin)結構域的基因Olfml1、Olfml3及嗅覺受體基因家族(olfactory receptor)，神經纖維瘤基因(neurofibromatosis 1，Nf1)都顯著上調，這些全新的發現對OA中神經相關的病理機制提供了全新的線索。

綜上所述，本研究利用小鼠基因表達芯片分析比較了TNF-α誘導的OA細胞模型基因表達的情況，并鑒定了一些DEGs和表達通路。此外，本研究還挖掘了一些新穎的OA相關的理論分子機制，后續這些病理假設還需要進一步探討驗證。