Identification and characterization of virus-encoded circular RNAs in host cells

Analysis workflow

The analysis workflow is outlined in Fig. 1. First, RNase R-treated RNA-Seq samples were collected, and raw reads were preprocessed by a quality control pipeline (Step 1). Then, the circRNA detection tool CIRI2 was used to identify viral circRNAs (Step 2). Finally, the viral circRNAs were annotated, analysed and stored in the Viral-circRNA database.

RNA-Seq

Total RNA-Seq or RNase R-treated RNA-Seq data are suitable for circRNA identification. To control for false positives, only RNase R-treated RNA-Seq samples were taken to identify viral circRNAs. Virus-infected and RNase R-treated RNA-Seq data were downloaded from the NCBI GEO [39] and covered six viruses, EBV, KSHV, rLCV, Middle East respiratory syndrome coronavirus (MERS-CoV), HPV and hepatitis B virus (HBV) (Table S1). Six poly(A)-selected tumour datasets infected by viruses included in our study (Table S2) were downloaded from the GEO database for false discovery rate (FDR) estimation.

Additionally, three RNase R-treated RNA-Seq samples, including two HPV-infected samples from patients with cervical cancer and one HBV-infected liver cancer cell line, Hep3B, were generated by our own laboratory. Human cervical cancer samples were collected from Tongji University Shanghai East Hospital. The raw and processed sequencing data were submitted to the SRA database (https://www.ncbi.nlm.nih.gov/sra) [40] under accession number PRJNA756819.

RNA quantification, sample preparation, library construction and Illumina sequencing

RNA degradation and contamination were monitored on 1 % agarose gels. RNA purity was checked using a NanoPhotometer spectrophotometer (IMPLEN). RNA concentration was measured using a Qubit RNA Assay Kit in a Qubit 2.0 Fluorometer (Life Technologies), and RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies). For circRNA sequencing, ribosome-depleted and RNase R-treated RNA samples were used as input material for the RNA sample preparations. Sequencing libraries were generated using the NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB) following the manufacturer’s recommendations. Clustering of the index-coded samples was performed on a cBot Cluster Generation System using a TruSeq PE Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina HiSeq platform, and 125 bp/150 bp paired-end reads were generated.

Ninety-six severe acute respiratory syndrome coronavirus-2 (SARS-CoV-2)-infected total or rRNA-depleted RNA-Seq samples across nine datasets (Table S3) were downloaded from the GEO database for additional analysis.

Processing of raw RNA-Seq data

The quality of sequencing reads was evaluated with FastQC v0.11.6 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) [41]. Raw reads from the RNA-Seq data were processed by removing adapters and low-quality bases using FASTQ v0.21.0 [42] with default parameters and by reducing sequence duplication using PRINSEQ-lite v0.20.4 [43] with derep=1.

Identification of viral circRNAs

The clean reads were then aligned to a reference genome of the virus with BWA-MEM using parameter -T 19 recommended by the CIRI2 manual. CIRI2 v2.0.6 [36] was applied with default parameters to identify viral circRNAs. The viral circRNA with at least two junction reads in its back-splicing site was initially identified. After basic statistics, to control for false positives, circRNAs with at least five junction reads were selected for further analysis. CircRNAs with the same back-splicing site identified from different samples were integrated as one circRNA.

Characterization of viral circRNAs

In the next step, these viral circRNAs were analysed in terms of their host genes, their length, transcription strand, the number of junction reads and their expression levels. The expression level of a viral circRNA was measured by the natural logarithm of counts per million mapped reads (CPM), which was calculated as follows: $CPM={\left(\raisebox{1ex}{$n$}\!\left/ \!\raisebox{-1ex}{$N$}\right.\right)\times 10}^6$ , where n and N represent the number of junction reads and the number of mapped reads, respectively.

Function analysis of viral circRNAs

CircRNAs with junction reads ≥5 and with lengths ≤1 kb and human miRNAs from miRBase [44] were collected to predict miRNA–circRNA interactions using miRanda v3.3a [45] and TarPmiR [46]. To control for false positives, strict parameter settings were used for the prediction, and only interactions identified by both tools were accepted for further analysis. The parameters for miRanda were Max Score ≥160 and Max Energy≤−20 kcal mol^–1, and for TarPmiR, they were Binding Probability=1 and Energy≤−20 kcal mol^–1.

To measure the quality of the predicted miRNA–circRNA interactions, the outputs of miRanda including two metrics – Max Score and Max Energy – as well as the output of TarPmiR – Energy – were integrated to calculate a combined quality score. Max Score, Max Energy and Energy were respectively scaled to between 0 and 1 by min-max normalization. Because lower energy means a more thermodynamically stable binding, Max Energy in miRanda and Energy in TarPmiR were inverted before scaling. The combined quality score was defined as the average of the scaled outputs of miRanda and TarPmiR, with higher scores indicating higher quality. Based on miRNA–circRNA interactions, the Cytoscape 3.9.1 software [47] was used to build an interaction network in which the width of edges was determined by the combined quality score.

Furthermore, miRNA–mRNA interactions were retrieved from the miRNA target databases TargetScan [48], miRDB [49] and miRTarBase [50]. Finally, Gene Ontology (GO) [51] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [52] pathway enrichment analyses were conducted on the overlapping miRNA targets shared by the three databases using the R package clusterProfiler v3.18.1 [53]. All the KEGG pathways and GO terms with q-values <0.05 were considered as significant enrichment.

Graphical visualization

The R package circlize v0.4.12 [54] were used to visualize viral circRNAs. Viral circRNAs with abundant isoforms were plotted by SpliceV [55]. The coverage curve of circRNAs on the viral reference genomes were generated by Gviz v1.34.1 [56].

Development of the Viral-circRNA Database

The Viral-circRNA Database (http://www.hywanglab.cn/vcRNAdb/) was built on the Django (version 2.2.7) framework and run on the Apache 2 web server with MySQL (version 5.7.24) as the database engine. The JavaScript libraries Echarts and Highcharts were used for drawing graphs. The database is available online without registration and was optimized for most prevalent browsers, such as Chrome.

Statistical analysis

The Benjamini and Hochberg procedure was used to calculate adjusted P-values in GO and KEGG enrichment analyses. The Wilcoxon rank-sum test was used to compare the ln(CPM) values between long and non-long circRNAs. All statistical analyses were executed in R.

Data collection and identification of viral circRNAs

To identify circRNAs of viral origin in host cells, we interrogated RNase R-treated RNA-Seq samples from virus-infected cell lines and patients with lymphoma, lung adenocarcinoma, cervical carcinoma and gastric carcinoma (Figs. 1 and 2a) from the NCBI GEO/SRA database. Among them, 69 % (27 of 39) of samples were from lymphoma patients, including the following subtypes: rhesus macaque lymphoma, primary effusion lymphoma, Burkitt’s lymphoma and B cell lymphoma. In addition, three supplementary RNase R-treated RNA-Seq samples were generated by our own laboratory, including two HPV-infected samples from patients with cervical cancer and one HBV-infected liver cancer cell line, Hep3B. In total, our analysis included 41 cancer samples (Table S1) infected with six types of viruses: KSHV, EBV, rLCV, MERS-CoV, HBV and HPV.

Raw reads from the RNA-Seq data were processed by removing low-quality reads and adaptor dimers and by controlling sequencing duplication. After clean reads were aligned to the viral genome, the circRNA detection tool CIRI2 was used to identify 3912 viral circRNAs (Fig. 1). Among them, 3198 were from MERS-CoV, followed by 590 from KSHV, 102 from EBV, 19 from rLCV and three from HPV18 (Fig. 2b). To estimate the FDR, we collected the poly(A)-selected RNA-Seq datasets with at least one dataset for each cancer type infected with the viruses in our study (Table S2) and identified the circRNAs with the tool CIRI2. CircRNAs detected in both poly(A)-selected and RNase R-treated samples were considered to be false positives. Therefore, the estimated FDR was 0.0043 (17 of 3912).

Lengths and expression levels of viral circRNAs

The number of junction reads of circRNAs in each sample varied from two to 34 232, and 70 % (4139 of 5924) of the circRNAs had at least five junction reads (Fig. 2c), suggesting their reliability as true positives. Based on the number of junction reads, the expression levels of viral circRNAs were measured as CPM. Of note, the expression levels of 45 % (2694 of 5924) of the circRNAs were greater than 1, and 12 % (714 of 5924) reached up to 10, showing that a small portion of circRNAs had relatively high expression levels (Fig. 2d). For circRNAs with the same back-splicing site identified from different samples, we integrated them as one circRNA. The statistics also showed that 50 % (1963 of 3912) of circRNAs ranged from 200 to 500 bp in length, which was in line with the characteristics of circRNAs. Interestingly, 8 % (311 of 3912) of circRNAs had almost head-to-tail joining of viral genomes with lengths ≥10 kb (Fig. 2e). Among 3912 viral circRNAs, 2602 were only identified in one sample, 738 across two samples, 515 across three samples, 36 across four samples and 21 across at least five samples, suggesting that the circRNAs identified across more samples are more likely to be true positives (Fig. 2f). However, the circRNAs only identified in one sample were not necessarily false positives due to the limited sample size.

To further control for false positives, we selected circRNAs with at least five junction reads covering the back-splicing sites as high-confidence viral circRNAs for the following analysis. As a result, the estimated FDR of these circRNAs declined to 0.0023 (six of 2612). We created circular diagrams to display back-splicing sites (inner circle) and coverage of junction reads (middle circle) collectively on the viral genome (outer circle) (Fig. 2g–j). In the inner circle, the starting position and ending position of a circRNA are connected by a curve, with red and blue representing the positive strand and negative strand, and darker colours implying higher expression. EBV circRNAs were ubiquitously highly expressed (Fig. 2g). In the middle circle, circRNAs expressed from the positive strand and negative strand are shown in purple and blue, respectively, with the y-axis representing the coverage of junction reads. In EBV (Fig. 2g), KSHV (Fig. 2h) and rLCV (Fig. 2j), the circRNAs were primarily expressed on the positive strand. However, in MERS-CoV (Fig. 2i), more circRNAs were expressed on the negative strand. Moreover, in EBV, MERS-CoV and rLCV (Fig. 2g, i, j), interestingly, there were some arrow-pointing circRNAs with back-splicing sites almost joining across the entire genome in a head-to-tail manner, demonstrating the long-range cyclization of genomic RNAs.

Host genes of viral circRNAs

To observe the host genes from which viral circRNAs were transcribed, we annotated the circRNAs (junction reads ≥5) by connecting their genomic positions to loci of viral genes (Fig. 3). The results showed that one viral gene could produce multiple circRNA isoforms, with the viral genes RPMS1/A73/BARF0 and EBNA in EBV (Fig. 3a), K7 and ORF68/ORF69 in KSHV (Fig. 3b), orf1ab in MERS-CoV (Fig. 3c) and RPMS1 in rLCV (Fig. 3d) expressing the most circRNAs. For example, in EBV (Fig. 3a, e), 21 circRNA isoforms were transcribed from RPMS1, followed by 13 from A73, 12 from BARF0 and 10 from EBNA. However, because A73 and BARF0 have overlapping genome loci with RPMS1, some circRNAs are shared by these genes. In KSHV, 124 circRNA isoforms were transcribed from K7 (Figs 3b and S1a). Similarly, we observed an obvious peak at K7 (Fig. 2h), indicating abundant expression of K7 circRNAs in KSHV. Moreover, 459 circRNAs were transcribed from orf1ab, a gene covering approximately 60 % of the genome in MERS-CoV (Figs 3c and S1b). Two types of gammaherpesviruses, RPMS1 and rLCV, share a highly conserved transcript repertoire [57]. Remarkably, RPMS1 in EBV produced the most circRNA isoforms, and its orthologue in rLCV produced 10 different isoforms (Figs 3d and S1c). The circRNAs expressed from EBV RPMS1 have been reported to be relevant to tumour proliferation, apoptosis, invasion and metastasis [58]. Our combined results were in line with previous reports on viral circRNAs, such as circRNAs transcribed from RPMS1 [22, 26, 28], EBNA [26], BHRF1 [26], LMP2 [22, 26] and BHLF1 [26] in EBV; K7 [22, 25] and vIRF-4 [22, 23, 25] in KSHV; and RPMS1 and EBNA in rLCV [23]. These circRNAs are worthy of further investigation to explore their potential functions in cancers.

Gammaherpesviruses such as EBV are involved in conversion of the infection state between latent and lytic phases. During much of its lifecycle, EBV remains silenced in the latent state to escape the host’s immune system, with only latent genes being expressed. During lytic infection, EBV switches to an active state, leading to extensive gene expression with the production of infectious viruses [59]. To evade immune recognition, EBV establishes different latent gene expression patterns consisting of EBNA, LMP, BART and an untranslated RNA called EBER [60]. As most of the EBV samples in our analysis were in a latent state [26], we observed that the major EBV circRNAs were transcribed from latent genes, such as EBNA, LMP and RPMS1/A73/BARF0 (Fig. 3a), indicating that these circRNAs might mediate immune escape similar to their linear transcripts. By comparison, most KSHV samples were in a lytic state, so KSHV circRNAs were primarily transcribed from lytic genes such as K7, ORF68 and ORF69, and only a small proportion of them were from latent genes, including ORF71, ORF72, ORF73 and K12 [61] (Fig. 3b). Significantly, K7, as the gene producing the most isoforms in KSHV, shares its locus with a long non-coding polyadenylated nuclear (PAN) RNA, which is the most abundant viral transcript during lytic replication [62]. Together, the expression of viral circRNAs is dependent on the infection state when host genes are activated. Therefore, as a supplement to linear transcripts, viral circRNAs may regulate the molecules of both viruses and infected cells because of their circular conformation and complementary sequence with linear RNAs to finally escape antiviral sensing of the host immune system in the latency state and manufacture more viruses in the lytic state.

Viral circRNAs as miRNA sponges

It is well known that endogenous circRNAs can function as miRNA sponges in eukaryotes [18]. To observe whether viral circRNAs can act as sponges to bind the miRNAs of the infected host cells, we used the miRNA target prediction tools miRanda (v3.3) [45] and TarPmiR [46] to predict the interactions between circRNAs and host miRNAs. We implemented this analysis on the circRNAs (length <1 kb and junction reads ≥5) and host miRNAs from the comprehensive miRNA database miRBase [44]. Based on the strict parameter settings (see Methods), our analysis identified 89 miRNA–circRNA pairs with 21 unique human miRNAs that were predicted by both tools.

To identify the important miRNAs with strong interaction to viral circRNAs, we integrated the outputs of miRanda including two metrics – Max Score and Max Energy – as well as the output of TarPmiR – Energy – to calculate a combined quality score of the interaction. Then we built a network in which wide edges represented a high quality score, and we focused on those hub miRNAs with multiple connections to circRNAs. The results showed that hsa-miR-6848–5p interacted with 21 KSHV circRNAs (Fig. S2a); hsa-miR-3085–5p (Fig. S2b), hsa-miR-6858–5p (Fig. S2c) and hsa-miR-8063 (Fig. S2d) respectively interacted with 18, 10 and nine MERS-CoV circRNAs; and hsa-miR-6747–5p interacted with four MERS-CoV circRNAs and one KSHV circRNA (Fig. S2e). More importantly, five MERS-CoV circRNAs interacting with hsa-miR-6747–5p had very high quality scores showing wide edges in the sub-network (Fig. S2e), providing more solid evidence for the potential role of hsa-miR-6747–5p in viral circRNA regulation. Other miRNAs with more than two interactions were shown in Fig. S2(f–k). According to the literature, hsa-miR-8063 was found to inhibit pancreatic cancer progression through inducing inhibition on ZFP91 [63], to suppress the epithelial–mesenchymal transition of lung cancer cells by regulating its target gene MAGAT3 [64], and to induce cell migration through inhibition of the tumour suppressor FOXN2 [65]. Hsa-mir-4745–5p was strongly dysregulated in Parkinson’s disease (PD) and can be directly regulated by HNF4a [66], which was identified as a longitudinal PD blood biomarker [67]. Hsa-mir-4745–5p was also observed to be upregulated in patients with colorectal cancer [68]. These all suggest there are potential interactions between host miRNAs and viral circRNAs, and some miRNAs deserve further investigation for their potential role in the regulation of viral pathology.

To further explore the functions of viral circRNAs, we used three miRNA target databases, TargetScan [48], miRDB [49] and miRTarBase [50], to retrieve the target genes of the 21 miRNAs above. This analysis led to the identification of 638 miRNA–mRNA interactions and 604 target genes. Enrichment of target genes based on KEGG pathway and GO analyses showed significantly enriched pathways and functions (Fig. 4). The top 10 enriched KEGG pathways showed a strong correlation with viruses and cancer (Fig. 4a), suggesting that viral circRNAs in tumour samples participate in viral infection and tumorigenesis. GO-enriched terms provided clues to how viral circRNAs play regulatory roles in the infected host cells. Viral circRNAs may regulate the cell cycle and protein structure and localization (Fig. 4b) through transcription repressor activity (Fig. 4c). In addition, cellular component terms, including nuclear envelope and nuclear speck, were enriched (Fig. 4d). Overall, the enriched functions were in line with the roles of viruses, indicating that some viral circRNAs function as host miRNA sponges to regulate infected host cells.

Negative strand bias of viral circRNAs in coronavirus

There were 3198 viral circRNAs identified in MERS-CoV, accounting for 82 % of the total in this study (Fig. 2b). These circRNAs were identified from one study in which the human lung cancer cell line Calu-3 was infected with MERS-CoV for 6 and 24 h in vitro, followed by RNase R-treated RNA-Seq. Our analysis identified dozens of circRNAs in cell lines at 6 h post-infection and thousands of circRNAs at 24 h post-infection. This result suggested that fast replication of the virus in the host cells was accompanied by the production of abundant viral circRNAs. Therefore, we considered that this cell line provided a good system to study viral circRNAs post-infection.

To gain more insight into the characteristics of viral circRNAs, we extracted viral circRNAs from MERS-CoV and observed the distributions of junction reads, expression levels and lengths (Fig. S3a–c), which were similar to those identified from the pooled viruses (Fig. 2c–e). Then, we focused on MERS-CoV circRNAs (junction reads ≥5) as strict false-positive controls. Interestingly, when mapping the junction reads of viral circRNAs onto positive-sense (+) and negative-sense (−) strands of the MERS-CoV genome (Fig. 4e), we found that circRNAs were widely expressed across the whole genome and were highly expressed in some regions of orf1ab, S, orf3/orf4a/orf4b and N/orf8b. Surprisingly, we observed that there were twice as many circRNAs on the negative-sense (−) strand as on the positive-sense (+) strand (Fig. 4e). As a positive-sense ssRNA virus, MERS-CoV utilizes full-length replication or discontinuous transcription to generate negative-sense RNA genomes or subgenomic negative-sense RNAs, which serve as a template for the synthesis of positive-strand RNAs rather than encoding proteins [69]. It is unclear whether these negative-sense circRNAs play regulatory roles in the process of positive-sense RNA synthesis.

Since these viral circRNAs were from in vitro experiments, technical factors and biological factors may also impact the process of RNA cyclization. To observe whether this negative-strand bias truly exists in positive-sense ssRNA viruses, we collected 96 SARS-CoV-2-infected total or rRNA-depleted RNA-Seq samples across nine datasets, instead of collecting the few RNase R-treated samples (Table S3), and identified 376 circRNAs. To control for false positives caused by RNA-Seq data without RNase R treatment, we narrowed down the circRNAs to those with junction reads ≥5, leading to 140 circRNAs for further analysis. The results showed that a negative-strand bias was still observed for viral circRNAs identified in SARS-CoV-2 (Fig. 4f). Together, these findings in MERS-CoV and SARS-CoV-2 reveal that positive-sense ssRNA viruses produce more abundant circRNAs encoded from the negative strand.

Long-range cyclization of the viral genome and subgenome

Although the lengths of circRNAs are hundreds of bases in most cases, in our analysis, after selecting circRNAs with at least five junction reads, there were still 9 % (191 of 2199) of the circRNAs in MERS-CoV and 5 % (two of 43) of those in EBV that were longer than half of the genome. This result suggested an interaction between 5′-terminal sequences and 3′-terminal sequences, causing a crosslink between two ends of the genome or subgenome, which leads to long-range cyclization.

The cyclization positions of all circRNAs in EBV and MERS-CoV are shown in Fig. 5(a, b). With the short-length grey circRNAs near the diagonal as the background, orange and red circRNAs located at the upper left corner were from long-range cyclization, covering ≥50 % and ≥90 % of the genome, respectively. We compared the expression levels between RNAs from long-range cyclization and short-length circRNAs in EBV (Wilcoxon rank-sum test P=0.018, Fig. 5c) and in MERS-CoV (Wilcoxon rank-sum test P<2×10⁻¹⁶, Fig. 5d) and observed a significantly higher expression of long circRNAs in MERS-CoV (Fig. 5d). Long-range cyclization was also observed in the SARS-CoV-2 samples described above (Fig. 5e and Wilcoxon rank-sum test P=0.049, Fig. 5f).

A database for viral circRNA in host cells

To enable knowledge sharing and reuse, we developed the Viral-circRNA Database (http://www.hywanglab.cn/vcRNAdb/, Fig. S4), which provides comprehensive knowledge of viral circRNAs identified in cancer samples to support further exploration of their physiological and pathological roles in host cells. The web interface of Viral-circRNA is user-friendly and allows users to search (Fig. S4a), browse (Fig. S4b) and download data. The database includes the following information on viral circRNAs: virus origin, reference genome, starting position, ending position, length, expression level, host gene, cancer type, RNA-Seq sample and visualization of the back-splicing position (Fig. S4c). In addition, an overall view of circRNAs produced from a given virus (Fig. S4a) and a data summary and basic statistics of the database are available (Fig. S4d). Links to external databases are also provided.

We performed a side-by-side comparison with VirusCircBase [70], a published viral circRNA database, from the perspective of sample type, RNA-Seq library, detection tool, etc. (Table S4). There were 3247 overlapped circRNAs which occupy 83 % (3247 of 3912) and 7 % (3247 of 46440) of the circRNA collection in our database and VirusCircBase database, respectively. Our analysis only focuses on tumour samples and carefully controls the false positives, by just taking advantage of RNase R-treated samples and employing the circRNA detection tool CIRI2, which has been shown to achieve excellent performance and is appropriate for viral genome detection. Therefore, our database provides less but more reliable circRNAs for promoting research in virus-infected tumours. More importantly, we provide several fields, including junction reads, supported samples, CPM, average junction reads across samples, etc. Users can filter the circRNAs based on their own criteria and download the data.