db/db小鼠糖尿病視網膜病變相關基因表達及通路的生物信息學分析

李 倩，何 飛，張 然

0引言

糖尿病視網膜病變(diabetic retinopathy， DR)是糖尿病的一種嚴重的微血管并發癥，也是導致失明的主要原因，通常主要影響年齡在30～70歲的人群[1]。糖尿病影響視網膜的整個神經血管單位，逐漸病變并最終纖維化，且發病率越來越高[2]。然而，DR目前的治療選擇仍然有限且有副作用，因此糖尿病患者最終失明的風險仍然存在[3]。雖然DR是糖尿病的常見并發癥，但我們對其潛在的分子機制知之甚少。近年來，C57BL/KsJ-db/db小鼠作為2型糖尿病的自發性糖尿病模型用于研究DR的發病機制[4]。本研究選取基因表達綜合數據庫(GEO)中的db/db小鼠芯片GSE55389，采用生物信息學方法進行全面分析，以期望從分子水平上得到新的生物標記物，靶點蛋白和通路，從而為DR的早期診斷和治療靶點研究提供新的突破。

1材料和方法

1.1材料基因芯片GSE55389數據集下載自GEO網站(http://www.ncbi.nlm.nih.gov/geo/)，該芯片數據集共有8個從視網膜提取的RNA樣本(其中DR組選用4例8周齡糖尿病db/db小鼠；對照組選取4例正常仔鼠)，然后與Affymetrix小鼠基因1.0ST陣列雜交。原始數據標準化后，在R語言(3.4.1版)(http://www.bioconductor.org)使用RMA(the robust multichip averaging)算法將其轉換為表達式值。

1.2方法

1.2.1差異表達基因分析采用RankProd軟件包進行差異表達基因(differentially expressed genes，DEGs)分析。RankProd是一種非參數Meta分析工具，其利用每個基因的折疊變化(FC)對每個區域內的基因進行排序和比較。通過對源數據集進行1000次置換計算顯著性值和錯誤發現率(FDR)值。僅假陽性率(PFP)值≤0.05的基因被認為是DR組和對照組之間的DEGs。

1.2.2 DEGs的基因本體分析和通路富集分析使用ClusterProfiler軟件包進一步對所識別的DEGs進行功能解釋，即基因本體(GO)分析和通路分析。GO分析使用0.05的閾值識別顯著豐富的GO術語。通路分析利用KEGG數據庫的超幾何檢驗進行富集分析，閾值為0.05。

1.2.3基因集富集分析使用基因集富集分析(GSEA)統計方法將基因組中基因的非平均強度值與其在生物通路中的發生聯系起來。對于富集研究，本研究允許最小基因集大小為10，以避免隨機注釋。加權富集統計用于計算1000個置換數的富集分數。FDR<0.05的富集基因集被認為是重要的通路。

1.2.4蛋白質相互作用分析使用在線數據庫STRING(10.5版)(https://string-db.org/)對篩選的DEGs進行蛋白質相互作用(PPI)分析。將綜合得分>0.4的配對用于構建PPI網絡，然后使用Cytoscape軟件構建網絡并分析DR中候選DEGs編碼蛋白之間的相互作用關系。為了探索更具體的蛋白質調控關系，本研究在Cytoscape軟件中使用了ClusterONE插件在P<0.1的條件下挖掘聚類模。

1.2.5轉錄因子結合分析基于調控基序和染色質免疫沉淀測序尋找候選轉錄因子(TFs)，使用cytoscapev3.5.1中的iRegulon應用程序搜索DR基因的調節因子。iRegulon通過掃描已知的TFs結合啟動子基序以及從DNA元素百科全書計劃(ENCODE)染色質免疫沉淀測序數據中發現的預測基序檢測TFs及其靶點。在“假定監管區域”“Motif排名數據庫”和“跟蹤排名數據庫”選項的上游設置20kb，其他選項被視為默認選項。DEGs被加載到Cytoscape[5]，并用于查詢iRegulon插件[6]。假定的調控區域選擇以轉錄起始位點(TSS)為中心的20kb為范圍，歸一化富集分數(NES)閾值設為4.0，模體相似性的最大FDR設定為0.001。運行iRegulon并以每個下調和上調的基因尋找TFs。

2結果

2.1 DR的DEGs在穩態條件下對肥胖db/db小鼠視網膜組織與正常小鼠視網膜組織進行全轉錄組分析，通過使用RankProd軟件包共得到66個DEGs，其中38個基因上調，28個基因下調，PFP≤0.05。

2.2 DEGs的GO分析如圖1、2，表1、2所示，在生物學過程組中，上調基因GO項無顯著性差異，下調基因主要集中在眼睛發育、相機型眼發育、視知覺、光刺激感覺知覺和相機型眼晶狀體發育。在細胞組分中，上調基因主要富集于細胞漿核糖體、細胞漿部分、呼吸鏈、核糖體亞單位、核糖體、細胞漿小核糖體亞單位、小核糖體亞單位、血液微粒、中間絲，下調基因主要富集于細胞體纖維。在分子功能中，上調基因主要富集于氧結合、過氧化物酶活性、氧化還原酶活性、過氧化物受體、抗氧化活性、NADH脫氫酶(泛醌)活性、NADH脫氫酶(醌)活性、血紅素結合、NADH脫氫酶活性、四吡咯結合，下調基因主要富集于晶狀體結構成分。上述結果表明，多數DEGs在細胞漿核糖體、眼睛發育和氧結合方面均顯著富集。

表1 DR中上調基因的顯著富集GO分析

表2 DR中下調基因的顯著富集GO分析

圖1 基因本體分析及DR中上調基因GO術語的顯著富集。

圖2 基因本體分析及DR中下調基因的GO術語顯著富集。

2.3通路分析上調基因主要在氧化磷酸化、帕金森病、心肌收縮、非洲錐蟲病、瘧疾和核糖體中富集，而在下調基因中沒有明顯的富集通路(表3)。

表3 db/db小鼠視網膜上調基因的顯著富集通路

2.4基因富集分析GSEA用于識別DR中的重要通路。FDR<0.05時，發現了12條上調通路和6條下調通路(表4、5)。上調通路主要包括脂肪消化吸收、脂肪細胞因子信號轉導、組氨酸代謝、谷胱甘肽代謝、氧化磷酸化。

表4 GSEA分析判定的上調通路

表5 GSEA分析判定的下調通路

2.5 PPI網絡構建利用STRING構建一個由40個基因和61個相互作用組成的PPI網絡(圖3)。Top2a、Cryaa、mt-Cytb、Mip、Rps3、Rpl27a、Rpl13a、Rpl23、Cryba1和Crygd是最顯著的10個節位基因，其可能在DR的病理進程中起重要作用(表6)。應用Cytoscape軟件中的ClusterONE插件對DEGs的PPI網絡進行聚類模塊分析，以預測蛋白質復合物，結果確定了5個DEGs模塊(圖3，表7)。模塊簇1富含參與氧化磷酸化的基因；模塊簇2富集于中間絲細胞骨架組織；模塊簇5在眼睛發育方面豐富。

表6 PPI網絡前10節點的度

表7 PPI網絡中標識的模塊簇

圖3 利用db/db小鼠視網膜DEGs構建PPI網絡 從PPI網絡中識別模塊簇，節點表示DEGs，線表示DEGs之間的度。

2.6 iRegulon預測調控基因集的TFs利用Cytoscape軟件中的iRegulon應用程序搜索DR相關基因集的TFs，對于每個上調和下調的基因集，基于上游20kb的motif預測了上調基因集的5個TFs和下調基因集的8個TFs。對上調基因的分析表明，Runx家族、Ifap2a和Ppara等靶基序的富集，表明這些TFs可能驅動差異基因表達(表8)。與下調基因相關的TFs包括Gata2、Hoxa13、Tbp和Runx3等(表9)。

表8 與上調基因相關的TFs預測

表9 與下調基因相關的TFs預測

3討論

DR是糖尿病最常見的微血管并發癥，也是眼科常見致盲性眼病之一，其發病的分子機制仍有待探究。db/db小鼠視網膜基因表達譜的應用，為本研究提供了從生物信息學方面研究DR發病機制的可能性，對基因芯片數據進行進一步的挖掘有助于從中發現新的生物學標記和診斷/治療靶點基因。

本研究對基因芯片GSE55389進行生物信息學軟件分析，共發現66個DEGs，并對這些DEGs進行GO富集分析、KEGG信號通路分析和PPI網絡分析。為避免只選擇差異明顯的基因而導致某些可能的信號通路被遺漏，本研究同時進行GSEA分析。利用iRegulon分析預測DR相關基因的TFs，其中，上調表達的有Runx2、Ifap2a、Ppara、MafA；下調表達的有Gata2、Hoxa13、Tbp和Runx3等，表明這些TFs可能驅動差異基因表達。Runx2是轉錄因子RUNX家族的一員，編碼1個帶有Runt-DNA結合域的核蛋白。研究發現，血管鈣化在終末期腎臟疾病中非常普遍，PARP1通過上調Runx2促進血管鈣化[7]；血管內皮生長因子(VEGF)直接作用于內皮細胞參與血管發生，是血管發生過程中最重要的血管生長因子，而Runx2是VEGF啟動子的調控轉錄因子[8]。過氧化物酶體增殖物激活受體(peroxisome proliferator-activated receptor，PPARs)是調節特定基因的配體依賴的核受體超家族中的成員，包括α、β、γ三個亞型。它們參與脂質和糖代謝的調節。研究證明PPARs在調節過氧化物酶體增殖、細胞分化、能量代謝、炎癥反應及血管保護作用等過程中扮演著重要角色[9-10]。

Maf家族蛋白是一類重要的轉錄因子，在調節細胞內各種蛋白的表達和細胞分化、細胞凋亡等發揮作用。目前已知有4種類型的Maf—MafA、MafB、c-Maf和NRL。其中，MafA在胰腺β細胞中特異表達，對胰島素的轉錄和分泌至關重要。作為一種葡萄糖調節的胰島β細胞特異性胰島素基因轉錄激活因子[11]，MafA能特異性結合胰島素增強元件RIPE3b并激活胰島素基因表達[12]，并且是體內葡萄糖刺激胰島素分泌的一種關鍵調節因子[13]。在db/db模型小鼠的研究表明，在2型糖尿病進展的早期階段，MafA蛋白在胰島β細胞中丟失[14]。Gata2是調節內皮細胞內皮素-1基因表達的轉錄激活劑，可結合共識序列5’-AGATAG-3’。該基因編碼的蛋白質在調節參與造血和內分泌細胞系發育和增殖的基因轉錄中起著至關重要的作用[15]。HOXA13在脊椎動物中編碼一類被稱為同源盒基因的轉錄因子的基因，被發現在4條不同的染色體上的A、B、C、D簇中。這些蛋白質的表達在胚胎發育過程中受到時空調控。有研究發現，HOXA13在發育過程中與胎盤血管形成模式和迷路內皮規范有關[16]。研究表明，Runx3的低表達與各類腫瘤異常血管增生、轉移[17]和肝臟異常分化[18]有關。糖尿病患者內皮祖細胞功能也與Runx3表達異常相關[19]。

綜上所述，本研究通過對db/db小鼠視網膜和正常小鼠視網膜組織樣本的基因表達譜進行全面的生物信息學分析，確定了以db/db小鼠為模型的DR的關鍵基因和通路，進一步分析了DR中轉錄因子與視網膜微血管增殖的可能分子機制，這些DEGs可能參與DR發生、發展的多種通路。除了上述生物信息學的探索外，還需要后續通過確鑿的實驗工作來確認所識別基因的功能。