衰老过程中骨骼肌microRNA-mRNA调控网络分析

窦媛媛 侯静雯 王枚

肌少症是衰老的常见特征，表现为增龄相关的进行性肌肉强度和质量的逐渐减弱。肌肉强度下降或肌肉生理功能减退伴随着老年人体育活动和能量消耗的减少[1]。近10年来肌少症的研究逐渐引起重视，并在基础和临床方面取得一些进展。尽管使用转录组和蛋白质组学方法表征了导致肌减少症的多种因素[2]，但对骨骼肌衰老相关变化的分子机制了解甚少。微小RNA(miRNA)是长度为20～22个核苷酸的内源性非编码的RNA，参与代谢、发育、癌症和衰老等生物过程[3]。最近，有研究证明了一些miRNA在衰老过程中差异表达，并以组织和细胞类型的特定方式参与衰老过程[4-7]。这些研究都证明了miRNA在肌肉衰老过程中的重要作用[8]。本研究对年轻和老年小鼠腓肠肌miRNA和mRNA数据库进行联合分析，通过下一代测序技术(next-generation sequencing，NGS)研究miRNA在骨骼肌衰老中的潜在作用，以精确鉴定成熟和新颖的miRNA，以期为解析肌少症的发生和寻找肌少症新的分子标记物及靶点提供理论依据。

1 材料和方法

1.1 基因表达谱数据下载 数据来源于美国国立生物技术信息中心 (national center for biotechnology information，NCBI)的基因表达综合数据库(gene expression omnibus, GEO)，包含5只青年小鼠和5只老年小鼠骨骼肌样本的RNAseq数据(GSE55163)及6只青年小鼠和6只老年小鼠骨骼肌样本的miRNAseq数据(GSE55164)。所有样本采用NextSeq 500 芯片进行表达基因的测序用于后续分析。

1.2 样本的聚类和差异基因分析 采用基于阵列间归一化功能对数据进行背景校正和标准化，筛选出对照探针和低表达探针，获得高质量的标准化数据。采用R语言中的limma 软件包的lmFit和eBayes函数分析差异表达基因。差异表达基因的筛选条件为P<0.05且log2 FC>1.5。

1.3 功能和途径富集分析 使用R语言Clusterprofiler包对差异表达基因进行基因本体(gene ontology，GO)和功能富集分析。利用Cytoscape的CluGO应用对差异表达基因进行京东基因与基因组百科全书(Kyoto encyclopedia of genes and genomes，KEGG)通路富集分析。此外，通过表达分析系统(GSEA)软件分析差异表达基因参与的GO功能和KEGG信号通路。P<0.05被定义为重要功能和通路分析的临界值。

1.4 miRNA-mRNA相互作用网络的构建 使用mirTarbase和TargetScan(6.2版，http：//www.targetscan.org)数据库来预测每个miRNA的潜在靶标mRNA(互作关系评分>0.5)。基于预测的miRNA-mRNA关系，选择显示表达相反的miRNA-mRNA,采用Cytoscape(3.0版，www.cytoscape.org/)构建miRNA和mRNA之间的相互作用网络。

1.5 样品收集 2020年4月从新疆医科大学实验室动物资源中心购买12只青年(6个月大)和13只老年(24个月大)雄性C57BL/6小鼠，每组小鼠体质量接近且健康。在处死之前，记录小鼠的体质量，分离腓肠肌。

1.6 实时荧光定量PCR(real-time fluorescence quantitative PCR，qRT-PCR) 使用mirVanaTMmiRNA分离试剂盒试剂(Thermo，USA)分离小鼠腓肠肌总RNA，将总RNA(1μg)反转录成cDNA。然后使用罗氏480II检测系统(Applied Biosystems, USA)进行实时定量PCR。

1.7 Western blot 检测蛋白表达 使用总蛋白提取试剂盒(碧云天，江苏)从腓肠肌中提取总蛋白。在10%十二烷基硫酸钠聚丙烯酰胺凝胶电泳凝胶(SDS-PAGE)上分离蛋白质(50μg)，并转移到聚偏氟乙烯膜(Millipore，USA)上。室温下用Tris缓冲盐水配置的5%奶粉将膜封闭1 h。用PBS洗涤3次后，将膜与一抗(Rybp抗体, 1∶600; GAPDH抗体，1∶2000; Abcam, USA)在37 ℃孵育2 h，然后加辣根过氧化物酶标记的二抗稀释液(1∶4000；sc-358914;Abcam, USA)，37 ℃摇床下孵育1 h。使用ECL检测试剂盒(Pierce，USA)检测信号。

1.8 数据统计与分析 使用SPSS 22.0统计软件，所有数据均以均数±标准差表示。采用两尾Student′st检验来分析2组之间的统计学差异。P<0.05为差异有统计学意义。

2 结果

2.1 测序数据处理结果 共鉴定到差异表达的基因877个，其中老年小鼠骨骼肌中上调基因747个，下调基因130个(图1 A和B)。共鉴定到差异表达的miRNA 84个，其中老年小鼠骨骼肌中上调miRNA 52个，下调miRNA 32个(图1 C和D)。

图1 差异miRNA和mRNA筛选

2.2 差异miRNA靶基因预测及miRNA-mRNA网络构建 共筛选到miRNA靶基因1765个，其中有112个基因与差异基因重复。构建了这112个基因和miRNA的miRNA-mRNA网络，结果包含21个miRNA。其中靶向基因超过7个的miRNA有9个，分别为：miR-122-5p，miR-127，miR-137-3p，miR-200a-3p，miR-375-3p，miR-449a-5p，miR-5101，miR-758-3p，miR-873a-5p。

2.3 代表性差异miRNA的验证 通过qRT-PCR分析了来自青年(6个月大)和老年(24个月大)小鼠腓肠肌的9种代表性miRNA表达变化。来自测序数据中代表性的4个下调的miRNA的qRT-PCR结果与测序结果一致(图2A)。来自测序数据的所有5个上调的miRNA在qRT-PCR中均上调，其中miR-375-3p在qRT-PCR结果中上调倍数高于测序结果(图2B)。

图2 代表性差异miRNA的验证

2.4 差异靶基因的GO分析 GO功能注释分类统计分析显示，差异靶基因显著富集的生物过程除了基因转录、表达和细胞迁移调节外，还有细胞周期和骨骼肌细胞分化过程。其中骨骼肌分化过程中包含了早期生长反应蛋白2(EGR2)、环1(Ring1)人转录因子 Yin-Yang1(YY1)和 Ring1 和 YY1 结合蛋白(Rybp)基因。而在细胞组成分中，差异表达基因与细胞外泌体、核及髓鞘形成中相关性较高。分子功能中，蛋白质及其复合物结合和RNA聚合酶Ⅱ转录共抑制活性过程中占比较高(图3)。

图3 差异靶基因的GO分析

2.5 EGR2、Rybp、Ring1和YY1基因的验证 通过qRT-PCR证实了骨骼肌分化过程中，老年组所有4个基因的mRNA水平与青年组相比均显著升高(P<0.05, 图4A)。其中，Rybp，Ring1和YY1结合蛋白在损伤性肌肉再生中负责调节骨骼肌的新生[9]。在老年小鼠的肌肉中，Rybp在蛋白质水平上也有所增加(图4B)。这些结果表明，Rybp及其推测的调控性miRNA miR-136可能促进了肌肉衰老过程。

图4 EGR2, Rybp, Ring 1和YY基因的验证

3 讨论

本研究通过系统地分析青年和老年小鼠骨骼肌的miRNA和mRNA基因表达谱，总共鉴定出84种差异表达miRNA，其靶基因与mRNA联合分析后发现21个miRNA与年龄相关。在下调的miRNA中，miR-127和miR-136a-5p属于骨骼肌中最丰富的miRNA。MiR-127被鉴定为富含心脏瓣膜的miRNA[10]；但是，以前没有关于其在骨骼肌发育中作用的报道。MiR-136a-5p是下调的miRNA之一，被鉴定为通过ROCK1基因促进成肌分化的成肌miRNA。MiR-137-3p是上调倍数最多的差异miRNA，它通过抑制MEF2A来抑制骨骼肌的分化[11]。因此，根据其表达和在肌发生中的作用，miR-127和miR-137-3p可能通过肌肉分化或再生参与了衰老过程。已有的研究报告了随着年龄的增长，骨骼肌中miRNA表达的全基因组表达谱[12-13]。最近，报道了基于恒河猴不同年龄骨骼肌miRNA表达谱的改变[14]。但在先前的研究中，还没有关于小鼠骨骼肌中差异表达miRNA的研究。因此，我们的数据对以前的研究进行了有效的补充。

大多数哺乳动物的miRNA基因散布在不同的基因组位置，但是一些miRNA基因聚集在相同的基因组中[15]。有些miRNA簇包含2～3个miRNA基因，而有些是包含超过50个miRNA基因的巨型簇。簇状的miRNA通常朝向相同的方向，并以类似表达模式的多顺反子方式转录。位于人14号染色体(小鼠12号染色体)上的印迹Dlk1-Dio3基因组区域包含父本表达的基因Dlk1、Rtl1和Dio3以及母本表达的非编码RNA基因Meg3(Gtl2)和Meg8(Rian)以及反义Rtl1(anti-Rtl1)[16]。此外，Dlk1-Dio3区包含54个miRNA，这是人类基因组中发现的最大miRNA簇[17]。本研究中鉴定出的4个下调的miRNA中有3个(75%)属于源自小鼠远端12号染色体上Dlk1- Dio3区的miRNA簇，并且在肌肉中的方向相同。

Dlk1-Dio3区域内的miRNA簇和基因在各种疾病(如癌症和精神分裂症)以及骨骼肌发育中受到差异调节[18]。绵羊Dlk1、Gtl2、Rtl1和Meg8的表达在骨骼肌中丰富，而Dlk1和Rtl1的过表达与Callipyge绵羊的肌肉肥大密切相关[19]。Gtl2-Dio3 miRNA簇介导的MEF2A对WNT信号的调节是骨骼肌再生所必需的[20]。通过对小RNA-Seq和mRNA-Seq数据的综合分析，我们使用TargetScan算法鉴定了由9个年龄调控的miRNA靶向的94个mRNA(去除重复)。GO分析显示，下调miRNA靶向的基因(例如Ccna2、Egr2、Klf4、Klf6、Mafb、Nfya和Rybp基因)主要参与转录调控。其中，有报道指出，在成肌过程中，miR-136下调了Ring1和Rybp，并在损伤诱导的肌肉再生中作为骨骼肌发生的负调节剂[9]。这些结果表明，Rybp及其推测的调控性miRNA miR-136可能有助于肌肉衰老。

综上，我们在骨骼肌中鉴定出9种与年龄相关的miRNA。通过分析miRNA与mRNA的相互作用，我们发现miRNA主要通过转录调节来促进肌肉衰老。本研究为miRNA在骨骼肌衰老中的作用机制研究奠定了基础。