基于生物信息學探討骨關節炎miRNA-mRNA調控網絡及其作用機制

陳鋒，陳濤，章曉云，陳躍平，柴源

(廣西中醫藥大學附屬瑞康醫院骨科，廣西 南寧 530011)

骨關節炎(osteoarthritis, OA)是一種慢性關節疾病，以軟骨退化、骨贅形成及滑膜炎癥為主要特征[1]。流行病學調查顯示，全球OA總患病率高達16%，且隨著人口老齡化的增加，OA的發病率也逐年上升[2]。OA引發的慢性疼痛和殘疾嚴重影響著人們的生活質量，并造成巨大的經濟負擔，但其發病機制仍未明確。因此，尋找OA發病機制中的關鍵靶點及相關信號通路，并采用相應的措施進行干預，對于延緩甚至逆轉OA進展具有重要意義。微小RNA(microRNA, miRNA)是一種內源性非編碼RNA，可通過抑制mRNA的翻譯或加速mRNA的降解來調控靶基因的表達[3]。近年研究證實，幾乎所有細胞的生命過程都依靠miRNA的調控作用進行修飾，并對基因的表達覆蓋了轉錄水平、轉錄后水平及表觀遺傳學水平的調控[4]，其在OA的發生、預防和治療中發揮著重要作用[5]。隨著高通量技術不斷應用于生物醫學研究中，通過適當的方法和工具可有效分析OA的基因組學和蛋白質組學。因此本研究從GEO數據庫檢索分析OA的miRNA和mRNA基因表達譜數據，探討OA中差異表達miRNAs的功能及調控的靶點，為靶向診斷、治療OA提供一定的科學依據。

1 材料和方法

1.1 數據收集

以“Osteoarthritis”、“mRNA”、“miRNA”、“Human”為關鍵詞，在GEO公共數據庫中檢索OA的相關芯片,獲得OA患者和健康對照的miRNA和mRNA表達譜數據集。miRNA表達譜數據集芯片GSE105027所處平臺為GPL21575，包含12例OA患者和12例健康對照的血清；mRNA表達譜數據集芯片GSE55235所處平臺為GPL96，包含12例OA患者和12例健康對照的關節滑膜組織。

1.2 miRNAs和mRNAs數據的差異分析

利用Perl語言對數據進行基因重注釋，獲得OA患者和健康對照的基因表達譜數據集，并利用R語言limma包對基因進行校正，再進行差異分析。以P<0.05和|log2(fold change，FC)|>1作為過濾條件，篩選差異表達的miRNAs和mRNAs。

1.3 差異表達miRNAs轉錄因子預測及功能富集分析

FunRich是一個用于蛋白基因功能富集和相互作用網絡分析的軟件。為了解差異表達miRNAs的作用機制，利用FunRich軟件對差異表達 miRNAs進行轉錄因子預測及功能富集分析。

1.4 miRNA-mRNA調控網絡的構建

利用FunRich軟件預測差異表達miRNAs的靶基因，將預測結果與步驟“1.2”所得的差異表達mRNAs取交集，整合得到差異表達miRNAs與差異表達mRNAs的調控網絡，并根據兩者的相關性篩選出mRNA交叉表達負調控的miRNA，最后將數據導入Cytoscape軟件對miRNA-mRNA調控網絡進行可視化處理。

1.5 蛋白互作網絡構建

STRING是一款用來分析蛋白與基因間相互作用的數據庫。為進一步探究差異表達mRNAs的作用機制，將基因導入STRING數據庫，限定研究物種為“Homo Sapiens”，并設置連接評分>0.98，獲得蛋白互作(protein-protein interaction, PPI)關系。將所得結果導入Cytoscape軟件中，利用“NetworkAnalyzer”工具進行可視化處理，構建PPI網絡，并利用“CytoHubba”工具根據度值篩選出Hub基因。

1.6 基因富集分析

利用R語言的clusterProfiler包對Hub基因進行GO功能富集分析，研究其主要生物功能；對Hub基因進行KEGG信號通路富集分析，研究其主要信號通路，P<0.05代表富集結果顯著。最后利用R語言的ggplot2包和GOplot包進行可視化處理。

2 結果

2.1 差異表達miRNAs和mRNAs的鑒定

利用Perl和R語言對芯片進行重注釋和數據校正后，再對其進行差異分析。結果顯示，與健康對照相比，OA患者共存在265個明顯改變的miRNAs，其中上調124個，下調141個；存在838個明顯改變的mRNAs，其中上調449個，下調389個。將所得結果分別繪制成火山圖,見圖1，圖中黑點代表無差異的基因,綠點代表下調的基因,紅點代表上調的基因。

圖1 差異表達miRNAs(A)和mRNAs(B)的火山圖

2.2 轉錄因子-miRNA調控網絡及富集分析結果

通過FunRich軟件對265個差異表達的miRNAs進行轉錄因子預測，共得到203個轉錄因子，根據其所調控的基因數與P值，選取富集程度最顯著的10個轉錄因子進行展示，見圖2。差異最顯著的10個轉錄因子分別為EGR1、SP1、POU2F1、SP4、MEF2A、NFIC、ZFP161、NKX6-1、FOXA1、RREB1，且其P值均<0.001，基本信息見表1；對265個差異表達的miRNAs進行功能富集分析，其生物過程主要涉及核堿基、核苷、核苷酸和核酸代謝的調控，調節細胞生長，轉錄等；細胞成分主要涉及細胞核、細胞質、內涵體等；分子功能主要涉及轉錄因子活性、轉錄調節活性、泛素特異性蛋白酶活性等，見圖3。

圖2 差異表達miRNAs的轉錄因子預測

圖3 差異表達miRNAs的GO功能富集分析

2.3 miRNA-mRNA網絡

利用FunRich軟件對265個差異表達miRNAs的下游靶標進行預測，共得到3 572個靶基因，與838個差異表達mRNAs取交集獲得194個候選靶標。然后，再根據差異表達miRNAs和差異表達mRNAs間的log2(FC)關系，篩選出交叉表達負調控的mRNAs 96個、miRNAs 25個。最后將數據導入Cytoscape軟件對miRNA-mRNA調控網絡進行可視化處理，見圖4。圖中紅色節點表示基因上調，綠色節點表示基因下調，節點間的連線表示兩者呈負調控關系。

圖4 miRNA-mRNA調控網絡

2.4 PPI網絡

借助STRING數據庫和Cytoscape軟件構建PPI網絡，見圖5。圖中共涉及46個節點、62條邊，節點代表mRNA，兩節點的連線則說明兩者間存在相互作用關系；節點越大、顏色越深則度值越大，邊的粗細反映連接評分，邊越粗，mRNA間的互作關系越緊密。根據度值篩選出排名前6的mRNAs為JUN、FOS、IL-6、TNF、IL-1β、CXCL8。這些度值較大的mRNAs在整個網絡中起著關鍵作用，可能是OA發生發展的Hub基因，基本信息見表2。

圖5 mRNA的蛋白互作網絡

表2 關鍵mRNA的基本信息

2.5 GO與KEGG富集分析

GO富集分析Hub基因的功能過程中，共確定了934個條目，其中包括908個生物過程、5個細胞成分和21個分子功能。根據P值，前10個富集的生物過程、細胞成分和分子功能如圖6所示。生物過程主要涉及對脂多糖的反應、對細菌來源分子的反應、DNA結合轉錄因子活性的調節等方面；細胞成分主要涉及RNA聚合酶Ⅱ轉錄調節復合物、轉錄調節復合物、吞噬杯等；分子功能主要涉及細胞因子活性、細胞因子受體結合、受體配體活性等。KEGG富集分析Hub基因共確定了70個條目，主要涉及IL-17信號通路、Toll樣受體信號通路、AGE-RAGE信號通路、TNF信號通路、NOD樣受體信號通路，如圖7所示。

圖6 Hub基因的GO富集分析氣泡圖

圖7 Hub基因的KEGG富集分析氣泡圖

3 討論

在OA的研究過程中，針對骨代謝、基質降解、氧化應激和免疫炎癥反應等方面的治療均已應用于臨床，但這些方法的療效并不理想[6]。越來越多的研究表明，OA的發生發展與miRNA密切相關，在最近的一項研究中發現[7]，miR-206可調控補體因子H(complement factor H,CFH)，干擾素誘導的含有解旋酶C域的蛋白1(interferon induced with helicase C domain 1,IFIH1)，神經損傷誘導蛋白1(nerve injury-induced protein 1,NINJ1)等基因的轉錄水平，并顯著抑制軟骨細胞的增殖，同時促進分解代謝酶和凋亡誘導因子的表達。這表明通過敲低miR-206表達可有效抑制OA的軟骨降解，并在介導OA發病機制中起著重要調控作用。通過基因測序技術更是發現在關節軟骨損傷過程中表達量發生明顯改變的miRNAs高達數百個。然而，針對OA尚缺乏完整的miRNA-mRNA調控網絡模式。因此，深入研究miRNA在OA發生發展中的作用及其調控機制有重要意義。

之前的研究已經證實miRNA的表達受轉錄因子的調控，且這些轉錄因子在OA發生發展過程中發揮重要的調控作用，主要涉及滑膜炎癥反應，軟骨與骨組織的破壞[8]。本研究預測了調控OA的差異表達miRNAs的上游轉錄因子共203個，根據P值選出排名前10的轉錄因子分別為EGR1、SP1、POU2F1、SP4、MEF2A、NFIC、ZFP161、NKX6-1、FOXA1、RREB1，其中EGR1差異最顯著。EGR1是一種核蛋白，具有轉錄調節功能，可潛在地調節miRNA的表達。研究表明，EGR1在miR-195上游啟動子區有結合位點，可抑制miR-195的表達。而EGR1又可通過與Krüppel樣因子5(Krüppel-like factor 5, KLF5)結合，減少其泛素化降解，從而上調基質金屬蛋白酶的轉錄水平，最終導致軟骨退變并誘發OA[9-10]。因此研究轉錄因子與miRNA組成的調控網絡可對疾病的發病機制在系統層次上提供重要的線索。但這些轉錄因子在OA中的具體機制仍有待進一步驗證。

對自身靶基因的調控是miRNAs在疾病發病過程中發揮作用的研究重點，本研究結果顯示265個差異表達miRNAs共調控了194個差異表達mRNAs。對265個差異表達miRNAs進行功能富集分析顯示其主要涉及核堿基、核苷、核苷酸和核酸代謝的調控，調節細胞生長，轉錄等過程，與既往研究完全符合[11]。而在miRNA-mRNA調控網絡中，同屬于原癌基因家族的JUN和FOS蛋白可在胞核中形成復合物，并與DNA內的轉錄因子激活蛋白1(activator protein 1, AP-1)調節位點結合，形成AP-1復合物來調節細胞生物過程。研究表明，TNF-α、IL-1、干擾素等炎癥因子通過MAPK信號介導AP-1表達，誘導基質金屬蛋白酶和膠原酶活化進而破壞軟骨基質[12-13]。Luo等[14]研究表明，miR-139-5p過表達會抑制c-Jun的水平進而減少內皮因子表達，而內皮因子過表達誘導的血管異常增生和炎性反應是介導OA軟骨退變、滑膜增厚及滑膜炎癥的關鍵途徑；miR-139-5p過表達可通過抑制c-Fos減少TNF-α、NF-κB的表達，而TNF-α在OA中是介導炎癥反應和軟骨降解的關鍵靶基因[15-16]；miR-222-3p可通過抑制c-Jun表達減少心肌纖維化，而在OA中纖維化是軟骨退變的病理變化之一[17-18]。因此，miRNA在OA的發病過程中與mRNA的互作對疾病進展影響重大，可為闡明OA的發病機制及后續研究提供依據。

為進一步探究差異表達mRNAs在OA中的作用機制，通過PPI網絡還發現IL-6、TNF、IL-1β、CXCL8等Hub基因。這些基因的生物過程主要涉及對脂多糖的反應、對細菌來源分子的反應、DNA結合轉錄因子活性的調節等方面，與OA病情的發生發展密切相關[19]。研究表明，IL-6、TNF-α、IL-1β等炎癥因子在OA中的表達顯著上調，且三者間的協同作用共同導致了巨噬細胞去極化，滑膜炎癥，軟骨細胞凋亡及基質降解等病理改變[20-21]。而巨噬細胞受上述炎癥因子刺激會產生趨化因子CXCL8，進而募集和激活相關炎癥因子，加劇OA的病變程度[22]。對Hub基因進行KEGG富集分析，發現其主要富集在IL-17、Toll樣受體、AGE-RAGE、TNF和NOD樣受體信號通路上。研究發現，OA患者外周血中IL-17水平明顯高于健康者,且與病情的嚴重程度呈正相關,通過降低血液中IL-17的水平能有效地控制OA患者的病情[23]；Toll樣受體和NOD樣受體在體內水平與IL-1β、IL-6、TNF-α等炎癥因子呈正相關，并會誘導破骨細胞分化并造成骨破壞，兩者參與的信號通路均與OA有著緊密的聯系[24-25]；AGE-RAGE信號通路活化后會導致炎癥因子的表達增多，進而引起細胞凋亡和組織損傷。通過抑制該通路可降低氧化應激與炎癥反應，有效防止細胞損傷[26]；TNF信號通路在炎癥反應早期會釋放大量TNF-α，進而誘導IL-1、IL-6等炎癥因子的產生，與滑膜炎癥及軟骨細胞凋亡等病理改變過程密切相關[27]。因此，我們認為通過作用上述靶點與通路可有效調控血液中的炎癥因子抑制機體炎癥反應，同時也可調控軟骨細胞和破骨細胞的增殖、分化與凋亡，從而緩解關節骨與軟骨破壞情況，達到治療OA的目的。

綜上所述，本研究借助生物信息學成功篩選并構建了OA的miRNA-mRNA調控網絡，探討了可能與OA發生發展相關的致病機制。雖然該調控網絡某些靶點的相關機制已闡明，但大多數miRNAs和mRNAs的致病機制仍不清楚，這對于挖掘新的機制和生物標志物很有意義，并可為今后闡明基因間的調控機制、協助預防和治療OA提供關鍵的數據支撐和參考依據。此外本研究存在一些缺點，首先，由于篩選條件的限制，只能對關鍵基因、關鍵信號通路進行分析，在一定程度上使得研究結果具有局限性；其次，缺乏進一步的實驗來驗證miRNA-mRNA調控網絡中基因間的互作關系。在后續的研究中，應通過相關實驗進一步討論和驗證本研究的預測。