PM2.5染毒HBE细胞前后差异表达的microRNAs及其生物信息学分析

蔡 颖，李闰冰 ，郑 凯，秦双建 ，李柏茹 ，张朝晖，徐新云

(1．深圳市疾病预防控制中心，广东 深圳518055；2．南华大学公共卫生学院，湖南衡阳 421001；3．中南大学湘雅公共卫生学院，湖南 长沙 410078)

随着中国城市化进程快速发展，雾霾作为一种灾害性天气现象，已经成为影响大众健康的公共卫生问题。大气细颗粒物(fine particulate matter，PM2.5)是雾霾的主要成分之一，PM2.5凭借粒径小和表面积大等优势，能较长时间悬浮在空气中，并能随着呼吸道进入肺泡，诱发许多呼吸系统疾病[1]。多项研究发现，哮喘发生率的增加与大气中PM2.5的浓度上升有很大关系[2-3]。呼吸系统是PM2.5暴露和发挥作用的主要靶器官，我国工业化发展的同时也伴随着PM2.5浓度的不断上升，对人类健康构成了重大威胁，因此，探索PM2.5对呼吸系统损害的分子机制对于呼吸系统疾病的预防和治疗具有重要意义[4-6]。

近年来，研究发现非编码RNA在生物发育过程中起着重要的基因调控作用，microRNA(miRNA)一直是非编码RNA中的研究热点，可以在转录水平及转录后水平调控基因表达。许多研究表明，miRNA表达失调在PM2.5诱导呼吸系统疾病中起着重要作用，包括哮喘[7]，肺癌[8]等。miRNA 还被证明与免疫系统调节有关[9]。上皮间质转化(epithelial-mesenchymal transition，EMT)是肿瘤进展中的一个关键形态学变化，研究报道了A549 细胞暴露于PM2.5后诱发EMT，当敲除miR-32后显著增强Smad1介导的信号通路活性，增强了细胞的迁移和侵袭能力，促进EMT 进程[10]。由此可见，miRNA在PM2.5介导机体发生损伤的过程中起着重要的调控作用。因此，本文选取HBE 细胞作为实验对象，分析PM2.5暴露后是否引起HBE 细胞中的miRNA 表达改变，以及利用生物信息学方法筛选差异miRNAs，并进行功能富集分析，筛选关键miRNAs 及靶基因，从而在非编码RNA 研究领域为PM2.5致HBE 细胞损伤的分子机制研究提供新思路。

1 材料与方法

1.1 主要试剂

HBE细胞购于中国上海细胞库，PM2.5样品来源于山西省太原市，DMEM 培养基购于Hyclone 公司。胎牛血清(FBS)、0.25% Trpsin-EDTA 购于美国Gibco公司。

1.2 PM2.5采集与处理

使用TH-150CⅢ型中流量大气采样器(武汉市天虹仪器有限责任公司)采集PM2.5样品，采样地点选择煤炭消耗大省山西太原市，采样流速100 L/min，每次采样时间24 h。将吸附有PM2.5的滤膜剪成2 cm×2 cm放于烧杯中，加入超纯水完全浸没滤膜，用一层锡箔纸使烧杯口密封，超声振荡30 min，去除烧杯中残留的滤膜，制成PM2.5悬浊液，于-80 ℃条件下冷冻24 h。然后将样品置于真空冷冻干燥机中干燥24 h，至PM2.5完全干燥后用紫外灯灭菌1～2 h，用超纯水溶解PM2.5颗粒后制成PM2.5储备液，分装后于-20 ℃保存待用，使用前振荡混匀。

1.3 细胞培养与PM2.5染毒

将HBE细胞培养于含10%胎牛血清的高糖DMEM培养基中，37 ℃、CO2体积分数为5%的温箱中传代培养。待细胞生长至对数期时，以未处理组作为对照组，根据前期实验研究，PM2.5染毒HBE 细胞的IC50为70.12 mg/mL，因此采用50 mg/mL 的PM2.5染毒浓度作为实验组(以下称PM2.5组)，于不含胎牛血清及双抗的高糖DMEM 培养基中进行PM2.5染毒培养，每组一式3份准备miRNA测序。

1.4 构建cDNA文库及miRNA测序

高通量测序(illumina HiSeqTM2500/MiSeq 测序平台)检测RNA 样品中miRNA 的表达情况，由天津诺禾致源生物信息科技有限公司完成。大致流程为：①RNA样品检测。琼脂糖凝胶电泳分析RNA降解程度以及是否有污染；Nanodrop 检测RNA 的纯度[D(260)/D(280)比 值]； Qubit 对 RNA 浓 度 进 行 精 确 定 量 ；Agilent 2100 精确检测RNA 的完整性。②使用Small RNA Sample Pre Kit 构建文库。③使用Qubit2.0 进行初步定量后，采用qPCR 方法对文库的有效浓度进行准确定量(文库有效浓度＞2 nmol/L)，以保证文库质量。④把不同浓度文库按照不同比例混合后进行HiSeq/MiSeq测序。

1.5 miRNA差异表达分析

使用热图进行层次聚类分析以显示样品之间差异miRNA表达量的聚类模式。

用火山图可以推断差异表达miRNA的整体分布情况，从差异倍数(fold change)和校正后的P值两个水平进行评估，按照调整后P＜0.05 的标准筛选差异表达miRNAs。

1.6 miRNA功能预测

由miRanda和RNAhybrid两个软件预测miRNAs的靶基因，对每组差异表达miRNAs 的靶基因的集合进行富集分析。Gene Ontology(http://www.geneontology.org/)是基因功能国际标准分类体系，从生物学过程、细胞组分和分子功能3 个方面对靶基因进行功能注释，KEGG(Kyoto Encyclopedia of Genes and Genomes)是有关信号通路的主要公共数据库，能确定候选靶基因参与的最主要生化代谢途径和信号转导途径。

1.7 相互作用网络图

miRNAs 与靶基因及靶基因之间的相互作用关系构建网络图为：①采用Cytoscape 软件对miRNAs 及其相应的靶基因间的相互作用关系绘制成可视化网络图。筛选与靶基因相互作用数量较多的miRNAs 作为核心miRNAs。②利用STRING在线数据库预测靶基因间的相互作用关系，再由Cytoscape绘制成网络图。筛选其中相互作用数量超过10个节点的靶基因作为核心靶基因。

2 结 果

2.1 PM2.5暴露导致HBE细胞内miRNAs表达谱改变

PM2.5组和对照组HBE 细胞miRNA 高通量测序结果见图1。两组分别检测到1 021和927个miRNA，两组取交集得811 个miRNAs。根据调整后P＜0.05 筛选标准，有27 个miRNAs 发生差异表达，其中7 个miRNAs 表达上调，20 个miRNAs 表达下调。miR-372-3p 上调倍数最大，FC 为 1.94；miR-223-3p 下调程度最大，FC 为0.11。火山图和热图展示了27 个差异miRNAs的分布情况，见图2。差异倍数较大的前15个miRNAs见表1。

图1 PM2.5染毒组和对照组HBE细胞中miRNA测序结果

图2 miRNAs表达谱分布图

2.2 miRNAs功能富集分析

应用GO和KEGG 富集分析来预测差异miRNAs 的功能，结果见图3 和图4。GO 分析结果显示，差异miRNAs 主要参与细胞代谢、蛋白质代谢等生物学过程，如单一生物代谢过程、代谢过程等。细胞组分主要集中在胞内如细胞质、细胞器等。对于分子功能，差异miRNAs主要体现在蛋白质结合、催化活性等。

表1 前15个差异表达miRNAs

为了分析PM2.5致HBE细胞损害的潜在机制，采用KEGG 数据库对靶基因参与调控的信号通路及途径进行分析。结果显示，代谢途径、PI3K-Akt 信号通路、Ras 信号通路、MAPK 信号通路等途径可能是PM2.5发挥毒性作用的潜在途径。

2.3 相互作用关系分析

图3 候选靶基因GO富集柱状图

图4 候选靶基因KEGG富集散点图

图5 miRNA-靶基因网络图

2.3.1 miRNAs-靶基因相互作用对miRNAs 进行功能预测得到miRNA参与调控的生物学过程及相关信号通路，再利用可视化软件展现miRNAs 及其靶基因间的对应关系，见图5。一个miRNA 可以调控多个靶基因，一个靶基因也可以由多个miRNA调控。将与靶基因相互作用数量较多的miRNA(上调5个，上调6个)作为核心miRNAs，绘制核心miRNAs-靶基因相互作用网络图。结果显示，上调核心miRNAs为：miR-371b-3p、miR-371a-5p、miR-27a-3p、miR-7-5p、miR-372-3p，其中miR-371b-3p调控靶基因最多。下调核心 miRNAs 为 miR-145-5p、miR-133a-3p、miR-204-5p、miR-125b-5p、miR-486-5p、miR-512-3p，其中miR-145-5p调控靶基因最多。

2.3.2 靶基因相互作用将上调和下调差异miRNA的靶基因分别输入STRING 在线数据库，利用Cytoscape 可视化软件对所得TSV 文件数据绘制网络图，见图6。并筛选由下调miRNAs 调控的3个相互作用节点大于10 的核心基因，为血管内皮生长因子(vascular endothelial growth factor A，VEGFA)、丝裂原激活的蛋白激酶3(mitogen-activated protein kinase 3，MAPK3)和 CC 趋化因子受体 5 型(C-C chemokine receptor type 5，CCR5)。

3 讨 论

图6 靶基因相互作用网络图

在过去的几十年里，miRNAs 已经越来越为人们所广泛研究[11-12]。miRNA 是一种高度保守的小型非编码RNA。通过与靶mRNA 的3′非编码区碱基互补配对，在转录后水平调控基因表达，发挥其生物学功能。miRNAs 能通过抑制翻译或者降解靶基因来维持细胞代谢、细胞凋亡、蛋白质合成等正常生理过程[13]，在异常表达情况下，miRNA 既能充当致癌基因也能作为抑癌基因[14-15]。有许多研究发现，PM2.5暴露能导致miRNAs 表达谱的改变[16]，为呼吸系统疾病的分子机制研究提供了新的方向。

利用Illumina 测序平台，快速、高效地读取高质量的测序数据。我们在PM2.5组和对照组分别检测到1 021 和 927 个 miRNA。按照调整后P＜0.05 筛选条件，我们共筛选到27 个差异miRNAs，包括7 个上调miRNAs和20个下调miRNAs。对靶基因进行富集分析显示，差异miRNAs 主要参与代谢途径、PI3K-Akt 信号通路、Ras 信号通路、MAPK 信号通路等关键信号转导途径。PI3K-AKT 信号通路是细胞内重要信号转导通路，广泛参与了细胞增殖、代谢及死亡等生命过程，在肿瘤发展进程中起着重要的调控作用[17]。胶质瘤细胞中MRPS16(线粒体核糖体蛋白S16)过表达通过激活PI3K-AKT 信号通路促进Snail 的表达，从而促进肿瘤的进展[18]。MAPK 信号途径是涉及细胞增殖、分化及凋亡等生物学过程的关键途径，已有许多证据表明了MAPK 信号途径失控与多种恶性肿瘤包括肺癌的发展、诊断及预后的相关性[19]。Soleimani 等[20]发现Ras/MAPK 信号通路中的miRNAs 在大肠癌细胞增殖、凋亡等过程中起着重要作用。

我们分别对上下调表达差异miRNAs 及其预测的靶基因构建相互作用网络图，结果显示，共鉴定11个核心miRNAs，含上调5 个：miR-371b-3p、miR-371a-5p、miR-27a-3p、miR-7-5p、miR-372-3p，其中miR-371b-3p 调控靶基因最多。下调6 个：miR-145-5p、 miR-133a-3p、 miR-204-5p、 miR-125b-5p、miR-486-5p、miR-512-3p，其中miR-145-5p 调控靶基因最多。Yue等[21]研究发现，miR-371a-5p过表达后能逆转肿瘤抑制因子lncRNA TUSC7 对胰腺癌细胞迁移、侵袭和EMT的抑制作用。有研究表明，miR-27a-3p 既具备肿瘤抑制作用[22]也可促进肿瘤发生发展[23]。miR-145-5p 是一种肿瘤抑制因子miRNA，在肺癌[24]、胶质母细胞瘤[25]等癌症中表达下调。miR-7-5p[26]、372-3p[27]、miR-133a-3p[28]、miR-204-5p[29]等也被证明在癌症的发生发展中起着关键作用。关于这些核心miRNAs 在PM2.5致HBE 细胞损伤中的分子机制尚不完全清楚，有待后续研究。

进一步对靶基因之间相互作用分析发现，VEGFA、MAPK3、CCR5位于核心位置。VEGFA在血管生成、细胞增殖及迁移等途径发挥重要作用，该基因被证明在许多癌症中表达上调[30]。MAPK3是哺乳动物中第一个被发现的丝裂原活化蛋白激酶，是ERK/MAPK 信号途径中重要组成部分。Cao 等[31]证明miR-129靶向MAPK3调控胃癌细胞增殖，当MAPK3基因表达抑制时能诱导胃癌细胞凋亡，降低顺铂耐药性。趋化因子受体是一类介导趋化因子行使功能的GTP-蛋白偶联的跨膜受体，已有许多研究表明CCR5 在激活癌症的生长及转移中具有重要作用[32]。本实验中这3个蛋白与下调的miRNAs 靶向相关，提示VEGFA、MAPK3、CCR5 在PM2.5染毒HBE 后有可能受到调控而发生表达改变，可以为下一步PM2.5的致病机制研究提供参考。

综上所述，本文采用miRNA测序技术测定了PM2.5染毒HBE 细胞后的差异miRNAs，主要富集在PI3KAkt 信号通路、肿瘤相关途径、Ras 信号通路、MAPK信号通路等，筛选了11 个核心miRNAs 和3 个核心基因，并且这些miRNAs 和蛋白与肿瘤的发生有着密切的关系，为进一步深入研究PM2.5毒性作用机制拓展了新思路。