Annotation and functional clustering of circRNA expression in rhesus macaque brain during aging

The abundance and function of circular RNAs (circRNAs) in mammalian brain have been reported, but their alterations in the biology of brain aging remain elusive. Here, using deep RNA profiling with linear RNA digestion by RNase R we explored a comprehensive map of changes in circRNA expression in rhesus macaque (macaca mulatta) brain in two age groups from adult (10 y) to aged (20 y) periods. Total 17,050 well expressed, stable circRNAs were identified. Cluster analysis reveals that dynamic changes in circRNA expression show the spatial-, sex- and age-biased specificities. On the basis of separate profiling of the RNAs, age-related circRNAs show differential correlation to host mRNA expression. Furthermore, two voltage-dependent L- and R-type calcium channel gene-derived circCACNA2D1 and circCACNA1E negatively regulate their host mRNA expression. Our results demonstrate the power of changes in circRNA expression to reveal insights into a potentially circRNA-mediated regulatory mechanism underlying the biology of brain aging.

Most circRNAs were low expressed with no more than 10 back splicing reads.
(b) Venn diagram showing the overlap between the published circRNAs and the predicted circRNAs by CIRI2.
(c) Box plot of the circRNA expression level by different groups: published specific, CIRI2 specific and theirs overlapped. The expression level of specific circRNAs was significant lower than that of overlapped circRNAs.
(d) Gels presentation of the circRNAs that were validated by Sanger sequencing. The arrows represent the location of back-splicing fragments. Nine circRNAs were selected as presentation.
(e) CIRI2 specifically detected circMARK4 was validated RT-qPCR and Sanger sequencing for sequences around the back splicing sites. The bases above the black line indicate the back splicing junction.
(f-k) Co-detected circRNAs by CRIR2 and find_circ were validated by RT-qPCR and Sanger sequencing. The presentation format was the same with (e).

Supplementary Fig. S2
Patterns and functional analysis of circRNA expression in rhesus macaque brain (c) Bar plot showed the top 10 enriched biological process terms from GO database for the circRNA host genes.

Supplementary Fig. S3
Analysis and characteristics of spatial-related circRNAs in rhesus macaque brain  (c) Ten-micron cryostat sections of 10-and 20-year-old sex-matched macaque PFC and hippocampus were immunostained with R-type voltage-dependent calcium channels CACNA1E antibody using HRP immunocytochemistry (brown). Scale bar, 50 µm.

Animals and Samples collection
Frozen postmortem brain tissue samples from rhesus macaque were provided by Kunming Primate Research Center of the Chinese Academy of Sciences (KPRC). Brain regions were systematically collected from well-characterized rhesus monkeys born and raised at the KPRC in outdoor, 6-acre enclosures that provide a naturalistic setting and normal social environment. For circRNA-seq analysis, two male and two female rhesus macaque brain specimens at each of stages representing adult (10-year-old) and old (20-year-old) were profiled. Extensive health, family lineage and dominance information were maintained on all animals.
According to a widely used macaque brain atlas and brain maps (http: //www.Brainmaps.org), tissues spanning eight anatomically distinct regions were selected and collected from each specimen. The detailed information was described as below: the prefrontal cortex was sampled at the main sulci, the posterior cingulate cortex was sampled at the Brodmann's area 23, the temporal cortex at the superior temporal gyrus, the parietal cortex at the middle sylvian fissure, the occipital cortex at the V1, and the cerebellar cortex was sampled at the cauda cerebellum.
The hippocampus (including CA1 and dentate gyrus) was also sampled. All the collected samples were washed with RNA later solution (AM7021, Ambion, USA) and put in freezing tubes to store at liquid nitrogen temperature.
All animal procedures were in strict accordance with the guidelines for the National Care and

CircRNA library preparation
Before the preparation of the circRNA-seq libraries, the total RNA samples for two males were mixed and those for females were mixed as well to eliminate the individual biases. Samples were collected from eight different brain regions, and consequently a total of 32 mixed samples from 10-year male and female, and from 20-year male and female were subjected to library preparation. For each sample, 10µg of total RNA was used for cirRNA-seq library preparation.

circRNA Prediction and filtering
After reads were obtained from sequencing, Cutadapt (1) was used to trim adaptors (with 10% error rate) and low quality bases (less than 16). Reads that were less than 16nt in both ends were discarded. Before prediction, we used bwa(2) method to align reads to the Rhesus Macaque (3) reference genome (with parameters: -k 19 -T 19). After alignment, CIRI2 software (4, 5) was used to predict circRNAs in macaque brain samples. Candidate circRNAs with no less than two head-to-tail junction reads were preserved.
To test the precision and sensitivity of the CIRI2 software, we downloaded the raw data of the SH-SY5Y neuronal cells (6). We aligned the reads to the human genome (GRCH38 version) and predicted the circRNAs with our prediction pipeline with the same parameters. Then we converted the published circRNAs locus from GRCH37 to GRCH38 by liftOver (7). CircRNAs which splice sites were detected in +/-5 nucleotide interval around the predicted circRNAs start or end were treated as overlap circRNAs. At the same time, we re-predicted circRNAs with our data by utilizing find_circ (version 1.2) software with default parameters (8). Overlapping analysis was performed for the prediction results between CIRI2 and find_circ.
CircRNAs are constituted by two parts by CIRI2: back splicing junction part and non-junction part. After obtaining the circRNAs, reads number across the back splicing junctions was treated as the raw expression level of each circRNAs. After that, we filtered the circRNAs expressed only in sporadic samples with low abundance. Based on the back-splicing reads number of each circRNA, we set the following requirement as the filtering threshold: 1) No less than 3 back-splicing reads was required to define an expressed circRNA.
2) One brain area is defined as expressed only if no less than 2 samples satisfied the 1st requirement.
3) One circRNA is expressed only if on less than 2 areas were expressed.

4)
If circRNAs can't meet the 3rd requirement, the circRNAs must be expressed in only one brain area and the total expressed samples must be on less than 6 (20% of total samples).

5)
If circRNAs can't meet the 3rd and 4th requirements, the total back-splice reads number of this circRNA must be no less than 30.
To obtain a normalized expression between different samples, we used the DESeq normalization method (9) to obtain the normalized expression level for circRNAs by dividing a size factor.

Conservation analysis for macaque circRNAs
To analyze the homologous feature between macaque circRNAs and other mammalian species, we downloaded the human and mouse circRNAs sequence from the circBase (10). Then we aligned the circRNA sequence to the human and mouse circRNAs by blastn with E-value 1e-3.
Meanwhile, we also adopted the conservation calculation method in Rybak-Wolf et al (6), and performed this analysis among three species: macaque, mouse, and human.

Characteristic analysis of circRNA expression
To find the characteristic of circRNA expression in macaque brain, we compared the circRNA expression pattern from multiple dimensions. First, we used WGCNA method (11) to classify circRNAs. Cubic power was chosen as the soft threshold to calculate block wise modules. The output is the circRNA modules according to their expression pattern. For each circRNA module, eigengenes was chosen to represent the expression pattern.
Due to the low and dynamic characteristic of circRNAs compared with host mRNAs, we used unpaired t-test model to filtrate the spatial, ageing and sexual-related circRNAs as described previously (12). For spatial-related circRNAs filtration, we compared one brain area with other areas by t-test and repeated this according the eight brain areas. For ageing-related circRNAs filtration, we compared 10-year-old samples with 20-year-old samples. For sexual-related circRNAs filtration, we compared female samples with male samples. After comparison, we set p-value 0.05 as the threshold for filtration.

Analysis of the relationships between circRNAs and host mRNAs
To explore the expression relationship of circRNAs and their host mRNAs, we used the poly(A) selected RNA-seq data from our another study (13). The biological samples of the two studies were the same and total RNAs were extracted from the brain tissue parallelly. So we could compare the expression level between circRNAs and host mature mRNAs. Based on the expression of each circRNA and host mRNA, Pearson correlation coefficient (PCC) and p-value were obtained for each mRNA-circRNA pair. Then we filter the result by a given threshold, absolute PCC no less than 0.3 and p-value no more than 0.1. Besides the positive correlation pairs, negative pairs with correlation coefficient less than 0 were also included.

RNA extraction and reverse transcription PCR
RNA was using PureLink micro-to-midi total RNA purification system (Invitrogen). RNA quantity and quality were evaluated by Nano Drop ND-1000 spectrophotometer (Nano Drop Thermo, Wilmington, DE) and RNA integrity was assessed by agarose gel electrophoresis.
Specific divergent primers were designed to amplify the circular and linear Rhesus macaque transcripts. Semi-quantitative RT-PCR was performed with Superscript III one-step RT-PCR system with platinum Taq High Fidelity (Invitrogen). For quantitative real-time PCR, cDNAs were prepared by using oligonucleotide (dT), random primers, and a Thermo Reverse Transcription kit (Signal way Biotechnology). qPCR was performed with SYBR Green I Master GAPDH was used as a reference. Both target and reference were amplified in triplicate wells.
And the relative level of each circular and linear transcript was calculated using 2 − △△ Ct method.
All PCR primer sequences can be found in Supplementary Table S3.

Immunohistochemistry
For DAB/bright field staining, ten micron cryostat sections of macaque brain were pretreated in 0.3% hydrogen peroxide in methanol for 30 min to remove endogenous peroxidase activity,

Protein extraction and western blot
Frozen macaque brain tissues or cultured fetal hippocampal neurons were lysed using 1× RIPA buffer (Sigma Aldrich), complemented with PIC (Protease Inhibitor Cocktail, Cat #11836153001, Roche) and Phosphatase Inhibitor Cocktail (Cat #11836153001, Roche). After protein quantification, samples were boiled for 10 min in SDS loading buffer (1:1 ratio). An aliquot (up to 50 µg) of the resulting sample was run on a SDS-PAGE gel and then transferred to a Hybond-PVDF membrane (Amersham Pharmacia). The membrane was blocked in 5% milk for 1 h at room temperature and incubated overnight with the appropriate primary antibody re-suspended in 5% milk in TBST at 4 °C. Following three time washes using TBST, the membrane was incubated for 2 h at room temperature with the appropriate secondary antibody re-suspended in 5% milk in TBST at room temperature, followed by 4 washes using TBST. Signal ECL (Pierce) amplification was detected by Tanon 5200 Multi chemiluminescent image system (Tanon, Shanghai, China). Signal intensity was quantified with Image J (National Institute of Health).

Fetal Rhesus macaque neural cell cultures
Fetal Rhesus macaque neuronal cultures were prepared with modification and maintained as described previously (14) .

siRNA knockdown in fetal macaque brain cultured neurons
Lentiviral siRNA preparation was prepared as described previously (15). piLenti TM -siRNA constructs including a set of circCacna2d1 and a set of circCacna1e were constructed using a standard protocol from Applied Biological Material Inc. Prism was used to perform t-tests and to visualize bar charts. Error bars represent SEM.

BASEscope Assays for circRNA detection
BASEscope assays were performed using BaseScope™ Detection Reagent Kit-RED (#322900-USM, Advanced Cell Diagnostics (ACD) Hayward, CA) in accord with the manufacturer′s protocol. CircCACNA2D1 and circCACNA1E junction site-targeting or -nontargeting label probes conjugated to HRP were ordered from ACD. Briefly, 10 µm cryostat macaque brain sections were pretreated with Hydrogen Peroxide followed by performing target retrieval using an Oster ® Steamer. Dried slides were placed on the slide rack, and incubated with RNAscope ® Protease III at 40°C for 30 minutes in the HybEZ™ system (ACD). Next, the slides were incubated at 40°C in order to hybridize probes (circCACNA2D1 and circCACNA1E) for 2 HRS BASEScope intensity was normalized to the control. Images were thresholded in Image J and, using the Image Calculator tool, a third image was generated showing those pixels which were positive in all input channels. Using the Particle Analysis tool, the size and number of the thresholded clusters were analyzed. Microsoft Excel was used to calculate the fraction of positive clusters. GraphPad Prism was used to perform ANOVAs and t-tests and to visualize bar charts.
Error bars represent SEM. The target genes, probed regions, and sequences of target probes are listed in Supplementary Table S4.

DIG-labeled anti-sense RNA probes synthesis by in vitro transcription
PCR primers were designed using standard primer designing tools (Primer Premier 5.0) to amplify 100-300 nt fragment corresponding to linear and circular RNA sequences or 100-150 nt fragment corresponding to head-to-tail junction (short circRNA specific probes and not overlapping with linear RNA sequence). T7 promoter sequence was added to the reverse primer to obtain an antisense probe in vitro transcription reaction. In vitro transcription was performed using T7 RNA polymerase (Roche) with DIG-RNA labeling mix (Roche) according to manufacturer's instruction. DNA templates were removed by DNAs I digestion and RNA probes purified by phenol chloroform extraction and subsequent precipitation. Probes were used at