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

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

树鼩（qú）脑后天发育及衰老过程中circRNAs特征分析

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

发表单位：中国医学科学院，北京协和医学院医学生物学研究所树鼩种质资源中心

1 摘要

    circular RNAs (circRNAs)是一类新的非编码RNA，分布在各种物种和组织中。目前，对树鼩脑中circRNAs的信息了解很少。本研究我们用RNAseq技术在树鼩幼儿（aged 47-52 days），青年（aged 15-18 months）,老年 (aged 78-86 months)期海马体、小脑共鉴定到35,007个circRNAs。我们发现在不同的脑区域及不同时间段，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（长时程增强）。我们发现树鼩29,087 (83.1%) 的circRNAs与人类是同源的。ceRNA网络分析发现novel_circRNA_007362可能海绵吸附24个miRNAs调控UBE4B表达。本研究我们得到了树鼩脑后天发育及衰老过程中circRNAs表达谱，有助于我们研究circRNAs在脑年龄变化及年龄相关疾病中的功能。

2 材料方法

1 收集动物组织

    所有的树鼩均是子一代，饲养在中国医学科学院，北京协和医学院医学生物学研究所，昆明，中国。所有树鼩健康状态相同，没有明显的肿瘤或疾病。3个发育时期（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℃保存。所有动物实验获得中国医学科学院，北京协和医学院医学生物学研究所动物研究伦理委员会的批准。

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纯化继续文库，PCR扩增。HiSeqTM2500测序。

3 circRNA测序分析

    原始下机数据FASTQ文件，去除含有接头reads，去除“N”超过10% reads，和Q-value ≤ 20碱基超过50%的低质量reads。使用短片段比对工具Bowtie2将reads比对到rRNA数据库，去除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对其类型，染色体分布，长度分布进行统计。Back spliced junction reads使用RPM标准化。RPM能够评估测序深度的影响。因此得到的RPM值可以直接用于不同样品的比较。BLAST比对circRNA到circBase进行注释。未能注释到的circRNAs认为是新的circRNAs。

4 无监督聚类和PCA

    应用无监督聚类和PCA分析样品间的相关性。无监督聚类使用gmodels分析，PCA用ggplot2作图。

5 分析差异表达circRNAs

    R包edgeR用于circRNAs差异表达分析。阈值为FC ≥ 2 and p < 0.05。gglpot2绘制火山图。

6 趋势分析

    为了分析差异表达circRNAs的表达模式，将表达数据进行0, log2 (v1/v0), and log2(v2/v0)转化，然后使用Short Time-series Expression Miner (STEM) software趋势分析。p ≤ 0.05认为有意义。

7 GO KEGG富集分析

    GO注释可以用于预测基因的功能。GO分类来源于GO 数据库(http://www.geneontology.org)，它是一个社区性的生信分析资源使用结构化的约定的词汇归类基因产物：包括Molecular Function（分子功能）, Biological Processes（生物学过程）, Cellular Components（细胞组分）。差异表达的circRNAs  (profile 0 and profile 7)比对到GO数据库，阈值为FDR ≤ 0.05。符合该阈值的GO术语认为是富集到的有意义GO术语。KEGG是一个通路数据库，KEGG分析有助于理解差异表达基因参与的角色。我们选择趋势分析得到的circRNA (profile 0 and profile 7)，进行GO 和KEGG富集分析(http://kobas.cbi.pku.edu.cn/home.do)。Pathways p< 0.05认为有意义。ggplot2绘制散点图。

8 共表达网络构建

    通过构建共表达网络鉴定circRNAs，miRNAs，mRNA的关系。我们选择与UMP（泛素介导的蛋白水解）和long-term depression(长期消退，p < 0.05,  profile 0 在海马体和小脑均下调)相关的circRNA，使用Cytoscape (V3.2.0)构建circRNA-miRNA-mRNA通路调控网络。

9 数据获取

    RNAseq数据已上传GEO数据库，编号为SRP132147。

10 同源分析

    将人、小鼠、猕猴circRNA数据与树鼩的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产物2%琼脂糖凝胶分离，Sanger测序。

3 结果

1 circRNAseq概述

    测序共产生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。基于find_circ预测，共得到35007个circRNAs。所有鉴定到的circRNAs都是第一次发现。如图Fig. 1A，大部分circRNAs的长度在400nt左右，与之前的研究一致。而且鉴定了不能类型的circRNAs, 大部分的类型是annot_exon，接着是exon_intron（Fig. 1B）。

2 不同脑区域不同时间circRNA表达谱

    对树鼩脑不同区域的不同时间的circRNAs进行无监督聚类分析和PCA分析以研究样品之间的相关性和时空动态变化。不同区域间circRNA的表达谱相对独立(Fig. S1)。箱形图用于比较样品间的表达量分布，以分析circRNAs在不同时间段的变化。如Fig. 2A，在不同脑区域，不同时间段，circRNAs整体上没有显著的改变。能在海马体检测到的大部分circRNAs也能在小脑鉴定到，反之亦然(Fig. 2B)。在海马体中独立表达的有983个circRNAs，在小脑独立表达的有1009个。在不同时间段独立表达的circRNAs比较少(Fig. 2C)。

3 随着年龄的增长，小脑区circRNAs表达量倾向于下调

    我们做了circRNAs表达谱分析以确定在海马体和小脑发育中是否有circRNAs差异表达。差异表达的计算基于back-spliced reads per million (RPM)。只有fold change (FC)≥ 2 and p-value < 0.05的认为是有差异的。我们比较了每个区域幼儿期（I）与青年期（Y），青年期（Y）与老年期（O），幼儿期（I）与老年期（O）。差异表达的circRNA在1101-3234不等(Fig. 3A)。火山图显示小脑区随着年龄增长circRNAs倾向于下调 (Fig. 3B-D)。而在海马体这种下调的趋势仅表现在I与Y比较(Fig. S2A-C)。小脑区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比较，上下调的数量几乎一致(Fig. S2C)。这些数据说明circRNAs的表达具有时空特异性。

4 脑不同时间circRNA的动态变化

    为了鉴定circRNA表达谱的动态变化，我们做了趋势分析以鉴定脑区域不同时间段的主要差异表达的circRNAs。在海马体和小脑中鉴定到7中主要的趋势(Fig. 4A-B)。海马体中4个趋势(profiles0, 1, 6, and 7)有显著意义，小脑3个趋势(profiles 0, 1, and 6)有显著意义。我们也分析了海马体和小脑区上调和下调的circRNAs。我们发现海马体有385个circRNAs上调，而小脑有197个；海马体有334个circRNAs下调，而小脑有439个(Fig. 4C)。这些发现说明年龄相关异常表达的circRNAs可能在脑发育和衰老中发挥重要作用。

5 差异表达的circRNAs富集分析

    为了推测脑不同区域上下调circRNAs发挥的作用，我们对profile 0（下调circRNAs）, profile7（上调circRNAs）进行GO，KEGG富集分析。对于profile 0，GO分析显示circRNAs在不同脑区域富集到相似的功能(Fig. S3A)。对于profile 7，GO富集到的条目差异很大(Fig. S3B)。在profile 0中海马体区富集到14个通路满足p < 0.05 (Tables S4-5)，而小脑区富集到30个通路。对top20个pathways绘制散点图(Fig. 5A-B)。在profile 7（上调circRNAs）中，海马体和小脑共鉴定到7和24个通路((p < 0.05，Tables S6-7)。柱状图呈现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。两个区域共有的有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与衰老相关。而且两个区域大部分的通路可以归类到environmental information processing, genetic information processing, and organismal systems (Tables S4-7)。

6 CircRNA-miRNA-mRNA共表达

    CircRNAs可以海绵吸附microRNAs进而调控基因表达。因此我们基于UMP通路和long-term depression构建共表达网络。在小脑，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)。基因UBE4B可以被这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个与小脑区一致）和5个circRNAs（其中两个与小脑区一样）调控，如Fig. 6B。这个复杂的ceRNA网络说明novel_circRNA_007362 and UBE4B可能通过miRNAs吸附结合作用靶基因参与UMP pathway调控脑发育和衰老。

7 树鼩circRNAs验证

    为了验证RNAseq鉴定到的circRNAs，我们自由选择了8个海马体circRNAs，6个小脑circRNAs。设计合适的引物(Table S2)扩增circRNA成环连接点。与测序结果一致，6个circRNAs (novel_circ_013859, novel_circ_020413, novel_circ_000954, novel_circ_012192, novel_circ_003099, and novel_circ_018752)随着年龄增长下调(Fig. 7A)。海马体中得到RNAseq倍数结果与RT-PCR定量结果一致(Fig.7B–D)。novel_circ_029333在both I vs O and I vs Y中更高。而且PCR产物进一步用琼脂糖凝胶验证。其中两个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在物种与组织的分布上存在特异性。在海马体和小脑区观察到丰富的circRNAs。之前有研究比较过人额皮质，肝，甲状腺和肌肉组织发现，circRNAs在脑中表达最丰富。我们的发现与此一致。Agnieszka et al.报道，circRNA在小鼠小脑区比其他脑部地区表达更丰富，而在小脑区表达的circRNAs与线性RNA的比例也很高。但在本研究，我们发现海马体与小脑circRNAs的整体表达没有明显的差异(Fig. 2A)。我们测序时去除了线性RNA，没有获得线性RNA的数据，不了解circRNAs与线性RNA表达的比例。

    我们发现小脑下调的circRNAs倾向于年老的组织。而在海马体，下调的circRNAs仅倾向于在I 与Y比较时。很明显在树鼩脑衰老时，有些circRNAs表达量下调，而有些表达量上调。不清楚为何会有这种现象。年老小鼠中，Gruner et al.发现一个明显的趋势，circRNAs在皮质和海马器上调，而心脏不是。这说明在树鼩中，骨骼肌相关的circRNAs并没有随着年龄改变，尽管19个circRNAs随着年龄增长下调，可能是由于组织或物种特异性。接下来，我们用qPCR鉴定了时间尺度上的circRNAs表达量，结果与RNAseq结果一致，因此可以确信我们的RNAseq结果是可靠地。但并清楚这些一致circRNAs的靶标为何依然有变化。

    为了得到与脑发育和衰老相关的circRNAs，我们比较了小脑区和海马体I, Y, and O三个时期的circRNA表达量。依据GO，KEGG分析，我们关注上调和下调的circRNAs。在profile 0（下调circRNAs）,GO富集分析显示，在这两个区域富集到相似的功能(Fig. S3A)。而且，UMP pathway, MAPK signaling pathway, phosphatidylinositol signaling system, and long-term depression在这两个区域均有鉴定到。其他通路如RNA degradation, the ErbB signaling pathway, and PI3K-Akt signaling pathway也伴随上述提到的4个通路出现(Fig. 4C–D)。在profile 7（上调circRNAs）,GO分析显示两个区域富集到的功能差异很大(Fig. S3B)。KEGG分析显示rap1 signaling pathway and long-term potentiation在两个区域均有鉴定到。其他通路如Fanconi anemia pathway, cAMP signaling pathway, and insulin signaling pathway也伴随着上述两个通路出现(Fig.S3C)。这些发现显示涉及到的通路与脑发育和衰老相关。之前研究报道UMP pathway在山羊初级毛囊，次级毛囊发育时发挥重要作用，UBE2O在控制初级毛囊，次级毛囊发育相关外显子表达时有重要作用。Ham et al.发现UMP涉及到正常脑与病理脑衰老。UMP可以降解有危害的聚集的蛋白增长寿命。另外，我们对明显参与UMP and long-term depression的circRNA-miRNA-mRNA构建ceRNA网络。我们发现novel_circRNA_007362可以结合24个miRNA进而影响UBE4B的表达。而且，novel_circRNA_032232通过结合9 miRNAs，影响BIRC6表达(Fig. 6B)，并且间接的影响UBE4B and HERC1的表达。根据KEGG分析，UBE4B, BIRC6, and HERC1是UMP通路中关键的基因。（此处讲了HERC1，BIRC6泛素化相关功能，此处省略）。ceRNA网络显示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调控小脑和海马体发育与衰老。总之大脑发育与衰老的调控是一个复杂的多种方面的，许多细节还未知。

    分析circRNAs同源性有助于我们理解共同的规律。之前研究显示circRNAs在人和小鼠脑中是高度保守的。比如，Agnieszka et al.发现15,849中有4,522小鼠circRNAs是与人一致的。Guo et al. 发现20%的小鼠circRNAs与人类是直系同源。而在树鼩中，约83.1%的circRNA与人类circRNA有同源关系。或许是树鼩与灵长类动物关系更近。而且，人，树鼩中共有的保守的circRNAs的功能和产生机制更相似。我们在树鼩海马体和小脑中鉴定到大量的circRNAs。我们的结果可以为鉴定影响脑发育，正常衰老，突触再生，重塑，Alzheimer's disease, and Parkinson’s disease相关的circRNAs提供一个借鉴。未来需要将更多研究放在与脑发育，衰老相关疾病circRNAs的功能上。

5 研究思路

材料：

    树鼩3个发育时期（I, Y, and O）海马体和3个发育时期（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 数据库进行GO富集分析，KOBAS KEGG富集分析。

    6 构建共表达网络，取与UMP和long-term depression相关的基因，miRNA, circRNA构建共表达网络。

    7 circRNA确认，qPCR, RT-PCR凝胶电泳，Sanger测序。

结果展示：

1 circRNA描述

    CircRNAseq产出序列信息（table S1），circRNA长度分析线图（Fig. 1A），circRNA类型柱状图（Fig. 1B）。

2 circRNA表达谱

    CircRNA分组PCA图（Fig. S1），箱形图比较不同样品circRNA整体表达量（Fig. 2A），样品circRNA鉴定数量venn图（Fig. 2B，Fig. 2C）。

3 不同发育时期circRNA比较

    各个分组分别比较差异基因数量统计（Fig. 3A），火山图（Fig. 3，Fig. S2）。

4 不同发育时期趋势分析

    趋势分析折线图（Fig. 4）。

5 差异表达的circRNAs功能富集分析

    GO柱状图（Fig. S3），通路散点图（Fig. 5），通路富集列表（Tables S4-7）。

6 CircRNA-miRNA-mRNA共表达

    共表达网络图(Fig. 6)。

7 树鼩circRNAs验证

    qPCR定量柱状图（Fig. 7A-D），电泳图（Fig. 7E）,Sanger测序峰图（Fig. 7F–G）。

8 同源分析

    Blastn同源分析结果列表（Table S8），不同物种同源结果venn图（Fig. 8）。

微信扫描二维码关注微信公众号了解更多文献信息