Research Paper Volume 10, Issue 4 pp 833—852
CircRNAs in the tree shrew (Tupaia belangeri) brain during postnatal development and aging
- 1 Center of Tree Shrew Germplasm Resources, Institute of Medical Biology, Chinese Academy of Medical Science and Peking Union Medical College, Kunming, China
- 2 Yunnan Key Laboratory of Vaccine Research and Development on Severe Infectious Diseases, Kunming, China
- 3 Yunnan Innovation Team of Standardization and Application Research in Tree Shrew, Kunming, China
Received: February 9, 2018 Accepted: April 23, 2018 Published: April 30, 2018
https://doi.org/10.18632/aging.101437How to Cite
Abstract
Circular RNAs (circRNAs) are a novel type of non-coding RNA expressed across different species and tissues. At present, little is known about the expression and function of circRNAs in the tree shrew brain. In this study, we used RNA-seq to identify 35,007 circRNAs in hippocampus and cerebellum samples from infant (aged 47–52 days), young (aged 15–18 months), and old (aged 78–86 months) tree shrews. We observed no significant changes in the total circRNA expression profiles in different brain regions over time. However, circRNA tended to be downregulated in the cerebellum over time. Real-time RT-PCR analysis verified the presence of circRNAs. KEGG analysis indicated the occurrence of ubiquitin-mediated proteolysis, the MAPK signaling pathway, phosphatidylinositol signaling system, long-term depression, the rap1 signaling pathway, and long-term potentiation in both brain regions. We also observed that 29,087 (83.1%) tree shrew circRNAs shared homology with human circRNAs. The competing endogenous RNA networks suggested novel_circRNA_007362 potential functions as a 24-miRNAs sponge to regulate UBE4B expression. Thus, we obtained comprehensive circRNA expression profiles in the tree shrew brain during postnatal development and aging, which might help to elucidate the functions of circRNAs during brain aging and in age-related diseases.
Introduction
Circular RNAs (circRNAs) are a novel type of non-coding RNA with a covalently closed-loop structure generated by back splicing. CircRNAs were first identified approximately 37 years ago [1], but they have only recently attracted the attention of researchers. CircRNAs are considerably more stable and more resistant to RNase R than linear RNAs [2]. However, the functions of most circRNAs remain unclear. Previous studies have shown circRNAs can act as “sponges” for microRNAs via the competing endogenous RNA (ceRNA) network in order to regulate target gene expression [3,4]. CircRNAs may also serve as sponges for RNA-binding proteins to post-transcriptionally modulate mRNA expression [4].
Recent studies have suggested that circRNAs can have critical effects on the modulation of cellular processes, such as apoptosis, proliferation, development, cancer, and aging [5–7]. Studies have also shown that circRNAs accumulate in the brain during aging, and thus they constitute a novel biomarker for aging [8]. For example, aged mice exhibit increased expression of circRNAs in the cortex and hippocampus [9]. It has been reported that FOXO3, a member of the forkhead family of transcription factors, plays an important role in cellular senescence and aging, and it is associated with hippocampal neuronal injury [5,10,11]. Furthermore, the circ-Foxo3 (the circRNA generated from Foxo3) is present in aged patients and murine cardiovascular disease models [5].
Many circRNAs have been detected in different tissues in humans, rhesus macaques, mice, sheep, pigs, and zebrafish [2,12–15]. It has also been reported that circRNAs are more enriched in neuronal tissues compared with other tissues [2,8,16,17]. In mice, researchers have identified 1,212 circRNAs in the hippocampus and 2,407 circRNAs in the cerebellum [2]. In the human brain, 65,731 circRNAs have been identified [2]. Ninety percent of mouse circRNAs and 73% of human circRNAs remain unannotated [18].
The tree shrew (Tupaia belangeri) is a novel and valuable laboratory animal due to its inexpensive maintenance, small body size, short reproductive cycle, and high brain-to body mass ratio [19,20]. This model has been used to study human diseases, such as herpes simplex viruses [21], hepatitis C viruses [22], and Coxsackie virus A16 [23], especially in the nervous system, as well as Alzheimer's disease [24], brain development and aging [25,26], social stress and depression [27,28]. Whole genome sequencing of the tree shrew has been completed, but only a limited number of transcriptome profiles are available for of various tissues in the tree shrew. Furthermore, very little information is available regarding circRNAs in the tree shrew.
In this study, we investigated the effects of age on the circRNA expression profiles in the tree shrew hippocampus and cerebellum using high-throughput sequencing. We also constructed the ceRNA networks based on target analysis. Numerous previously unreported circRNAs were identified in the tree shrew brain with increasing age, which may improve our understanding of the roles of circRNAs during brain aging and in age-related brain diseases.
Results
Overview of circRNA-Seq
In total, 465,791,754 clean reads (235,543,658 for hippocampus and 230,248,096 for cerebellum) were generated. After removing the low quality, poly-N-containing, and adapter-containing reads from the clean reads, 433,887,754 high quality clean reads (216,130,360 from the hippocampus and 217,757,394 from the cerebellum) were obtained from the two tissues (Table S1). The high quality clean reads from which the mapped rRNA reads had been removed were mapped to the tree shrew genome using TopHat2, and the unmapped reads were selected. Based on theoretical predictions, 35,007 circRNAs were detected using find_circ. All of the identified circRNAs were novel. As shown in Fig. 1A, most of the circRNAs were approximately 400 nt in length, as found in previous studies [29]. In addition, different types of circRNAs were identified (summarized in Fig. 1B), where the most common type was annot_exon followed by exon_intron.
Figure 1. Characterization of circRNAs in the tree shrew brain. (A) Length distribution of circRNAs. (B) Type distribution of circRNAs. Red indicates annot_exon, brown represents antisense, green indicates exon_intron, light blue shows intergenic, blue represents intronic, and pink indicates one_exon.
CircRNA expression profiles in different brain regions over time
Unsupervised hierarchical clustering and principal component analysis (PCA) were performed to understand the relationships between circRNA expression and spatiotemporal dynamics in the tree shrew brain. The circRNA expression patterns in the two brain regions comprised relatively independent clusters, and the first two principal coordinates corresponded to the brain region and age (Fig. S1). A box view was used to compare the expression distributions in the samples in order characterize the changes in circRNA expression over time. As shown in Fig. 2A, there were no significant changes in total circRNA expression in different brain regions over time. The majority of the circRNAs annotated in the hippocampus were also detected in the cerebellum, and vice versa (Fig. 2B). In total, 983 circRNAs were expressed exclusively in the hippocampus and 1009 were expressed exclusively in the cerebellum. A small number of circRNAs were expressed exclusively in each age group (Fig. 2C).
Figure 2. Characteristics of circRNA expression. (A) Box plots comparing circRNAs measured based on the log of the expression level (base 2) in the two brain regions over time. CER: cerebellum, HIP: hippocampus. (B) Venn diagram showing the overlapping circRNAs identified in the hippocampus and cerebellum. (C) Venn diagram showing the overlapping circRNAs detected in different brain regions over time. Different colors and shapes represent different brain regions and time points. HIP-I, HIP-Y, and HIP-O represent the three ages in the hippocampus. CER-I, CER-Y, and CER-O represent the three ages in the cerebellum.
CircRNAs expressed in the cerebellum tended to be downregulated with increasing age
We performed expression profiling to determine whether circRNAs were differentially expressed in the hippocampus and cerebellum during development. Differences in circRNA expression between samples were calculated based on the back-spliced reads per million (RPM) mapped reads method. Only differentially expressed circRNAs with a fold change (FC) ≥ 2 and p-value < 0.05 were considered. We compared infants (I) versus young (Y) animals, Y versus old (O) animals, and I versus O for each brain region. The numbers of differentially expressed circRNAs were high in the group comparison, where they ranged from 1,101‒3,234 with increasing age in the different brain regions (Fig. 3A). Volcano plots indicated that circRNAs tended to be downregulated with age in cerebellum samples (Fig. 3B–D). However, in the hippocampus, the tendency to be downregulated with increasing age was only found when I was compared with Y (Fig. S2A–C). In the cerebellum, we found that 1,339 (4.1%) circRNAs were significantly downregulated and 935 (2.9%) were significantly upregulated in I compared with O animals (Fig. 3D). In I vs Y, 1,272 (3.9%) were significantly downregulated and 1,142 (3.5%) were significantly upregulated (Fig. 3B). In Y vs O, 970 (2.9%) were significantly downregulated, and 739 (2.2%) were significantly upregulated (Fig. 3C).
Figure 3. Differential circRNA expression in the cerebellum with aging. (A) Bar plot representation of differentially expressed circRNAs in different age groups. (B–G) Volcano plots showing the dysregulated circRNAs in the cerebellum in different age groups. The log of FC (base 2) is plotted on the x-axis and the negative log of FDR (base 10) is plotted on the y-axis (FC >2 and p < 0.05). The green and red points on the graph represent downregulated and upregulated circRNAs, respectively.
In the hippocampus, we observed that 1,636 (5.0%) circRNAs were significantly downregulated and 1,386 (4.3%) were significantly upregulated in I compared with Y animals (Fig. S2A). However, in Y vs O, 512 (1.7%) were significantly downregulated, and 588 (1.9%) were significantly upregulated (Fig. S2B). In I vs O, the numbers of significantly downregulated and upregulated circRNAs were almost unanimously (Fig. S2C). These results indicated that circRNA expression could be spatially and temporally specific.
Dynamics of circRNA expression in the brain over time
To characterize the changes in circRNA expression, we performed trend analysis to identify the predominant differentially expressed circRNAs in each brain region over the course of development. We identified seven main circRNA profiles in the hippocampus and cerebellum, where each represented a characteristic expression pattern (Fig. 4A–B). Four patterns (profiles 0, 1, 6, and 7) were significant in the hippocampus and three patterns were significant in the cerebellum (profiles 0, 1, and 6). Among these patterns, profile 0 was significantly downregulated in both the hippocampus and cerebellum. We also analyzed the up- and downregulated circRNAs in both the hippocampus and cerebellum. We found that 385 circRNAs were increased in the hippocampus and 197 circRNAs in the cerebellum, whereas 334 circRNAs were decreased in the hippocampus and 439 circRNAs in the cerebellum (Fig. 4C). These findings suggested that the age-associated dysregulation of circRNAs in the tree shrew brain may play roles in brain development and aging.
Figure 4. Overall expression patterns of circRNAs in different regions of the brain. (A) Expression patterns of circRNAs in the hippocampus. (B) Expression patterns of circRNAs in the cerebellum. The colored profile represents significant enrichment (p < 0.05). The uncolored profiles denote no significant enrichment. The profile ID and the number of genes in the profile are shown above the chart. (C–D) Venn diagram showing the up- and downregulated circRNAs in both the hippocampus and cerebellum. C: upregulated circRNAs. D: downregulated circRNAs.
Enrichment of differentially expressed circRNAs
To elucidate the possible roles of the up- and downregulated circRNAs in different brain regions, enriched Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were obtained for the circRNAs in profile 0 (downregulated circRNAs) and profile 7 (upregulated circRNAs). For profile 0, GO analyses showed that the circRNAs in the two brain regions were enriched in terms of similar functions (Fig. S3A). However, for profile 7, GO analyses showed that the circRNAs in the two brain regions were enriched in terms of quite divergent functions (Fig. S3B). In profile 0, 14 pathways in the hippocampus and 30 pathways in the cerebellum satisfied the requirement of p < 0.05 (Table S4, Table S5). We used bubble charts to represent the top 20 pathways in profile 0 in each brain region (Fig. 5A–B). In profile 7 (upregulated circRNAs), seven and 24 KEGG pathways (p < 0.05) were identified in the hippocampus and cerebellum, respectively (Table S6, Table S7). We used bar plots to represent the KEGG pathways in profile 7 for each brain region (Fig. S4). In profile 0, RNA degradation, the neurotrophin signaling pathway, and synaptic vesicle cycle were identified in the hippocampus. The phosphatidylinositol signaling system, phospholipase D signaling pathway, PI3K-Akt signaling pathway, and ErbB signaling pathway were identified in the cerebellum. The ubiquitin-mediated proteolysis (UMP) pathway, MAPK signaling pathway, phosphatidylinositol signaling system, and long-term depression were identified in both brain regions (Fig. 5A–B). In profile 7, the synaptic vesicle cycle, UMP, and the Fanconi anemia pathway were identified in the hippocampus. The cAMP signaling pathway, phosphatidylinositol signaling system, and insulin signaling pathway were identified in the cerebellum. The rap1 signaling pathway and long-term potentiation were identified in both brain regions (Fig. S3C). These results suggested that the circRNAs involved in these pathways may be related to aging in the hippocampus and cerebellum of the tree shrew. Furthermore, most of the pathways in the two brain regions were classified into environmental information processing, genetic information processing, and organismal systems (Tables S4–7).
Figure 5. KEGG analysis of profile 0 in the hippocampus and cerebellum. The bubble chart shows enriched differentially expressed genes in signaling pathways. (A) Bubble chart of the top 20 pathways in the hippocampus. (B) Bubble chart of the top 20 pathways in the cerebellum. The Y-axis label represents the pathway and the X-axis label represents the rich factor (rich factor = amount of differentially expressed genes enriched in the pathway/amount of all genes in background gene set). The color and size of the bubble represent enrichment significance and the amount of differentially expressed genes enriched in the pathway, respectively.
CircRNA-miRNA-mRNA co-expression
CircRNAs can serve as “sponges” for microRNAs to regulate host gene expression, so we constructed a co-expression network based on the UMP pathway and long-term depression. In the cerebellum, novel_circRNA_007362 was predicted to combine with 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, tch-miR-135b-5p, tch-miR-15b-5p, tch-miR-195-5p, tch-miR-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 was predicted to be targeted by the same 24 miRNAs and six circRNAs (novel_circRNA_007356, novel_circRNA_007357, novel_circRNA_007362, novel_circRNA_007364, novel_circRNA_007373, and novel_circRNA_007379) (Fig. 6A). In the hippocampus, UBE4B was predicted to be targeted by seven miRNAs (six of these miRNAs were the same as those in the cerebellum) and five circRNAs (two circRNAs were the same as those in the cerebellum), as shown in Fig. 6B. This complicated ceRNA network suggested that novel_circRNA_007362 and UBE4B might play regulatory roles in the UMP pathway via the miRNAs mentioned above and their targets during brain development and aging.
Figure 6. Competing endogenous RNA network based on circRNA–miRNA–mRNA interactions. (A) Network in the cerebellum. (B) Network in the hippocampus. Red circles represent miRNAs. Blackish green circles represent genes. Pink blue circles represent circRNAs.
Validation of tree shrew circRNAs
To validate the circRNAs generated by RNA-seq, we randomly selected eight circRNAs from the hippocampus and six circRNAs from the cerebellum. We then designed appropriate primers (Table S2) and amplified the circular junctions. Similar to our RNA-seq data, the expression levels of the six circRNAs (novel_circ_013859, novel_circ_020413, novel_circ_000954, novel_circ_012192, novel_circ_003099, and novel_circ_018752) decreased with age (Fig. 7A). The FCs detected in the hippocampus by RNA-seq were in agreement with the RT-PCR quantification results (Fig. 7B–D). The expression of novel_circ_029333 was higher in both I vs O and I vs Y. In addition, the PCR products were analyzed using agarose gels to validate that a single DNA was amplified (Fig. 7E). In particular, two of the selected circRNAs, i.e., novel_circ_013859 and novel_circ_007625, were further confirmed by Sanger sequencing (Fig. 7F–G), where the real-time RT-PCR results were in agreement with the RNA-seq data, thereby demonstrating that the identified circRNAs had different expression levels in vivo.
Figure 7. RT-PCR and Sanger sequencing validation of the selected circRNAs. (A) Six circRNAs from the cerebellum were quantified by RT-qPCR analysis. (B–D) Expression level changes for eight circRNAs from the hippocampus using RT-qPCR between I vs O, I vs Y, and Y vs O. (E) RT-qPCR products were visualized using electrophoresis. (F–G) Two representative examples of qPCR products confirmed by Sanger sequencing. F: novel_circ_013859. G: novel_circ_007625.
Homology of tree shrew circRNAs
We compared tree shrew circRNAs with human, mouse, and rhesus macaque circRNAs. The results indicated that 29,087 (83.1%) tree shrew circRNAs shared homology with human circRNAs (Table S8), 3,783 (10.8%) with mouse circRNAs, and 392 (1.1%) with rhesus macaque circRNAs. Only 56 tree shrew circRNAs (0.2%) shared homology with all three species and 5,695 (16.3%) of the tree shrew circRNAs did not share homology with any of the other three species (Fig. 8).
Figure 8. Venn diagram representation showing the homology among four species, i.e., human, mouse, rhesus macaque, and tree shrew. The alignment threshold was e-value < 1e-10 by BLASTN (Table S8). Rhesus: rhesus macaque, has: human, mmu: mouse, Tupaia_Belangeri: tree shrew.
Discussion
In the present study, we identified 35,007 circRNAs in the tree shrew brain using high-throughput sequencing. We found that most of the circRNAs measured approximately 400 nt in length, as also reported in previous studies [29]. In the tree shrew, the major circRNAs were annot_exon, whereas in the sheep pituitary, the major circRNAs were intergenic [13]. This finding may suggest that the distribution of circRNAs can exhibit species and tissue specificity.
We observed that abundant circRNAs were expressed in the tree shrew hippocampus and cerebellum. A previous study that compared circRNA expression in the human frontal cortex, liver, thyroid gland, and muscle found that circRNAs were most abundant in the brain [2]. Our results were consistent with these findings. Agnieszka et al. [2] also reported enriched circRNA expression in the mouse cerebellum compared with other brain regions, where the ratio of circRNA expression to linear RNA expression was also higher in the cerebellum. However, in our study, we found no significant difference in total circRNA expression levels in the hippocampus and cerebellum (Fig. 2A). We removed the linear RNA when sequencing, so there were no data concerning the ratio of circRNA relative to linear RNA in the tree shrew.
We found that the downregulation of circRNAs tended to increase with age in the cerebellum. However, in the hippocampus, the downregulation of circRNAs only tended to occur in I versus Y. Clearly, some circRNAs were downregulated whereas others were upregulated with increasing age in the tree shrew brain. It is unclear why this phenomenon occurred. In aged mice, Gruner et al. [30] found a striking trend where circRNAs were upregulated in the cortex and hippocampus, but not in the heart. It has also been shown that in the rhesus macaque, the abundance of most skeletal muscle circRNAs did not change with increasing age, although 19 circRNAs were downregulated with age [12], possibly due to tissue or species specificity. Furthermore, we employed reverse transcription followed by real-time quantitative PCR (qPCR) to identify the expression levels in the brain over time, and the qPCR results were in agreement with the RNA-sequencing, thereby confirming that our RNA-sequencing results were high quality. However, it was unclear whether the corresponding targets of these circRNAs also changed.
To determine the circRNAs related to brain development and aging, we compared the circRNA expression levels in the cerebellum and hippocampus of I, Y, and O animals. We focused on the upregulated and downregulated circRNAs in the hippocampus and cerebellum based on GO and KEGG analyses. In profile 0 (downregulated circRNAs), GO analyses showed that circRNAs in the two brain region were enriched in terms of similar functions (Fig. S3A). Furthermore, the UMP pathway, MAPK signaling pathway, phosphatidylinositol signaling system, and long-term depression were identified in the hippocampus and cerebellum. Other pathways such as RNA degradation, the ErbB signaling pathway, and PI3K-Akt signaling pathway were also incorporated in the four pathways mentioned above (Fig. 4C–D). In profile 7 (upregulated circRNAs), GO analyses showed that circRNAs in the two brain regions were enriched in terms of quite divergent functions (Fig. S3B), and KEGG analyses identified the rap1 signaling pathway and long-term potentiation in both brain regions. Other pathways such as the Fanconi anemia pathway, cAMP signaling pathway, and insulin signaling pathway were also incorporated in the two pathways mentioned above (Fig. S3C). These findings demonstrated the involvement of pathways related to brain development and aging. Previous studies have shown that the UMP pathway was important for primary and secondary hair follicle development in cashmere goats [31], and that UBE2O played an important role in controlling development-associated exon expression in the primary and secondary hair follicles [31]. Ham et al. also found that UMP was involved with normal and pathological brain aging [32]. UMP can degrade damaged and aggregated proteins that increase with age [33]. Furthermore, we constructed ceRNA networks comprising circRNA-miRNA-mRNA that significantly participated in UMP and long-term depression. We found that a ceRNA, novel_circRNA_007362, competed to bind 24 miRNAs, thereby affecting UBE4B expression. In addition, novel_circRNA_032232 affected BIRC6 expression by competing to bind nine miRNAs (Fig. 6B), and it indirectly influenced the expression levels of UBE4B and HERC1. According to KEGG analysis, UBE4B, BIRC6, and HERC1 are important elements of UMP. UBE4B is an E3/E4 ligase that triggers substrate polyubiquitination and it regulates neuronal death and survival [34]. HERC1 (HECT domain and RCC1 domain) is an E3 ubiquitin ligase protein and previous studies have shown that the HERC1 E3 ubiquitin ligase-tbl mutation affects neurons in the cerebral neocortex, the CA3 region of the hippocampus, and the ventral horn of the spinal cord [35]. BIRC6 (Baculoviral IAP Repeat Containing 6), a member of the inhibitors of apoptosis proteins (IAPs) family, encodes an inhibitor of apoptosis and a chimeric E2/E3 ubiquitin ligase in mammals, which regulates p53 and the mitochondrial pathway of apoptosis and it is essential for mouse embryonic development [36,37]. In summary, the ceRNA network suggested that novel_circRNA_007362 and novel_circRNA_032232 might interact with 24 miRNAs (Fig. 6A) as well as tch-let-7a-3p, tch-miR-136-5p, tch-miR-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 (Fig. 6B) to regulate development and aging of the cerebellum and hippocampus. In general, the regulation of brain development and aging is complex and multifaceted, and many of the details remain unknown.
Analyzing the homology of circRNAs can help to understand common characteristics. Previous studies indicated that circRNAs were highly conserved in the mouse and human brains. For example, Agnieszka et al. [2] showed that 4,522 out of 15,849 mouse circRNAs were conserved in humans. Guo et al. found that 20% of mouse circRNAs were orthologous with human circRNAs [38]. However, in the tree shrew, 83.1% of the circRNAs shared homology with human circRNAs, probably because tree shrews are more closely related to primates than rodents. Furthermore, the conserved expression of circRNAs in humans and tree shrews suggests that their functions and mechanisms of biogenesis might also be similar.
In summary, we identified abundantly expressed circRNAs in the tree shrew hippocampus and cerebellum. Our results may provide a valuable platform for identifying circRNAs that influence brain development, normal aging, synapse regeneration and remodeling, Alzheimer's disease, and Parkinson’s disease. Further studies should focus on the functions of these circRNAs in brain development and age-related disease.
Materials and Methods
Animals and tissue collection
All of the tree shrews were the first filial generation, and they were raised at the Institute of Medical Biology, Chinese Academy of Medical Science and Peking Union Medical College in Kunming, China. The tree shrews were healthy and consistent with the group standards for tree shrews (T/CALAS 08-2017 and T/CALAS 09-2017), without visible signs of tumorigenesis or disease. Nine male tree shrews in three different stages of development (I, Y, and O) were used. The I group comprised three tree shrews aged 47–52 days and weighing 100–110 g, the Y group comprised three animals aged 15–18 months and weighing 120–150 g, and the O group comprised three animals aged 78–86 months and weighing 120–150 g (Table S3). The life span of tree shrews in captivity is up to 10 years, and their lifespan in approximately one-eighth of the human lifespan [26,39]. Thus, the ages of our youngest and oldest animals corresponded to human ages of about 1 and 57 years, respectively. Tree shrews were sacrificed and the hippocampus and cerebellum were dissected immediately, and stored at -80°C until use. All of the animal procedures were approved by the ethical committee for Animal Research in the Institute of Medical Biology, Chinese Academy of Medical Science and Peking Union Medical College.
RNA extraction and RNA-Seq
Total RNA was isolated from each sample using TRIzol reagent (Invitrogen) according to the manufacturer’s protocol. The quality of the RNA was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies).
For the circRNA-seq, after isolating the total RNA, we used RNase R to degrade the linear RNA, and purification was performed using the RNeasy MinElute Cleanup Kit (Qiagen). Next, we removed the ribosomal RNA. We then added fragmentation buffer to fragment the circRNA into short pieces, and random primers were used to synthesize the first-strand cDNA. Subsequently, the second-strand cDNA was synthesized using dNTPs, RNase H, DNA polymerase I, and buffer. The library fragments were purified with VAHTSTM DNA Clean Beads, before repairing the ends, adding poly-(A), and ligating to Illumina sequencing adapters. Next, we used uracil-N-glycosylase to digest the second-strand cDNA. The digested products were purified with VAHTSTM DNA Clean Beads and PCR amplified, thereby completing the preparation of the entire library. CircRNA sequencing was performed using the Illumina HiSeqTM2500 platform by Gene Denovo Biotechnology Co. (Guangzhou, China).
CircRNA sequencing analysis
Using the raw FASTQ files, we removed adapter contamination, reads containing more than 10% unknown nucleotides, and low quality reads containing more than 50% low quality (Q-value ≤ 20) bases. Next, we used the short reads alignment tool Bowtie2 to map reads to the rRNA database and removed the rRNA mapped reads. The sequences were then aligned with the tree shrew genome using TopHat2 [40] (version 2.0.3.12). Mapped reads were discarded and unmapped reads were collected for circRNA identification. The two ends of unmapped reads were intercepted (default = 20 bp) and aligned with the tree shrew genome to find unique anchor positions within a splice site. Anchor reads that aligned in the reverse orientation (head-to-tail) indicated that circRNA splicing had occurred. The anchor reads were mapped to the tree shrew genome again and the results were submitted to find_circ [3] software to identify circRNAs.
The identified circRNAs were subjected to statistical analysis based on their type, chromosome distribution, and length distribution. To quantify circRNAs, back-spliced junction reads were scaled to RPM mapped reads. The RPM method can eliminate the effects of different amounts of sequencing data on the calculation of circRNA expression levels. Therefore, the calculated expression levels could be used directly to compare differential expression among samples. circRNAs were compared with circBase [18] using BLAST for annotation. The circRNAs that could not be annotated were defined as novel circRNAs. Box plots and Venn diagrams were generated using online tools (www.omicshare.com/tools/Home/Soft/box and www.omicshare.com/tools/Home/Soft/venn).
Unsupervised hierarchical clustering and PCA
Unsupervised hierarchical clustering and PCA were performed to explore the relationships between samples. Unsupervised hierarchical clustering analyses were performed using gmodels in R (https://cran.r-project.org/web/packages/gmodels/index.html) and PCA was conducted using ggplot2 in R (https://cran.r-project.org/web/packages/ggplot2/index.html).
Analysis of differentially expressed circRNAs
The edgeR package (http://www.r-project.org/) was used to identify differentially expressed circRNAs between samples. We identified circRNAs with FC ≥ 2 and p < 0.05 in comparisons of samples as significant differentially expressed circRNAs. Volcano plots were generated using gglpot2 in R.
Trend analysis
In order to examine the patterns of differentially expressed circRNAs, the expression data for each sample were normalized to 0, log2 (v1/v0), and log2 (v2/v0), and then clustered using Short Time-series Expression Miner (STEM) software [40]. The clustered profiles with p ≤ 0.05 were considered to be significant profiles.
GO and KEGG pathway analyses
GO annotations are useful for predicting the functions of gene products across numerous species [41]. The GO categories were derived from the GO database (http://www.geneontology.org), which is a community-based bioinformatics resource for classifying gene product functions using structured, controlled vocabularies, and it includes three domains: Molecular Function, Biological Processes, and Cellular Components [42,43]. The differentially expressed circRNAs (profile 0 and profile 7) were mapped to GO terms in the GO database, where FDR ≤ 0.05 was taken as the threshold. GO terms that satisfied this condition were defined as significantly enriched GO terms. KEGG is a major public pathway-related database [44]. KEGG analysis can help to understand the biological roles of differentially expressed genes. We selected the differentially expressed circRNAs in profile 0 and profile 7 obtained by trend analysis for GO and KEGG analyses (http://kobas.cbi.pku.edu.cn/home.do), respectively. Pathways with p < 0.05 were considered significant pathways. Bubble charts were generated using gglpot2 in R.
Construction of co-expression network
We identified the interactions among circRNAs, miRNAs, and mRNAs by constructing co-expression networks. We selected genes involved with UMP and long-term depression (p < 0.05, and downregulated in profile 0 in both the hippocampus and cerebellum), and used Cytoscape (V3.2.0) to establish a circRNA–miRNA–mRNA pathway regulation network.
Data access
RNA-seq data were submitted to the NCBI Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/) under accession number SRP 132147.
Homology analysis
The numbers of noncoding circRNAs in the human, mouse, and rhesus macaque genomes that shared homology with tree shrew brain circRNAs were identified based on a threshold (e-value < 1e-10) using BLASTN (Table S8). The human, mouse, and rhesus macaque data were downloaded from: http://www.ncbi.nlm.nih.gov.
Quantitative real-time RT-PCR
To confirm the individual circRNAs, total RNA was extracted from samples using TRIzol reagent (Invitrogen) and reverse transcribed into cDNA with random primers by a RevertAid First Stand cDNA Synthesis Kit (GeneCopoeia, Cat AORT-0020), according to the manufacturer’s instructions. qPCR analysis was performed using the SYBR Green Kit (GeneCopoeia, Cat AORT-0020) with the PikoReal PCR system (Thermo, USA) under the following conditions: 95°C for 10 min, followed by 45 cycles at 95°C for 15 s and 60°C for 60 s. All of the primers used are listed in Table S2. Relative expression levels were calculated using the 2–△△Ct method. GAPDH was used as a reference. PCR products were tested by electrophoresis on 2% agarose gels, then gel extracted and Sanger sequenced.
Acknowledgements
We thank Dr Feng Yue from the Yunnan Provincial Center for Molecular Medicine, Faculty of Life Science and Technology, Kunming University of Science and Technology, Kunming, China for his useful suggestions.
Conflicts of Interest
The authors declare no conflicts of interest.
Funding
This study was funded by research grants from Yunnan Science and Technology Talent and Platform Program (2017HC019), the National Natural Science Foundation of China (31601907), the Yunnan Province Major Science and Technology Project (2017ZF007), and the Key Laboratory of Yunnan Province for Ophthalmic Research and Disease Control (2017DG008).
References
- 1. Hsu MT, Coca-Prados M. Electron microscopic evidence for the circular form of RNA in the cytoplasm of eukaryotic cells. Nature. 1979; 280:339–40. https://doi.org/10.1038/280339a0 [PubMed]
- 2. Rybak-Wolf A, Stottmeister C, Glažar P, Jens M, Pino N, Giusti S, Hanan M, Behm M, Bartok O, Ashwal-Fluss R, Herzog M, Schreyer L, Papavasileiou P, et al. Circular RNAs in the Mammalian Brain Are Highly Abundant, Conserved, and Dynamically Expressed. Mol Cell. 2015; 58:870–85. https://doi.org/10.1016/j.molcel.2015.03.027 [PubMed]
- 3. Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M, Loewer A, Ziebold U, Landthaler M, et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature. 2013; 495:333–38. https://doi.org/10.1038/nature11928 [PubMed]
- 4. Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK, Kjems J. Natural RNA circles function as efficient microRNA sponges. Nature. 2013; 495:384–88. https://doi.org/10.1038/nature11993 [PubMed]
- 5. Maiese K. Disease onset and aging in the world of circular RNAs. J Transl Sci. 2016; 2:327–29. https://doi.org/10.15761/JTS.1000158 [PubMed]
- 6. Lukiw WJ. Circular RNA (circRNA) in Alzheimer’s disease (AD). Front Genet. 2013; 4:307. https://doi.org/10.3389/fgene.2013.00307 [PubMed]
- 7. Bachmayr-Heyda A, Reiner AT, Auer K, Sukhbaatar N, Aust S, Bachleitner-Hofmann T, Mesteri I, Grunt TW, Zeillinger R, Pils D. Correlation of circular RNA abundance with proliferation--exemplified with colorectal and ovarian cancer, idiopathic lung fibrosis, and normal human tissues. Sci Rep. 2015; 5:8057. https://doi.org/10.1038/srep08057 [PubMed]
- 8. Westholm JO, Miura P, Olson S, Shenker S, Joseph B, Sanfilippo P, Celniker SE, Graveley BR, Lai EC. Genome-wide analysis of drosophila circular RNAs reveals their structural and sequence properties and age-dependent neural accumulation. Cell Reports. 2014; 9:1966–80. https://doi.org/10.1016/j.celrep.2014.10.062 [PubMed]
- 9. Hannah Gruner MC. Matthew Baue & Pedro Miura. CircRNA accumulation in the aging mouse brain. Sci Rep. 2016; 6:38907. https://doi.org/10.1038/srep38907
- 10. Park SH, Sim YB, Lee JK, Lee JY, Suh HW. Characterization of temporal expressions of FOXO and pFOXO proteins in the hippocampus by kainic acid in mice: involvement of NMDA and non-NMDA receptors. Arch Pharm Res. 2016; 39:660–67. https://doi.org/10.1007/s12272-016-0733-9 [PubMed]
- 11. Wang S, Chong ZZ, Shang YC, Maiese K. WISP1 neuroprotection requires FoxO3a post-translational modulation with autoregulatory control of SIRT1. Curr Neurovasc Res. 2013; 10:54–69. https://doi.org/10.2174/156720213804805945 [PubMed]
- 12. Abdelmohsen K, Panda AC, De S, Grammatikakis I, Kim J, Ding J, Noh JH, Kim KM, Mattison JA, de Cabo R, Gorospe M. Circular RNAs in monkey muscle: age-dependent changes. Aging (Albany NY). 2015; 7:903–10. https://doi.org/10.18632/aging.100834 [PubMed]
- 13. Li C, Li X, Ma Q, Zhang X, Cao Y, Yao Y, You S, Wang D, Quan R, Hou X, Liu Z, Zhan Q, Liu L, et al. Genome-wide analysis of circular RNAs in prenatal and postnatal pituitary glands of sheep. Sci Rep. 2017; 7:16143. https://doi.org/10.1038/s41598-017-16344-y [PubMed]
- 14. Shen Y, Guo X, Wang W. Identification and characterization of circular RNAs in zebrafish. FEBS Lett. 2017; 591:213–20. https://doi.org/10.1002/1873-3468.12500 [PubMed]
- 15. Venø MT, Hansen TB, Venø ST, Clausen BH, Grebing M, Finsen B, Holm IE, Kjems J. Spatio-temporal regulation of circular RNA expression during porcine embryonic brain development. Genome Biol. 2015; 16:245. https://doi.org/10.1186/s13059-015-0801-3 [PubMed]
- 16. Ashwal-Fluss R, Meyer M, Pamudurti NR, Ivanov A, Bartok O, Hanan M, Evantal N, Memczak S, Rajewsky N, Kadener S. circRNA biogenesis competes with pre-mRNA splicing. Mol Cell. 2014; 56:55–66. https://doi.org/10.1016/j.molcel.2014.08.019 [PubMed]
- 17. You X, Vlatkovic I, Babic A, Will T, Epstein I, Tushev G, Akbalik G, Wang M, Glock C, Quedenau C, Wang X, Hou J, Liu H, et al. Neural circular RNAs are derived from synaptic genes and regulated by development and plasticity. Nat Neurosci. 2015; 18:603–10. https://doi.org/10.1038/nn.3975 [PubMed]
- 18. Glažar P, Papavasileiou P, Rajewsky N. circBase: a database for circular RNAs. RNA. 2014; 20:1666–70. https://doi.org/10.1261/rna.043687.113 [PubMed]
- 19. Cao J, Yang EB, Su JJ, Li Y, Chow P. The tree shrews: adjuncts and alternatives to primates as models for biomedical research. J Med Primatol. 2003; 32:123–30. https://doi.org/10.1034/j.1600-0684.2003.00022.x [PubMed]
- 20. Li CH, Yan LZ, Ban WZ, Tu Q, Wu Y, Wang L, Bi R, Ji S, Ma YH, Nie WH, Lv LB, Yao YG, Zhao XD, Zheng P. Long-term propagation of tree shrew spermatogonial stem cells in culture and successful generation of transgenic offspring. Cell Res. 2017; 27:241–52. https://doi.org/10.1038/cr.2016.156 [PubMed]
- 21. Li L, Li Z, Wang E, Yang R, Xiao Y, Han H, Lang F, Li X, Xia Y, Gao F, Li Q, Fraser NW, Zhou J. Herpes Simplex Virus 1 Infection of Tree Shrews Differs from That of Mice in the Severity of Acute Infection and Viral Transcription in the Peripheral Nervous System. J Virol. 2015; 90:790–804. https://doi.org/10.1128/JVI.02258-15 [PubMed]
- 22. Feng Y, Feng YM, Lu C, Han Y, Liu L, Sun X, Dai J, Xia X. Tree shrew, a potential animal model for hepatitis C, supports the infection and replication of HCV in vitro and in vivo. J Gen Virol. 2017; 98:2069–78. https://doi.org/10.1099/jgv.0.000869 [PubMed]
- 23. Li JP, Liao Y, Zhang Y, Wang JJ, Wang LC, Feng K, Li QH, Liu LD. Experimental infection of tree shrews (Tupaia belangeri) with Coxsackie virus A16. Dongwuxue Yanjiu. 2014; 35:485–91. https://doi.org/10.13918/j.issn.2095-137.2014.6.485 [PubMed]
- 24. Fan Y, Luo R, Su LY, Xiang Q, Yu D, Xu L, Chen JQ, Bi R, Wu DD, Zheng P, Yao YG. Does the Genetic Feature of the Chinese Tree Shrew (Tupaia belangeri chinensis) Support Its Potential as a Viable Model for Alzheimer’s Disease Research? J Alzheimers Dis. 2018; 61:1015–28. https://doi.org/10.3233/JAD-170594 [PubMed]
- 25. Wei S, Hua HR, Chen QQ, Zhang Y, Chen F, Li SQ, Li F, Li JL. Dynamic changes in DNA demethylation in the tree shrew (Tupaia belangeri chinensis) brain during postnatal development and aging. Zool Res. 2017; 38:96–102. https://doi.org/10.24272/j.issn.2095-8137.2017.013 [PubMed]
- 26. Keuker JI, de Biurrun G, Luiten PG, Fuchs E. Preservation of hippocampal neuron numbers and hippocampal subfield volumes in behaviorally characterized aged tree shrews. J Comp Neurol. 2004; 468:509–17. https://doi.org/10.1002/cne.10996 [PubMed]
- 27. Meyer U, van Kampen M, Isovich E, Flügge G, Fuchs E. Chronic psychosocial stress regulates the expression of both GR and MR mRNA in the hippocampal formation of tree shrews. Hippocampus. 2001; 11:329–36. https://doi.org/10.1002/hipo.1047 [PubMed]
- 28. Fuchs E. Social stress in tree shrews as an animal model of depression: an example of a behavioral model of a CNS disorder. CNS Spectr. 2005; 10:182–90. https://doi.org/10.1017/S1092852900010038 [PubMed]
- 29. Dong WW, Li HM, Qing XR, Huang DH, Li HG. Identification and characterization of human testis derived circular RNAs and their existence in seminal plasma. Sci Rep. 2016; 6:39080. https://doi.org/10.1038/srep39080 [PubMed]
- 30. Gruner H, Cortés-López M, Cooper DA, Bauer M, Miura P. CircRNA accumulation in the aging mouse brain. Sci Rep. 2016; 6:38907. https://doi.org/10.1038/srep38907 [PubMed]
- 31. Ji XY, Wang JX, Liu B, Zheng ZQ, Fu SY, Tarekegn GM, Bai X, Bai YS, Li H, Zhang WG. Comparative Transcriptome Analysis Reveals that a Ubiquitin-Mediated Proteolysis Pathway Is Important for Primary and Secondary Hair Follicle Development in Cashmere Goats. PLoS One. 2016; 11:e0156124. https://doi.org/10.1371/journal.pone.0156124 [PubMed]
- 32. Ham S, Kim TK, Lee S, Tang YP, Im HI. MicroRNA Profiling in Aging Brain of PSEN1/PSEN2 Double Knockout Mice. Mol Neurobiol. 2017. https://doi.org/10.1007/s12035-017-0753-6 [PubMed]
- 33. Vilchez D, Saez I, Dillin A. The role of protein clearance mechanisms in organismal ageing and age-related diseases. Nat Commun. 2014; 5:5659. https://doi.org/10.1038/ncomms6659 [PubMed]
- 34. Zeinab RA, Wu H, Sergi C, Leng R. UBE4B: a promising regulatory molecule in neuronal death and survival. Int J Mol Sci. 2012; 13:16865–79. https://doi.org/10.3390/ijms131216865 [PubMed]
- 35. Ruiz R, Pérez-Villegas EM, Bachiller S, Rosa JL, Armengol JA. HERC 1 Ubiquitin Ligase Mutation Affects Neocortical, CA3 Hippocampal and Spinal Cord Projection Neurons: An Ultrastructural Study. Front Neuroanat. 2016; 10:42. https://doi.org/10.3389/fnana.2016.00042 [PubMed]
- 36. Okumu DO, East MP, Levine M, Herring LE, Zhang R, Gilbert TS, Litchfield DW, Zhang Y, Graves LM. BIRC6 mediates imatinib resistance independently of Mcl-1. PLoS One. 2017; 12:e0177871. https://doi.org/10.1371/journal.pone.0177871 [PubMed]
- 37. Ren J, Shi M, Liu R, Yang QH, Johnson T, Skarnes WC, Du C. The Birc6 (Bruce) gene regulates p53 and the mitochondrial pathway of apoptosis and is essential for mouse embryonic development. Proc Natl Acad Sci USA. 2005; 102:565–70. https://doi.org/10.1073/pnas.0408744102 [PubMed]
- 38. Guo JU, Agarwal V, Guo H, Bartel DP. Expanded identification and characterization of mammalian circular RNAs. Genome Biol. 2014; 15:409. https://doi.org/10.1186/s13059-014-0409-z [PubMed]
- 39. Michaelis T, de Biurrun G, Watanabe T, Frahm J, Ohl F, Fuchs E. Gender-specific alterations of cerebral metabolites with aging and cortisol treatment. J Psychiatr Res. 2001; 35:231–37. https://doi.org/10.1016/S0022-3956(01)00025-5 [PubMed]
- 40. Ernst J, Bar-Joseph Z. STEM: a tool for the analysis of short time series gene expression data. BMC Bioinformatics. 2006; 7:191. https://doi.org/10.1186/1471-2105-7-191 [PubMed]
- 41. Inglis DO, Skrzypek MS, Arnaud MB, Binkley J, Shah P, Wymore F, Sherlock G. Improved gene ontology annotation for biofilm formation, filamentous growth, and phenotypic switching in Candida albicans. Eukaryot Cell. 2013; 12:101–08. https://doi.org/10.1128/EC.00238-12 [PubMed]
- 42. Blake JA, Dolan M, Drabkin H, Hill DP, Li N, Sitnikov D, Bridges S, Burgess S, Buza T, McCarthy F, Peddinti D, Pillai L, Carbon S, et al, and Gene Ontology Consortium. Gene Ontology annotations and resources. Nucleic Acids Res. 2013; 41:D530–35. https://doi.org/10.1093/nar/gks1050 [PubMed]
- 43. Balakrishnan R, Harris MA, Huntley R, Van Auken K, Cherry JM. A guide to best practices for Gene Ontology (GO) manual annotation. Database (Oxford). 2013; 2013:bat054. https://doi.org/10.1093/database/bat054 [PubMed]
- 44. Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T, Yamanishi Y. KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008; 36:D480–84. https://doi.org/10.1093/nar/gkm882 [PubMed]