梅花鹿鹿茸快速生长期环状RNA鉴定及生物信息学1)

韩若冰 李和平

(东北林业大学，哈尔滨，150040)

鹿茸的快速生长及其每年的周期性再生现象在哺乳动物中极为罕见。在梅花鹿鹿茸的快速生长期，其生长速度可达到12.5 mm·d-1[1]，不仅如此，鹿茸的生长速度堪比癌细胞，却极少发生癌变，这些独特的生物学特性使得鹿茸受到研究者的广泛关注。了解鹿茸生长发育的分子调控机制对于鹿茸生物学的研究甚至是癌症治疗的进一步研究均具有重要的科学意义。

迄今为止，已有众多研究关注鹿茸生长发育的分子调控机制。随着高通量测序技术的迅速发展，选用梅花鹿鹿茸作为试验材料，从转录组水平上探究鹿茸生长发育的研究应运而生。通过对不同生长时期的梅花鹿鹿茸进行转录组测序，得到了几种可能对鹿茸生长发挥重要调控作用的转录因子：PER1、EGFR1、GAS1[2]。通过对骨化期梅花鹿鹿茸的转录组测序得到了Sp7、SOX4、MMP9、MMP13等30种与鹿茸软骨内骨化相关的基因[3]。然而关于非编码RNA在鹿茸生长发育过程中的研究还较为匮乏，主要集中于miRNA的研究，Yao et al.[4]通过高通量测序技术鉴定得到一些参与鹿茸软骨发生和软骨内骨化的miRNA(miR-3072-5p、miR-1600、miR-34-5p、miR-6889-5p、miR-6729-5p)。

环状RNA是一类在转录后水平具有重要基因表达调控作用的内源性非编码RNA[5-6]，不具有5’末端及3’末端。近些年来已经成为1个新的研究热点。然而，对于鹿科动物中环状RNA的研究还较少。本研究运用高通量测序技术鉴定鹿茸2个生长发育时期(鹿茸生长前期、快速生长期)的环状RNA并进行生物信息学分析，试图挖掘对鹿茸生长发育发挥调控作用的环状RNA，旨在为梅花鹿的多组学研究及其生长发育机理的探究提供参考。

1 材料与方法

1.1 样品采集

试验材料采集自哈尔滨金地鹿业有限公司的3只健康成年雄性梅花鹿鹿茸的顶端间充质组织。在其鹿茸生长至小鞍子茸型期，取左茸顶端组织的间充质组织为鹿茸生长前期样品；在其生长至二杠茸型时，采集右茸顶端组织的间充质组织为快速生长期样品。样品采集后，将间充质组织切成小块，迅速放入液氮中保存备用。

1.2 RNA提取及文库构建、测序

按照TRIzol试剂盒的操作说明书提取6个间充质组织样品的总RNA，使用分光光度计检测提取总RNA的浓度、纯度并使用Agilent Bioanalyzer检测总RNA的完整性。提取得到的总RNA放置于冰箱-80 ℃保存以防止RNA降解。得到质检合格的RNA后，去除其中的核糖体RNA，合成cDNA的第一链和第二链并且进行PCR扩增，利用Illumina HiseqTM 4000对构建的文库进行测序。

1.3 环状RNA的鉴定与特征分析

利用TopHat[7]将过滤后的数据与马鹿参考基因组进行比对，在比对结果中提取未比对上的数据，截取每一条未比对上的数据的两端以得到短序列读段，再次比对到参考基因组，并利用find_circ[5]软件鉴定得到环状RNA。鉴定条件为：只保留有且只有1个清晰断点的环状RNA；每条读长的2个短序列读段比对到基因组上的位置重叠不能超过2 bp且只允许2 bp错配；独特的读长需大于2条；环状RNA的长度小于100 kbp等。鉴定得到环状RNA后，对环状RNA类型、长度及环状RNA在参考基因组上的分布进行统计。

1.4 表达量计算与差异分析

采用每千个碱基的转录每百万映射读取的转录本(RPM)计算环状RNA的表达量。根据得到的环状RNA的表达量对鹿茸2个生长发育时期的环状RNA进行差异分析，从而得到显著差异表达的环状RNA，显著差异表达的环状RNA的筛选条件设置为P<0.05，|log2FC|>1。

1.5 富集分析

为探究环状RNA在鹿茸生长发育过程中发挥的作用，对环状RNA的来源基因进行GO富集分析和KEGG功能注释分类。在差异表达基因中显著性富集的通路以及差异表达基因中显著富集的GO条目筛选阈值均设置为q≤0.05。

1.6 circRNA-miRNA-mRNA调控网络的构建

运用RNAhybrid (v2.1.2) 、svm_light (v6.01)，Miranda (v3.3)、TargetScan (v7.0)对最显著的10个环状RNA进行靶向miRNA的预测，结合本研究前期测得的miRNA数据[8]，从中筛选差异表达的miRNA。除此之外，利用mirTarBase预测得到与这些差异表达的miRNA具有靶向关系的mRNA，并从中筛选出差异表达的mRNA，进而得到miRNA-mRNA靶向作用关系，最终利用cytoscape[9]构建与鹿茸生长发育相关的circRNA-miRNA-mRNA调控网络。

2 结果与分析

2.1 环状RNA的鉴定

通过对测序得到的原始数据进行过滤，得到的Q20为98.93%～99.05%，Q30为96.25%～96.58%，GC含量为51.35%～53.34%。从比对结果中提取未比对上的数据并截取其两端，得到的短序列读段为45865210～60425688，比对到参考基因组上的比对率为63.83%～67.32%。

在鹿茸发育的前期和快速生长期共鉴定得到7 790个环状RNA。环状RNA主要来源于蛋白质编码基因的已注释的外显子、单一外显子、基因反义链、内含子、外显子和内含子和基因间区[5，10]。对鉴定得到的环状RNA类型进行了统计(表1)，其中，大多数环状RNA来自于蛋白质编码基因的基因间区及外显子、内含子区，来自于内含子区的环状RNA数目最少，而环状RNA的长度主要分布于1～1 000之间(表2)。对环状RNA在染色体上的分布进行统计(表3)，结果表明7 790个环状RNA分布在参考基因组的35条染色体上，其中分布在CM008042.1上的环状RNA数量最少，而分布在CM008012.1上的环状RNA数量最多。

表1 环状RNA的类型分布

表2 环状RNA的长度分布

表3 环状RNA的染色体分布

2.2 环状RNA的差异表达分析

通过对鉴定得到的2个发育时期的环状RNA表达量进行差异基因的筛选，共得到208个显著差异表达的环状RNA，包括显著上调的39个差异表达环状RNA以及显著下调的169个差异表达环状RNA(图1)。差异最显著的10个环状RNA的详细信息见表4。

红色表示表达量上调的环状RNA；绿色表示表达量下调的环状RNA；黑色表示表达量没有差异。

表4 差异最显著的前10个环状RNA

2.3 差异表达基因的富集分析

对差异环状RNA的来源基因进行富集分析可用于环状RNA的功能预测。GO富集分析结果表明，环状RNA的来源基因富集到166条与分子功能相关的GO条目上，显著富集到生物调控，细胞过程、代谢过程等GO条目上；富集到156条与细胞组分相关的GO条目上，显著富集到细胞、细胞部分等GO条目上；富集到980条与生物过程相关的GO条目上，包括结合与催化活性等相关的GO条目(表5)。

表5 差异表达的circRNA的来源基因的GO富集分析

环状RNA来源基因的KEGG结果显示共富集到81条信号通路上，包括粘着斑、PI3K-Akt信号通路、MAPK信号通路、TNF信号通路、FoxO信号通路(图2)。其中，有13个差异显著的环状RNA的6个来源基因显著富集到与细胞增殖、细胞周期等相关的PI3K-Akt信号通路上，分别是novel_circ_004050、novel_circ_004052、novel_circ_004059、novel_circ_004060、novel_circ_004062、novel_circ_004064、novel_circ_004069、novel_circ_004074、novel_circ_004863、novel_circ_002034、novel_circ_004646、novel_circ_003420、novel_circ_006430。

图2 差异表达环状RNA的来源基因的KEGG富集分析

2.4 控网络的构建

通过预测得到的环状RNA、miRNA、mRNA之间的靶向互作关系构建最终的调控网络，我们发现在10个差异最显著的环状RNA中，有2个环状RNA预测得到的miRNA在鹿茸发育的前期和快速生长期的表达量差异不显著。因此，最终的调控网络由8个差异表达的环状RNA和与其具有靶向互作关系的25个差异表达的miRNA所构成，通过进一步预测，得到这些差异表达的miRNA有120个存在靶向关系的差异表达的mRNA(图3)。

图3 差异表达的circRNA -miRNA-mRNA调控网络

3 讨论

自从环状RNA在植物类RNA病毒中首次被发现后[11]，在多种生物体内均发现有环状RNA的存在[12]，其表达具有时空特异性。根据竞争性内源RNA机制(ceRNA)，环状RNA通过与miRNA竞争性结合进而调控下游靶基因的表达[13]，在多种肿瘤发生及组织器官的生长发育过程中均发挥着重要的调控作用，但在动物体中关于环状RNA的研究还较少，且主要集中于环状RNA的鉴定。

本研究首次鉴定了鹿茸生长发育过程中的环状RNA，得到鹿茸发育前期及快速生长期的间充质组织中共有7 790个环状RNA，其中1 901个环状RNA仅在鹿茸发育前期具有特异性表达，1 478个环状RNA仅在鹿茸的快速生长期有特异性表达，有208个环状RNA为显著差异表达。我们猜测这些在不同生长发育时期有特异性表达的环状RNA在鹿茸的某个生长发育时期发挥着独特的调控作用。除此以外，我们发现一些环状RNAs在鹿茸发育的2个时期均为显著高表达(表达量高于其他环状RNA 100倍以上，例如novel_circ_002441、novel_circ_007423、novel_circ_001328、novel_circ_007778、novel_circ_001447、novel_circ_005051、novel_circ_000884、novel_circ_005048、novel_circ_005975)，这些高表达的环状RNA通过与miRNA及mRNA之间的靶向作用关系对鹿茸的生长发育发挥着某些重要的调控作用，但具体的作用机制还需进一步的功能试验来探究。

环状RNA的功能分析主要通过其来源基因富集分析及与其具有靶向作用关系的miRNA和mRNA来进行，显著差异表达的环状RNA来源基因的GO富集分析主要集中于生物调控、细胞过程、代谢过程、细胞代谢过程，而KEGG富集到代谢通路、粘着斑、PI3K-Akt信号通路、MAPK信号通路、FoxO信号通路上，说明差异显著的环状RNA主要涉及鹿茸发育的生物调控以及代谢等过程。在构建的调控网络中，1个环状RNA可能会同时靶向多个miRNA，并且1个miRNA也可能会与多个环状RNA具有靶向互作关系。其中，novel_circ_004646的靶向能力最强，靶向了13个差异表达的miRNA，而novel_circ_005407、novel_circ_000593只靶向了1个差异表达的miRNA。在与环状RNA具有靶向关系的miRNA中，novel_circ_004646的其中1个靶向miR-457-x的靶向能力最强，与29个差异表达的mRNA均具有互作关系。但环状RNA、miRNA、mRNA三者之间是否具有真实的靶向关系，其对鹿茸发育发挥着怎样的调控作用以及如何发挥这些调控作用都需要进一步的研究。