葡萄幼叶环状RNA（CircRNA）的鉴定与分析

陈迎春，吴新颖，蒋锡龙，张倩倩，李益，慕茜，杨立英，王咏梅，张加魁*，王鹏飞*

（山东省葡萄研究院/农业部华东都市农业重点实验室，山东济南 250100）

在生物细胞内，除了含有信使RNA（messenger RNA，mRNA）之外，还含有各种类型的非编码RNA。例如microRNA（miRNA），小干扰RNA（small interfering RNA，siRNA），反式作用RNA（trans-acting siRNA，tasiRNA）以及长链非编码RNA（long noncoding RNA，lncRNA）等。而环状RNA（circRNAs）则是一类最新发现的独特的非编码RNA[1-3]。而在其起初发现后的十几年内，环状RNA曾被认为是一种RNA剪切拼接的错误[4]。目前，随着高通量深度测序技术的飞速发展，大量的环状RNA在动物中被发现。研究显示环状RNA在动物细胞内稳定并高水平表达，而这些被发现的环状RNA也被证实在很多生物学过程中发挥重要作用[1,2,6-8]。

环状RNA可以起源于外显子（此类为外显子起源环状RNA），内含子（此类为内含子环状RNA）以及基因间区[9-11]。也可以是部分源自内含子而部分源自外显子（此类为外显子-内含子起源RNA，exon-intron circRNA，EIciRNA）[12-15]，甚至可以源自转运RNA

（tRNA）的内含子（此类为tricRNA）[16]。环状RNA是由于RNA头与尾的反向剪切而形成的。与线性RNA不同，环状RNA不含有5'端的帽子结构和3'端的尾部，其形态是上游RNA片段和下游RNA片段相连形成封闭的环状[10]。因此，环状RNA不会被RNA酶R降解[11]。然而大部分的环状RNA被发现是源自蛋白编码基因的外显子[17]。研究表明序列的互补及外显子跳跃是环状RNA形成的原因[18-22]。而有的研究显示，环状RNA也可以通过含外显子的套索前体产生[23]。此外，RNA绑定蛋白也涉及环状RNA的形成过程。例如MBNL1、ADAR1及Quaking也可以在环状RNA合成中起到重要作用[24-26]。而环状RNA一旦被合成，由于其缺少起始密码子和终止密码子，因而不能编码蛋白序列[24,26]。

环状RNA目前被发现的功能主要是作为miRNA海绵。例如已经发现的ciRS-7（也被称作CDR1as），就被发现包含超过70个保守的miR7结合位点[27]。由于ciRS-7对miR7的结合，使得miR-7的活性大大降低，从而增加miR-7靶基因的表达水平。而环状NRAciRS-7一旦被降解，miR7则会被释放[28]。Sry则是另一种被证实为miRNA海绵的睾丸特异表达的环状RNA[29]。该环状RNA含有16个miR138的绑定位点。这些发现增加了对miRNA调控网络的理解，并增加了我们对竞争性内源RNA（ceRNA）网络机制的认识[30]。环状RNA也可以作为一种蛋白海绵。例如果蝇和人中的环状RNA circMbl。这种环状RNA上存在许多的muscle blind蛋白结合位点[26]。环状RNA circMbl可以清除多余的muscle blind蛋白，从而调控该蛋白的表达水平[26]。有的环状RNA（ecircRNA）可作为“mRNA陷阱”，隔离翻译起始站点，从而导致线性mRNA无法翻译。例如，小鼠formin（FMN）基因可产生作为mRNA陷阱的ecircRNA[31]。此外，环状RNA可以通过与Pol II的互作从而正调控Pol II转录[14]。最近研究还显示，EIciRNAs可以通过与U1 snRNA的互作微调控父母本基因的表达[16]。这一发现不仅揭示了环状RNA在转录调控中的作用，但也揭示了特异EIciRNA和U1 snRNA相互作用的调控机制。总的来说，环状RNA与转录机制相互作用细胞中基因表达调控机制提供了新的观点。

在哺乳动物中，环状RNA在突触中高表达，并且在神经元分化过程中差异表达[26]。环状RNA在许多肿瘤中也有广泛的表达，其表达水平与人类肿瘤的临床特征密切相关。因此，环状RNA在癌症中将可能被作为一种疾病的生物标志物[32]。而与对动物环状RNA的研究相比，对植物中环状RNA的研究相对较少[33]。在水稻中，一些环状RNA在磷充足和磷饥饿的条件下差异表达，显示了环状RNA可能在对磷饥饿的应激反应中起到作用[34]。这些结果表明，环状RNA也在植物中大量存在，并且可能在非生物胁迫响应中起到重要的作用。

葡萄是一种重要的果树。而目前对葡萄环状RNA的研究尚未开展。本研究将利用高通量测序及计算机预测的方法鉴定葡萄叶中环状RNA的数量、种类以及在叶中表达水平。并探索其来源的基因及其功能。并初步预测这些葡萄环状RNA可能靶向的miRNA。本研究将丰富对葡萄中环状RNA的了解，并为葡萄miRNA调控网络的研究打下基础。

1 材料和方法

1.1 试验材料

采集酿酒葡萄品种‘赤霞珠’一年生自根苗的幼叶，用于总RNA的提取。

1.2 方法

1.2.1 RNA提取

利用TRIZOL试剂盒（购自Invitrogen公司,USA）提取幼叶的总RNA，操作步骤按照试剂盒说明书。

1.2.2 环状RNA文库构建及高通量测序

用Ribo-ZeroTM Magnetic kit植物叶片专用试剂盒去除总RNA中rRNA。委托华大基因公司构建环状RNA测序文库并进行高通量测序，测序平台为Illumina Hiseq 2000系统。

1.2.3 生物信息学法鉴定环状RNA

葡萄基因组序列被下载自葡萄基因组网站（http：//genomes.cribi.unipd.it/grape/index.php）。得到的clean read用Bowtie2软件比对到葡萄基因组上，去除不能比对上的read，留下能比对到基因组上的read进行下一步分析。利用CIRI和find circ软件分析比对到基因组上的read，找到接合位点测序读段对。CIRI通过两次扫描比对生成的SAM（Sequence Alignment/Map）文件来检测circRNA。最后，经过一系列过滤得到候选的环状RNA。其本质就是找到正确的接合位点测序读段对，而依据接合位点测序读段对判断出这是环状RNA的一部分，从而鉴定发现环状RNA。

1.2.4 环状RNA表达量的分析

根据比对环状RNA的接合位点测序读段对数来计算环状RNA的表达量，由于使用了CIRI、fnd circ这两个软件来预测，取两者最终的接合位点测序读段对数结果的平均值。本文采用RPB作为环状RNA的均一化表达量数值。RPB＝比对上基因组的所有reads标准化到十亿后跨过back-spliced位点的junction reads数目。

1.2.5 Nr和GO注释与分类，KEGG注释及KEGG pathway通路分析

利用blastP软件在NCBI Nr数据库检索涉及基因编码蛋白的功能注释。利用在线软件Blast2Go（https：//www.blast2go.com/）对该研究涉及基因编码的蛋白进行分析，搜索其对应的GO功能注释。然后利用在线软件BGI WEGO（http：//wego.genomics.org.cn/cgi-bin/wego/index.pl）对这些注释过的基因进行分类。利用KEGG在线数据库（http：//www.kegg.jp/kegg/ko.html）中的线软件BlastKOALA（http：//www.kegg.jp/blastkoala/）对涉及的基因进行比对分析，搜索其对应的KEGG功能注释及KO号。利用这些基因KO号在KO数据库（https：//www.kegg.jp/kegg/ko.html）中进行搜索，从而比对到这些基因所处的KEGG 通路。

1.2.6 靶向环状RNA的miRNA初步预测及miRNA靶基因的预测

将利用psRobot（http：//omicslab.genetics.ac.cn/psRobot/）和psRNATarget（http：//plantgrn.noble.org/psRNATarget/）在线软件预测miRNA的靶基因以及与miRNA可以结合的环状RNA。葡萄miRNA序列下载自miRbase（http：//www.mirbase.org/）。

2 结果与分析

2.1 环状RNA的鉴定及其表达分析

综合两个软件的鉴定结果，在葡萄幼叶中共鉴定出1172个环状RNA。根据其在基因组上的起始、终止位置定位，将其分为基因来源的环状RNA和基因间区来源的环状RNA。其中基因来源的环状RNA有1147个，这些环状RNA来源于1478个蛋白编码基因。基因间区来源的环状RNA有25个。这些环状RNA在基因组上的位置显示，其在葡萄的各个染色体上均有分布。暂时以这些葡萄环状RNA的染色体号，及在染色体上的起始、终止位置命名，作为其ID。

通过分析这些葡萄环状R N A表达量认为，表达量最高的环状RNA为chr12∶3260512|3271102，其表达量（R P B）为5 0 3 1 3；表达量第二高的为c h r 6∶5 6 2 2 4 1 3|5 6 2 9 7 5 6；第三高的为chr4∶3368165|3380481。在这些环状RNA中，表达量最低的环状RNA表达量（RPB）为81。

2.2 环状RNA来源基因的Nr，GO和KEGG注释分析

我们分析了这些葡萄环状R N A来源基因的功能。Nr注释结果显示，这些环状RNA来源基因的功能注释为942种。这些基因中，有的具有抗病相关功能，例如抗TMV蛋白N，抗病基因座类受体蛋白激酶，抗病蛋白（VIT_201s0011g01110，CircRNA-chr1∶968621|979574）。有的与发育相关，例如类细胞分裂周期蛋白5。有的与激素信号相关，例如AFR。有的与植物抗逆相关，例如hsp70。有的与表观修饰相关，例如DNA（胞嘧啶-5）甲基转移酶1。一些环状RNA来源于转录因子家族，例如NAC转录因子家族，GATA转录因子家族。此外，这些环状RNA的来源基因还包括查尔酮合成酶、鲨烯合酶、细胞色素 P450 CYP72A219、LEAF RUST 10、NAD(P)H脱氢酶、质磷酸磷酸酶2同种型X2、肌动蛋白相关蛋白3同种型X1。

GO注释结果显示，这些环状RNA来源基因可以被分为3个大类：生物学过程、分子功能以及细胞组分。这些基因被注释为细胞组分的有358个，被注释为分子功能的为431个，而被注释为生物学过程的有341个。这些基因在细胞组分大类下又被分为14个条目。其中大多数基因被注释为细胞、细胞部分以及细胞器。这些基因在分子功能大类下又可被分为8个条目。其中大多数基因被注释为催化活性，绑定及刺激反应。这些基因在生物学过程大类下又可被分为16个条目，其中大多数基因被注释为代谢过程、细胞过程及定位（图1）。

图1 环状RNA来源基因的GO分类Figure 1 GO classification of circular RNA-derived genes

KEGG注释结果显示，这些环状RNA来源基因可以被比对到276个KEGG通路上，涉及的生物学过程包括：代谢途径、次级代谢途径、碳代谢、RNA降解途径、剪切体及核糖体等生物学过程（表1）。

表1 环状RNA来源基因富集在KEGG通路中的前20种通路Table 1 Top 20 pathways in which the circular RNA-derived genes are enriched in the KEGG pathway

2.3 与环状RNA结合miRNA的鉴定与分析

根据生物信息学预测，共鉴定出可以靶向这些环状RNA的miRNA73个，包括vvi-miR393a、vvi-miR393b、vvi-miR396a、vvi-miR535b、vvi-miR845c、vvi-miR156a等。涉及最大的miRNA家族为miR156、miR166家族。在这些环状RNA中，能被miRNA靶向的仅有74个。我们又利用生物信息学方法预测了这73个miRNA的靶基因，共有710个靶基因被鉴定。很多涉及的miRNA都对应了10以上靶基因，例如vvi-miR3630-3p，其靶基因包括13个，分别是VIT_208s0040g00990、VIT_207s0031g02270、VIT_207s0005g05420、VIT_204s0023g00310、VIT_202s0033g00870、VIT_202s0033g00850、VIT_202s0033g00840、VIT_202s0033g00800、VIT_202s0033g00790、VIT_202s0033g00700、VIT_202s0033g00670、VIT_202s0033g00660及VIT_216s0098g00970。vvi-miR3630-3p靶向的环状RNA为CircRNA- chr2：15496105|15589650。该环状RNA的来源基因为VIT_202s0033g00850，也是vvi-miR3630-3p靶基因之一。这些靶基因涉及许多生物学功能，有的具有抗逆相关功能，例如HSP90.1，HSP83及LEA2。有的具有抗病相关功能，例如类RPP13蛋白1及抗TMV resistance类 N-蛋白。有的与激素信号相关，例如类生长素蛋白、AFR18、AFR23、类生长四转运蛋白、生长素诱导单笔及乙烯不敏感蛋白 2。有的与发育相关，例如细胞周期检查点控制蛋白、细胞分裂周期20.2、类APC复合物辅因子、细胞分裂蛋白FtsZ同系物1、生长调节因子1。很多靶基因为转录因子，例如类AP2、TOE3、ERF038、ERF084、类GATA24-like、myb、NAC25、bHLH77及GAMYB。此外，这些靶基因还包括类黄酮3',5'-羟化酶2、SPL16、SPL6和SPL7，及跨膜蛋白45B。而SPL基因和未知功能的基因占据最多比例。

GO注释结果显示，这些靶基因也可以被分为3个大类：分别为其中被注释为细胞组分的有516个，被注释为分子功能的为559个，而被注释为生物学过程的有528个。在细胞组分大类下又可被分为16个条目，其中大多数靶基因被注释为细胞、细胞部分以及细胞器；在分子功能大类下又可被分为12个条目，其中大多数基因被注释为催化活性、绑定及生物调控；在生物学过程大类下又可被分为30个条目，其中大多数靶基因被注释为代谢过程、细胞过程及刺激反应，图2。

图2 靶基因的GO分类Figure 2 GO classification of target genes

KEGG注释结果显示，这些靶基因可以被比对到184个KEGG通路上，涉及的生物学过程包括：代谢途径、次级代谢途径，植物激素信号转导途径，植物与病原菌互作途径，泛素介导的蛋白水解及碳代谢等生物学过程（表2）。

表2 靶基因富集在KEGG通路中的前20种通路Table 2 Top 20 pathways in which the target genes are enriched in the KEGG pathway

3 讨论与结论

在本研究中，共鉴定了1172个葡萄幼叶中的环状RNA。之前的研究显示，Ye等在水稻和拟南芥中鉴定出12037个和6012个环状RNA[34]。Lu等报道了2354个水稻中的环状RNA[35]。Wang等在小麦中分离出88个环状RNA[36]。Zuo等在番茄中发现854个环状RNA，其中163个环状RNA显示出了对低温的响应[38]。Zhao等在大豆中发现了5372个环状RNA[37]。我们在葡萄中鉴定的环状RNA数目与番茄相似，比其他物种少。这可能是由于我们只选取了一个组织进行鉴定。

被鉴定环状RNA来源的基因涉及很多功能，包括抗逆、抗病、发育等。目前已知的环状RNA一个重要功能就是作为CeRNA机制的一部分，通过参与调控miRNA来调控miRNA靶基因的表达[28]。因此分析这些环状RNA可能涉及的miRNA，发现这些环状RNA可以结合73种葡萄miRNA。这1172个环状RNA中有74个可以结合miRNA，大部分不能结合。说明可能很多环状RNA不涉及ceRNA机制。这些能结合miRNA的环状RNA证明在葡萄中也存在ceRNA机制，即“环状RNA-miRNA-靶基因”的三联单元。例如“CircRNA-chr2：15496105|15589650-vvi-miR3630-3p-靶基因”单元比较复杂，涉及的靶基因较多。

我们比较了靶基因和环状RNA来源基因，发现有很多功能类似，但是也有很多功能不同。GO分析显示二者最富集的条目都是细胞、细胞部分、细胞器、催化活性及绑定。但是二者被注释的条目种类差别很大。KEGG分析显示二者所包含成员被比对到的KEGG通路中种类也有所不同。例如只有环状RNA来源基因中有成员可以被定位到剪切体及核糖体通路，而只有miRNA靶基因中有成员可以被比对到植物激素信号转导通路，植物与病原菌互作通路及泛素介导的蛋白水解通路。

本结果显示，葡萄中有限的或较少的环状RNA可以通过调控miRNA来调控更多不同种类靶基因，从而发挥多种功能。而预测环状RNA结合miRNA方法，主要是检测本测序得到的接合位点的测序读段对序列上的miRNA结合位点。如果将来能够通过实验确定全部环状RNA的完整序列，则可能会发现更多的miRNA结合位点。