CircRNAs in the tree shrew (Tupaia belangeri) brain during postnatal development and aging

樹鼩（qú）腦后天發(fā)育及衰老過程中circRNAs特征分析

期刊：AGING-US 影響因子：4.867

發(fā)表單位：中國醫(yī)學科學院，北京協(xié)和醫(yī)學院醫(yī)學生物學研究所樹鼩種質資源中心

1 摘要

    circular RNAs (circRNAs)是一類新的非編碼RNA，分布在各種物種和組織中。目前，對樹鼩腦中circRNAs的信息了解很少。本研究我們用RNAseq技術在樹鼩幼兒（aged 47-52 days），青年（aged 15-18 months）,老年 (aged 78-86 months)期海馬體、小腦共鑒定到35,007個circRNAs。我們發(fā)現(xiàn)在不同的腦區(qū)域及不同時間段，circRNAs整體均沒有顯著變化。但小腦傾向于下調表達。RT-PCR驗證了circRNAs的存在。KEGG主要富集到ubiquitin-mediated proteolysis（泛素介導的蛋白水解）, the MAPK signaling pathway, phosphatidylinositol signaling system（磷脂酰肌醇信號通路）, long-term depression（長期消退）, the rap1 signaling pathway, and long-term potentiation（長時程增強）。我們發(fā)現(xiàn)樹鼩29,087 (83.1%) 的circRNAs與人類是同源的。ceRNA網(wǎng)絡分析發(fā)現(xiàn)novel_circRNA_007362可能海綿吸附24個miRNAs調控UBE4B表達。本研究我們得到了樹鼩腦后天發(fā)育及衰老過程中circRNAs表達譜，有助于我們研究circRNAs在腦年齡變化及年齡相關疾病中的功能。

2 材料方法

1 收集動物組織

    所有的樹鼩均是子一代，飼養(yǎng)在中國醫(yī)學科學院，北京協(xié)和醫(yī)學院醫(yī)學生物學研究所，昆明，中國。所有樹鼩健康狀態(tài)相同，沒有明顯的腫瘤或疾病。3個發(fā)育時期（I, Y, and O）9只用于研究。I組三只年齡在47-52天，體重100-110g；Y組三只年齡15-18月，體重120-150g；O組三只年齡78-86月，體重120-150g（Table S3）。樹鼩的壽命最長10年，其壽命大約是人類的1/8。因此本研究最年輕與最老的動物約相當于人類1和57歲。手術獲得其海馬體、小腦，-80℃保存。所有動物實驗獲得中國醫(yī)學科學院，北京協(xié)和醫(yī)學院醫(yī)學生物學研究所動物研究倫理委員會的批準。

2 RNA提取和RNAseq

    TRIzol reagent (Invitrogen)提取total RNA。Agilent 2100 Bioanalyzer (Agilent Technologies)評估RNA質量。

    RNase R消化線性RNA，RNeasy MinElute Cleanup Kit (Qiagen)純化。而后去除核糖體RNA。加入fragmentation buffer將RNA片段化，random primers合成第一鏈cDNA。而后，dNTPs, RNase H, DNA polymerase I, and buffer合成第二條鏈。VAHTSTM DNA Clean Beads純化文庫片段，末端修復，添加poly-(A), 連接測序接頭。而后用uracil-N-glycosylase(UNG酶)消化第二條鏈cDNA。VAHTSTM DNA Clean Beads純化繼續(xù)文庫，PCR擴增。HiSeqTM2500測序。

3 circRNA測序分析

    原始下機數(shù)據(jù)FASTQ文件，去除含有接頭reads，去除“N”超過10% reads，和Q-value ≤ 20堿基超過50%的低質量reads。使用短片段比對工具Bowtie2將reads比對到rRNA數(shù)據(jù)庫，去除rRNA reads。而后使用TopHat2(version2.0.3.12)將reads比對到樹鼩基因組，去掉能夠匹配到基因組的reads，收集未能比對的reads用于circRNA鑒定。未能比對的reads，兩端各取20bp，然后比對到參考基因組，尋找可以連接的位點。anchor reads比對方向相反（head-to-tail）暗示可能有circRNA形成。anchor reads比對到參考基因組的結果用find_circ軟件進一步鑒定circRNAs。鑒定到的circRNAs對其類型，染色體分布，長度分布進行統(tǒng)計。Back spliced junction reads使用RPM標準化。RPM能夠評估測序深度的影響。因此得到的RPM值可以直接用于不同樣品的比較。BLAST比對circRNA到circBase進行注釋。未能注釋到的circRNAs認為是新的circRNAs。

4 無監(jiān)督聚類和PCA

    應用無監(jiān)督聚類和PCA分析樣品間的相關性。無監(jiān)督聚類使用gmodels分析，PCA用ggplot2作圖。

5 分析差異表達circRNAs

    R包edgeR用于circRNAs差異表達分析。閾值為FC ≥ 2 and p < 0.05。gglpot2繪制火山圖。

6 趨勢分析

    為了分析差異表達circRNAs的表達模式，將表達數(shù)據(jù)進行0, log2 (v1/v0), and log2(v2/v0)轉化，然后使用Short Time-series Expression Miner (STEM) software趨勢分析。p ≤ 0.05認為有意義。

7 GO KEGG富集分析

    GO注釋可以用于預測基因的功能。GO分類來源于GO 數(shù)據(jù)庫(http://www.geneontology.org)，它是一個社區(qū)性的生信分析資源使用結構化的約定的詞匯歸類基因產(chǎn)物：包括Molecular Function（分子功能）, Biological Processes（生物學過程）, Cellular Components（細胞組分）。差異表達的circRNAs  (profile 0 and profile 7)比對到GO數(shù)據(jù)庫，閾值為FDR ≤ 0.05。符合該閾值的GO術語認為是富集到的有意義GO術語。KEGG是一個通路數(shù)據(jù)庫，KEGG分析有助于理解差異表達基因參與的角色。我們選擇趨勢分析得到的circRNA (profile 0 and profile 7)，進行GO 和KEGG富集分析(http://kobas.cbi.pku.edu.cn/home.do)。Pathways p< 0.05認為有意義。ggplot2繪制散點圖。

8 共表達網(wǎng)絡構建

    通過構建共表達網(wǎng)絡鑒定circRNAs，miRNAs，mRNA的關系。我們選擇與UMP（泛素介導的蛋白水解）和long-term depression(長期消退，p < 0.05,  profile 0 在海馬體和小腦均下調)相關的circRNA，使用Cytoscape (V3.2.0)構建circRNA-miRNA-mRNA通路調控網(wǎng)絡。

9 數(shù)據(jù)獲取

    RNAseq數(shù)據(jù)已上傳GEO數(shù)據(jù)庫，編號為SRP132147。

10 同源分析

    將人、小鼠、獼猴circRNA數(shù)據(jù)與樹鼩的circRNAs序列進行BLASTN同源比對(e-value < 1e-10)。

11 qPCR驗證

    為了確認個別的circRNAs,TRIzol reagent (Invitrogen)提取total RNA, random primers by a RevertAid First Stand cDNA Synthesis Kit (GeneCopoeia, Cat AORT-0020)反轉錄為cDNA。SYBR Green Kit (GeneCopoeia, Cat AORT-0020)在PikoReal PCR system (Thermo, USA)進行qPCR，條件為：95°C for 10 min, followed by 45 cycles at 95°C for 15 s and 60°C for 60 s。引物在列表Table S2。2^-△△Ct計算相對表達量，GAPDH作為內參。PCR產(chǎn)物2%瓊脂糖凝膠分離，Sanger測序。

3 結果

1 circRNAseq概述

    測序共產(chǎn)生465,791,754 raw reads (海馬體：235,543,658，小腦：230,248,096)。去除低質量的，高“N”，含有接頭的reads后，共得到433,887,754 高質量clean reads (海馬體：216,130,360，小腦：217,757,394)(Table S1)。去除比對到rRNA的reads，而后用TopHat2比對到樹鼩基因組，選擇未能比對的reads?；趂ind_circ預測，共得到35007個circRNAs。所有鑒定到的circRNAs都是第一次發(fā)現(xiàn)。如圖Fig. 1A，大部分circRNAs的長度在400nt左右，與之前的研究一致。而且鑒定了不能類型的circRNAs, 大部分的類型是annot_exon，接著是exon_intron（Fig. 1B）。

2 不同腦區(qū)域不同時間circRNA表達譜

    對樹鼩腦不同區(qū)域的不同時間的circRNAs進行無監(jiān)督聚類分析和PCA分析以研究樣品之間的相關性和時空動態(tài)變化。不同區(qū)域間circRNA的表達譜相對獨立(Fig. S1)。箱形圖用于比較樣品間的表達量分布，以分析circRNAs在不同時間段的變化。如Fig. 2A，在不同腦區(qū)域，不同時間段，circRNAs整體上沒有顯著的改變。能在海馬體檢測到的大部分circRNAs也能在小腦鑒定到，反之亦然(Fig. 2B)。在海馬體中獨立表達的有983個circRNAs，在小腦獨立表達的有1009個。在不同時間段獨立表達的circRNAs比較少(Fig. 2C)。

3 隨著年齡的增長，小腦區(qū)circRNAs表達量傾向于下調

    我們做了circRNAs表達譜分析以確定在海馬體和小腦發(fā)育中是否有circRNAs差異表達。差異表達的計算基于back-spliced reads per million (RPM)。只有fold change (FC)≥ 2 and p-value < 0.05的認為是有差異的。我們比較了每個區(qū)域幼兒期（I）與青年期（Y），青年期（Y）與老年期（O），幼兒期（I）與老年期（O）。差異表達的circRNA在1101-3234不等(Fig. 3A)。火山圖顯示小腦區(qū)隨著年齡增長circRNAs傾向于下調 (Fig. 3B-D)。而在海馬體這種下調的趨勢僅表現(xiàn)在I與Y比較(Fig. S2A-C)。小腦區(qū)I與O比較，1,339 (4.1%) circRNAs顯著下調，935 (2.9%)顯著下調(Fig. 3D)。I 與Y比較，1,272 (3.9%)顯著下調，1,142 (3.5%)顯著上調(Fig. 3B)。Y與O比較，970 (2.9%)顯著是下調，739 (2.2%)顯著上調。

在海馬體，I 與Y比較，1,636 (5.0%)顯著下調，1,386 (4.3%)顯著上調(Fig. S2A)。而Y與O比較，512(1.7%)顯著下調，588(1.9%)顯著上調(Fig. S2B)。I與O比較，上下調的數(shù)量幾乎一致(Fig. S2C)。這些數(shù)據(jù)說明circRNAs的表達具有時空特異性。

4 腦不同時間circRNA的動態(tài)變化

    為了鑒定circRNA表達譜的動態(tài)變化，我們做了趨勢分析以鑒定腦區(qū)域不同時間段的主要差異表達的circRNAs。在海馬體和小腦中鑒定到7中主要的趨勢(Fig. 4A-B)。海馬體中4個趨勢(profiles0, 1, 6, and 7)有顯著意義，小腦3個趨勢(profiles 0, 1, and 6)有顯著意義。我們也分析了海馬體和小腦區(qū)上調和下調的circRNAs。我們發(fā)現(xiàn)海馬體有385個circRNAs上調，而小腦有197個；海馬體有334個circRNAs下調，而小腦有439個(Fig. 4C)。這些發(fā)現(xiàn)說明年齡相關異常表達的circRNAs可能在腦發(fā)育和衰老中發(fā)揮重要作用。

5 差異表達的circRNAs富集分析

    為了推測腦不同區(qū)域上下調circRNAs發(fā)揮的作用，我們對profile 0（下調circRNAs）, profile7（上調circRNAs）進行GO，KEGG富集分析。對于profile 0，GO分析顯示circRNAs在不同腦區(qū)域富集到相似的功能(Fig. S3A)。對于profile 7，GO富集到的條目差異很大(Fig. S3B)。在profile 0中海馬體區(qū)富集到14個通路滿足p < 0.05 (Tables S4-5)，而小腦區(qū)富集到30個通路。對top20個pathways繪制散點圖(Fig. 5A-B)。在profile 7（上調circRNAs）中，海馬體和小腦共鑒定到7和24個通路((p < 0.05，Tables S6-7)。柱狀圖呈現(xiàn)profile 7的KEGG通路(Fig. S4)。對于profile 0，海馬體富集到RNA degradation, the neurotrophin signaling pathway, and synaptic vesicle cycle。小腦富集到phosphatidylinositol signaling system, phospholipase D signaling pathway, PI3K-Akt signaling pathway, and ErbB signaling pathway。兩個區(qū)域共有的有ubiquitin-mediated proteolysis (UMP) pathway, MAPK signaling pathway, phosphatidylinositol signaling system, and long-term depression(Fig. 5A-B)。對于profile 7，海馬體富集到 synaptic vesicle cycle, UMP, and Fanconi anemia pathway，而小腦富集到cAMP signaling pathway, phosphatidylinositol signaling system, and insulin signaling pathway。共有富集是rap1 signaling pathway and long-term potentiation(Fig. S3C)。這些結果說明海馬體和小腦這些circRNAs與衰老相關。而且兩個區(qū)域大部分的通路可以歸類到environmental information processing, genetic information processing, and organismal systems (Tables S4-7)。

6 CircRNA-miRNA-mRNA共表達

    CircRNAs可以海綿吸附microRNAs進而調控基因表達。因此我們基于UMP通路和long-term depression構建共表達網(wǎng)絡。在小腦，novel_circRNA_007362可以結合24  miRNAs(tch-let-7e-5p, tch-let-7i-5p, tch-let-7f-5p, tch-miR-125a-5p, tch-miR-1301, tch-miR-135a-5p, tchmiR-135b-5p, tch-miR-15b-5p, tch-miR-195-5p, tchmiR-1-5p, tch-miR-218-5p, tch-miR-22-3p, tch-miR-26a-5p, tch-miR-26b-5p, tch-miR-296-3p, tch-miR-335-5p, tch-miR-34a-5p, tch-miR424-5p,tch-miR-491-5p, tch-miR-455-3p,tch-miR-503, tch-miR-592, tch-miR-9771e, and tch-miR-98-5p)?；騏BE4B可以被這24個miRNA及6個circRNAs調控(novel_circRNA_007356, novel_circRNA_007357, novel_circRNA_007362, novel_circRNA_007364, novel_circRNA_007373, and novel_circRNA_007379) (Fig. 6A)。海馬體中基因UBE4B可以被7個miRNA（6個與小腦區(qū)一致）和5個circRNAs（其中兩個與小腦區(qū)一樣）調控，如Fig. 6B。這個復雜的ceRNA網(wǎng)絡說明novel_circRNA_007362 and UBE4B可能通過miRNAs吸附結合作用靶基因參與UMP pathway調控腦發(fā)育和衰老。

7 樹鼩circRNAs驗證

    為了驗證RNAseq鑒定到的circRNAs，我們自由選擇了8個海馬體circRNAs，6個小腦circRNAs。設計合適的引物(Table S2)擴增circRNA成環(huán)連接點。與測序結果一致，6個circRNAs (novel_circ_013859, novel_circ_020413, novel_circ_000954, novel_circ_012192, novel_circ_003099, and novel_circ_018752)隨著年齡增長下調(Fig. 7A)。海馬體中得到RNAseq倍數(shù)結果與RT-PCR定量結果一致(Fig.7B–D)。novel_circ_029333在both I vs O and I vs Y中更高。而且PCR產(chǎn)物進一步用瓊脂糖凝膠驗證。其中兩個circRNAs, novel_circ_013859 and novel_circ_ 007625進一步用Sanger測序確認(Fig. 7F–G)。RT-PCR結果與RNAseq結果一致，說明鑒定到的circRNAs在體內確實有不同的表達譜。

8 樹鼩circRNAs同源性分析

    我們將樹鼩circRNAs與人，小鼠，獼猴circRNAs進行同源性分析。結果說明29,087 (83.1%)樹鼩circRNAs與人類的circRNAs是同源的(Table S8)，3,783(10.8%)與小鼠同源，392 (1.1%) 與獼猴同源。只有56個circRNAs(0.2%)在三個物種中均有鑒定到。5,695 (16.3%)未能與任一物種比對到同源序列(Fig. 8)。

4 討論

    本研究我們用RNAseq技術在樹鼩腦中鑒定到35007個circRNAs。大部分circRNAs的長度在400nt左右，與之前的報道的一致。大部分circRNAs是annot_exon，而在綿羊腦垂體中，大部分的circRNAs位于intergenic。這說明circRNAs在物種與組織的分布上存在特異性。在海馬體和小腦區(qū)觀察到豐富的circRNAs。之前有研究比較過人額皮質，肝，甲狀腺和肌肉組織發(fā)現(xiàn)，circRNAs在腦中表達最豐富。我們的發(fā)現(xiàn)與此一致。Agnieszka et al.報道，circRNA在小鼠小腦區(qū)比其他腦部地區(qū)表達更豐富，而在小腦區(qū)表達的circRNAs與線性RNA的比例也很高。但在本研究，我們發(fā)現(xiàn)海馬體與小腦circRNAs的整體表達沒有明顯的差異(Fig. 2A)。我們測序時去除了線性RNA，沒有獲得線性RNA的數(shù)據(jù)，不了解circRNAs與線性RNA表達的比例。

    我們發(fā)現(xiàn)小腦下調的circRNAs傾向于年老的組織。而在海馬體，下調的circRNAs僅傾向于在I 與Y比較時。很明顯在樹鼩腦衰老時，有些circRNAs表達量下調，而有些表達量上調。不清楚為何會有這種現(xiàn)象。年老小鼠中，Gruner et al.發(fā)現(xiàn)一個明顯的趨勢，circRNAs在皮質和海馬器上調，而心臟不是。這說明在樹鼩中，骨骼肌相關的circRNAs并沒有隨著年齡改變，盡管19個circRNAs隨著年齡增長下調，可能是由于組織或物種特異性。接下來，我們用qPCR鑒定了時間尺度上的circRNAs表達量，結果與RNAseq結果一致，因此可以確信我們的RNAseq結果是可靠地。但并清楚這些一致circRNAs的靶標為何依然有變化。

    為了得到與腦發(fā)育和衰老相關的circRNAs，我們比較了小腦區(qū)和海馬體I, Y, and O三個時期的circRNA表達量。依據(jù)GO，KEGG分析，我們關注上調和下調的circRNAs。在profile 0（下調circRNAs）,GO富集分析顯示，在這兩個區(qū)域富集到相似的功能(Fig. S3A)。而且，UMP pathway, MAPK signaling pathway, phosphatidylinositol signaling system, and long-term depression在這兩個區(qū)域均有鑒定到。其他通路如RNA degradation, the ErbB signaling pathway, and PI3K-Akt signaling pathway也伴隨上述提到的4個通路出現(xiàn)(Fig. 4C–D)。在profile 7（上調circRNAs）,GO分析顯示兩個區(qū)域富集到的功能差異很大(Fig. S3B)。KEGG分析顯示rap1 signaling pathway and long-term potentiation在兩個區(qū)域均有鑒定到。其他通路如Fanconi anemia pathway, cAMP signaling pathway, and insulin signaling pathway也伴隨著上述兩個通路出現(xiàn)(Fig.S3C)。這些發(fā)現(xiàn)顯示涉及到的通路與腦發(fā)育和衰老相關。之前研究報道UMP pathway在山羊初級毛囊，次級毛囊發(fā)育時發(fā)揮重要作用，UBE2O在控制初級毛囊，次級毛囊發(fā)育相關外顯子表達時有重要作用。Ham et al.發(fā)現(xiàn)UMP涉及到正常腦與病理腦衰老。UMP可以降解有危害的聚集的蛋白增長壽命。另外，我們對明顯參與UMP and long-term depression的circRNA-miRNA-mRNA構建ceRNA網(wǎng)絡。我們發(fā)現(xiàn)novel_circRNA_007362可以結合24個miRNA進而影響UBE4B的表達。而且，novel_circRNA_032232通過結合9 miRNAs，影響B(tài)IRC6表達(Fig. 6B)，并且間接的影響UBE4B and HERC1的表達。根據(jù)KEGG分析，UBE4B, BIRC6, and HERC1是UMP通路中關鍵的基因。（此處講了HERC1，BIRC6泛素化相關功能，此處省略）。ceRNA網(wǎng)絡顯示novel_circRNA_007362 and novel_circRNA_032232可能與24 miRNAs結合包括tch-let-7a-3p, tch-miR-136-5p, tchmiR-146a-5p, tch-miR-17-5p, tch-miR-183-5p, tch-miR-188-5p, tch-miR-23b-3p, tch-miR-326, and tch-miR-93-5p調控小腦和海馬體發(fā)育與衰老。總之大腦發(fā)育與衰老的調控是一個復雜的多種方面的，許多細節(jié)還未知。

    分析circRNAs同源性有助于我們理解共同的規(guī)律。之前研究顯示circRNAs在人和小鼠腦中是高度保守的。比如，Agnieszka et al.發(fā)現(xiàn)15,849中有4,522小鼠circRNAs是與人一致的。Guo et al. 發(fā)現(xiàn)20%的小鼠circRNAs與人類是直系同源。而在樹鼩中，約83.1%的circRNA與人類circRNA有同源關系?；蛟S是樹鼩與靈長類動物關系更近。而且，人，樹鼩中共有的保守的circRNAs的功能和產(chǎn)生機制更相似。我們在樹鼩海馬體和小腦中鑒定到大量的circRNAs。我們的結果可以為鑒定影響腦發(fā)育，正常衰老，突觸再生，重塑，Alzheimer's disease, and Parkinson’s disease相關的circRNAs提供一個借鑒。未來需要將更多研究放在與腦發(fā)育，衰老相關疾病circRNAs的功能上。

5 研究思路

材料：

    樹鼩3個發(fā)育時期（I, Y, and O）海馬體和3個發(fā)育時期（I, Y, and O）小腦用于RNAseq。

研究方法：

    1 RNAseq：TRIzol提取total RNA，Agilent 2100質檢，RNase R消化線性RNA，去核糖體RNA，鏈特異性建庫，HiSeqTM2500測序。

    2 find_circ軟件鑒定circRNAs，RPM表達量標準化，BLAST同源circRNA比對。

    3 PCA分析樣品分組情況。

    4 R包edgeR用于circRNAs差異表達分析，閾值為FC ≥ 2 and p < 0.05。STEM進行趨勢分析。

    5 GO 數(shù)據(jù)庫進行GO富集分析，KOBAS KEGG富集分析。

    6 構建共表達網(wǎng)絡，取與UMP和long-term depression相關的基因，miRNA, circRNA構建共表達網(wǎng)絡。

    7 circRNA確認，qPCR, RT-PCR凝膠電泳，Sanger測序。

結果展示：

1 circRNA描述

    CircRNAseq產(chǎn)出序列信息（table S1），circRNA長度分析線圖（Fig. 1A），circRNA類型柱狀圖（Fig. 1B）。

2 circRNA表達譜

    CircRNA分組PCA圖（Fig. S1），箱形圖比較不同樣品circRNA整體表達量（Fig. 2A），樣品circRNA鑒定數(shù)量venn圖（Fig. 2B，F(xiàn)ig. 2C）。

3 不同發(fā)育時期circRNA比較

    各個分組分別比較差異基因數(shù)量統(tǒng)計（Fig. 3A），火山圖（Fig. 3，F(xiàn)ig. S2）。

4 不同發(fā)育時期趨勢分析

    趨勢分析折線圖（Fig. 4）。

5 差異表達的circRNAs功能富集分析

    GO柱狀圖（Fig. S3），通路散點圖（Fig. 5），通路富集列表（Tables S4-7）。

6 CircRNA-miRNA-mRNA共表達

    共表達網(wǎng)絡圖(Fig. 6)。

7 樹鼩circRNAs驗證

    qPCR定量柱狀圖（Fig. 7A-D），電泳圖（Fig. 7E）,Sanger測序峰圖（Fig. 7F–G）。

8 同源分析

    Blastn同源分析結果列表（Table S8），不同物種同源結果venn圖（Fig. 8）。

微信掃描二維碼關注微信公眾號了解更多文獻信息