hsa-miR-130家族生物信息分析及在SLE中差异表达分析①

栾鹏飞 周少岚 党 洁 贾 维 马占兵 霍正浩 （宁夏医科大学基础医学院医学遗传学系与细胞生物学系，宁夏医科大学生育力保持教育部重点实验室，银川750004）

系统性红斑狼疮（systemic lupus erythematosus，SLE）是一种病因未明，以血清中产生自身抗体为主要特征、可累及多系统多器官的慢性自身免疫性疾病，好发于育龄妇女［1］。其病因学包括环境和遗传易感等多种因素［2］。越来越多的证据表明，微小RNA（miRNA）在SLE 的发病机制中起关键作用。miRNA 是一类长度在18~36 nt，短小而高度进化保守的非编码小RNA 分子，miRNA 的种子序列（seed）通过直接与靶标mRNA 3'末端非翻译区（3'UTR）的miRNA 响应原件（miRNA response elements，MREs）完全互补或不完全互补结合，发挥转录后基因沉默或表达抑制效应，抑制转录后水平基因表达［3］。miRNA 在众多生物学过程中均扮演十分重要的角色，包括调节细胞生长、凋亡，机体的新陈代谢及信号转导等。在众多miRNA 当中，hsa-miR-130 家族是一簇被广泛研究的miRNA 分子。有研究表明，hsa-miR-130家族与免疫疾病密切相关，该家族成员包括 miR-130a/b、miR-301 和 miR-454［4-5］。据估计，哺乳动物基因组都包含2 000个以上的miRNA［6］，这些miRNAs 调节每个基因组中约60%的基因，甚至有多个miRNA 调控某一靶标基因的表达，作用呈现网络化和系统化［7-9］。有效预测miRNA 靶基因并结合系统生物学手段对其功能进行分析和研究是研究其作用机制的有效途径，而基于生物信息学，使用多种算法工具联合可有效降低靶标预测假阳性率［10］。目前，已有研究表明，hsa-miR-130 家族分子与细胞的增殖、免疫代谢过程等相关，但其在SLE发病中的具体分子作用机制和临床表达特征尚未见报道。因此，本研究拟通过miRNA 家族保性分析，实时荧光定量PCR 技术明确hsa-miR-130a-3p 在SLE 中的表达特性，通过生物信息学预测hsa-miR-130 家族靶基因，并对其靶基因集合进行基因本体（gene ontology，GO）富集分析、信号通路（KEGG）富集分析及靶标蛋白质相互作用网络分析（PPI）等研究，以期为深入研究hsa-miR-130家族在SLE中作用机制提供理论和实验依据。

1 资料与方法

1.1 资料 本研究经宁夏医科大学伦理委员会批准（批准号：2019001），所有参与者均提供书面知情同意书。选取2018 年1 月至2018 年11 月宁夏医科大学附属风湿科医院收治的18 例SLE 患者作为病例组（男、女各9 例），样本纳入标准严格遵守1982年和1997 年修订的美国风湿病学会SLE 分类标准（Hochberg，1997）。另按性别和年龄与SLE 患者组相匹配，选取同期18 例没有类风湿性关节炎、肾脏疾病、结缔组织疾病或原发性干燥综合征病史的健康个体作为对照组。

1.2 方法

1.2.1 荧光定量表达分析 RT-qPCR 法检测临床SLE 患者及对照组外周血单个核细胞（PBMCs）中hsa-miR-130a-3p 的表达水平，所用引物见表1。采用天根总RNA 柱法提取试剂盒（DP118）提取18 例SLE 患者及健康志愿者PBMCs 中的总RNA。使用茎环引物逆转录获取cDNA，逆转录条件：37℃15 min，85℃5 s，4℃保存。目标基因的表达水平检测使用SYBR-Green 法，以cDNA 为模板，扩增条件：95℃ 预变性30 s；95℃ 变性5 s；62℃ 退火30 s，72℃延伸30 s，重复40 个循环。数据处理以GAPDH 为内参基因，2-ΔΔCt方法计算目标基因的相对表达水平。

表1 实时荧光定量PCR引物序列Tab.1 RT-qPCR primer sequences

1.2.2 进化保守性分析使用 NCBI（https：//www.ncbi.nlm.nih.gov/）、UCSC（http：//genome.ucsc. edu/）、miRbase（http：//www. mirbase. org/）等在线数据库查找hsa-miR-130 家族miRNA 前体的碱基序列及染色体定位，经DNAMAN 软件多序列连配分析，比较其在物种间的进化保守性。

1.2.3 靶基因预测和富集分析 分别使用不同算法的靶基因预测软件miRanda（http：//www. microrna.org/）、targetscan7.2（http：//www.targetscan.org/vert_72/）和 PicTar（https：//pictar. mdc-berlin. de/）联合预测的基因作为进一步分析的靶基因数据集。应用软件Cytoscape3.7.2 的CluoGO 插件进行靶基因GO 功能富集分析，分别投射细胞组分（cellular component，CC）、分子功能（molecular function，MF）和生物学过程（biological process，BP）三大基因本体术语。使用R 包ClusterProfilter 进行KEGG 通路富集分析。

1.2.4 靶基因互作网络（PPI）分析 将联合预测靶基因名称列表提交至 STRING（https：//string-db.org/）蛋白互作网络分析数据库，选择物种为Human species 执行蛋白质相互作用网络分析。在Cytoscape 中将STRING 输出的网络图源文件打开，使用MCODE插件从网络图中找到关键基因。

1.3 统计学分析 采用SPSS25.0 软件分析数据，所有计量资料以表示，采用单因素方差分析（ANOVA）进行统计学分析，组间比较采用t检验，RT-qPCR 实验结果采用 Mann-WhitneyU检验，ClusterProfilter 富集性分析采用Fisher 确切概率法检验，P＜0.05为差异有统计学意义。

2 结果

2.1 SLE 生化指标 临床SLE 患者年龄、性别、系统性狼疮疾病活动度指数（SLEDAI）评分、血沉（ESR）、抗ds-DNA 抗体水平（Anti-dsDNA）、抗核抗体（ANA）、超敏反应蛋白（CRP）等临床指标如表2所示。统计分析表明，男、女性患者免疫学检查结果ANA 均为阳性、抗dsDNA 均大于100.0，且二者同时存在，特异性高；不同性别间SLEDAI 评分差异有统计学意义（P＜0.01）。

表2 SLE患者的临床特征（）Tab.2 Clinical characteristics of SLE patients（）

表2 SLE患者的临床特征（）Tab.2 Clinical characteristics of SLE patients（）

Note：ESR reference range：adult males 0~15 mm/ h，adult females 0~20 mm/ h；anti-dsDNA reference range：0.00~7.00；ANA reference range：negative（-）；CRP reference range：0.00~2.87.

CRP（mg/ L）12.46±4.28 4.55±16.43 35.00 0.627 SLE Samples Female Male t/ T P Age（year）47.67±4.78 37.33±4.26 1.613 0.126 SLEDAI 12.89±0.78 14.33±1.00 3.414 0.004 ESR（mm/ h）29.33±5.47 26.44±26.10 28.50 0.289 Anti-dsDNA（U/ ml）134.99±96.12 121.73±52.08 40.00 0.965 ANA++--

2.2 hsa-miR-130a-3p 在SLE 患者组与对照组中的表达 hsa-miR-130a-3p 的RT-qPCR 曲线和溶解曲线如图1A、B，结果表明引物特异性较高，内参基因和目的基因扩增效率均大于85%，Ct 值范围在18~30个循环以内，能够用于miRNA差异表达的有效检测。与健康对照组相比，SLE 患者外周血PBMCs 中hsa-miR-130a-3p 表达降低，差异有统计学意义（图1C，P＜0.05）。

图1 hsa-miR-130a-3p RT-qPCR检测结果Fig.1 Detection results of hsa-miR-130-3p RT-qPCR

2.3 hsa-miR-130 家族进化保守性分析 miRbase分别检索人（hsa）、黑猩猩（ptr）、大鼠（rno）、小鼠（mmu）等其他物种miRNA 前体序列，多序列比对结果显示，hsa-miR-130 种子区序列AGTGCAAT 在脊椎动物中高度保守，见图2。

图2 不同物种中hsa-miR-130前体序列中的种子序列Fig.2 Seed sequences of precursor of hsa-miR-130 in different species

2.4 hsa-miR-130的靶基因预测 采用不同算法软件交叉预测miRNA 靶标，可明显提高预测作用靶标的阳性率。检索 PicTar、Targetscan 及 miRanda 数据库，参数设置默认，结果分别得到483、1 028和4 470个靶基因，取三者预测靶基因的交集后共得到285 个靶基因集合用于后续功能富集分析（图3）。

图3 hsa-miR-130家族靶基因预测情况Fig.3 Prediction of hsa-miR-130 family target genes

2.5 hsa-miR-130 靶基因功能注释及信号通路分析 通过Cytoscape CluGO 插件对285 个靶基因进行GO 富集分析，结果表明，靶基因涉及的细胞成分（CC）主要包括Ccr4-Not 复合物、膜细胞器及cul3 环泛素连接酶复合物等；在生物学过程（BP）上则与细胞通讯、信号转导及细胞因子和趋化因子介导的信号通路等密切相关；分子功能（MF）包括多种蛋白酶的活性、转录因子的活动、离子通道的活动等；采用R 包ClusterProfilter 对靶基因集合进行生物通路富集分析，其中45个靶基因具有相关生物通路。以人的所有基因为背景基因，发现hsa-miR-130的靶基因富集于细胞VEGF/VEGFR信号通路、雌激素膜受体信号转导以及与细胞周期调控密切相关的其他信号通路中（P＜0.05）。见图4。

图4 hsa-miR-130靶基因的GO功能分析（Top10）Fig.4 GO function analysis of hsa-miR-130 target gene（Top10）

2.6 靶基因转录因子结合位点（TFBS）分析 对hsa-miR-130 预测的共同靶标mRNA 启动子区转录因子进行富集分析，结果如图5 所示，有众多与SLE发生发展密切关联的转录因子。基因表达的调控是疾病发生发展的核心，而基因表达的特异性调节主要通过转录因子（transcription factors，TF）实现，miRNA 作为重要的转录后基因表达调控因子，在各种生物学过程中发挥关键作用，而hsa-miR-130如何影响众多转录因子发挥作用，有待更深层次的探究。

图5 hsa-miR-130共同靶标mRNA 启动子区转录因子富集分析Fig.5 Enrichment analysis of transcription factors in mRNA promoter region of common target of hsamiR-130

2.7 靶基因互作网络 联合预测得到hsa-miR-130靶基因共285个，根据STRING蛋白质相互作用数据信息可构建PPI网络，进一步从该PPI网络中进行聚类从而构建功能模块，使用Cytoscape 插件MCODE从网路图中找到排名前20 的关键基因，结果如图6所示。通过互作网络分析及聚类分析以发现整个网络中的关键基因，如ESR1、SKP1、KIT 等与SLE 发病密切相关。其中女性性类固醇——雌二醇通过其受体（ESR）发挥作用，致使SLE 出现性别二态性的机制至今尚不明确，而ESR1 作为互作网络中的关键基因受hsa-miR-130的调控，能否导致疾病的进展有待深入研究；另外KIT 基因位于B 细胞受体信号通路（BCR）上游，而作为hsa-miR-130靶基因之一的KIT 或将受其调控进而影响下游STAT 转录因子的表达，或有可能激活炎症通路，导致SLE。

图6 PPI网络关键基因（Top20）Fig.6 PPI network key genes（Top20）

3 讨论

自 1993 年 LEE 等［11］在研究线虫发育过程中发现首个miRNA 基因以来，越来越多的miRNA 被陆续发现，其功能及在基因表达调控中发挥的作用也受到越来越多的关注。目前为止，共有超过48 000个miRNAs 被识别，包括2 693 个人类起源的mi-RNAs［12-13］。由于单个 miRNA 可靶向多达 400 种不同的mRNA，因此预测超过一半的人类基因直接受miRNA调控［8］。有研究表明，miRNAs表达异常是许多病理过程的特征，包括癌症、代谢紊乱、炎症、心血管、神经发育及自身免疫性疾病［14］。

研究表明，SLE 患者中存在大量异常表达的miRNAs，参与了SLE 免疫异常和器官损伤，而应用miRNA 治疗 SLE 具有重要临床意义［1］。既往文献报道hsa-miR-130 在急性心肌梗死、非小细胞肺癌、膜性肾病及其他肿瘤中下调发挥作用（表3），且已有研究报道hsa-miR-130a 参与众多重要的生物学过程［15-18］。对 hsa-miR-130 预测的共同靶标 mRNA 启动子区转录因子进行富集分析，结果显示出众多受其影响的转录因子，其中参与了细胞生长和分化等重要生命过程的转录因子Sp1，是目前研究极为广泛的系统转录因子，该分子可影响细胞周期调控分子的表达和生长相关信号转导通路及细胞的凋亡过程，从而对机体免疫系统产生作用［19］；FOS则在研究中被证明与其他基因存在通过不同途径间的相互作用，以及作为关键转录因子与SLE 发病过程密切 相 关［20-21］；EGR1 曾 被 报 道 与 IGF-2、APEX、SRD5A1、CD44 和 EGFR 等基因相互作用，并通过转录调控激活这些基因在B 细胞中的表达，加速SLE病情的发展，也被确定与SLE的发病密切相关［22］。

表3 hsa-miR-130调控靶基因参与人类的部分疾病Tab.3 hsa-miR-130 regulatory target genes are involved in some human diseases

此外，有文献指出，hsa-miR-130a可与miR-326、miR-155、miR-210 共同作为多形性成胶质细胞瘤患者生存的预后和预测标志物［4］；也有报道表明hsamiR-130a-3p可能参与狼疮性肾炎发病，但其具体分子机制尚不明确［23］。此外有研究表明，在SLE 患者中，hsa-miR-130a-3p 可通过上调Klotho（抗衰老基因）保护脂多糖诱导的肾小球细胞损伤［23］；与hsamiR-130a-3p 同家族的另一分子hsa-miR-130b-3p 曾在WANG 等［24］的研究中被发现与蛋白尿和狼疮性肾炎的慢性肾损伤指数密切相关，虽然没有直接证据证明循环中的hsa-miR-130b-3p 是否会被肾小管上皮细胞吸收并影响其功能，但有研究指出，从间充质干细胞释放的微泡可通过改变mRNA 和mi-RNA 水平从而保护由缺血再灌注损伤所引起的急性肾损伤和随后出现的慢性肾损伤，因此增加体内循环hsa-miR-130b-3p 的水平能否有效改善病情值得进一步研究［25］。

本研究通过生物信息学方法预测hsa-miR-130的靶基因，并对其靶基因进行GO注释及KEGG通路分析。预测得到的285 个hsa-miR-130 的靶基因存在于细胞的各个组分中，具有细胞周期调控、DNA表达调控和自噬调节等分子功能，并富集于细胞VEGF/VEGFR 信号通路等其他信号通路，而众多文献表明SLE 患者皮损中VEGF 蛋白表达水平增高，提示VEGF 蛋白在SLE 的发生、发展过程中可能发挥重要作用，但SLE 患者低表达的hsa-miR-130a 与VEGF 信号通路中的关键基因Flt-1 和KDR 是否存在直接靶向关系有待进一步验证［26-27］；此外，KEGG通路分析显示雌激素膜受体信号转导通路在细胞增殖、分化和免疫系统中具有重要作用，而这也在相关研究中得到证实［28］。此外，据报道，异常的mi-RNA 表达模式是由雌激素引起的，且越来越多的miRNA 被鉴定出受雌激素调控，其中包括miR-146a和 miR-155，都与 SLE 密切相关［28-30］。有研究表明SLE 患者血浆中雌激素水平与 IL-6、IL-21、IL-23 密切正相关，女性SLE 患者雌激素水平较健康女性显著升高。异常雌激素可能通过作用于miR-21、miR-25及miR-186等分子参与干扰素途径在SLE 中发挥重要作用［31］，而在 TargetScan Human 7.2 中可查到ESR1基因与hsa-miR-130家族种子区域的配对类型均为7mer-m8，并且二者与hsa-miR-130 家族非种子区域3'端存在不同长度的互补位型（表4）。hsamiR-130 家族种子区域与ESR1 基因的3'UTR 完全互补配对可增强hsa-miR-130对ESR1基因的沉默效果，故可利用hsa-miR-130 与ESR1 基因的高特异性来研究hsa-miR-130对ESR1基因富集的雌激素代谢通路（如SLE 的激素治疗）。此外，本文对hsa-miR-130 家族靶标PPI 网络度量分析后，雌激素受体（ESR1）作为PPI 核心关键基因，生物信息学分析显示其可能受hsa-miR-130调控从而参与SLE发病，但其通过哪种机制发挥作用有待深入验证。

表4 ESR1基因3'-UTR与hsa-miR-130家族结合位点Tab.4 Binding sites of ESR1 gene 3'-UTR with the hsa-miR-130 family

目前有关hsa-miR-130 在SLE 发病中的分子机制尚未揭示。本文检测了SLE患者及健康人群外周血 PBMCs 中 hsa-miR-130a-3p 的表达，发现 hsa-miR-130a-3p在SLE患者中明显下调，结合生物信息学分析结果，提示该分子可能参与SLE的发病，为进一步研究miRNA参与SLE的发病机制研究奠定基础。