葡萄幼葉環(huán)狀RNA（CircRNA）的鑒定與分析

陳迎春，吳新穎，蔣錫龍，張倩倩，李益，慕茜，楊立英，王詠梅，張加魁*，王鵬飛*

（山東省葡萄研究院/農業(yè)部華東都市農業(yè)重點實驗室，山東濟南 250100）

在生物細胞內，除了含有信使RNA（messenger RNA，mRNA）之外，還含有各種類型的非編碼RNA。例如microRNA（miRNA），小干擾RNA（small interfering RNA，siRNA），反式作用RNA（trans-acting siRNA，tasiRNA）以及長鏈非編碼RNA（long noncoding RNA，lncRNA）等。而環(huán)狀RNA（circRNAs）則是一類最新發(fā)現(xiàn)的獨特的非編碼RNA[1-3]。而在其起初發(fā)現(xiàn)后的十幾年內，環(huán)狀RNA曾被認為是一種RNA剪切拼接的錯誤[4]。目前，隨著高通量深度測序技術的飛速發(fā)展，大量的環(huán)狀RNA在動物中被發(fā)現(xiàn)。研究顯示環(huán)狀RNA在動物細胞內穩(wěn)定并高水平表達，而這些被發(fā)現(xiàn)的環(huán)狀RNA也被證實在很多生物學過程中發(fā)揮重要作用[1,2,6-8]。

環(huán)狀RNA可以起源于外顯子（此類為外顯子起源環(huán)狀RNA），內含子（此類為內含子環(huán)狀RNA）以及基因間區(qū)[9-11]。也可以是部分源自內含子而部分源自外顯子（此類為外顯子-內含子起源RNA，exon-intron circRNA，EIciRNA）[12-15]，甚至可以源自轉運RNA

（tRNA）的內含子（此類為tricRNA）[16]。環(huán)狀RNA是由于RNA頭與尾的反向剪切而形成的。與線性RNA不同，環(huán)狀RNA不含有5'端的帽子結構和3'端的尾部，其形態(tài)是上游RNA片段和下游RNA片段相連形成封閉的環(huán)狀[10]。因此，環(huán)狀RNA不會被RNA酶R降解[11]。然而大部分的環(huán)狀RNA被發(fā)現(xiàn)是源自蛋白編碼基因的外顯子[17]。研究表明序列的互補及外顯子跳躍是環(huán)狀RNA形成的原因[18-22]。而有的研究顯示，環(huán)狀RNA也可以通過含外顯子的套索前體產生[23]。此外，RNA綁定蛋白也涉及環(huán)狀RNA的形成過程。例如MBNL1、ADAR1及Quaking也可以在環(huán)狀RNA合成中起到重要作用[24-26]。而環(huán)狀RNA一旦被合成，由于其缺少起始密碼子和終止密碼子，因而不能編碼蛋白序列[24,26]。

環(huán)狀RNA目前被發(fā)現(xiàn)的功能主要是作為miRNA海綿。例如已經發(fā)現(xiàn)的ciRS-7（也被稱作CDR1as），就被發(fā)現(xiàn)包含超過70個保守的miR7結合位點[27]。由于ciRS-7對miR7的結合，使得miR-7的活性大大降低，從而增加miR-7靶基因的表達水平。而環(huán)狀NRAciRS-7一旦被降解，miR7則會被釋放[28]。Sry則是另一種被證實為miRNA海綿的睪丸特異表達的環(huán)狀RNA[29]。該環(huán)狀RNA含有16個miR138的綁定位點。這些發(fā)現(xiàn)增加了對miRNA調控網(wǎng)絡的理解，并增加了我們對競爭性內源RNA（ceRNA）網(wǎng)絡機制的認識[30]。環(huán)狀RNA也可以作為一種蛋白海綿。例如果蠅和人中的環(huán)狀RNA circMbl。這種環(huán)狀RNA上存在許多的muscle blind蛋白結合位點[26]。環(huán)狀RNA circMbl可以清除多余的muscle blind蛋白，從而調控該蛋白的表達水平[26]。有的環(huán)狀RNA（ecircRNA）可作為“mRNA陷阱”，隔離翻譯起始站點，從而導致線性mRNA無法翻譯。例如，小鼠formin（FMN）基因可產生作為mRNA陷阱的ecircRNA[31]。此外，環(huán)狀RNA可以通過與Pol II的互作從而正調控Pol II轉錄[14]。最近研究還顯示，EIciRNAs可以通過與U1 snRNA的互作微調控父母本基因的表達[16]。這一發(fā)現(xiàn)不僅揭示了環(huán)狀RNA在轉錄調控中的作用，但也揭示了特異EIciRNA和U1 snRNA相互作用的調控機制?？偟膩碚f，環(huán)狀RNA與轉錄機制相互作用細胞中基因表達調控機制提供了新的觀點。

在哺乳動物中，環(huán)狀RNA在突觸中高表達，并且在神經元分化過程中差異表達[26]。環(huán)狀RNA在許多腫瘤中也有廣泛的表達，其表達水平與人類腫瘤的臨床特征密切相關。因此，環(huán)狀RNA在癌癥中將可能被作為一種疾病的生物標志物[32]。而與對動物環(huán)狀RNA的研究相比，對植物中環(huán)狀RNA的研究相對較少[33]。在水稻中，一些環(huán)狀RNA在磷充足和磷饑餓的條件下差異表達，顯示了環(huán)狀RNA可能在對磷饑餓的應激反應中起到作用[34]。這些結果表明，環(huán)狀RNA也在植物中大量存在，并且可能在非生物脅迫響應中起到重要的作用。

葡萄是一種重要的果樹。而目前對葡萄環(huán)狀RNA的研究尚未開展。本研究將利用高通量測序及計算機預測的方法鑒定葡萄葉中環(huán)狀RNA的數(shù)量、種類以及在葉中表達水平。并探索其來源的基因及其功能。并初步預測這些葡萄環(huán)狀RNA可能靶向的miRNA。本研究將豐富對葡萄中環(huán)狀RNA的了解，并為葡萄miRNA調控網(wǎng)絡的研究打下基礎。

1 材料和方法

1.1 試驗材料

采集釀酒葡萄品種‘赤霞珠’一年生自根苗的幼葉，用于總RNA的提取。

1.2 方法

1.2.1 RNA提取

利用TRIZOL試劑盒（購自Invitrogen公司,USA）提取幼葉的總RNA，操作步驟按照試劑盒說明書。

1.2.2 環(huán)狀RNA文庫構建及高通量測序

用Ribo-ZeroTM Magnetic kit植物葉片專用試劑盒去除總RNA中rRNA。委托華大基因公司構建環(huán)狀RNA測序文庫并進行高通量測序，測序平臺為Illumina Hiseq 2000系統(tǒng)。

1.2.3 生物信息學法鑒定環(huán)狀RNA

葡萄基因組序列被下載自葡萄基因組網(wǎng)站（http：//genomes.cribi.unipd.it/grape/index.php）。得到的clean read用Bowtie2軟件比對到葡萄基因組上，去除不能比對上的read，留下能比對到基因組上的read進行下一步分析。利用CIRI和find circ軟件分析比對到基因組上的read，找到接合位點測序讀段對。CIRI通過兩次掃描比對生成的SAM（Sequence Alignment/Map）文件來檢測circRNA。最后，經過一系列過濾得到候選的環(huán)狀RNA。其本質就是找到正確的接合位點測序讀段對，而依據(jù)接合位點測序讀段對判斷出這是環(huán)狀RNA的一部分，從而鑒定發(fā)現(xiàn)環(huán)狀RNA。

1.2.4 環(huán)狀RNA表達量的分析

根據(jù)比對環(huán)狀RNA的接合位點測序讀段對數(shù)來計算環(huán)狀RNA的表達量，由于使用了CIRI、fnd circ這兩個軟件來預測，取兩者最終的接合位點測序讀段對數(shù)結果的平均值。本文采用RPB作為環(huán)狀RNA的均一化表達量數(shù)值。RPB＝比對上基因組的所有reads標準化到十億后跨過back-spliced位點的junction reads數(shù)目。

1.2.5 Nr和GO注釋與分類，KEGG注釋及KEGG pathway通路分析

利用blastP軟件在NCBI Nr數(shù)據(jù)庫檢索涉及基因編碼蛋白的功能注釋。利用在線軟件Blast2Go（https：//www.blast2go.com/）對該研究涉及基因編碼的蛋白進行分析，搜索其對應的GO功能注釋。然后利用在線軟件BGI WEGO（http：//wego.genomics.org.cn/cgi-bin/wego/index.pl）對這些注釋過的基因進行分類。利用KEGG在線數(shù)據(jù)庫（http：//www.kegg.jp/kegg/ko.html）中的線軟件BlastKOALA（http：//www.kegg.jp/blastkoala/）對涉及的基因進行比對分析，搜索其對應的KEGG功能注釋及KO號。利用這些基因KO號在KO數(shù)據(jù)庫（https：//www.kegg.jp/kegg/ko.html）中進行搜索，從而比對到這些基因所處的KEGG 通路。

1.2.6 靶向環(huán)狀RNA的miRNA初步預測及miRNA靶基因的預測

將利用psRobot（http：//omicslab.genetics.ac.cn/psRobot/）和psRNATarget（http：//plantgrn.noble.org/psRNATarget/）在線軟件預測miRNA的靶基因以及與miRNA可以結合的環(huán)狀RNA。葡萄miRNA序列下載自miRbase（http：//www.mirbase.org/）。

2 結果與分析

2.1 環(huán)狀RNA的鑒定及其表達分析

綜合兩個軟件的鑒定結果，在葡萄幼葉中共鑒定出1172個環(huán)狀RNA。根據(jù)其在基因組上的起始、終止位置定位，將其分為基因來源的環(huán)狀RNA和基因間區(qū)來源的環(huán)狀RNA。其中基因來源的環(huán)狀RNA有1147個，這些環(huán)狀RNA來源于1478個蛋白編碼基因。基因間區(qū)來源的環(huán)狀RNA有25個。這些環(huán)狀RNA在基因組上的位置顯示，其在葡萄的各個染色體上均有分布。暫時以這些葡萄環(huán)狀RNA的染色體號，及在染色體上的起始、終止位置命名，作為其ID。

通過分析這些葡萄環(huán)狀R N A表達量認為，表達量最高的環(huán)狀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。在這些環(huán)狀RNA中，表達量最低的環(huán)狀RNA表達量（RPB）為81。

2.2 環(huán)狀RNA來源基因的Nr，GO和KEGG注釋分析

我們分析了這些葡萄環(huán)狀R N A來源基因的功能。Nr注釋結果顯示，這些環(huán)狀RNA來源基因的功能注釋為942種。這些基因中，有的具有抗病相關功能，例如抗TMV蛋白N，抗病基因座類受體蛋白激酶，抗病蛋白（VIT_201s0011g01110，CircRNA-chr1∶968621|979574）。有的與發(fā)育相關，例如類細胞分裂周期蛋白5。有的與激素信號相關，例如AFR。有的與植物抗逆相關，例如hsp70。有的與表觀修飾相關，例如DNA（胞嘧啶-5）甲基轉移酶1。一些環(huán)狀RNA來源于轉錄因子家族，例如NAC轉錄因子家族，GATA轉錄因子家族。此外，這些環(huán)狀RNA的來源基因還包括查爾酮合成酶、鯊烯合酶、細胞色素 P450 CYP72A219、LEAF RUST 10、NAD(P)H脫氫酶、質磷酸磷酸酶2同種型X2、肌動蛋白相關蛋白3同種型X1。

GO注釋結果顯示，這些環(huán)狀RNA來源基因可以被分為3個大類：生物學過程、分子功能以及細胞組分。這些基因被注釋為細胞組分的有358個，被注釋為分子功能的為431個，而被注釋為生物學過程的有341個。這些基因在細胞組分大類下又被分為14個條目。其中大多數(shù)基因被注釋為細胞、細胞部分以及細胞器。這些基因在分子功能大類下又可被分為8個條目。其中大多數(shù)基因被注釋為催化活性，綁定及刺激反應。這些基因在生物學過程大類下又可被分為16個條目，其中大多數(shù)基因被注釋為代謝過程、細胞過程及定位（圖1）。

圖1 環(huán)狀RNA來源基因的GO分類Figure 1 GO classification of circular RNA-derived genes

KEGG注釋結果顯示，這些環(huán)狀RNA來源基因可以被比對到276個KEGG通路上，涉及的生物學過程包括：代謝途徑、次級代謝途徑、碳代謝、RNA降解途徑、剪切體及核糖體等生物學過程（表1）。

表1 環(huán)狀RNA來源基因富集在KEGG通路中的前20種通路Table 1 Top 20 pathways in which the circular RNA-derived genes are enriched in the KEGG pathway

2.3 與環(huán)狀RNA結合miRNA的鑒定與分析

根據(jù)生物信息學預測，共鑒定出可以靶向這些環(huán)狀RNA的miRNA73個，包括vvi-miR393a、vvi-miR393b、vvi-miR396a、vvi-miR535b、vvi-miR845c、vvi-miR156a等。涉及最大的miRNA家族為miR156、miR166家族。在這些環(huán)狀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靶向的環(huán)狀RNA為CircRNA- chr2：15496105|15589650。該環(huán)狀RNA的來源基因為VIT_202s0033g00850，也是vvi-miR3630-3p靶基因之一。這些靶基因涉及許多生物學功能，有的具有抗逆相關功能，例如HSP90.1，HSP83及LEA2。有的具有抗病相關功能，例如類RPP13蛋白1及抗TMV resistance類 N-蛋白。有的與激素信號相關，例如類生長素蛋白、AFR18、AFR23、類生長四轉運蛋白、生長素誘導單筆及乙烯不敏感蛋白 2。有的與發(fā)育相關，例如細胞周期檢查點控制蛋白、細胞分裂周期20.2、類APC復合物輔因子、細胞分裂蛋白FtsZ同系物1、生長調節(jié)因子1。很多靶基因為轉錄因子，例如類AP2、TOE3、ERF038、ERF084、類GATA24-like、myb、NAC25、bHLH77及GAMYB。此外，這些靶基因還包括類黃酮3',5'-羥化酶2、SPL16、SPL6和SPL7，及跨膜蛋白45B。而SPL基因和未知功能的基因占據(jù)最多比例。

GO注釋結果顯示，這些靶基因也可以被分為3個大類：分別為其中被注釋為細胞組分的有516個，被注釋為分子功能的為559個，而被注釋為生物學過程的有528個。在細胞組分大類下又可被分為16個條目，其中大多數(shù)靶基因被注釋為細胞、細胞部分以及細胞器；在分子功能大類下又可被分為12個條目，其中大多數(shù)基因被注釋為催化活性、綁定及生物調控；在生物學過程大類下又可被分為30個條目，其中大多數(shù)靶基因被注釋為代謝過程、細胞過程及刺激反應，圖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個葡萄幼葉中的環(huán)狀RNA。之前的研究顯示，Ye等在水稻和擬南芥中鑒定出12037個和6012個環(huán)狀RNA[34]。Lu等報道了2354個水稻中的環(huán)狀RNA[35]。Wang等在小麥中分離出88個環(huán)狀RNA[36]。Zuo等在番茄中發(fā)現(xiàn)854個環(huán)狀RNA，其中163個環(huán)狀RNA顯示出了對低溫的響應[38]。Zhao等在大豆中發(fā)現(xiàn)了5372個環(huán)狀RNA[37]。我們在葡萄中鑒定的環(huán)狀RNA數(shù)目與番茄相似，比其他物種少。這可能是由于我們只選取了一個組織進行鑒定。

被鑒定環(huán)狀RNA來源的基因涉及很多功能，包括抗逆、抗病、發(fā)育等。目前已知的環(huán)狀RNA一個重要功能就是作為CeRNA機制的一部分，通過參與調控miRNA來調控miRNA靶基因的表達[28]。因此分析這些環(huán)狀RNA可能涉及的miRNA，發(fā)現(xiàn)這些環(huán)狀RNA可以結合73種葡萄miRNA。這1172個環(huán)狀RNA中有74個可以結合miRNA，大部分不能結合。說明可能很多環(huán)狀RNA不涉及ceRNA機制。這些能結合miRNA的環(huán)狀RNA證明在葡萄中也存在ceRNA機制，即“環(huán)狀RNA-miRNA-靶基因”的三聯(lián)單元。例如“CircRNA-chr2：15496105|15589650-vvi-miR3630-3p-靶基因”單元比較復雜，涉及的靶基因較多。

我們比較了靶基因和環(huán)狀RNA來源基因，發(fā)現(xiàn)有很多功能類似，但是也有很多功能不同。GO分析顯示二者最富集的條目都是細胞、細胞部分、細胞器、催化活性及綁定。但是二者被注釋的條目種類差別很大。KEGG分析顯示二者所包含成員被比對到的KEGG通路中種類也有所不同。例如只有環(huán)狀RNA來源基因中有成員可以被定位到剪切體及核糖體通路，而只有miRNA靶基因中有成員可以被比對到植物激素信號轉導通路，植物與病原菌互作通路及泛素介導的蛋白水解通路。

本結果顯示，葡萄中有限的或較少的環(huán)狀RNA可以通過調控miRNA來調控更多不同種類靶基因，從而發(fā)揮多種功能。而預測環(huán)狀RNA結合miRNA方法，主要是檢測本測序得到的接合位點的測序讀段對序列上的miRNA結合位點。如果將來能夠通過實驗確定全部環(huán)狀RNA的完整序列，則可能會發(fā)現(xiàn)更多的miRNA結合位點。