Identification of Circular RNAs using RNA Sequencing

Circular RNAs (CircRNAs) are endogenous, non-coding RNAs that have gained attention given their pervasive expression in the eukaryotic transcriptome^1,2,3. They are formed when exons back-splice to each other and hence were initially considered to be splicing artifacts^4,5. However, recent studies have demonstrated that circRNAs exhibit cell type, tissue, and developmental stage specific expression^3,6 and are evolutionarily conserved^2,3. Furthermore, they are involved in mediation of protein-protein interactions⁷, micro-RNA (miRNA) binding^3,8,9,10, and regulation of parental gene transcription¹¹.

In classical RNA sequencing (RNA-seq), circRNAs may be completely lost during library construction as a result of poly-A selection for mRNAs or may be difficult to isolate given their low abundance. However, recent circRNA characterization studies have incorporated a pre-treatment step using RNase R in order to enrich for circRNAs^2,12,13. RNase R is an exoribonuclease that digests linear RNAs, leaving behind circular RNA structures. CircRNA enrichment protocols were optimized by generating and comparing data from two commercially available whole transcriptome library construction kits, with and without an RNase R pre-treatment step, and using varying amounts of total RNA input (1 to 4 µg). The optimized protocol was next used to evaluate the abundance of circRNAs across five different brain regions (cerebellum [BC], inferior parietal lobe [IP], middle temporal gyrus [MG], occipital cortex [OC] and superior frontal gyrus [SF]) and four other tissue types (liver [LV], lung [LU], lymph node [LN] and pancreas [PA]). RNA-seq libraries were paired end sequenced and data was analyzed using six different circRNA prediction algorithms: find_circ³, CIRI¹⁴, Mapsplice¹⁵, KNIFE¹⁶, DCC¹⁷, and CIRCexplorer¹⁸. Based on our analysis, the highest number of unique circRNAs was detected when using a stranded total RNA library preparation kit with RNase R pre-treatment and 4 µg total input RNA. The optimized protocol is described here. As previously reported^19,20, the highest enrichment of circRNAs was observed in the brain compared to other tissue types.

This research has been performed in compliance with all institutional, national and international guidelines for human welfare. Brain tissues were obtained from the Banner Sun Health Research Institute Brain and Body Donation Program in Sun City, AZ. The operations of the Brain and Body Donation Program are approved by the Western Institutional Review Board (WIRB protocol #20120821). All subjects or their legal representatives signed the informed consent. Commercial (non-brain) biospecimens were purchased from Proteogenex.

1. RNase R Treatment

NOTE: In the following steps, the reaction volume is adjusted to a total volume of 50 µL. This is the minimum sample volume to be used in the RNA cleanup & concentrator kit (see Table of Materials). Additionally, the optimized protocol described here is for an input amount of 4 µg of total RNA. A longer incubation time for RNase R treatment is recommended for an input amount >4 µg.

1. Dilute total RNA to 4 µg in 39 µL RNase-free water in a microcentrifuge tube and mix well by pipetting.
2. In a separate tube, dilute the RNase R to a working concentration of 2 U/µL with 1x RNase R Reaction Buffer. Make only enough for immediate use.
3. Pipette 39 µL of total RNA and 5 µL of 10x RNase R Reaction Buffer into a 1.5 mL reaction tube and mix well by pipetting (50 µL will be the total reaction volume). Next, add 6 µL of RNase R (2 U/µL).
4. Adjust the pipette to the full reaction volume (50 µL) and mix well by pipetting up and down 10 times.
5. Place the tube in a 37 °C water bath for 10 min. Make sure that the full reaction volume is immersed in the water bath.
6. Place the tube on ice and immediately proceed with RNA cleanup & concentration (section 2).

2. Purifying RNA Using an RNA Cleanup and Concentrator Kit

NOTE: When using high quality RNA (RIN>8, DV200>80%), RNase R treatment may result in loss of approximately 60% of RNA. Using a 4 µg input, it is estimated that 2–2.5 µg of treated RNA is left after section 1.

1. Before starting, prepare the RNA Wash Buffer by adding 48 mL of 100% ethanol to the buffer concentrate and mix well by pipetting. Place purification columns into collection tubes (see Table of Materials) and place in a tube rack.
   NOTE: Use the following centrifugation settings for all following steps: 10,000–16,000 x g. If DNase I treatment has already been performed, skip DNase I treatment at this stage.
2. Add 2 volumes of RNA Binding Buffer to the RNase R treated sample, and mix well by pipetting (total volume: 150 µL).
3. Add 1 volume of 100% ethanol to the RNA Binding Buffer and RNase R treated sample mixture, and mix well by pipetting (total volume: 300 µL).
4. Transfer the entire volume to the column and centrifuge the column for 30 s. Discard flow through.
5. Add 400 µL of RNA Prep Buffer directly to the column, centrifuge the column for 30 s, and discard flow through.
6. Add 700 µL of RNA Wash Buffer directly to the column, centrifuge the column for 30 s, and discard flow through.
7. Add 400 µL of RNA Wash Buffer directly to the column, centrifuge the column for 2 min, and transfer the column to a fresh RNase-free 1.5 mL tube.
8. Add 11 µL of RNase-free water directly to the column by holding the pipette tip right above the column filter and ensuring that water lands only on the column filter.
9. Incubate the column for 1 min at room temperature and centrifuge for 1 min.
10. Before discarding the column, check for flow through in the RNase-free tube. If elution was successful, store sample at -80 °C or immediately proceed with library preparation. The final total elution volume of approximately 10 µL is used for library construction.
    NOTE: Stopping point: Leave RNA at -80 °C for up to 7 days before continuing with library preparation.

3. circRNA Library Prep

NOTE: See Table of Materials for kit, which contains most reagents used in this section.

1. rRNA depletion and fragmentation
   1. Transfer 10 µL of purified RNA from step 2.10 to a clean well in a new 96-well 0.3 mL PCR plate. To the well, add 5 µL of rRNA Binding Buffer followed by 5 µL rRNA Removal Mix. Gently pipette up and down 10 times to mix.
   2. Seal plate and incubate for 5 min at 68 °C on a pre-programmed, pre-heated thermocycler block. After completion of the 5 min incubation, place plate on bench and incubate at room temperature for 1 min.
   3. Remove seal from plate. Add 35 µL of vortexed room temperature rRNA removal beads to sample. Adjust the pipette to 45 µL and pipette up and down 10–20x to mix thoroughly. Incubate plate for 1 min at room temperature.
   4. Transfer plate to a magnetic stand and incubate on the stand for 1 min or until the solution clears. Transfer all of the supernatant (~45 µL) to new well on the same plate, or new plate (depending on how many samples you are working with).
   5. Vortex the RNA cleanup beads (see Table of Materials) until well dispersed, and add 99 µL of beads to each sample. Pipette up and down 10x to mix. Incubate the plate at room temperature for 10 min.
   6. Transfer the plate to the magnetic stand and incubate an additional 5 min or until the solution clears. Remove and discard all of the supernatant from the well.
   7. With the plate still on the magnetic stand, add 200 µL of freshly prepared 80% EtOH to the well without disrupting the beads. Incubate for 30 s, then remove and discard ethanol. Repeat for a total of 2 washes.
   8. Add 11 µL of Elution Buffer to each well and pipette up and down 10 times to mix. Incubate at room temperature for 2 min, and then transfer to the magnetic stand until the solution clears (1–5 min).
   9. Transfer 8.5 µL of the supernatant from the well to a new well on the same plate or to a new plate. Add 8.5 µL of the Elute, Primer, Fragment High mix to each well containing sample. Pipette up and down 10 times to mix thoroughly.
   10. Seal plate and incubate for 8 min at 94 °C on a pre-programmed, pre-heated thermocycler block. Remove from thermocycler when it reaches 4 °C and centrifuge briefly.
       NOTE: Proceed immediately to the Synthesize First Strand cDNA protocol.
2. Synthesize cDNA
   1. For each sample being prepared, mix 9 µL First Strand Synthesis Mix with 1 µL of reverse transcriptase (see Table of Materials). Add 8 µL of the mixture to the sample. Pipette up and down 6 times to mix.
      1. Seal plate and incubate on a pre-programmed, pre-heated thermocycler block using the following parameters: 25 °C for 10 min, 42 °C for 15 min, 70 °C for 15 min, 4 °C hold. Proceed immediately to second strand synthesis.
   2. Add 5 µL of Resuspension buffer to each sample followed by 20 µL of Second Strand Marking master mix. Pipette the entire volume up and down 6 times.
   3. Seal plate and incubate on a pre-programmed, pre-heated thermocycler block set to 16 °C for 1 h. Following incubation, remove the plate from the thermocycler and let it equilibrate to room temperature.
   4. Vortex PCR purification beads (see Table of Materials) and add 90 µL beads to each well of sample. Pipette up and down 10 times to mix thoroughly. Incubate at room temperature for 10 min.
   5. Transfer the bead/sample mix to the magnetic stand and incubate for 5 min or until liquid clears. Remove and discard supernatant.
   6. Add 200 µL of 80% EtOH to each sample. Incubate samples on the magnetic stand at room temperature for 30 s. Discard supernatant. Repeat 1x.
   7. Allow beads to dry at room temperature for 6 min, and then remove from magnetic stand.
   8. Resuspend beads in 19.5 µL of Resuspension buffer. Pipette up and down 10 times to mix thoroughly. Incubate at room temperature for 2 min, then transfer to the magnetic stand and incubate for an additional 1 min or until the liquid clears.
   9. Transfer 17.5 µL of supernatant to new well/new plate.
      NOTE: If not proceeding immediately, samples can be stored at -20 °C for up to 7 days.
3. Library preparation
   1. Add 12.5 µL of A-Tailing Mix to each well containing supernatant. Pipette the entire volume up and down 10 times to mix.
   2. Incubate the reaction on a pre-programmed, pre-heated thermocycler block set to 37 °C using the following parameters: 37 °C for 30 min, 70 °C for 5 min, 4 °C hold. When samples reach 4 °C, proceed immediately to adapter ligation.
   3. To each sample add 2.5 µL of Resuspension buffer, 2.5 µL of a unique RNA adapter, and 2.5 µL of Ligation Mix. Pipette up and down 10x to mix.
   4. Incubate samples on a pre-programmed, pre-heated thermocycler block at 30 °C for 10 min.
   5. Add 5 µL of Stop Ligation buffer to each sample and pipette up and down to mix.
   6. Add 42 µL of mixed PCR purification beads to each sample and mix thoroughly. Follow steps 3.2.6 through 3.2.10, but change the resuspension volume to 52 µL and final elution volume to 50 µL.
   7. Repeat the PCR purification bead protocol again with the 50 µL elution from step 3.3.6, but change the resuspension volume to 22 µL and final elution volume to 20 µL.
      NOTE: If not proceeding immediately, samples can be stored at -20 °C for up to 7 days.
   8. Add 5 µL of PCR Primer Cocktail and 25 µL of PCR Master Mix to each sample. Mix by pipetting up and down 10 times. Incubate the reaction on a pre-programmed, pre-heated thermocycler block using the following parameters: 98 °C for 30 s; then 8 cycles of 98 °C for 10 s, 60 °C for 30 s, and 72 °C for 30 s; then 72 °C for 5 min, then 4 °C hold.
      NOTE: Optimization of the total number of PCR cycles may be needed to generate sufficient amounts of library for sequencing.
   9. Follow the protocol for PCR bead purification (steps 3.2.4 through 3.2.9), except add 50 µL of well-mixed PCR purification beads and change the resuspension volume to 32.5 µL with a final elution volume of 30 µL.
      NOTE: Samples should be stored at -20 °C.
4. Quantification and Quality Control using a nucleic acid analyzer
   1. Allow tapes and reagents to equilibrate at room temperature for 30 min.
   2. Mix 2 µL of library with 2 µL of HS D1000 buffer, and add to a compatible well plate.
   3. Seal tightly with compatible foil seal, and vortex for 1 min at 2,000 rpm.
   4. Spin down and load plate on the analyzer following software prompts.
      NOTE: Libraries should be approximately 260 bp in size.

4. Data Analysis Workflow

1. Sequence RNA-seq libraries (see Table of Materials) to generate 82 bp paired-end reads. Convert raw sequencing data in the form of basecall files (.bcl) to FASTQs using the bcl2fastq tool (v0.2.19).
2. Detect circRNAs.
   NOTE: Based on previously reported evidence that an ensemble circRNA detection approach performs better compared to using a single detection tool^21,22, we suggest using multiple tools for circRNA detection. Here, circRNAs were identified using six existing circRNA prediction algorithms: find_circ, CIRI, CIRCexplorer, Mapsplice, KNIFE, and DCC, applying the recommended parameter settings for each algorithm.
   1. Download and install each circRNA detection algorithm on a Linux high performance computing cluster using the instructions provided by the developers.
   2. Align RNA-seq FASTQs against the reference genome (GRCh37), utilizing the aligner recommended for each tool.
   3. Following alignment, execute circRNA detection algorithms by applying their respective recommended parameter settings.
   4. Each tool will output a multi-column results file with the list of detected circRNAs, extract the circRNA co-ordinates and the number of supporting reads from this in order to quantify the number of candidates detected in each sample/test condition.
3. Convert the circRNA coordinates output by CIRI, Mapsplice, and DCC to 0-based coordinates to be consistent with the other three algorithms.
4. Select circRNAs with two or more supporting reads or downstream analyses and comparisons. Table 1 summarizes all the parameters evaluated in our study along with the total number of sequencing reads generated for each sample.
5. For each sample/test condition, count the number of detected circRNAs normalized to the number of mapped reads generated for that library, per million. Summarize the results across the various tools/samples in box plots, as detailed in the Representative Results.

Data generated using a commercially available universal control RNA (UC) and using two library preparation kits, both of which include a ribo-depletion step in their protocols, was first assessed. Using an analytical workflow (Data analysis workflow, section 4), overall, a higher number of circRNAs was detected in the TruSeq datasets compared to the Kapa ones (Figure 1). Although the ribosomal RNA (rRNA) percentages were below 5% in datasets from both kits for lower input amounts (1, 2 ug), Kapa datasets had higher rRNA content for 4, 5, and 10 ug inputs (Table 2). Hence, based on the number of detected circRNAs and rRNA depletion efficiency, further experiments were performed using the TruSeq kit.

Next, the significance of RNase R pre-treatment was tested by comparing the data generated from RNase R pre-treated and non-pre-treated libraries. To this end, total RNA was extracted from the MG of healthy elderly individuals and sequencing data generated from libraries with (N = 3) and without (N = 3) pre-treatment using RNase R²³ was compared. A higher number of circRNAs was consistently identified in the pre-treated libraries compared to the non-pre-treated ones (Figure 2). This is expected since pre-treatment removes linear RNAs, thus enriching for circRNA species.

Thirdly, the amount of input RNA that would be optimal for detecting a higher diversity of circRNAs was tested. Libraries were prepared using 1, 2, and 4 µg of total input RNA which was extracted from MG, OC, and SF brain regions, and as well as UC RNA. Comparing the abundance of circRNAs detected from each library, the highest diversity of circRNA species was observed when using 4 µg input RNA compared to 2 and 1 µg (Figure 3), as reflected by the number of unique circRNAs identified. One caveat to note is that although various incubation times during RNase R treatment were not tested, a trend whereby an increasing number of circRNAs were detected across total RNA inputs of 1 to 4 µg was observed when controlling for all other parameters.

This optimized protocol was then applied across multiple tissue types to compare circRNA abundances. Five brain regions, including BC, MG, OC, IP, and SF, from four healthy elderly individuals, were tested, along with four other tissue types, including LV, LU, LN and PA, from six healthy donors. Overall, a higher abundance of circRNAs was observed in the brain compared to other tissue types (Figure 4), as has been previously reported^19,20.

Figure 1: CircRNA detection using TruSeq vs. Kapa total RNA kits. Sequencing data was generated for UC RNA using two separate total RNA library preparation kits, each with 1, 2, 4, 5, and 10 µg input RNA and ribonuclease R (RNase R) pre-treatment. The number of circRNAs detected by the tools in each sample was normalized to the number of mapped reads, per million (Y-axis). Please click here to view a larger version of this figure.

Figure 2: CircRNA detection with and without RNase R pre-treatment. Sequencing data generated using the TruSeq kit was used to compare the impact of RNase R pre-treatment. RNA was extracted from the middle temporal gyrus (MG) of healthy elderly controls for this analysis. The normalized number of circRNAs detected (Y-axis) was calculated similarly to Figure 1. RNase R+ = pre-treated with RNase R, RNase R- = not pre-treated with RNase R. Please click here to view a larger version of this figure.

Figure 3: CircRNA detection using varying amount of input RNA. Using RNA extracted from MG, occipital cortex (OC), and superior frontal gyrus (SF), as well as UC RNA, the number of unique circRNAs detected when using 1, 2, and 4 µg of input RNA, each whose library was constructed using the TruSeq kit and RNase R pre-treatment, was compared. The normalized number of circRNAs detected (Y-axis) was calculated similarly to Figure 1. Please click here to view a larger version of this figure.

Figure 4: CircRNA detection in brain vs. other tissue types. CircRNA enriched datasets using RNA extracted from various brain regions including cerebellum (BC), inferior parietal lobe (IP), MG, OC, and SF, as well as four other tissue types including liver (LV), lung (LU), lymph node (LN), and pancreas (PA) were generated. CircRNA enrichment was carried out using the Illumina TruSeq kit with RNase R pre-treatment and 4 µg of total input RNA. Box plots represent the number of circRNAs detected by at least three of the six tools across the samples from each brain region/tissue type. Please click here to view a larger version of this figure.

Table 1: Sample and test conditions summary. Summary of all parameters and test conditions evaluated in this study, along with the total number of sequencing reads generated for each sample. UC = universal control, MG = middle temporal gyrus, OC = occipital cortex, BC = cerebellum, IP = inferior parietal lobe, SF = superior frontal gyrus, LU = lung, LV = liver, LN = lymph node, PA = pancreas, RNase R = ribonuclease R, RNase R+ = pre-treated with RNase R, RNase R- = not pre-treated with RNase R.

Table 2: rRNA percentages in TruSeq vs. Kapa libraries.

In this study, two commercially available library preparation kits, pre-treatment options, and input RNA amounts were tested in order to optimize a circRNA enrichment protocol for construction of circRNA sequencing libraries. Based on this study’s assessments, a number of key aspects and critical steps in creating circRNA sequencing libraries are apparent. Our evaluation confirms the utility of RNase R pre-treatment, as reflected by the increased number of circRNAs detected. Overall, a higher diversity of circRNAs when using the Illumina TruSeq library kit with RNase R pre-treatment and 4 µg of input RNA was observed. These results align with previous findings that the RNase R enrichment step is beneficial for detection of circRNAs².

Additional key aspects of circRNA library construction include the amount of total RNA that is available for sequencing as well as the type of tissue that the RNA extracted from. Although a 4 µg input of total RNA was found to yield the highest number of detected circRNAs, the majority of RNAseq studies utilize <=1 µg of total RNA such that obtaining higher amounts may be challenging, particularly for analysis of human specimens. Identification of circRNAs remains feasible for lower input amounts, but it is relevant to acknowledge that the specificity of the analysis may be impacted. This study further highlights the higher number of circRNAs that are detected in human brain compared to other tissues, as previously reported^19,20. It is thus critical to acknowledge the differential expression of circRNAs across different tissue types. Furthermore, additional research in the context of disease will be important for shedding light into how circRNAs may be involved in pathogenic processes.

The performance of the two assessed RNA library preparation kits also highlights that although different commercially-available kits may demonstrate significant similarities, differences are still observed when analyzing circRNAs. Two major findings from this comparison include decreased rRNA depletion and a lower number of circRNAs identified using one approach. While one possibility is that a higher abundance of rRNA in a sample may interfere with creating sequence-able circRNA library molecules, this finding emphasizes the need to assess seemingly similar kits, particularly when reagents are proprietary.

Although the data presented here provides insights into the existence and abundance of circRNAs in various tissue types, this study has a few technical limitations. Firstly, while RNase R treatment reduces the population of linear RNAs in a sample, it is not well understood if this depletion step introduces any biases in circRNA detection and whether it may deplete circRNAs. Previous studies have reported that in some cases, circRNAs are sensitive to RNase R^2,24,25. Secondly, it is unclear if increasing the total RNA input above 4 µg will result in a linear increase in the number of identified circRNAs. As previously mentioned, available total RNA is often limited in research studies so lower input amounts were considered here. Of note, circRNAs can still be detected when using lower inputs but it is important to acknowledge that lower inputs are associated with detection of a lower number of circRNAs. Thirdly, the optimized protocol presented here utilizes RNAs extracted from a specific set of tissues. Given the variable distribution of circRNA expression across different tissue types, the association between total RNA input amounts and the number of identified circRNAs may differ across tissues.

With increasing interests in understanding the biological role of circRNAs, new strategies are also being developed to better enable characterization and identification of circRNAs. One new bioinformatics approach enables identification of circRNAs that may be lowly expressed through reconstruction of full-length circRNAs, and also enables quantification of expression of specific circRNA isoforms²⁶. This approach takes advantage of features described as reverse overlap (RO) reads that may occur on the 3’ or 5’ ends of circRNA library molecules. Development of new strategies for identifying circRNAs, encompassing both laboratory approaches and bioinformatics tools, will contribute to the field’s understanding of the function and impact of circRNAs.