胰島素分泌細胞誘導過程中miRNA-mRNA 調控網絡的生物信息學分析

王 濤，潘 鑫

（錦州醫科大學醫療學院醫學系，遼寧 錦州 121013）

miRNA 作為一類小的非編碼RNA，其長度大約在22 個核苷酸左右，能夠與靶mRNA 3'端非翻譯區互補，進而降解或阻止該mRNA 翻譯出蛋白質，實現轉錄后水平的基因表達調控[1]。糖尿病作為一種慢性疾病，部分是由胰島β 細胞受損而引起，胰腺或胰島移植對于該類疾病治療效果較好，但是由于供體缺乏以及免疫排斥導致臨床應用受限[2，3]，而干細胞分化為胰島素分泌細胞（insulin-producing cells，IPCs）再行移植，有望徹底根治[4，5]。生物信息學通過整合多個學科領域信息和知識，能夠從整體層面揭示疾病或生物過程的復雜分子機制[6-8]。本研究運用生物信息學方法分析人胚胎干細胞（human embryonic stem cells，hESCs）向IPCs 誘導分化過程中mRNA 的表達譜，預測靶向作用于該mRNA 的miRNA，構建miRNA-mRNA 調控網絡。通過對調控網絡的分析，擬找到IPCs 誘導分化的關鍵節點，為進一步開展分子靶向干預提供新思路。

1 資料與方法

1.1 資料來源 選擇GEO（https://www.ncbi.nlm.nih.gov/geo/）數據庫內的GSE42094 數據集作為研究對象，該數據集包含23 個樣本數據（GSM1032319-41），涉及未分化hESCs、IPCs 分化的5 個時期、胰腺內胰島和胎胰。選擇未分化hESCs（hESCs 組）和IPCs 分化的5 個時期（Diff1、Diff2、Diff3、Diff4 和IPCs 組）研究mRNA 的表達譜。應用數據庫自帶軟件GEO2R對6 組數據進行分析，獲取差異表達的mRNA。

1.2 GO 功能和KEGG 路徑分析 選用DAVID（https://david.ncifcrf.gov/tools.jsp）數據庫進行功能富集分析，主要分析GO 功能以及KEGG 路徑，以P＜0.05為差異有統計學意義。

1.3 蛋白互作網絡分析 選用STRING（https://stringdb.org/）數據庫進行蛋白互作網絡（protein-protein interaction，PPI）分析，選擇Homo sapiens 作為“Organism”的條件，其余參數均為系統默認設置。

1.4 miRNA 預測 針對差異表達基因預測其靶miRNA，選用miRTarBase（https://mirtarbase.cuhk.edu.cn/）數據庫，選擇reporter assay、western blot、qPCR 和microarray4 種實驗驗證的miRNA-mRNA靶向匹配作為預測miRNA 的入選標準，并應用Cytoscape 3.9.1 軟件對miRNA-mRNA 的調控網絡進行可視化。

1.5 預測miRNA 的驗證 將預測的miRNA 與PMID：20735361[9]的Table S4 數據取交集，分析miRNA 在人胚胎T3 干細胞向胰腺胰島樣細胞團分化過程中的表達，其Table S4 數據顯示出胚胎T3 干細胞（hES-T3 cells grown on mouse embryonic fibroblast feeder，T3ES）、胚胎內胚層（embryoid bodies differentiated from T3 cells，T3EB）和胰腺胰島樣細胞團（pancreatic islet-like cell clusters derived from T3 cells，T3pi）3 組標準化的miRNA 表達量。

1.6 轉錄因子-靶基因調控網絡的構建 針對驗證的miRNA 以及顯著下調基因查找調控其表達的轉錄因子，選用FunRich 3.1.3 軟件分析轉錄因子，以P＜0.05 為差異有統計學意義，并作為轉錄因子的入選標準，同時應用Cytoscape 3.9.1 插件iRegulon 預測轉錄因子，將二者分別與Mahmoud H 等[10]報道的β細胞分化相關的重要轉錄因子取交集，對轉錄因子-靶基因調控網絡可視化。

2 結果

2.1 差異及因及其GO 功能和KEGG 路徑分析 共得到差異表達的基因188 個，選取誘導后顯著下調的基因共21 個（表1），對其進行GO 功能和KEGG路徑分析，發現下調基因中POU5F1、DNMT3B、NANOG 和UCA1 集中在基因表達調控，POU5F1 和NANOG 與內胚層的決定以及成體干細胞種群維持有關，PPP1R16B 和PPP2R2C 與蛋白磷酸化調節活性有關。DAVID 數據庫進行KEGG 路徑顯著性分析并沒有找到有統計學意義的信號通路。

表1 表達下調的差異基因

2.2 蛋白互作網絡分析 通過PPI 分析21 個誘導后顯著下調的基因，發現LOC729860、UCA1 和HLADPB2 未被STRING 數據庫識別，因此該網絡有18個節點，14 條邊，其PPI 富集P值為1.0e-09，差異有統計學意義。進一步分析發現POU5F1、UTF1、DNMT3B、NANOG、GABRB3、DPPA2 和TCL1B 形成互作網絡，剩余11 個蛋白相互之間沒有關聯，見圖1。

圖1 蛋白互作網絡分析

2.3 miRNA-mRNA 調控網絡 通過miRTarBase 數據庫查找靶向作用于顯著下調基因的miRNA，選擇reporter assay、western blot、qPCR 和microarray 4 種實驗驗證的miRNA-mRNA 靶向匹配作為預測miRNA 的入選標準，共得到25 個預測的miRNA，應用Cytoscape 3.9.1 軟件對miRNA-mRNA 調控網絡進行可視化（圖2），發現PPI 網絡的7 個關聯蛋白也處于該調控網絡的核心區域，預測的大部分miRNA 圍繞其展開，其中，miR-335-5p 能夠調控POU5F1、DNMT3B、NANOG、FAM124B 和LECT1 的表達，進一步在GenomeNet 數據庫（https://www.genome.jp/pathway/hsa04550+5460）中對這5 個基因進行KEGG 路徑分析，發現POU5F1 和NANOG 處于干細胞多能性調控信號通路（hsa04550）之中，調控干細胞的自我更新。

圖2 miRNA-mRNA 調控網絡

2.4 預測miRNA 的驗證 將預測的miRNA 與PMID：20735361 的Table S4 數據取交集，得到24 個驗證的miRNA，只有1 個預測的miRNA（miR-200b-3p）未被驗證，大部分預測的miRNA 得到驗證，即存在于干細胞及其誘導分化的過程中，表明這種基于4 種實驗預測miRNA 的方法具備一定的可信度（表2）。表2 顯示miR-26b-5p、miR-106b-5p、miR-200c-3p、miR-148a-3p、miR-375、miR-204-5p、miR-7-5p、miR-302a-3p 和miR-124-3p 表達量均在胚胎內胚層階段（T3EB）顯著升高而又在胰腺胰島樣細胞團階段（T3pi）顯著降低，miR-26a-5p、miR-335-5p 和miR-134-5p 表達量隨著誘導進程逐漸升高。這兩類miRNA 的表達特點與顯著下調基因表達量在T3EB階段后顯著降低形成對照，提示二者之間存在負性調控。結合圖2 分析發現miR-26b-5p 和miR-335-5p 將miRNA-mRNA 調控網絡核心區域的7 個基因往外延伸，其中，LRRC2 通過miR-26b-5p、FAM124B 和LECT1 通過miR-335-5p 與核心區域關聯。進一步將僅通過microarray 實驗預測的miRNA 剔除后作調控網絡圖（圖3），發現大部分miRNA 主要調控POU5F1、DNMT3B 和NANOG 的表達，該部分miRNA 是通過reporter assay、western blot 和qPCR 實驗驗證了其與靶基因的匹配關系，其中，miR-200c-3p、miR-148a-3p、miR-204-5p、miR-302a-3p、miR-26a-5p、miR-134-5p 和miR-335-5p與其靶基因表達模式存在負性調控，可得到進一步驗證，miR-369-5p 也被上述3 種實驗所驗證，但是其在整個誘導過程中表達量未發生變化。結合表2發現，僅通過microarray 實驗預測的miRNA，如miR-26b-5p、miR-106b-5p、miR-375、miR-7-5p 和miR-124-3p 與其靶基因的表達模式存在負性關系，或可佐證該部分miRNA 與其靶基因的靶向調控。

圖3 3 種方法預測的調控網絡

表2 PMID：20735361 結果驗證預測的miRNA

2.5 顯著下調基因的轉錄因子 針對21 個顯著下調基因查找調控其表達的轉錄因子，FunRich 軟件沒有得到有統計學意義的轉錄因子，通過iRegulon 分析得到49 個轉錄因子，將其與188 個差異表達基因取交集，得到POU5F1、NANOG 和TFAP2A 3 個轉錄因子，并應用Cytoscape 可視化（圖4A），發現POU5F1 和NANOG 本身即為顯著下調基因，其表達量見圖4B 和4C，其中NANOG 調控9 個基因表達，包括其自身亦即存在順式作用調控，POU5F1 調控6個基因表達，也存在順式作用調控，且二者之間互為對方轉錄因子，存在相互轉錄調控現象，而TFAP2A調控6 個基因表達，其中包含POU5F1 和NANOG基因。ZDHHC22、PPP1R16B 和GABRB3 3 個基因同時受上述3 個轉錄因子的調控，結合POU5F1 和NANOG 基因表達量顯著下降的變化，可以理解該部分顯著下調基因主要存在轉錄水平的調控。進一步將iRegulon 分析得到的49 個轉錄因子與Mahmoud H 等[10]報道的β 細胞分化相關的重要轉錄因子取交集，得到ONECUT1（HNF6）、NEUROD1和NKX2-2 3 個轉錄因子（圖4D），發現ONECUT1能夠調控8 個基因表達，其中包括NANOG 基因，NEUROD1 與 ONECUT1 共同調控 PPP2R2C、ZDHHC22 和LECT1 3 個基因，而LECT1 基因也受NKX2-2 轉錄因子調控，可見該部分基因除受轉錄水平調控外還存在較強的轉錄后水平調控，以至于其表達量在分化后期顯著下降。

圖4 轉錄因子-靶基因調控網絡

2.6 驗證miRNA 的轉錄因子 針對24 個驗證的miRNA 查找調控其表達的轉錄因子，FunRich 軟件得到158 個有統計學意義的轉錄因子，將其與Mahmoud H 等[10]報道的β 細胞分化相關的重要轉錄因子取交集，得到NKX6-1、ONECUT1、PAX6、MAFB、HNF4A、PDX1、PAX4 和MNX1（HB9）8 個轉錄因子（圖5A）。iRegulon 插件分析得到52 個轉錄因子，將其與188 個差異表達基因取交集，得到TFAP2A、NR2F2 和FOXA2 3 個轉錄因子，進一步與Mahmoud H 等報道的β 細胞分化相關的重要轉錄因子取交集，得到ONECUT1、FOXA2 和PAX4 3 個轉錄因子，將上述兩部分轉錄因子融合分析，發現FOXA2 既是β 細胞分化相關的重要轉錄因子也是188 個差異表達基因之一，其表達量見圖5B，可見其在整個誘導分化過程中一直處于高表達狀態，對整個誘導分化進程發揮轉錄水平的調控。應用Cytoscape 對轉錄因子-miRNA 可視化（圖5C），發現轉錄因子TFAP2A 能夠調控7 個miRNA 表達，其中miR-375、miR-148a-3p、miR-124-3p 和miR-200c-3p 表達量均在T3EB 階段顯著升高而又在T3pi 階段顯著降低，這與TFAP2A 的表達模式一致（圖5D），表明存在較強轉錄水平的調控。而miR-142-3p、miR-9-5p 和miR-34a-5p 表達量在整個誘導分化過程中一直下降，提示這3 個miRNA 還受其他轉錄因子的調控或者越到分化后期其轉錄后水平的調控越強以至于檢測到的miRNA 表達量越低。

圖5 miRNA 的轉錄因子

3 討論

本研究通過對GSE42094 的6 組樣本數據分析，得到188 個差異表達的基因，選取誘導后顯著下調的21 個基因進行后續研究。選用DAVID 數據庫對其進行GO 功能和KEGG 路徑顯著性分析，發現其中POU5F1 和NANOG 基因與基因表達調控、內胚層的決定和成體干細胞種群維持3 個功能均有相關性，但沒有找到有統計學意義的KEGG 路徑。進一步在GenomeNet 數據庫中檢索POU5F1 和NANOG 基因的KEGG 路徑，發現二者同處于干細胞多能性調控信號通路，能夠調控干細胞的自我更新，這與其GO 功能相呼應。另外，該通路與TGF-β、MAPK、PI3K-Akt 以及Wnt 4 條信號通路相關聯，這與Mahmoud H 等[10]報道的體外β 細胞分化的關鍵信號通路相一致。

通過構建miRNA-mRNA 調控網絡發現，LOC729860、UTF1、KLKB1、DPPA2、NAP1L3、HLA -DPB2 等6 個基因沒有對應miRNA 調控，而PPP1R16B、RPL36A、VSIG10、ZDHHC22 等4 個 基因雖有相應miRNA 調控，但不滿足設定的入選標準，因此，這10 個基因在應用Cytoscape 3.9.1 軟件進行miRNA-mRNA 調控網絡可視化過程中未有miRNA 呈現。網絡顯示miR-335-5p 能夠調控POU5F1、NANOG、FAM124B、DNMT3B 和LECT1 5個基因的表達。其中，POU5F1 和NANOG 基因較為熟知[11-13]，而FAM124B 在人HeLa 細胞內能夠與染色質解旋酶DNA 結合蛋白7（chromodomain helicase DNA binding protein 7，CHD7）相互作用，而CHD7基因在小鼠ESCs 內與OCT4、NANOG 以及SOX2同處于一個增強子元件，提示FAM124B 作為一個關聯因子參與增強子介導的轉錄[14，15]。DNA 甲基轉移酶3B（DNA methyltransferases 3B，DNMT3B）是DNA 甲基轉移酶蛋白家族成員之一，DNMT1 主要負責維持DNA 甲基化，DNMT3A 和DNMT3B 負責從頭合成甲基轉移酶[16，17]。在hESCs 內，DNMT3B 的表達量是已分化細胞的近30 倍，其缺失能夠影響線粒體的生成以及提高細胞內α-酮戊二酸的水平，進而引起hESCs 的多向分化[18]。白細胞源性趨化蛋白1（leukocyte cell-derived chemotaxin 1，LECT1）又稱軟骨源性抑制因子-1（chondromodulin-1），是一個25 kDa 的糖蛋白，由跨膜蛋白經轉錄后修飾生成，其有助于蛋白聚糖的合成以及體外軟骨細胞的生長[19]。在hESCs 內，LECT1 表達量是已分化細胞的近40 倍，其和POU5F1、NANOG 一樣被作為干性維持基因[20]。因此，miR-335-5p 所靶向調控的這5 個基因均與干性維持有關，提示miR-335-5p 通過負性調控該類基因的表達，抑制hESCs 的干細胞多能性調控信號通路，促使其轉向分化。另外，miR-335-5p 表達量在誘導分化過程中逐漸升高已由文獻數據所驗證，這與上述干性維持基因的顯著下調形成對照。

2021 年，Sabouri E 等[21]報道miR-375、miR-7、miR-186、miR-26a、miR-690 以及miR-302a 6 個IPCs 誘導分化過程中的重要miRNA。本研究同樣發現miR-375 和miR-26a 能夠靶向調控DNMT3B 基因表達、miR-7 能夠調控TCL1B 以及miR-302a 能夠調控NANOG 基因表達。可見該miRNA-mRNA調控網絡具備一定可信度，通過進一步查找雙方的轉錄因子也可部分地佐證，因此，基于顯著下調基因構建的該miRNA-mRNA 調控網絡在干細胞誘導分化過程中發揮重要作用。

綜上所述，本研究基于顯著下調基因構建的miRNA-mRNA 調控網絡參與干細胞誘導分化進程并發揮重要作用。