基于生物信息学分析胰岛素分泌细胞诱导过程中miRNA-mRNA 调控网络

王 涛，侯润博，潘 鑫

(锦州医科大学医疗学院医学系1，医学技术系2，辽宁 锦州 121013)

Ⅰ型糖尿病作为一种慢性疾病主要是由胰岛β细胞受损引起，胰腺或胰岛移植可以彻底根治，但是由于供体缺乏以及免疫排斥大大限制该治疗方案的临床应用[1，2]，干细胞分化为胰岛素分泌细胞（insulinproducing cells，IPCs）再行移植，有望成为Ⅰ型糖尿病治疗的新策略[3，4]。miRNA 作为一类小的非编码RNA，其长度大约在22 个核苷酸左右，能够与靶mRNA 3’端非翻译区互补，进而降解或阻止该mRNA 翻译出蛋白质，实现转录后水平的基因表达调控[5，6]。生物信息学通过整合多个学科领域信息和知识，能够整体层面揭示疾病或生物过程的复杂分子机制[7]，而且其通过对已有实验数据的分析，能够减少重复实验所带来的资源浪费。GEO 数据库是生物信息学研究重要的中心资源之一，其存储大量高通量数据，但是大部分数据并未被充分利用[8]。目前，运用生物信息学分析干细胞向IPCs 诱导分化过程中miRNA-mRNA 的调控网络报道较少见。因此，本研究运用生物信息学方法分析干细胞向IPCs 诱导分化过程中的mRNA 表达谱，并构建miRNAmRNA 调控网络。

1 资料与方法

1.1 一般资料 通过NCBI 的GEO 数据库（https://www.ncbi.nlm.nih.gov/geo/）检索人胚胎干细胞（human embryonic stem cells，hESCs）向IPCs 诱导分化过程中mRNA 表达谱的相关数据，最终只有GSE42094 数据集符合要求，该数据集基于GPL10558 平台，包含23 个样本数据（GSM1032319-41），涉及未分化hESCs（WA09）、IPCs 分化的5 个时期、胰腺内胰岛和胎胰。

1.2 数据处理与差异表达分析 选择GSE42094 数据集内的未分化hESCs 和IPCs 分化的5 个时期进行研究，进入“Analyze with GEO2R”项目，点击“Define groups”设定分组，其中，hESCs 组包含3 个样本，Diff1 组3 个，Diff2 组3 个，Diff3 组2 个，Diff4 组2个，IPCs 组3 个，对6 组数据进行差异表达分析，筛选标准设定为adj.P＜0.05。

1.3 GO 功能和KEGG 路径显著性分析 为分析差异表达基因的潜在生物学功能，选用DAVID（https://david.ncifcrf.gov/tools.jsp）数据库进行功能富集分析，主要分析GO 功能以及KEGG 路径2 个方面，以P＜0.05 为差异有统计学意义。

1.4 蛋白互作网络分析 为分析差异表达基因的相互作用及其调控网络，选用STRING（https://stringdb.org/）数据库进行蛋白互作网络PPI 分析，选择Homo sapiens 作为“Organism”的条件，其余参数均为系统默认设置。

1.5 miRNA 预测 针对差异表达基因预测其靶miRNA，选用miRTarBase（https://mirtarbase.cuhk.edu.cn/）数据库，选择reporter assay、western blot、qPCR 和microarray 四种实验验证的miRNA-mRNA 靶向匹配作为预测miRNA 的入选标准，并应用Cytoscape 9.1.0 软件对miRNA-mRNA 的调控网络进行可视化。

1.6 预测miRNA 的验证 将预测的miRNA 与PMID:20 735 361[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）三组标准化的miRNA 表达量。

2 结果

2.1 差异表达基因 通过对GSE42094 的6 组样本数据分析，差异表达基因的前250 个被显示，除去没有名称和重复的基因，共得到188 个。进一步对hESCs和IPCs 两组数据分析得到差异表达基因，并与上述188 个基因取交集后，选取诱导后上调和下调的各前10 个差异表达基因进行后续研究，见表1。该20个基因代表整个hESCs 向IPCs 诱导过程中表达量差异最具有统计学意义的基因。

表1 诱导后上调和下调的前10 个差异表达基因

2.2 GO 功能和KEGG 路径显著性分析 通过对诱导后上调和下调的各前10 个差异表达基因进行GO 功能和KEGG 路径显著性分析，发现上调基因在“外泌体”和“胞外”比较集中，主要涉及“谷胱甘肽代谢”和“PPAR 信号通路”，而下调基因在“基因表达调控”“内胚层的决定”以及“成体干细胞种群维持”比较集中。

2.3 蛋白互作网络分析 通过PPI 分析发现10 个上调基因中有6 个基因所表达的蛋白质具有相互作用，有2 个基因显示二者间直接作用（GSTA1 和GSTA2），剩余两个没有关联（图1A）。而下调的10 个基因中LOC729860 和UCA1 在STRING 数据库中没有收录，其余8 个基因中有5 个表达的蛋白质具有相互作用，剩余3 个没有关联（图1B）。将上调和下调的基因进行融合分析，发现下调基因KLKB1 能够与上调基因FGG 直接进行蛋白相互作用，且上调基因蛋白互作网络通过VTN 与下调基因网络的POU5F1 直接作用，将两个网络相关联（图1C）。

图1 蛋白互作网络分析

2.4 miRNA-mRNA 调控网络 通过miRTarBase 数据库查找靶向作用于差异表达基因的miRNA，选择reporter assay、western blot、qPCR 和microarray 4 种实验验证的miRNA-mRNA 靶向匹配作为预测miRNA 的入选标准，发 现 CDH17、APOA2、FAM162B、GSTA2、LOC729860、UTF1、KLKB1 等7个基因没有对应miRNA 调控，而APOA1、TTR、PPP1R16B、RPL36A 4 个基因虽有相应miRNA 调控，但不满足设定的入选标准，因此，该11 个基因在应用Cytoscape 9.1.0 软件进行miRNA-mRNA 调控网络可视化过程中，未有miRNA 呈现，其余9 个基因共对应23 个预测的miRNA。其中，DNMT3B 基因显示有10 个miRNA 调控其表达，NANOG 基因有9个，而SPINK1 和UCA1 基因仅各有1 个miRNA调控其表达，且该2 个基因与miRNA-mRNA 调控网络没有直接关联。进一步分析发现miR-335-5p能够调控VTN、NANOG、POU5F1 和DNMT3B 4 个基因的表达（图2），其中VTN 为上调基因（图3A），NANOG、POU5F1 和DNMT3B 为下调基因（图3B～图3D）。

图2 miRNA-mRNA 调控网络

图3 GEO2R 所示基因的表达量

2.5 预测miRNA 的验证 将预测的miRNA 与PMID:20 735 361 的Table S4 数据取交集，得到21 个验证的miRNA，剩余2 个预测的miRNA（miR-200b-3p和miR-144-3p）未被验证，可见绝大部分预测的miRNA 得到验证，亦即存在于干细胞及其诱导分化的过程中（表2），同时表明基于4 种实验预测miRNA 的方法具备一定的可信度。表2 显示miR-302a-3p、miR-34a-5p、miR-34c-5p、miR-29b-3p、miR-375 和miR-29a-3p 表达量诱导分化后显著降低，而miR-134-5p、miR-335-5p 和miR-26a-5p 表达量随着诱导分化进程逐渐升高，其中miR-335-5p表达量的逐渐升高与NANOG、POU5F1 和DNMT3B表达量的逐渐降低（图3B～图3D）形成对照，提示二者之间存在靶向关系且为负性调控。结合图2 进一步分析发现上调基因与下调基因网络通过多个miRNA 相关联，其中上调基因VTN 通过miR-26b-5p和miR-335-5p 与下调基因网络多个基因关联，并且VTN 通过miR-26b-5p 与GSTA1 和GSTA2 的二者间直接作用相关联。上调基因FGG 通过miR-29家族（miR-29a-3p、29b-3p、29c-3p）与下调基因DNMT3B 相关联，提示miRNA 通过靶向不同表达类型的基因进而调控干细胞的诱导分化进程。

表2 PMID:20 735 361 数据验证预测的miRNA

表2 （续）

2.6 miR-335-5p 的结构和功能信息 运用miRBase数据库（https://www.mirbase.org/）检索miR-335-5p的结构信息，其前体hsa-mir-335 的序列及其茎环结构如图4A 所示，标红部分为成熟miR-335-5p 的序列（16-UCAAGAGCAAUAACGAAAAAUGU-38）。进一步检索发现前体hsa-mir-335 在169 篇文献中被提及，共涉及864 个词句，其词汇云如图4B 所示，可见其在癌症和肿瘤方向研究最广泛，其次就是诱导和分化等方向，与本研究的方向相符。

图4 miR-335-5p 的结构和功能信息

3 讨论

本研究为分析hESCs 向IPCs 诱导分化过程中的mRNA 表达谱，采用生物信息学方法筛选GEO数据库内的数据集得到GSE42094，该数据集包含未分化hESCs、IPCs 分化的5 个时期、胰腺内胰岛和胎胰等样本数据。本研究选择与研究目的直接相关的未分化hESCs 和IPCs 分化的5 个时期来研究，通过数据库自带软件GEO2R 对数据进行分析，得出差异表达的mRNA 共188 个，分别选取诱导后上调和下调的各前10 个差异表达基因作为后续研究。

通过GO 功能和KEGG 路径显著性分析，发现下调基因多涉及基因表达调控以及多能干细胞调控信号通路等与干细胞诱导分化机制相关的内容，而上调基因涉及内容则比较宽泛。进一步分析差异表达基因的相互作用，发现上调基因蛋白互作网络与下调基因网络通过上调基因VTN 和下调基因POU5F1 的直接作用而相关联，提示VTN 和POU5F1 基因可能为该蛋白互作网络的核心基因。玻连蛋白（Vitronectin，VTN）是一种细胞外基质糖蛋白，属于整合素家族成员之一，其在伤口愈合和凝血过程中起着重要作用，进一步研究发现其参与多种生物学过程包括细胞粘附和迁移[10-12]，也能够通过JNK 和ERK 信号通路诱导INS-1 细胞的增殖[13]。本研究发现玻连蛋白在Diff3 组表达量开始升高，在IPCs 组达到峰值，表明其参与干细胞的诱导分化过程。POU5F1 又称OCT4，是一种维持干细胞多能性的转录因子[14，15]，其与NANOG 在胚胎早期发育中起着关键作用，并可作为胚胎干细胞的标志物以及干性维持[16，17]。本研究发现hESCs 组POU5F1 表达量最高，随着诱导进程逐渐下降，在Diff3 组表达量就已经处于较低水平。

构建miRNA-mRNA 调控网络，发现miR-335-5p能够调控VTN、NANOG、POU5F1 和DNMT3B 基因的表达，其中VTN 为上调基因，NANOG、POU5F1 和DNMT3B 为下调基因。有学者研究发现在糖尿病模型GK 大鼠的胰岛内miR-335-5p 表达量显著增高，其能够直接靶向抑制Stxbp1 mRNA 的表达，进一步研究发现miR-335-5p 在INS-1 832/13 细胞内的过表达能够抑制葡萄糖刺激的胰岛素分泌[18]。Tang XW 等[19]发现在妊娠期糖尿病小鼠中过表达miR-335-5p 能够抑制VASH1 的表达进而降低胰岛素释放的水平，并通过激活TGF-β 信号通路诱发胰岛素抵抗。最近，Li G 等[20]研究发现miR-335-5p 能够下调SLC2A4 基因表达进而抑制胰腺细胞的生长，与二型糖尿病的进展密切相关。但是，miR-335-5p是否参与干细胞诱导分化为IPCs 未见报道。

通过将预测的miRNA 与文献PMID:20 735 361数据取交集，得到18 个验证的miRNA，miR-335-5p包含其中，结果显示miR-335-5p 在内胚层阶段（T3EB）其表达量升高2 倍多，至胰腺胰岛样细胞团阶段（T3pi）升高接近5 倍，可见其在整个干细胞及其诱导分化过程中均有表达，且随着诱导分化进程逐渐升高，表明miR-335-5p 不但参与糖尿病的发生进展，而且参与干细胞诱导分化为IPCs 的整个过程。

综上所述，本研究共筛选得到差异表达的基因188 个，分别选取诱导后上调和下调的各前10 个进行分析，发现上调基因蛋白互作网络通过VTN 与下调基因网络的POU5F1 直接作用，将两个网络相关联。miRNA-mRNA 调控网络发现miR-335-5p 能够调控VTN、NANOG、POU5F1 和DNMT3B 4 个基因的表达，进一步验证发现miR-335-5p 表达量随着诱导分化进程逐渐升高，与NANOG、POU5F1 和DNMT3B 表达量逐渐降低形成对照，提示二者之间存在靶向关系且为负性调控。因此，该miRNAmRNA 调控网络显示miR-335-5p 靶向多个基因参与hESCs 向IPCs 诱导分化过程。