key: cord-1044749-kq75oikl
authors: Yang, Shaomin; Zhou, Hong; Liu, Mingde; Jaijyan, Dabbu; Cruz‐Cosme, Ruth; Ramasamy, Santhamani; Subbian, Selvakumar; Liu, Dongxiao; Xu, Jiayu; Niu, Xiaoyu; Li, Yaolan; Xiao, Lizu; Tyagi, Sanjay; Wang, Qiuhong; Zhu, Hua; Tang, Qiyi
title: SARS‐CoV‐2, SARS‐CoV, and MERS‐CoV encode circular RNAs of spliceosome‐independent origin
date: 2022-03-31
journal: J Med Virol
DOI: 10.1002/jmv.27734
sha: 59c49152122faf9465fd530f05f9cd61e3508c8d
doc_id: 1044749
cord_uid: kq75oikl

Circular RNAs (circRNAs) are a newly recognized component of the transcriptome with critical roles in autoimmune diseases and viral pathogenesis. To address the importance of circRNA in RNA viral transcriptome, we systematically identified and characterized circRNAs encoded by the RNA genomes of betacoronaviruses using both bioinformatical and experimental approaches. We predicted 351, 224, and 2764 circRNAs derived from severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2), SARS‐CoV, and Middle East respiratory syndrome coronavirus, respectively. We experimentally identified 75 potential SARS‐CoV‐2 circRNAs from RNA samples extracted from SARS‐CoV‐2‐infected Vero E6 cells. A systematic comparison of viral and host circRNA features, including abundance, strand preference, length distribution, circular exon numbers, and breakpoint sequences, demonstrated that coronavirus‐derived circRNAs had a spliceosome‐independent origin. We further showed that back‐splice junctions (BSJs) captured by inverse reverse‐transcription polymerase chain reaction have different level of resistance to RNase R. Through northern blotting with a BSJ‐spanning probe targeting N gene, we identified three RNase R‐resistant bands that represent SARS‐CoV‐2 circRNAs that are detected cytoplasmic by single‐molecule and amplified fluorescence in situ hybridization assays. Lastly, analyses of 169 sequenced BSJs showed that both back‐splice and forward‐splice junctions were flanked by homologous and reverse complementary sequences, including but not limited to the canonical transcriptional regulatory sequences. Our findings highlight circRNAs as an important component of the coronavirus transcriptome, offer important evaluation of bioinformatic tools in the analysis of circRNAs from an RNA genome, and shed light on the mechanism of discontinuous RNA synthesis.

Circular RNAs (circRNAs) are a class of single-stranded ncRNA species with a covalently closed circular configuration. 5 The lack of 5′-and 3′-ends makes circRNAs resistant to exonuclease-mediated degradation and thus more stable than linear RNAs. 6 CircRNAs encoded by DNA genomes are produced during gene transcription and by the spliceosome either through back-splicing of exons or from intron lariats by escaping debranching. 7 They can encode proteins 8 or function as decoys for microRNAs (miRNAs) and proteins. 9 Accumulating data show that circRNAs are important pathological biomarkers for cancers, neurological diseases, and autoimmune diseases. [10] [11] [12] Viral-derived circRNAs have also been identified from several DNA viruses and are implicated in viral pathogenesis. [13] [14] [15] [16] 

The RNA-sequencing (RNA-Seq) data sets of total RNAs harvested from SARS-CoV-or SARS-CoV-2-infected Vero E6 (African green monkey kidney) cells at 24 h postinfection (GSE153940 and GSE56193) 18 The interactions between human/African green monkey miRNAs with SARS-CoV-2 full-length circRNAs and host DEGs were predicted using miRanda (-sc 150 -en −7). 30 

For each experimentally identified BSJ or forward-splice junction (FSJ), 60 bp around the 5′-and 3′-breakpoints were compared using blastall (BLAST) to identify homologous and reverse complementary sequences of 6 bp or longer. 

Inverse RT-PCR products were separated on 2% agarose gels.

Candidate BSJ sequences were gel purified (Zymoclean Gel DNA Recovery kit, Zymo) and TA cloned (ThermoFisher) according to the manufacturers' instructions. At least eight colonies were checked for YANG ET AL. Hybridization was performed at 62°C overnight. Washing, staining, and imaging was performed as instructed by the manufacturer.

Pairs of AmpFISH probes for circRNA 29122 | 28295 and smFISH probes for ORF1 were designed and purified as previously described. 36 To specifically amplify the circRNA but not sgRNAs carrying the N gene sequence, one probe targeted the donor sequence and the other targeted the acceptor sequence of the circRNA back-splice junction. Hybridization chain reaction (HCR) hairpin sequences were added to the probes as previously described. 36 Probe sequences were provided in Table 2 . RNA FISH was performed as described previously. 36 Briefly, Vero E6 cells were seeded on the glass coverslips (thickness 0.1 mm) coated with 0.1% gelatin and cultured in DMEM media with 10% FBS and antibiotics penicillin (100 IU/ml)-streptomycin (100 µg/ml). Cells were infected with SARS-CoV-2 at an MOI of 0.01 and coverslips were fixed at 24 h postinfection. For fixation, cells were washed with 1× phosphatebuffered saline (PBS; pH 7.4) and fixed with 4% paraformaldehyde/ PBS for 10 min at room temperature (RT). After being washed with 1× PBS, the cells were incubated with 70% ethanol for 10 min at RT.

Then the cells were permeabilized with 0.05% Triton X-100 at 4°C for 10 min. For the cells to be treated with RNase R, the cells were first equilibrated with 1x RNase R reaction buffer (Lucigen, Cat#RNR07250) for 30 min followed by RNAse R (Lucigen, cat#RNR07250) treatment at 37°C for 3 h. Cells were equilibrated with hybridization wash buffer (10% formamide/2× solution of sodium citrate [SSC]) for 5 min at RT, incubated with probes specific for positive-and negative-stranded linear or circRNA at 37°C overnight in a hybridization buffer containing 10% formamide/2× SSC. 37 Coverslips were washed twice with hybridization wash buffer.

Amplification reactions of hybridized probe was performed in HCR buffer containing 5× SSC and mounted as described previously. 36 RNase R treatment was performed at 37°C for 3 h followed by AmpFISH. Z-stacks images were acquired using Axiovert 200M microscope using the same setting for all samples and were processed by maximum projection. We applied find_circ in parallel with CIRI2 and found that AG | GT signal-biased algorithms, find_circ, 25 resulted in false positive reads (Figure S1F,G and Data S2). To determine the proportion of total BSJs to total CoV-mapping reads, which were identified by CIRI2 and find_circ, we performed a statistic analysis of CoV BSJs. CIRI2 identified 0.01477%, 0.00296%, and 0.43997% of total CoV reads that were BSJs of SARS-CoV-2, SARS-CoV, and MERS-CoV circRNAs, respectively (Table S1) .

Strikingly, find_circ found 0.36121%, 0.00733%, and 0.43962% of total CoV reads that were BSJs of SARS-CoV-2, SARS-CoV, and MERS-CoV circRNAs, respectively. Therefore, CIRI2 is more reliable. and local back-splicing in regions corresponding to the N gene and the 3′-untranslated region (UTR) of CoVs. Additionally, we noticed a general enrichment of local BSJ-spanning reads along the diagonal of the graph across three data sets. Although the local BSJ reads were low outside the 3′-end of SARS-CoV-2 and SARS-CoV genome ( Figure 1C ,D), a few local BSJs located in ORF1a/b, S, and between ORF3a and M of the MERS-CoV-2 genome had moderate to high read counts ( Figure 1E ). Interestingly, similar distribution patterns have been reported for SARS-CoV-2 FSJs in sgRNAs. 3, 38 These results suggest that circRNAs may be conserved in betacoronaviruses.

To better characterize CoV circRNAs, we performed de novo reconstruction of full-length SARS-CoV-2, SARS-CoV, and MERS-CoV circRNAs using CIRI-full. 26 CoV circRNAs tend to be negative-stranded ( Figure 2B ), which is opposite to the prediction of strand preference by Cai et al. 17 Third, CoVs tend to produce single-exon circRNAs, whereas host circRNAs undergo further intron excision, 6 resulting in multiple FSJs in the circRNAs ( Figure 2C ). Our analysis predicts that 12%-35% of CoV circRNAs contain an FSJ, further supporting the existence of gaps in CoV circRNAs. Alternative intron inclusion of DNA genome-encoded circRNAs gives rise to diverse circular isoforms, 39 

Next, we systematically validated SARS-CoV-2 BSJ hotspots predicted by our bioinformatic analyses. We extracted total RNA from Vero E6 cells mock-treated or infected with SARS-CoV-2 at 24 hpi.

Forward and reverse inverse PCR primers were designed in such a way that all donor or acceptor sequences in each hotspot will be picked up ( Figure 3A,B) . To validate the two major back-splicing events, we performed RT-PCR with divergent primer sets targeting the distant BSJ hotspot 29001~29903 | 1~500 ( Figure S3A ) and the local BSJ hotspot 28501~29500 | 27501~28500 ( Figure S3B ). We also performed inverse RT-PCR with five sets of primers targeting abundant SARS-CoV-2 circRNAs predicted by CIRI2 ( Figure 3C ). Some predicted full-length CoV circRNAs contain ORFs and thus potentially encode proteins ( Figure S2-4) . We used divergent primers flanking individual ORFs to validate these circRNAs ( Figure S3C ). To determine whether the inverse PCR products were BSJs rather than nonspecific PCR product, we gel-purified candidate BSJ amplicons based on the molecular weight ( Figure 3C and S3A-C, red arrowheads), subcloned, and Sanger-sequenced at least eight colonies for each candidate. Using this pipeline, we identified 75 BSJs from 169 clones (Table 1 ORF6, N, ORF10, and the 3′-UTR to the 5′-UTR, as well as local fusion within N, from N to ORF7a, ORF7b, and ORF8, and from ORF6 to M (Data S4). Third, we found the breakpoints of a given circRNA unexpectedly flexible. For example, circRNA 29761 | 13 had seven variants with the 3′-breakpoints ranging from genomic location 29759 to 29767 and the 5′-breakpoints from 5 to 19 (Table 3) .

Insertion of additional nucleotides between the breakpoints was also observed ( Figure 3D ,E and Table 3 ). This is in sharp contrast to the accurate GU | AG back-splice breakpoints seen in circRNAs reported thus far. It further suggests the mechanism driving CoV RNA backsplicing is error-prone.

Taken together, we experimentally confirmed that the SARS- 

Next, we tried to validate our bioinformatic predictions of SARS-CoV-2 circRNA length and composition with experimental data.

About 30% of sequenced BSJ-spanning amplicons were over 500 nt, indicating that full-length circRNAs were even longer (Table 4 ). This result supports our prediction that SARS-CoV-2 circRNAs are longer than host circRNAs. Although 81.7% of sequenced amplicons carried two fragments separated by a BSJ, we found 16% amplicons with three fragments and 2%, which were even more complex (Table 4) . inverse PCR result also showed the existence of BSJs that were exclusively negative-stranded ( Figure 3J , Arrowheads Ⅰ, Ⅲ

, and Ⅳ).

We concluded that SARS-CoV-2 BSJs can either be strand-specific or exist as both sense and antisense RNA. Our results indicate that bioinformatic prediction of strandness could be unreliable with circRNAs from an RNA genome.

RNase R is a 3′-5′ exoribonuclease that digests all linear RNAs except lariat or circRNA structures. As our experimental data set suggests BSJ-spanning transcripts may not be circularized, we performed RNase R sensitivity assays to determine whether BSJcontaining sgRNAs were truly circular. We first examined the genome-wide resistance of SARS-CoV-2 RNA to RNase R treatment.

Agarose gel electrophoresis of total RNA extracted from SARS-CoV-

showed that ribosomal RNAs were completely degraded after 45 min of RNase R treatment ( Figure 4A ). Northern blotting with a DIGlabeled host β-actin probe confirmed the degradation of actin messenger RNA (mRNA) by RNase R ( Figure 4B ). We further showed with a BSJ-spanning probe targeting SARS-CoV-2 N gene that gRNA and canonical sgRNAs containing the N gene sequence were efficiently removed by RNase R ( Figure 4C ). 41 As no signal was detected in the mock sample with the N BSJ probe ( Figure S4A ), the smear signal around and below canonical sgRNAs bands were likely noncanonical sgRNAs detected by the N BSJ probe ( Figure S4A ). RT-PCR with convergent primers showed that SARS-CoV-2 RNAs were not completely degraded by RNase R, suggesting that some RNA components are resistant to RNase R ( Figure 4E) Figure 1C ) and experimentally identified BSJs (Table 1 ).

Next, we performed inverse RT-PCR on total RNA with or without RNase R treatment. Our result showed that some bands are resistant to RNase R whereas others are susceptible ( Figure 4F ,G and S4B,C). Specifically, bands corresponding to BSJ 29195 | 27789 ( Figure 4E ,F) were resistant to RNase R whereas bands corresponding to BSJ 29122 | 28295 and BSJ 29761 | 13 were more susceptible.

These results suggest that not all BSJ-containing RNAs were circularized.

Lastly, using the BSJ-spanning probe targeting 29122 | 28295, we identified three distinct bands at 0.3, 1.0, and 1.5 kb after RNase R treatment ( Figure 4D ). This is consistent with our bioinformatic prediction that the length of SARS-CoV-2 circRNAs fall into three groups, one below 0.5 kb, one near 1 kb, and one around 1.5 kb ( Figure 2A ). As our inverse Figure 5A , only when the pair of target sequences are juxtaposed, an amplified signal will be produced in AmpFISH. 36 We found that that in addition to the we designed specific juxtaposed primers to amplify only circRNA that is resistant to RNase R.

Repetitive intronic elements, including the primate-specific Alu elements, enables intramolecular RNA looping, thereby promoting cellular circRNA biogenesis in cis. 35, 42, 43 As for RNA recombination in CoVs, the prevailing model predicts that discontinuous RNA synthesis is mediated by homologous motifs called transcription-regulatory sequences (TRSs).

There is a leader TRS (TRS-L) located in the 5′-UTR that is identical or highly homologous to the different body TRSs (TRS-Bs) located in front of The GO Biological Processes analysis showed that SARS-COV-2-circRNAs-miRNAs-upregulated host genes were mainly associated with "muscle tissue development," "ossification," and "response to virus" ( Figure 7B ). GO Cellular Component analysis revealed the enrichment of "apical part of cell" and "apical plasma membrane" genes ( Figure 7C ). Molecular function of the candidate genes fell into the classifications of "cytokine receptor binding" and "growth factor activity" (Figure 7D ). A few cellular genes were downregulated by the viral circRNAs ( Figure 5E ) probably indirectly. In addition, KEGG pathways analysis showed that SARS-COV-2-circRNAs-miRNAsupregulated genes were involved in "Tumor necrosis factor signaling pathway," "Cytokine−cytokine receptor interaction," and "Mitogen- RNA transcription than the current model, which predicts that stemloop structures in TRSs induce RdRp stalling yet lacks details in the mechanism of template-switching. We propose that reverse complementarity bring distant genomic loci to physical proximity, and that local homologous sequences enable template-switching. When the 5′-and the 3′-end of nascent transcript are close enough, RNA circularization occurs. We found that the recently published RNA-RNA interactome of SARS-CoV-2 gRNA and sgRNA highly correlated with our identified short-and long-range BSJ hotspots. 44 The identified local (<1000 nt) RNA-RNA interaction also explained the enrichment of BSJ-and FSJ-spanning reads along the diagonal of junction plots ( Figure 1C-E) . 3, 38 Our work provides insights into the understanding of CoV gene function during viral propagation, immune evasion, and pathogenesis. circ-FBXW7 ( Figure S5A-C) . The SARS-CoV-2 genome exhibits strong affinity to host miRNAs 17 and RNA-binding proteins. 50 CoV circRNAs may act as decoys to indirectly regulate host gene expression. Recent studies have showed that foreign circRNAs activate innate immunity through the nucleic acid senor RIG-I, 51 and that RNA circularization diminishes immunogenicity compared to the linear form. 52 Our future plan is to investigate the biological functions of SARS-CoV-2 circRNAs. Moreover, we will perform a deep circRNA sequencing for SARS-CoV-2 variants, including Alpha, Beta, Delta, Gamma, and Omicron to see whether the circRNAs in SARS-CoV-2 are conserved. 

The authors declare noconflicts of interest.

All data are available in the manuscript or the Supporting information.

The following reagent was deposited by the Centers for Disease Control and Prevention, and was obtained through BEI Resources, NIAID, NIH: SARS-CoV-2, Isolate USA-WA1/2020, NR-52281.

ORCID Shaomin Yang https://orcid.org/0000-0001-7932-4777

COVID-19, SARS and MERS: are they closely related?

Emerging coronaviruses: genome structure, replication, and pathogenesis

The architecture of SARS-CoV-2 transcriptome

SARS-CoV-encoded small RNAs contribute to infection-associated lung pathology

Circular RNA expression: its potential regulation and function

Circular RNAs are abundant, conserved, and associated with ALU repeats

Circular RNA and its mechanisms in disease: from the bench to the clinic

The biogenesis, functions, and challenges of circular RNAs

CircHIPK3 sponges miR-558 to suppress heparanase expression in bladder cancer cells

The landscape of circular

Intrathecal circHIPK3 shRNA alleviates neuropathic pain in diabetic rats

Huachansu suppresses TRPV1 upregulation and spinal astrocyte activation to prevent oxaliplatininduced peripheral neuropathic pain in rats

The Epstein Barr virus circRNAome

Kaposi's sarcoma-associated herpesvirus-encoded circRNAs are expressed in infected tumor tissues and are incorporated into virions

Human papillomavirus 16 E7 oncoprotein alters the expression profiles of circular RNAs in Caski cells

Circular RNAs: new epigenetic signatures in viral infections

Identification and characterization of circRNAs encoded by MERS-CoV, SARS-CoV-1 and SARS-CoV-2

Discovery of SARS-CoV-2 antiviral drugs through large-scale compound repurposing

Competing endogenous RNA network profiling reveals novel host dependency factors required for MERS-CoV propagation

Transcriptomic analysis reveals novel mechanisms of SARS-CoV-2 infection in human lung cells

Fast and accurate short read alignment with Burrows-Wheeler transform

Fast gapped-read alignment with Bowtie 2

Qualimap 2: advanced multi-sample quality control for high-throughput sequencing data

Circular RNA identification based on multiple seed matching

Circular RNAs are a large class of animal RNAs with regulatory potency

Reconstruction of full-length circular RNAs enables isoform-level quantification

STAR: ultrafast universal RNA-seq aligner

featureCounts: an efficient general purpose program for assigning sequence reads to genomic features

Moderated estimation of fold change and dispersion for RNA-seq data with DESeq. 2

MicroRNA targets in Drosophila

clusterProfiler: an R package for comparing biological themes among gene clusters

The Vienna RNA websuite

IRESfinder: identifying RNA internal ribosome entry site in eukaryotic cell using framed kmer features

Targeting mitochondria-located circRNA SCAR alleviates NASH via reducing mROS output

Short intronic repeat sequences facilitate circular RNA production

High-fidelity amplified FISH for the detection and allelic discrimination of single mRNA molecules

Imaging individual mRNA molecules using multiple singly labeled probes

The SARS-CoV-2 subgenome landscape and its novel regulatory features

Accurate quantification of circular RNAs identifies extensive circular isoform switching events

Circular RNA profiling reveals an abundant circHIPK3 that regulates cell growth by sponging multiple miRNAs

SARS-CoV-2 reverse genetics reveals a variable infection gradient in the respiratory tract

Complementary sequence-mediated exon circularization

Analysis of intron sequences reveals hallmarks of circular RNA biogenesis in animals

The short-and long-range RNA-RNA interactome of SARS-CoV-2

Endogenous microRNA sponges: evidence and controversy

The hepatitis delta (delta) virus possesses a circular RNA

New insights about the regulation of Nidovirus subgenomic mRNA synthesis

The natural history of group I introns

Why do RNA viruses recombine?

SARS-CoV-2 contributes to altering the post-transcriptional regulatory networks across human tissues by sponging RNA binding proteins and micro-RNAs

Sensing Self and foreign circular RNAs by intron identity

RNA circularization diminishes immunogenicity and can extend translation duration in vivo

SARS-CoV-2, SARS-CoV, and MERS-CoV encode circular RNAs of spliceosome-independent origin