類風濕關節炎滑膜組織線粒體功能相關差異表達基因的分析

巴燕娜，袁 晴，賀 彤，穆 楠，吳振彪，袁志棟,3※

(1.第四軍醫大學西京醫院臨床免疫科，西安 710032； 2.第四軍醫大學基礎醫學院生理與病理生理學教研室，西安 710032；3.贛南醫學院基礎醫學院，江西 贛州 341000)

類風濕關節炎(rheumatoid arthritis，RA)是一種以關節滑膜炎性病變為主的系統性自身免疫性疾病，主要表現為慢性、對稱性和進行性的多關節炎[1]。該病影響了近1%的人群，女性發生率是男性的2～3倍[2]。滑膜組織中成纖維細胞樣滑膜細胞在RA中顯示出異常的行為，被認為與關節損害有關[3-4]，這種細胞對凋亡的抵抗性是RA的一個主要特征[5]。因此，系統分析RA滑膜組織細胞的基因表達水平，不但可加深對RA病理性過程及臨床癥狀的理解，而且能更好地理解藥物治療的作用機制[6]。

線粒體是真核細胞中的一種重要細胞器，擁有獨立的基因組DNA，共37個基因。這些基因和核基因組編碼的定位于線粒體內蛋白的基因對線粒體的氧化磷酸化、能量轉換等功能活動是必需的。RA患者高度失調的微血管結構造成了缺氧微環境，異常的滑膜細胞代謝和線粒體功能障礙隨之而來，增加活性氧類，誘導炎癥[7]。另外，在細胞應激或損傷時，線粒體成分可被釋放到胞質或細胞外，可作為危險信號發揮作用[8]。因而，研究線粒體相關基因的異常表達對探討RA的發病機制非常關鍵。

轉錄組芯片分析已被用來研究RA病理發生機制、治療效果和生物標志物[9-13]。雖然線粒體與炎癥、免疫密切相關，但線粒體在RA病理發生機制中的功能一直未受到研究者的重視。為探索這些問題，本研究通過系統分析RA與正常人滑膜組織的基因表達芯片數據，以期確定RA病理發生與線粒體之間的關系及鑒定關鍵的影響滑膜組織病變的線粒體相關基因。

1 材料與方法

1.1滑膜活檢組織的公共芯片數據 RA滑膜活檢組織的基因表達數據下載于美國國家生物技術信息中心的GEO數據庫的芯片數據，ID為GSE77298。該數據集包含16例末期RA患者和7例無關節疾病的受試者的滑膜活檢組織，提取總RNA后在Affymetrix芯片上檢測，然后獲得滑膜組織的基因表達數據[10]。但由于該芯片設計上未考慮線粒體基因組編碼的基因，因此只能檢測線粒體相關基因而不能檢測線粒體編碼基因的表達變化。

1.2分析方法

1.2.1基因表達芯片數據的分析 下載芯片數據后，用R軟件安裝simpleaffy、affyPLM和limma等程序包后對芯片基因表達數據依次處理，利用主成分分析、熱圖聚類分析評估樣本間的總體相似度，剔除錯誤分組的樣本，鑒定差異表達基因，繪制差異表達基因的火山圖和熱圖。

1.2.2線粒體相關基因功能注釋與蛋白網絡分析 本研究定義核基因組編碼但其蛋白產物定位于線粒體，進而與線粒體結構、功能相關的基因為線粒體功能相關基因。線粒體功能相關基因信息來源于MitoMiner(MitoMiner 4.0v2018JUN)數據庫[14]，該數據庫收錄了1 626個基因，包括人的1 184個已知其蛋白產物定位于線粒體的基因和442個預測其蛋白產物定位于線粒體的基因。將其與鑒定的全轉錄組差異表達基因列表取交集，即提取顯著性差異表達的線粒體功能相關基因。差異表達的線粒體相關基因的功能注釋與蛋白質網絡分析采用在線分析工具Metascape[15]，而Metascape應用了圖論聚類算法——分子復合物檢測MCODE(Molecular Complex Detection)[16]。

2 結 果

2.1滑膜組織差異表達的基因 首先對芯片數據進行標準化，并根據主成分分析和熱圖聚類分析結果評估樣本間的總體相似度，發現RA1、RA2、RA9和RA15樣本數據與總體的差異較大，進而剔除這4例樣本，并分析基因差異表達，獲得1 986個在RA滑膜組織中顯著差異表達的基因(校正后的P<0.05且表達值變化2倍以上)，其中表達水平2倍以上變化的上調基因有1 490個，下調的基因有496個，見圖1、圖2。這些差異表達基因中的IGLC1、SPP1、MMP1、JCHAIN和CHI3L2為表達水平上調倍數變化6以上的基因；PLIN1和ADIPOQ(脂聯素)是表達水平下調倍數變化6以上的基因。

Up：上調(紅色)；Down：下調(綠色)；logFC：RA患者與正常人標準化后的基因表達值比值，即以2為底倍數變化的對數；-log10(adj.P.Val)：以10為底校正后P值的對數的相反數；NotSig：無顯著變化

圖1 基因表達的火山圖

HC1-7：正常對照；RA：類風濕關節炎患者；藍色：基因表達下調；紅色：基因表達上調

圖2 差異表達基因的層級聚類熱圖

2.2差異表達的線粒體功能相關基因 本研究結果顯示MitoMiner 4.0數據庫收錄的1 355個線粒體功能相關基因在芯片表達數據中有表達，但由于芯片設計的局限性，線粒體基因組編碼的基因未檢測到表達。在前面鑒定的1 986個顯著差異表達基因中發現有169個線粒體功能相關基因，占目前所知的蛋白產物定位于線粒體的基因總數量的10.39%，其中表達上調的基因數為131個，下調的基因數為38個。

2.3差異表達線粒體相關基因的功能注釋 為了解169個差異表達線粒體相關基因的功能，對這些基因進行注釋后獲得了以條形圖形式列出了P值排列前20的富集通路，離散色階代表不同的統計學差異。P值最顯著的通路是核苷酸代謝過程通路(GO:0009117)，其次是蛋白質在內質網中的蛋白質加工(KEGG通路：HSA04141)和線粒體結構組織通路(GO:0007005)。另外，還包括凋亡信號的調控(GO:2001233)、內質網應激反應(GO:0034976)、氧化應激反應(GO:0006979)和活性氧類代謝過程(GO:0072593)通路等(圖3A)。差異表達線粒體相關基因的功能通路網絡中顏色代表了聚類簇的ID，聚類簇ID與圖3A富集的通路一致，每個節點代表了一個富集的術語，但功能通路網絡聚類簇更體現了富集的術語間的功能相關性(圖3B)。這些異常表達線粒體功能相關基因的功能富集通路體現了RA相關功能通路出現了擾動。

2.4差異表達的線粒體相關基因的蛋白相互作用網絡 利用Metascape生成了差異表達的線粒體相關基因的蛋白質互作網絡(圖4A)。Metascape檢測的蛋白-蛋白相互作用網絡的密集連接區，它們代表了分子復合物。每一個MCODE成分分別進行了通路和過程富集分析，這些MCODE成分中根據P值打分最高的術語是KEGG通路HSA04141[內質網蛋白質加工，log10(P)=-17.1]、GO:0009117[核苷酸代謝過程，log10(P)=-14.1]和GO:0006753[核苷磷酸鹽代謝過程，log10(P)=-14.0]。蛋白質互作網絡中4種不同的顏色表示識別到的4個MCODE成分(圖4A、圖4B)，分別是MCODE_1：Reactome 通路R-HSA-1799339[靶向膜的SRP依賴的共翻譯蛋白，log10(P)=-11.1]、Reactome通路R-HSA-72766[翻譯，log10(P)=-8.6]和HSA04141[內質網蛋白質加工，log10(P)=-4.2]；MCODE_2：Reactome通路R-HSA-446203 [天冬酰胺N-連接糖基化，log10(P)=-9.5]、KEGG通路HSA04141[內質網蛋白質加工，log10(P)=-8.8]和KEGG通路HSA00510[N-聚糖生物合成，log10(P)=-6.2]；MCODE_3：Reactome通路R-HSA-382551[小分子轉運，log10(P)=-3.1]和MCODE_4：GO:0050878[體液水平的調控，log10(P)=-3.8]。

4A：線粒體相關基因的蛋白質互作網絡圖； 4B：蛋白質網絡中密集連接的蛋白質組

圖4 差異表達的線粒體相關基因編碼的蛋白質互作用網絡

3 討 論

線粒體在細胞內不僅參與物質的代謝、能量轉化和鈣離子存儲等生命活動，還在病理狀態下由損壞的線粒體導致先天性免疫系統異常激活而引起自身性免疫疾病[17]。然而，在RA中的線粒體功能一直未引起研究者的重視。為了探究線粒體功能相關的基因在RA中的表達變化和功能，本研究分析了GEO數據庫基因表達芯片數據GSE77298，獲得了1 986 個在RA滑膜組織中顯著差異表達的基因，這些基因與DisGeNET數據庫[18]收錄的1 832個RA相關基因的交集為322個，與Liu等[19]報道的RA患者中顯著差異表達上調的228個(實際是227個)基因和下調的36個基因相比，相一致的上調基因110個，相一致的下調基因10個；另一項研究報道在滑膜組織中RA組比正常組高表達的基因1 018個，下調的基因893個[20]，與本研究結果相同的分別為686個和257個。另外，在RA滑膜組織中顯著差異表達、倍數變化6以上的上調基因SPP1、MMP1在以前的研究中[21-23]也得到了驗證。而ADIPOQ是倍數變化6以上的表達下調基因，Harshan等[20]的研究結果顯示脂聯素在RA患者滑膜組織中表達水平較正常人顯著降低，本研究結果與上述文獻一致。在這1 986個顯著差異表達基因中，本研究鑒定了169個與正常對照存在統計學差異表達的線粒體功能相關基因。而這些線粒體功能相關基因部分也得到了DisGeNET數據庫[18]及RA相關研究報道[19-20,24-28]的數據支持與驗證。169個線粒體相關基因的25個與DisGeNET數據庫中的一致，其中FHL1、G0S2和SLIT3為表達下調的3個基因，其他的22個為上調基因。SLIT3能夠抑制RA中Robo3誘導的滑膜成纖維細胞的入侵[24]，因此SLIT3的表達下調可能促進RA的病情惡化。另外，G0S2可作為生物標志物用來預測響應抗腫瘤壞死因子治療RA患者[27]。而表達上調倍數變化最大的為TYMS，其遺傳多態性與RA患者的甲氨蝶呤臨床治療的效果相關[25-26]。

本研究發現的部分差異表達基因不僅與已往研究報道[18-28]一致，而且通過功能通路注釋和互作網絡分析，結果顯示GO:0009117(核苷酸代謝過程)和KEGG通路HSA04141(內質網蛋白質加工)不僅分別是差異表達的線粒體相關基因功能富集通路中P值最顯著的通路之一，也分別是蛋白互作網絡MCODE成分中根據P值得分最高的生物學過程術語之一。將蛋白質互作數據與通路和生物學過程富集數據相結合顯示了差異表達的線粒體功能相關基因的蛋白質互作網絡可能在RA中存在異常。而有報道發現青蒿素能通過抑制PDK1(GO:0009117共有34個差異表達的線粒體功能相關基因之一)誘導的蛋白激酶B和核糖體S6激酶家族成員2磷酸化激活來抑制RA成纖維樣滑膜細胞的遷移和入侵[28]。

另外，富集的通路還包括線粒體結構組織通路(GO:0007005)、凋亡信號的調控(GO:2001233)、內質網應激反應(GO:0034976)、氧化應激反應(GO:0006979)和活性氧類代謝過程(GO:0072593)通路，以及蛋白質互作網絡中識別到的MCODE成分R-HSA-1799339 (靶向膜的SRP依賴的共翻譯蛋白)等。以上功能與通路的富集注釋結果以及已有相關研究報道[18-28]的驗證顯示了這些線粒體功能相關基因的異常表達可能與RA滑膜組織的病理變化有關。

需要指出的是，本研究分析的是芯片數據，由于該芯片探針未覆蓋線粒體基因組編碼的基因，所以分析結果僅鑒定了核基因組編碼的線粒體功能相關基因，并且這些基因均為編碼蛋白基因。因此，后續更完整的研究應囊括1 626個定位線粒體內的相關編碼基因、線粒體基因組編碼的基因，以及非編碼RNA基因，以期能夠全面系統地揭示線粒體對RA的發生、發展及滑膜組織病變所起的作用，找到有效的RA疾病診斷分子標志物和治療的藥物靶標。