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

Abstract

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.

Keywords: RNA biology; SARS-CoV-2; circular RNA; coronavirus; virology.

Figure 1

Predicted severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2), SARS‐CoV, and Middle East respiratory syndrome coronavirus (MERS‐CoV) circular RNA (circRNA) abundance and landscapes. (A) Illustration of CIRI2‐based identification of circRNAs. Gapped reads are aligned to the donor and acceptor sequences. If the genomic location of the donor is downstream of the acceptor, it is considered as a back‐splice junction (BSJ). The 5′‐ and 3′‐breakpoints are then determined. (B) Ranked expression level of de novo identified circRNAs from CoV and host genomes. (C–E) Frequency of circularization events in SARS‐CoV‐2 (C), SARS‐CoV (D), and MERS‐CoV (E). Counts of BSJ‐spanning reads (starting from a coordinate in the x axis and ending in a coordinate in the y axis) indicated by color. The counts were aggregated into 500 nt bins for both axes. Distribution of start/end position was shown as histograms on the x and y axis.

Figure 2

Comparison of predicted full‐length coronavirus (CoV) circular RNAs (circRNAs) and host circRNAs. (A) Length distribution of circRNAs derived from host genomes and CoVs. Average length indicated by dashed lines. (B) Strand distribution of host and viral circRNAs. (C) Statistics of forward‐splice junctions (FSJs) in host and viral circRNAs. Number of full‐length circRNAs used in the analysis: severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2), 127; SARS‐CoV, 122; Middle East respiratory syndrome coronavirus (MERS‐CoV), 836; Monkey, 4815; Human, 31,807.

Figure 3

Experimental identification of severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2) circular RNAs (circRNAs) in Vero E6 cells. (A) Illustration of inverse reverse‐transcription polymerase chain reaction (RT‐PCR) with divergent primers would selectively amplify different regions of circRNAs but not linear RNAs. (B) Schematic diagram showing divergent primers were designed to amplify as much predicted back‐splice junctions (BSJs) in a given hotspot as possible. (C–I) Inverse RT‐PCR result with selected primer sets were shown in C. Bands indicated by red arrows were sequenced. Representative Sanger sequencing results were shown in D–I. BSJ breakpoints were indicated by dashed lines. Donor (green) and acceptor (red) sequences and downstream/upstream sequences (gray) flanking the junction were aligned to the BSJ sequence. (J) Strand‐specific RT inverse PCR result. Bands corresponding to BSJ 29122 | 28295 (225 bp), BSJ 29195 | 27789 (804 bp), and BSJ 29761 | 13 (431 and 119 bp) were indicated by red arrowheads.

Figure 4

Severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2) back‐splice junctions (BSJs) are resistant to RNase R treatment. (A) Gel electrophoresis of total RNA with and without RNase R treatment. (B) Northern blotting of host actin mRNA. (C,D) Northern blotting with a BSJ‐spanning N probe (28809–29122 | 28295–28494) under short (C) and long (D) exposure conditions. The size of three circular RNA (circRNA) bands were marked. (E) Reverse‐transcription polymerase chain reaction (RT‐PCR) with convergent primers targeting action and different regions of SARS‐CoV‐2 genome. Divergent primers targeting circHIPK3 was used as positive control for RNase R treatment. (F,G) Inverse RT‐PCR on RNA treated with and without RNase R. Red arrowheads indicate RNase‐R resistant bands.

Figure 5

Localization of severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2) circular RNAs (circRNAs) in infected Vero E6 cells. (A) Design of amplified fluorescence in situ hybridaization (AmpFISH) probes for the detection of a juxtaposed pair of sequences in the circRNA 29122 | 28295. Only when the two sequences are next to each other, which occurs in the circular RNA but not in linear RNA, the donor and acceptor probes can interact with each other and give rise to an amplified hybridization chain reaction (HCR) signal. (B) Multiplex imaging of the linear RNA corresponding to ORF1 obtained by single‐moleclule FISH (smFISH; green) and circRNA 29122 | 28295 (green) obtained from AmpFISH in the same set of cells. Nuclei were stained with 4′,6‐diamidino‐2‐phenylindole (DAPI). RNAse R treatment is expected to degrade linear but not the circular forms of RNAs. (C) Mean fluorescence intensity in single infected cells that were either treated with RNAse R (green) or left untreated (gray). Linear RNA was detected using ORF1‐specific smFISH probe sets (top) and the circRNA was detected by back‐splice junction (BSJ)‐specific AmpFISH probes pair (bottom).

Figure 6

Predicted secondary structure in N gene and the possible formation of back‐splice junction (BSJ) 29122 | 28295. (A) RNA structure prediction for genomic region 28260–29141 of the severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2) N gene. (B) Stem structure is predicted near the BSJ of circular RNA (circRNA) 29122 | 28295. Short homology (CCCAA) immediately upstream of the 5′‐ and 3′‐breakpoints (black arrowheads) are highlighted. (C,D) Proposed model of SARS‐CoV‐2 back‐splicing using BSJ 29122 | 28295 as an example. (C) Strong RNA–RNA basepairing causes RNA‐dependent RNA polymerase (RdRp) to stall. (D) Inverse homology near the base of the stem allows nascent transcript to bypass the stem and switch template.

Figure 7

Potential severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2) circular RNAs (circRNAs) function in the competitive endogenous RNA (ceRNA) coregulatory network. (A) Schematic diagram of potential SARS‐CoV‐2 circRNAs ceRNA coregulatory network analysis method. (B–D) The Gene Ontology (GO) enrichment of SARS‐CoV‐2 circRNAs competitive differentially expressed genes (DEGs). It shows the top 20 significantly enriched GO terms including Biological Process of SARS‐CoV‐2 circRNAs associated upregulated genes (B), cellular component of SARS‐CoV‐2 circRNAs associated upregulated genes (C), and molecular function of SARS‐CoV‐2 circRNAs associated upregulated genes (D), and cellular component of SARS‐CoV‐2 circRNAs associated downregulated genes (E). (F) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of SARS‐CoV‐2 circRNAs associated upregulated genes. (G) Schematic of hsa‐miR‐3194‐5p‐binding sites and the sequences of SARS‐CoV‐2 circRNA 29122 | 28295 and GPR115. (H) Inverse reverse‐transcription polymerase chain reaction (RT‐PCR) with divergent primers targeting SARS‐CoV‐2 circRNA 29122 | 28295 and convergent primers and targeting both linear RNA 28295–29122and circRNA 29122 | 28295. (I) The messenger RNA (mRNA) expression levels of GPR115 using quantitative real‐time PCR. Vector1, pcDNA3.1 (linear RNA‐overexpressed Vector); Vector2, pcD‐ciR (circRNA‐overexpressed Vector); Vector3, pcDNA3.1 CircRNA mini (circRNA‐overexpressed Vector); empty, empty Vector. Two‐tailed unpaired t test. **p < 0.01, ***p < 0.001 circRNA 29122 | 28295 versus empty Vector; ##p < 0.01 pcD‐ciR or pcDNA3.1 CircRNA mini versus pcDNA3.1.