胰島素分泌細(xì)胞誘導(dǎo)過(guò)程中miRNA-mRNA 調(diào)控網(wǎng)絡(luò)的生物信息學(xué)分析

王 濤，潘 鑫

（錦州醫(yī)科大學(xué)醫(yī)療學(xué)院醫(yī)學(xué)系，遼寧 錦州 121013）

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

1 資料與方法

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

1.2 GO 功能和KEGG 路徑分析 選用DAVID（https://david.ncifcrf.gov/tools.jsp）數(shù)據(jù)庫(kù)進(jìn)行功能富集分析，主要分析GO 功能以及KEGG 路徑，以P＜0.05為差異有統(tǒng)計(jì)學(xué)意義。

1.3 蛋白互作網(wǎng)絡(luò)分析 選用STRING（https://stringdb.org/）數(shù)據(jù)庫(kù)進(jìn)行蛋白互作網(wǎng)絡(luò)（protein-protein interaction，PPI）分析，選擇Homo sapiens 作為“Organism”的條件，其余參數(shù)均為系統(tǒng)默認(rèn)設(shè)置。

1.4 miRNA 預(yù)測(cè) 針對(duì)差異表達(dá)基因預(yù)測(cè)其靶miRNA，選用miRTarBase（https://mirtarbase.cuhk.edu.cn/）數(shù)據(jù)庫(kù)，選擇reporter assay、western blot、qPCR 和microarray4 種實(shí)驗(yàn)驗(yàn)證的miRNA-mRNA靶向匹配作為預(yù)測(cè)miRNA 的入選標(biāo)準(zhǔn)，并應(yīng)用Cytoscape 3.9.1 軟件對(duì)miRNA-mRNA 的調(diào)控網(wǎng)絡(luò)進(jìn)行可視化。

1.5 預(yù)測(cè)miRNA 的驗(yàn)證 將預(yù)測(cè)的miRNA 與PMID：20735361[9]的Table S4 數(shù)據(jù)取交集，分析miRNA 在人胚胎T3 干細(xì)胞向胰腺胰島樣細(xì)胞團(tuán)分化過(guò)程中的表達(dá)，其Table S4 數(shù)據(jù)顯示出胚胎T3 干細(xì)胞（hES-T3 cells grown on mouse embryonic fibroblast feeder，T3ES）、胚胎內(nèi)胚層（embryoid bodies differentiated from T3 cells，T3EB）和胰腺胰島樣細(xì)胞團(tuán)（pancreatic islet-like cell clusters derived from T3 cells，T3pi）3 組標(biāo)準(zhǔn)化的miRNA 表達(dá)量。

1.6 轉(zhuǎn)錄因子-靶基因調(diào)控網(wǎng)絡(luò)的構(gòu)建 針對(duì)驗(yàn)證的miRNA 以及顯著下調(diào)基因查找調(diào)控其表達(dá)的轉(zhuǎn)錄因子，選用FunRich 3.1.3 軟件分析轉(zhuǎn)錄因子，以P＜0.05 為差異有統(tǒng)計(jì)學(xué)意義，并作為轉(zhuǎn)錄因子的入選標(biāo)準(zhǔn)，同時(shí)應(yīng)用Cytoscape 3.9.1 插件iRegulon 預(yù)測(cè)轉(zhuǎn)錄因子，將二者分別與Mahmoud H 等[10]報(bào)道的β細(xì)胞分化相關(guān)的重要轉(zhuǎn)錄因子取交集，對(duì)轉(zhuǎn)錄因子-靶基因調(diào)控網(wǎng)絡(luò)可視化。

2 結(jié)果

2.1 差異及因及其GO 功能和KEGG 路徑分析 共得到差異表達(dá)的基因188 個(gè)，選取誘導(dǎo)后顯著下調(diào)的基因共21 個(gè)（表1），對(duì)其進(jìn)行GO 功能和KEGG路徑分析，發(fā)現(xiàn)下調(diào)基因中POU5F1、DNMT3B、NANOG 和UCA1 集中在基因表達(dá)調(diào)控，POU5F1 和NANOG 與內(nèi)胚層的決定以及成體干細(xì)胞種群維持有關(guān)，PPP1R16B 和PPP2R2C 與蛋白磷酸化調(diào)節(jié)活性有關(guān)。DAVID 數(shù)據(jù)庫(kù)進(jìn)行KEGG 路徑顯著性分析并沒(méi)有找到有統(tǒng)計(jì)學(xué)意義的信號(hào)通路。

表1 表達(dá)下調(diào)的差異基因

2.2 蛋白互作網(wǎng)絡(luò)分析 通過(guò)PPI 分析21 個(gè)誘導(dǎo)后顯著下調(diào)的基因，發(fā)現(xiàn)LOC729860、UCA1 和HLADPB2 未被STRING 數(shù)據(jù)庫(kù)識(shí)別，因此該網(wǎng)絡(luò)有18個(gè)節(jié)點(diǎn)，14 條邊，其PPI 富集P值為1.0e-09，差異有統(tǒng)計(jì)學(xué)意義。進(jìn)一步分析發(fā)現(xiàn)POU5F1、UTF1、DNMT3B、NANOG、GABRB3、DPPA2 和TCL1B 形成互作網(wǎng)絡(luò)，剩余11 個(gè)蛋白相互之間沒(méi)有關(guān)聯(lián)，見(jiàn)圖1。

圖1 蛋白互作網(wǎng)絡(luò)分析

2.3 miRNA-mRNA 調(diào)控網(wǎng)絡(luò) 通過(guò)miRTarBase 數(shù)據(jù)庫(kù)查找靶向作用于顯著下調(diào)基因的miRNA，選擇reporter assay、western blot、qPCR 和microarray 4 種實(shí)驗(yàn)驗(yàn)證的miRNA-mRNA 靶向匹配作為預(yù)測(cè)miRNA 的入選標(biāo)準(zhǔn)，共得到25 個(gè)預(yù)測(cè)的miRNA，應(yīng)用Cytoscape 3.9.1 軟件對(duì)miRNA-mRNA 調(diào)控網(wǎng)絡(luò)進(jìn)行可視化（圖2），發(fā)現(xiàn)PPI 網(wǎng)絡(luò)的7 個(gè)關(guān)聯(lián)蛋白也處于該調(diào)控網(wǎng)絡(luò)的核心區(qū)域，預(yù)測(cè)的大部分miRNA 圍繞其展開(kāi)，其中，miR-335-5p 能夠調(diào)控POU5F1、DNMT3B、NANOG、FAM124B 和LECT1 的表達(dá)，進(jìn)一步在GenomeNet 數(shù)據(jù)庫(kù)（https://www.genome.jp/pathway/hsa04550+5460）中對(duì)這5 個(gè)基因進(jìn)行KEGG 路徑分析，發(fā)現(xiàn)POU5F1 和NANOG 處于干細(xì)胞多能性調(diào)控信號(hào)通路（hsa04550）之中，調(diào)控干細(xì)胞的自我更新。

圖2 miRNA-mRNA 調(diào)控網(wǎng)絡(luò)

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

圖3 3 種方法預(yù)測(cè)的調(diào)控網(wǎng)絡(luò)

表2 PMID：20735361 結(jié)果驗(yàn)證預(yù)測(cè)的miRNA

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

圖4 轉(zhuǎn)錄因子-靶基因調(diào)控網(wǎng)絡(luò)

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

圖5 miRNA 的轉(zhuǎn)錄因子

3 討論

本研究通過(guò)對(duì)GSE42094 的6 組樣本數(shù)據(jù)分析，得到188 個(gè)差異表達(dá)的基因，選取誘導(dǎo)后顯著下調(diào)的21 個(gè)基因進(jìn)行后續(xù)研究。選用DAVID 數(shù)據(jù)庫(kù)對(duì)其進(jìn)行GO 功能和KEGG 路徑顯著性分析，發(fā)現(xiàn)其中POU5F1 和NANOG 基因與基因表達(dá)調(diào)控、內(nèi)胚層的決定和成體干細(xì)胞種群維持3 個(gè)功能均有相關(guān)性，但沒(méi)有找到有統(tǒng)計(jì)學(xué)意義的KEGG 路徑。進(jìn)一步在GenomeNet 數(shù)據(jù)庫(kù)中檢索POU5F1 和NANOG 基因的KEGG 路徑，發(fā)現(xiàn)二者同處于干細(xì)胞多能性調(diào)控信號(hào)通路，能夠調(diào)控干細(xì)胞的自我更新，這與其GO 功能相呼應(yīng)。另外，該通路與TGF-β、MAPK、PI3K-Akt 以及Wnt 4 條信號(hào)通路相關(guān)聯(lián)，這與Mahmoud H 等[10]報(bào)道的體外β 細(xì)胞分化的關(guān)鍵信號(hào)通路相一致。

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

2021 年，Sabouri E 等[21]報(bào)道m(xù)iR-375、miR-7、miR-186、miR-26a、miR-690 以及miR-302a 6 個(gè)IPCs 誘導(dǎo)分化過(guò)程中的重要miRNA。本研究同樣發(fā)現(xiàn)miR-375 和miR-26a 能夠靶向調(diào)控DNMT3B 基因表達(dá)、miR-7 能夠調(diào)控TCL1B 以及miR-302a 能夠調(diào)控NANOG 基因表達(dá)。可見(jiàn)該miRNA-mRNA調(diào)控網(wǎng)絡(luò)具備一定可信度，通過(guò)進(jìn)一步查找雙方的轉(zhuǎn)錄因子也可部分地佐證，因此，基于顯著下調(diào)基因構(gòu)建的該miRNA-mRNA 調(diào)控網(wǎng)絡(luò)在干細(xì)胞誘導(dǎo)分化過(guò)程中發(fā)揮重要作用。

綜上所述，本研究基于顯著下調(diào)基因構(gòu)建的miRNA-mRNA 調(diào)控網(wǎng)絡(luò)參與干細(xì)胞誘導(dǎo)分化進(jìn)程并發(fā)揮重要作用。