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

李 倩，何 飛，張 然

0引言

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

1材料和方法

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

1.2方法

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

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

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

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

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

2結(jié)果

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

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

表1 DR中上調(diào)基因的顯著富集GO分析

表2 DR中下調(diào)基因的顯著富集GO分析

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

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

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

表3 db/db小鼠視網(wǎng)膜上調(diào)基因的顯著富集通路

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

表4 GSEA分析判定的上調(diào)通路

表5 GSEA分析判定的下調(diào)通路

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

表6 PPI網(wǎng)絡前10節(jié)點的度

表7 PPI網(wǎng)絡中標識的模塊簇

圖3 利用db/db小鼠視網(wǎng)膜DEGs構(gòu)建PPI網(wǎng)絡 從PPI網(wǎng)絡中識別模塊簇，節(jié)點表示DEGs，線表示DEGs之間的度。

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

表8 與上調(diào)基因相關的TFs預測

表9 與下調(diào)基因相關的TFs預測

3討論

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

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

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

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