中老年猕猴单个核细胞miRNA表达差异及其靶基因预测

朱向情，王 强，蔡学敏，李自安，庞荣清，潘兴华

衰老是机体细胞内基因结构、表达谱及表达水平、基因表达调控机制发生改变，导致生物活性分子、组织细胞及其组成成分、组织结构以及器官功能逐渐退化的生理过程。关于衰老发生学说，有体细胞突变、自由基、差错灾难、脂褐素累积、内分泌功能减退等[1]，这些学说从不同角度解释了衰老发生的原因，但均存在片面性。最近研究发现，机体免疫功能降低或紊乱、慢性炎症是机体衰老的重要原因[2]，但调控机制尚不清楚。外周血单个核细胞是机体免疫应答和防御的主要细胞组分，miRNA是一种含有约22个核苷酸的非编码小分子RNA，不具编码功能，但能够与靶mRNA结合并调控基因表达和细胞生物进程[3]。miRNA对衰老相关基因表达也具有调控作用，已发现了mir-499、mir-34c等衰老调节miRNA，但数据十分有限[4]。本研究采用二代测序生物信息分析技术，研究了中年与老年猕猴外周血中单个核细胞的miRNA的表达差异及其靶基因，为揭示衰老发生机制提供新的科学数据。

1 材料与方法

1.1 实验动物 选取8周岁和20～21周岁雌性健康猕猴各3只，分别设为中年组、老年组。每只采集外周血3.0 ml，用于单个核细胞分离和miRNA提取。所有猕猴均由中国科学院昆明动物学研究所动物中心提供，均饲养于医院实验动物中心，动物使用经过医院实验动物伦理委员会的批准，所有动物在普通级环境下、按照常规条件饲养。

1.2 miRNA测序及数据预处理 miRNA测序委托上海欧易生物医学科技有限公司完成，使用Illumina公司Hiseq 4000测序仪测序。对测序产生的原始Fastq格式文件数据进行引物与载体序列去除、测序片段碱基的质量检验和长度筛选等，筛选高质量的测序片段。数据预处理通过perl脚本完成。

1.3 miRNA的识别 经过数据预处理后，根据Read长度统计片段长度分布情况，长度在21～25 nt的Read视为预测的miRNA。采用Bowtie软件[5]将miRNA序列与miRBase数据库中的miRNA进行比对，比对上的序列为已知miRNA。将未比对上的序列与Rfam数据库中的序列进行比对，去除可能存在的其他类型的小RNA，将过滤后的序列与猴转录组序列进行比对，去除由降解mRNA测序产生的序列。将未同转录组序比对上的序列与RepBase数据库进行比对，去除重复序列。剩余序列采用miRDeep2软件[6]结合猴同源miRNA序列以及RNAfold软件预测的二级结构，识别新的miRNA。

1.4 miRNA差异分析 比对已知的miRNA序列，进行表达量统计以及差异分析，miRNA表达量计算采用TPM计算度量指标，TPM是指以每百万比配成对序列做miRNA表达量指标，利用DESeq2软件包[7]中的负二项分布检验计算基因差异表达量。采用NB（负二项分布检验的方式）对reads数进行差异显著性检验，根据P＜0.05筛选出差异表达miRNA。

1.5 差异miRNA靶基因预测 对差异表达分析得到的miRNA，利用miranda算法预测miRNA靶基因。miranda算法是基于miRNA-3'UTR序列匹配与能量稳定性评估计算步骤综合预测miRNA靶基因，该算法采用动态规划算法搜索miRNA与3'UTR互补同时稳定形成双链的区域，预测miRNA靶基因过程中所使用到的阈值参数为：S＞=150，ΔG＜=-30 kcal/mol和 Demand strict 5'seed pairing。其中S指匹配区域内single-residue-pair match scores；ΔG是双链形成时的自由能。

1.6 预测靶基因GO和KEGG富集分析 利用Bioconductor（3.4 版本）（http://www.bioconductor.org/）项目中的GO和KEGG富集分析软件包clusterProfiler（3.2.10版本），对预测到的差异miRNA靶基因进行GO和KEGG富集分析。

2 结果

2.1 差异表达miRNA 共发现差异表达miRNA 55个，其中21为已知miRNA（表1），在这些已知差异表达miRNA中，表达上调的有4个，分别为miR-449a-5p、miR-135a-1-3p、miR-9-5p 和 miR-615-3p；表达下调的有17个。发现新的差异表达的miRNA 34个，其中上调和下调miRNA分别为11个和23个，这些新miRNA的作用靶基因上不清楚。根据P值大小，前10个差异表达miRNA的信息见表2。

2.2 差异miRNA的靶基因预测 对21个已知差异miRNA进行靶基因预测，共找到1634个靶基因，其中值最小的miR-1262-3p预测到3个靶基因，分别为：VPS18 （基因登录号：NM_001257653.1）、FGF14（基因登录号：XM_015121402.1）和 FAM167A（基因登录号：XM_015144755.1）。

2.3 靶基因GO分析 对P≤0.05的已知miRNA进行靶基因预测，发现其作用靶基因有1492个。进一步对靶基因进行GO聚类分析，证明这些基因主要设计生物进程、细胞组份和分子功能方面。与生物进程相关的靶基因主要涉及转录及调控、细胞凋亡、信号转导、细胞分化等（图1A），与细胞组份相关的靶基因主要涉及细胞膜、核、质膜、细胞液等（图1B），与分子功能相关的靶基因主要涉及与ATP结合、DNA结合、RNA结合分子等（图1C）。

2.4 miRNA靶基因KEGG通路富集分析 对靶基因进行KEGG通路富集分析，发现这些靶基因主要富集在35条信号通路上，其中前20个显著富集的KEGG通路见图2，所涉及的KEGG通路主要包括胰岛素分泌、MAPK、PI3K-Akt、wnt、mTOR 信号通路及与细胞凋亡信号通路和免疫相关的toll样受体信号通路等。

图1GO聚类分析显示显著富集的前20个GO term柱状图

图2 富集的前20个KEGG通路柱状图

3 讨论

miRNA是机体中参与衰老调控的一类重要分子，其调控方式主要是通过碱基配对结合到衰老相关基因的mRNA的3'末端非编码区，影响mRNA的翻译，进而影响衰老进程[3，8]。有研究报道，mir-499、mir-34c、mir-1285、mir-34a 等 miRNA 能够通过与靶基因结合影响细胞衰老进程[8]。本研究共发现55个在中年与老年猕猴单个核细胞中有显著差异表达的miRNA，其中21个miRNA已被注释，上调的 4个 miRNA分别为 miR-135a-1-3p、miR-449a-5p、miR-615-3p和miR-9-5p，提示这些下调或上调miRNA的差异表达参与了机体衰老进程。进一步55个差异表达的miRNA共调控1634个靶基因。对这些靶基因进行GO分析结果显示，这些靶基因的作用主要涉及生物进程、细胞组份及ATP、DNA、RNA结合功能等；KEGG富集结果显示，这些靶基因主要被显著富集在胰岛素分泌、wnt以及mTOR等与衰老相关的通路上，主要调节细胞增殖、分化和衰老等[9]。本研究发现的miRNA可以通过调控这些通路中的关键基因表达来促进或延缓细胞衰老，可为临床上衰老相关疾病的诊断提供新的候选因子，也为衰老相关疾病的治疗提供新的靶点。

表1 识别的差异表达的已知miRNA

表2 新发现的差异表达miRNA