Circular RNA Profiling by Illumina Sequencing via Template-Dependent Multiple Displacement Amplification

Circular RNAs (circRNAs) are newly discovered incipient non-coding RNAs with potential roles in disease progression in living organisms. Significant reports, since their inception, highlight the abundance and putative functional roles of circRNAs in every organism checked for, like O. sativa, Arabidopsis, human, and mouse. CircRNA expression is generally less than their linear mRNA counterparts which fairly explains the competitive edge of canonical splicing over non-canonical splicing. However, existing methods may not be sensitive enough for the discovery of low-level expressed circRNAs. By combining template-dependent multiple displacement amplification (tdMDA), Illumina sequencing, and bioinformatics tools, we have developed an experimental protocol that is able to detect 1,875 novel and known circRNAs from O. sativa. The same method also revealed 9,242 putative circRNAs in less than 40 million reads for the first time from the Nicotiana benthamiana whose genome has not been fully annotated. Supported by the PCR-based validation and Sanger sequencing of selective circRNAs, our method represents a valuable tool in profiling circRNAs from the organisms with or without genome annotation.

With the assistance from computational algorithms [7,8,[27][28][29][30][31][32], numerous approaches have been developed to detect and validate the presence of circRNAs in different species across the kingdoms [12,15,33,34]. A major challenge in circRNA discovery might be attributed to its extremely low abundance in samples. The cutting-edge method involves the enrichment of circRNA by enzymatic digestion, RNA-Seq, and bioinformatics identification, followed by PCRbased validation [13,35,36]. A limitation in this approach is the sensitivity because library preparation in next-generation sequencing (NGS) is often associated with the loss of lowabundant molecules [12,37]. Thereby significant sequencing depth is required in order to identify putative circRNAs [12].
In the present study, we have introduced a step of templatedependent multiple displacement amplification (tdMDA) prior to library preparation. Together with a newly developed computer program, we have built an experimental pipeline that shows an enhanced sensitivity to identify circRNAs from the plants. . Extracted RNA was treated with 4U of Turbo DNase (2U/ L, Ambion, Austin, TX, USA) at 37 ∘ C for 30 minutes and then inactivated at 72 ∘ C for 10 minutes, followed by phenol/chloroform purification. Ten microgram of DNase-treated RNA was subjected to 10U of RNase R (20U/ L, Epicentre, Madison, WI, USA) digestion at 37 ∘ C for 15 minutes. Both integrity and concentration were determined respectively on 1.2% agarose gel and Nanodrop ND-1000 (Thermo Scientific, Waltham, MA, USA).

Materials and Methods
RNA amplification was achieved by RT-tdMDA protocol [38]. In that protocol, background amplification in MDA was eliminated by using exo-resistant random pentamer primers with their 5 ends blocked by C18 spacer [38]. Its efficiency was first evaluated by using extracted RNA and plasmid APTR9 [39,40] with the final primer concentration ranging from 50 to 200 M. About 1 g of RNase Rtreated RNA was converted into cDNA using RevertAid or RevertAid H minus first strand cDNA synthesis kit according to manufacturer's instructions (Thermo Scientific, Waltham, MA, USA). Approximately 50 ng of converted cDNA was directly used for tdMDA in a 20-L reaction consisting of 2 l of 10 mM dNTP mix, 2 l of 10X reaction buffer, 2 l of 500 M 5 end-blocked exo-resistant random pentamer primers, 0.6 l of Phi29 DNA polymerase (10U/ l, Thermo Scientific, Waltham, MA, USA), and 2 l of pyrophosphatase (0.01U/ l) (Thermo Scientific, Waltham, MA, USA). The reaction mixture was incubated at 28 ∘ C for 18 hours and terminated by heating at 65 ∘ C for 10 minutes. An aliquot of 3 l reaction was loaded on 1% TAE agarose gel to check for tdMDA performance.
. . Identification of circRNA from the tdMDA Amplicons by Cloning and Sequencing. The tdMDA amplicon was randomly digested with the restriction enzymes, SacI, HindIII, SspI, BamHI, EcoRV, and EcoRI (10U/ l, Thermo Scientific, Waltham, MA, USA) for 3 hours at 37 ∘ C. The HindIII, SacI digested fragments were purified by GeNei PCR purification kit (Bangalore, Karnataka, India) and cloned in pOK12 or pBluescript II KS (+) vector at the corresponding site at 16 ∘ C for 12 hours. After the confirmation by restriction digestion, a total of nine clones were sequenced, seven for N. benthamiana and two for O. sativa. The mapped clone sequences such as HindIII 10, HindIII 33, HindIII 38, and SacI 11 for N. benthamiana and HindIII 1 and HindIII 2 for O. sativa were subjected to prediction for their possibility of forming circRNA in Plant-circBase [41] (http://ibi.zju.edu.cn/plantcircbase/index.php). Predicted putative circRNA sequences were validated by RT-PCR with divergent primers (Table 1) [42].
. . Identification of circRNA from the tdMDA Amplicons by Illumina Sequencing. About 200 ng of tdMDA products was used for library construction using Illuminacompatible NEXTflex Rapid DNA sequencing kit (BIOO Scientific, Austin, Texas, USA) according to manufacturer's instructions, followed by sequencing at the Illumina NextSeq 500 platform (150-nt paired end) at Genotypic Technology, Bangalore as previously described [40,43]. Under genomic annotations from Ensembl plant release 29 [44], trimmed reads at Phred 23 were aligned with O. sativa Indica genome and N. benthamiana draft genome for subsequent circRNA identification using DCC software (v 0.4.4) [45]. In addition to the consideration of non-canonical splice junction, other parameters were also included for circRNA identification as postulated in the DCC [45]. All the analysis was carried out using Biolinux 8 OS [46].

. . Validation of circRNAs Derived from tdMDA-Illumina
Sequencing. Divergent primers were designed from the cir-cRNA derived from NGS-tdMDA containing the junction site (Table 1). Most primers designed for O. sativa and N. benthamiana were tested for the validation of corresponding circRNAs using genomic DNA and cDNA by the standardised annealing temperature (T A ). Divergent PCR products were subjected to sequencing or digestion with restriction enzymes.
. . Northern Hybridization. Non-radioactive northern hybridization was performed with the purified PCR fragment (>200 nt) as the probe, which spanned the corresponding cir-cRNA junction site. Probe preparation was followed with the DIG DNA labelling kit (Roche, Basel, Switzerland) according to manufacturer's instructions.

. . Analysis on circRNA Conservation and miRNA Binding
Sites. NCBI-BLASTN was used to examine the conservation of circRNAs with other reported plant species, such as S. bicolour [47], S. italica [47], B. distachyon [47], and those included in plant circular RNA database [41]. The psRNATarget, a plant small RNA target analysis server (2017 release) [48], was applied to annotate the possible role of reported O. sativa and N. tabacum miRNAs on predicted circRNAs. Table 1: List of divergent primers designed for use in confirmation of putative circRNA.

Results
. . e Elimination of Background Amplification by tdMDA. Template independent amplification (TIA) in MDA is a major concern owing to high concentrations of random hexamers and an extended incubation period [38]. To eliminate TIA, we followed the protocol proposed by Wang et al. [38]. Total RNA extracted from O. sativa was mixed with the plasmid pAPTR9 that harboured Bhendi yellow vein mosaic virus (BYVMV). After DNase treatment, BYVMV-specific PCR yielded no amplification, suggesting a complete DNA digestion in the template ( Figure S1a). This template was subsequently used to test the efficiency of RT-tdMDA protocol [38]. In four primer concentrations (50, 100, 150, and 200 M), no amplicon was found in the controls (no template) ( Figure S1b). In contrast, 50 ng of template along with 50 M final primer concentration showed an apparent amplification ( Figure S1b). Therefore, the use of exo-resistant random pentamer primers with blocked 5 ends efficiently eliminates TIA.
. . Novel circRNAs Identified by RT-tdMDA, Cloning, and Sanger Sequencing. After DNase and RNase R treatment ( Figure S2), total RNA from N. benthamiana and O. sativa plants was successfully amplified (Figure 1). Again, no amplification was observed from the negative controls (no template) (Figure 1). Of seven sequenced clones derived from N. benthamiana, four sequences, named SacI 11, HindIII 10, HindIII 33, and HindIII 38, showed 100% sequence identity in BLAST analysis against N. benthamiana genome [49] (https://solgenomics.net/organism/Nicotiana benthamiana/ genome). Two clones from O. sativa, HindIII 1 and HindIII 2, were aligned onto the intron region in chromosomes 7 and 1 of O. sativa with 100% and 99% identity respectively. These sequences were then analyzed in PlantcircBase for the potential of circRNA formation. As a result, three sequences from N. benthamiana, HindIII 10, HindIII 33, and SacI 11, were predicted to be putative circRNAs. The HindIII 10 sequence was partially mapped onto the intron domain of N3 disease resistance protein gene of Nicotiana paniculata with 96% sequence identity ( Figure S3,  Therefore, HindIII 33 sequence might be an intronic-exonic circRNA ( Figure S3, Figure S5, and Figure S6). The SacI 11 sequence (Acc. No. MF066173) was predicted as a circRNA. However, this sequence was mapped onto Nicotiana sylvestris mitochondrial genome and thus not included for further experimentation. Analysis in PlantcircBase also suggested that two clones from O. sativa, HindIII 1 (osa circ 032545) and HindIII 2 (osa circ 000547), were existing intronic and exon-intronic circRNAs respectively. For putative circRNA sequences HindIII 10 and HindIII 33, PCR amplification with divergent primers provided a positive result when cDNA was used as template whereas no amplification was observed at the same size when various concentration of genomic DNA was used as template ( Figure S4).

. . Complete circRNA Profiles Revealed by Illumina Sequencing.
Encouraged by the positive outcome from the cloning and Sanger sequencing, the amplicons from RT-tdMDA were subjected to Illumina sequencing for possible capture of entire circRNA repertoire. The total number of 150-nt paired end reads obtained from O. sativa and N. benthamiana were 21,818,956 and 38,060,238, respectively. Using the raw reads from the O. sativa, the DCC computational pipeline discovered thousands of circular splicing events that yielded 1,875 circRNAs (Figure 2(a)). These putative circRNAs are predominantly distributed on the chromosomes 1 and 5 (Figure 2 came from the genes without any particular chromosome assignment (Figure 2(b)). With respect to circRNA types, the intergenic-intergenic type (n=1,359) was the most abundant type followed by the intronic-intronic type (n=182) and the exonic-exonic type (n=123) (Figure 2(c)). Furthermore, 79% of putative circRNAs had the length between 100 and 999 nt whereas ∼1% and ∼20% had the size below 100 nt and larger than 1000 nt respectively (Figure 2(d)). The smallest circRNA was found to be of 32 nt between positions 5,187-5,219 in the genome. This putative circRNA is surprisingly an intergenic-intergenic type with CT/AC splice junction on scaffold ID AAAA02040137.1. On the other hand, the largest size of circRNA identified was 737,782 nt on the chromosome 11 between positions 18852424 and19590206, which harbours many functional genes like MIR genes, tRNA genes, and the genes encoding for hypothetical protein. The largest circRNA is assumingly formed in a non-canonical manner and categorised as an intron-intergenic type. Finally, all putative circRNAs were associated with a total of 578 genes in which ∼72% had the translation of hypothetical proteins. The gene ID BGIOSGA000405 on the chromosome 1 contributed the maximum number of circRNAs (n=35) while most genes gave only one or two circRNAs (Figure 3(a)). Individually, less than 10% of the genes could produce more than two circRNAs. Similar analysis in N. benthamiana yielded 9,242 cir-cRNAs, including 6,080 intergenic-intergenic, 1,257 intronintron, 1,009 intron-intergenic, and 896 intergenic-intronic circRNAs (Figures 4(a) and 4(b)). Interestingly, no exonic cir-cRNAs were identified probably because of unavailability of complete genome annotation. In comparison with O. sativa, circRNAs from N. benthamiana were larger in size with 69% of total identified circRNAs above 1,000 nt (Figure 4(c)). The Niben101Scf02816, an intergenic-intergenic circRNA formed by the non-canonical splice junction, had the smallest size with about 35 nt located between positions 85,132, and 85,167 on the genome. The longest circRNA was Niben101Scf03154 with 299,801 nt, also an intergenic-intergenic type between 589 and 300,390 genome positions.

. . Validation of circRNAs Identified by RT-tdMDA and
Illumina Sequencing. PCR and northern hybridization were used for the validation of selective circRNAs. For cir-cRNA Niben101Scf27324 (Nb circ7 primer, Table 1), PCR produced two distinct DNA bands with the sizes at ∼ 150 and ∼250 bp. There was no DNA amplification at similar sizes upon the use of the genomic DNA as template (Figures 5(a) and 5(b)). Sanger sequencing of PCR product mapped the larger fragment to the circRNA with extra sequence, suggesting a potential alternative splicing event involved for the biogenesis of this particular circRNA ( Figure S7). Further evidence came from the northern hybridization that signalled an apparent presence of cir-cRNA Niben101Scf27324 ( Figure S8). Divergent PCR also gave positive amplification for the osi circ10 ( Figure 5(c)) as well as other putative circRNAs including Nb circ3, Nb circ6, osi circ2, osi circ4, osi circ6, and osi circ8 (data not shown). Their authenticity was supported by restriction digestion of purified DNA bands from the gel (data not shown).
. . CircRNAs Are Conserved across the Species. Several reports claimed conservative nature of circRNA across species [8]. Therefore, we compared our circRNAs with all circRNAs either reported [15] or deposited in the plant circular RNA database [41]. Of 1,875 cicrRNAs from the O. sativa, significant similarity was found for 1,120 (60%) to O. Sativa ssp. Nipponbare, 549 (29.2%) to A. thaliana, and 145 (7.75%) to T. aestivum (Figure 6(a)). For N. benthamiana, the sequence similarity was also shared for 55 circRNAs with S. tuberosum, 60 circRNAs with A. thaliana, and 44 circRNAs with O. Sativa (Figure 7(a)). There was no or little conservation between our putative circRNAs and the circRNAs discovered in the plants like S. bicolour, S. italica, B. distachyon, H. vulgare, and P. trifoliata. This is probably due to a rare number of the circRNAs identified from these plants that could not provide a full scenario to explore the conservation (Table 2).

Discussion
CircRNAs encompass a transcript family with distinctive structures. Various methods are used to detect the circRNAs in both plants and animals [10,13]. The difficulty in circRNA identification lies in the inability to separate the circRNAs from other RNA species based on their size or electrophoretic mobility. Molecular techniques that involve amplification or fragmentation may destroy their circular nature since circRNAs lack a free 5 or 3 end. Likewise, methods that use polyadenylation ends, such as rapid amplification of cDNA ends (RACE) or poly (A) enrichment, cannot be employed for circRNA identification. These hindrances have been overcome by the emergence of exonuclease based enrichment procedures and high throughput sequencing techniques [12]. For instance, RNA sequencing has been used for the identification of circRNAs in Arabidopsis and O. sativa [14]. However, owing to an extremely low expression of circRNAs comparing to their linear mRNA counterparts, a high sequencing depth is demanded for productive capture of circRNAs. In order to improve the sensitivity, we exploited the use of tdMDA to identify circRNAs from the N. benthamiana and O. sativa plants.
Our data have demonstrated the feasibility of tdMDA in the discovery of circRNAs from small amount of RNA samples. First, both novel and known circRNAs were identified from RT-tdMDA product by random enzymatic digestion, cloning, and Sanger sequencing. Two novel circRNAs from the N. benthamiana were validated by divergent PCR and had no significant similarity with Arabidopsis and other plant candidates. The function of newly identified circRNAs from N. benthamiana has to be deciphered. Second, Illumina sequencing of the RT-tdMDA product and bioinformatics analysis captured 1,875 and 9,242 putative circRNAs from O. sativa and N. benthamiana respectively. The authenticity of selective cirRNAs was confirmed by PCR and northern hybridization. Using RNA sequencing, Jakobi et al. [50] reported the prediction of 575 circRNAs from 33.5 million reads in adult mice heart. Assuming a similar abundance among the samples, we could be able to identify much higher number of circRNAs from almost equal number of reads using the same computational pipeline. Earlier, Jeck and Sharpless (2013) stressed on the need of having 300,000-300,000,000 reads using traditional sequencing to get a single circRNA event whereas exonic circRNA is thought to present  roughly 1% of poly(A) RNA [12,31]. Again, Wang et al. in 2017 analyzed over 90 million raw reads and could obtain only 88 circRNAs. Analyzing say, more than 500 million reads will surely increase the chances of getting low abundance circRNAs but it will spike up the overall cost tremendously. Our method reduces the cost significantly without compromising on findings of lowly expressed circRNAs. Finally, the conservative nature of most predicted circRNAs across the plants further suggests the methodological reliability. Besides tdMDA, our experimental pipeline also takes the power of the bioinformatics tool. We have applied the DCC that gives the expression count of putative circRNAs as well as the linear RNAs expected from the same genome positions. Interestingly, 19% of O. sativa circRNAs showed overexpression with respect to their linear counterparts (Figure 3(b)). Functional aspects of circRNAs have not been fully understood in spite of the reports for their roles in miRNA sponging, transcriptional inhibition and protein formation [12]. Our analysis revealed there are 33 circRNAs that bind to 156 miRNAs in O. sativa Nipponbare (Table S1). This number is translated into 473 miRNA binding sites as a single circRNA can bind to more than one miRNA or vice versa (Table S2). Of the 473 miRNA binding sites, 391 sites (∼83%) are cleavage specific while the remaining 82 sites (∼17%) are possibly getting sponged by their targets (Figure 6(b)). In N. benthamiana, 2,099 circRNAs could have 8,149 miRNA binding sites on 163 published N. tabacum miRNAs (Table S3 and S4). Approximately 85% (n=6,916) of miRNA binding sites are cleavage specific and 15% (n=1,233) are target inhibitory in action that need to be deciphered (Figure 7(b)).

Conclusion
In summary, through the combination of tdMDA and bioinformatics tools, we have established an experimental protocol

Supplementary Materials
Table S1: list of 156 miRNAs binding with 33 predicted Indica circRNAs. Table S2: number of miRNA binding sites for each of the 33 Indica circRNAs with their position in the genome. Table S3: number of miRNA binding sites for each of the 2099 N. benthamiana circRNAs with their position in the genome. Table S4: list of N. tabacum miRNAs with their number of binding sites on 2099 N. benthamiana circRNAs. Figure S1: absence of TIA in MDA. Figure S2: RNA extraction from N. benthamiana and O. sativa plants. Figure S3: linear and predicted circular maps of the putative circRNAs of N. benthamiana. Figure S4: divergent PCR for circRNA confirmation. Figure S5: BLAST analysis of N. benthamiana cloned sequences in Sol Genomics Network. Figure S6: BLAST analysis of N. benthamiana cloned sequences in NCBI. Figure S7: mapping and validation of N. benthamiana circRNA. Figure S8: northern blotting for confirmation of circRNA. (Supplementary Materials)