Epstein-Barr virus circRNAome as host miRNA sponge regulates virus infection, cell cycle, and oncogenesis

ABSTRACT

Epstein-Barr virus (EBV) is an oncogenic virus that infects more than 90% of the world’s population. The proteins and miRNAs encoded by EBV are involved in multiple human malignancies. Recently R-resistance RNA-seq demonstrated that EBV-encoded circular RNAs. The current research aims to explore their functions in EBV-associated malignancies. Total 56 miRNAs were sponged by circRNAome. 24 and 9 in EBV host B and epithelial cells out of 56 miRNAs were detectable by miRNA-seq. 18 and 5 miRNAs were down-regulated in both types of host cells, respectively, after EBV infection. The network between five miRNAs and their targets included 1414 genes, 1419 nodes, and 2423 edges. These targets were enriched in multiple categories, and most of them were up-regulated in EBV-infected cells. These data represented the first report that EBV circRNAs could sponge the miRNAs to promote the up-regulated expression of their targets, involving in malignancies associated with EBV.

Graphical abstract

KEYWORDS:

• Circrna
• EBV
• cancer
• miRNA
• herpes virus
• cell cycle
• oncogenesis

Introduction

Epstein-Barr virus (human herpesvirus 4, EBV/HHV4) is a human oncogenic gamma herpesvirus, which causes a variety of human malignancies including lymphoid and nonlymphoid malignancies [1–3]. The lymphoid malignancies associated with EBV include Burkitt’s lymphoma, Hodgkin’s disease and non-Hodgkin’s lymphoma, and T-cell lymphoma [3]. In addition, EBV is also associated with nasopharyngeal carcinoma, leiomyocarcomas, and gastric carcinoma [4–7]. EBV may also cause lymphoproliferative disorders, such as posttransplant lymphoproliferative disease (PTLD) in solid organ or hematopoietic transplant patients [8,9]. Although studies indicated that EBV-encoded proteins and miRNAs played critical roles in the EBV associated with oncogenesis [10], the mechanism is still not fully understood. Recently the R-resistance RNA-sequencing identified the circRNAs encoded by EBV genome in EBV positive PTLD, gastric cancer, and different latently infected B cell lines [11,12]. These circRNAs were related to EBV-encoded oncogenic genes such as LMP2, RPMS1, and BART [11,12]. Up to now, the functions of these circRNAs encoded by EBV genome have not been studied yet. Whether these circRNAs are involved in the oncogenesis mediated by EBV infection needs to be investigated.

The circRNAs are generated by covalently backsplicing of a downstream 5′ splice donor site onto an upstream 3′ splice acceptor site [13,14]. Currently, circRNAs have been demonstrated to exert their functions by acting as miRNA sponges [15], competitively inhibiting their corresponding linear mRNAs [16], and binding to the RNA binding proteins [17]. Additionally, it has been found that some of circRNAs could be translated into proteins using cap-independent mechanism [18–20]. Considering miRNAs regulate multiple aspects of cancers including cell growth, invasion, angiogenesis, metastasis, and immune evasion [21,22], it is hypothesized that circRNAome composed of circRNAs encoded by EBV genome would sponge the miRNAs of host cells to regulate their behaviors, and promote the progression of EBV-associated malignancies. By creatively combining bioinformatics prediction, RNA-seq analysis, and database searching, it has been found that EBV infection caused the decrease in expression of several host miRNAs, and the increase in expression of their targets. These targets are enriched in multiple categories associated with EBV related cancers such as herpes virus infection, cell cycle, and cancer pathway. This first evidence indicated that EBV circRNAome was involved in the EBV-associated malignancies by sponging miRNAs of host cells.

Methods and materials

Identification of miRNAs which interacted with circrnaome of EBV

According to the genome sequence of Mutu-EBV, the circRNAs sequences were acquired, which were identified by R-resistance RNA sequencing [11]. The sequences were inputted to the miRDB database to identify the binding sites on circRNAs for miRNAs using MirTarget algorithm [23,24]. The binding sites for all human miRNAs were predicted and scored. The scores were weighed according to the numbers of miRNA binding sites and the abundance of circRNAs in circRNAome.

Database extraction

miRNA expression data were extracted from miRmine database which collected the expression profiles from different publicly available miRNA-seq datasets [25]. These included miRNAs from LCLd3 cells infected with or without EBV infection, and normal HEK293T and HEK293 cells. The original data of the RNA sequencing were downloaded from SRA database which stores raw sequencing data for research community to enhance reproducibility and allows for new discoveries by comparing data sets [26,27].

Rna-seq data analysis

The RNA-seq data for mock (SRR5309459, https://www.ncbi.nlm.nih.gov/sra?term=SRX2609768) and EBV infection (SRR5309460, https://www.ncbi.nlm.nih.gov/sra?term=SRX2609769) were downloaded from SRA database. The original reads were trimmed using FASTX-toolkit (V 0.0.13, http://hannonlab.cshl.edu/fastx_toolkit/index.html) and aligned to reference genome using TopHat2 (V2.1.1, https://ccb.jhu.edu/software/tophat/index.shtml). The differential expression analysis was performed using Cuffdiff (V2.2.1, http://cole-trapnell-lab.github.io/cufflinks/cuffdiff/). The different data were used to evaluate the alteration of the target expressions of miRNAs.

Network establishment

The network between miRNAs, of which the expression were down-regulated by EBV infection, suggesting these miRNAs were sponged by circRNAome of EBV, and their targets were established by using miRNet database, which integrated 11 miRNA databases including miRTarBase, TarBase, miRecords, miRDB, and miRanda, etc [28]. These targets were determined by high throughput sequencing methods such as CLASH, PAR-CLIP, and HITS-CLIP which were designed to map the miRNA interactome. The targets were then employed to analyze the enrichment in the pathways based on GO or KEGG database by hypergeometric test algorithm. The statistical value was calculated by the embedded statistical method of miRNet.

Statistical analysis

The up-regulated targets of miRNAs after EBV infection were acquired from the RNA-seq data with or without EBV infection. The data were analyzed using Fisher exact test in R (V3.6.1, https://www.r-project.org/). When P value was less than 0.05, it was statistically significant.

Results

Circrnaome of EBV as a miRNA sponge

Recently two independent research groups identified the circRNAs from EBV genome by R-resistance RNA sequencing [11,12], and their location was shown in Figure 1a, but the functions of this circRNAs were not known yet. In the current research, we explored the possibility that the circRNAome constituted by all circRNAs from EBV transcripts as a sponge negatively regulated the functions of host cell miRNAs. Using miRDB screening, we identified 56 potential miRNAs in human which could interact with cicrRNAome of EBV (Figure 1b). They were ranked using the weighed score according to their binding ability and binding sites, and the abundance of circRNAs which interacted with the corresponding miRNAs (Figure 1b). 24 out of 56 miRNAs could be detected in LCLd3 cells, a B cell line (Figure 2a). 18 out of 24 miRNAs have decreased expression after EBV infection (Figure 2a and Table 1). Table 1 lists all the decreased miRNA and their scores and decreased folds. The decreased fold of miRNAs was positively correlated with the weighed score (Figure 2b). This data indicated that circRNAome of EBV could inhibit the miRNA functions of host cells as a sponge.

Table 1. Down-regulated miRNAs after EBV infection in LCLd3 cells.

miR-Name	Mature-miR-sequence	Total score	Decrease (fold)
hsa-miR-4532	CCCCGGGGAGCCCGGCG	218.02	18.00
hsa-miR-1226-3p	UCACCAGCCCUGUGUUCCCUAG	140.76	11.00
hsa-miR-3195	CGCGCCGGGCCCGGGUU	125.81	10.50
hsa-miR-4749-5p	UGCGGGGACAGGCCAGGGCAUC	117.18	8.00
hsa-miR-4706	AGCGGGGAGGAAGUGGGCGCUGCUU	117.18	7.00
hsa-miR-4673	UCCAGGCAGGAGCCGGACUGGA	48.38	5.00
hsa-miR-4286	ACCCCACUCCUGGUACC	37.2	4.91
hsa-miR-4690-3p	GCAGCCCAGCUGAGGCCUCUG	35.34	4.00
hsa-miR-3667-5p	AAAGACCCAUUGAGGAGAAGGU	8.8	3.00
hsa-miR-4632-5p	GAGGGCAGCGUGGGUGUGGCGGA	25.01	3.00
hsa-miR-4447	GGUGGGGGCUGUUGUUU	10.58	2.00
hsa-miR-4505	AGGCUGGGCUGGGACGGA	10.26	2.00
hsa-miR-15b-5p	UAGCAGCACAUCAUGGUUUACA	3.72	1.89
hsa-miR-424-5p	CAGCAGCAAUUCAUGUUUUGAA	3.72	1.85
hsa-miR-6836-3p	AAUUCCCUUGUAGAUAACCCGG	5.2	1.75
hsa-miR-30c-2-3p	CUGGGAGAAGGCUGUUUACUCU	3	1.61
hsa-miR-4521	GCUAAGGAAGUCCUGUGCUCAG	1.86	1.53
hsa-miR-30c-1-3p	CUGGGAGAGGGUUGUUUACUCC	3	1.47

The scores of miRNAs were dependent on the expression enrichment and binding numbers with circRNAs. The decrease fold was calculated according to the counts of miRNAs in infected or un-infected LCLd3 cells.

Figure 1. The potential miRNAs interacted with circRNAome of EBV. a. The circRNA location in the EBV genome. b. The potential miRNAs that interacted with circRNAome of EBV. The miRNAs were identified using miRDB database screening. The score was the weighed score according to the binding ability, binding sites of miRNAs and expression abundance of circRNAs in R-resistance RNA-sequencing data.

Figure 2. EBV infection decreased expression of miRNAs which were targeted by circRNAome of EBV. a. EBV infection decreased the expression of miRNAs that interacted with circRNAome of EBV. The miRNA expression in LCLd3 cell with or without EBV infection was from miRmine database. P-value was calculated using Fisher exact test in R. b. The positive correlation was between the decreased fold of miRNAs upon EBV infection and their scores targeted by EBV circRNAome.

The down-regulated miRNAs by EBV infection in lcld3 cell are highly expressed in epithelial cells, another type of EBV host cells

In addition to B cells, EBV could also infect epithelial cells [29]. HEK293T and HEK293 cells, the epithelial cells derived from human kidney, were usually used as host cells to investigate EBV entry [30,31]. The miRNA expression data of HEK293T and HEK293 were acquired from miRmine [25], showing that 9 out of 56 predicted miRNAs were detected in HEK293T or HEK293, respectively (Figure 3). Interestingly, the down-regulated miRNAs such as has-miR-15b-5P, has-miR-4286 by EBV infection had the highest expression (red bars, Figure 3). This suggested that the sponge function of EBV circRNAome would essentially impact on the heavers of host cells.

Figure 3. The expression of miRNAs which potentially interacted with EBV circRNAome in HEK293T (a) and HEK293 (b) cells. The miRNA expression in HEK293T and HEK293 was from miRmine database. The red, purple and blue bars indicated miRNAs expression level that was decreased, increased or not altered after EBV infection.

The network between miRNAs sponged by EBV circrnaome and their targets

Considering the high expression level of miR-15b-5p, miR-30c-1-3p, miR-30c-2-3p, miR-424-5p, and miR-4286 in both HEK 293T and HEK293 cells (Figure 3), and their expression is decreased after EBV infection (Figure 2), they were used to establish the network with their targets. Using miRNA database [28], five miRNAs with their targets composed a comprehensive network with 1414 targets, 1419 nodes and 2423 edges (Figure 4a). The enrichment analysis indicated that their targets are enriched in cell cycle, cancer pathway and virus infection, etc., in KEGG pathway (Figure 4b). This was consistent with the fact that EBV infection was associated with multiple human cancers [1–3,5,9].

Figure 4. The network of the top five-decreased miRNAs associated with EBV circRNAome and their targets (a) and the enrichment of their targets (b). The network was generated using miRNet. The network was composed of 1414 targets, 1419 nodes, and 2423 edges. The enrichment assay indicated that the targets were enriched at herpes virus infection, cell cycle, cancer pathway, etc.

Up-regulated expression of miRNAs’ targets is associated with herpes infection, cell cycle and cancer pathway in host cells after EBV infection

The subset networks (Figure 5a–7a) associated with herpes virus infection, cell cycle and cancer pathway were extracted from the entire network (Figure 4a). The expression of 12 out of 17 targets related to herpes virus infection was increased in EBV-infected cells, while 4 was decreased, and 1 was not altered (Figure 5b). Fisher exact test indicated that there was significant difference between the increased targets and others (p = 0.0381). The targets with twofold increase in expression were indicated in Figure 5c, including TAB2, IFNGR2, etc. The expression of 21 out of 31 targets linked to cell cycle was increased (Figure 6b), while 7 was decreased, 2 was not altered, and 1 was not detectable. Fisher exact test indicated there was significant difference between the increased targets and others (p = 0.0092). The expression of 15 out of 21 were more than twofold (Figure 6c) in EBV-infected cells, including CCNE1/2, CDCD25A, E2F3, etc. There were 51 targets associated with cancer pathway and the expression of 31 targets was increased (Figure 7b) after EBV infection, while 12 was decreased, 3 was not altered, and 5 was not detectable. Fisher exact test indicated that the increased targets were significantly different from others (p = 0.0016). The expression of 19 targets was greater than twofold (Figure 7c), including the well-known important genes such as CDK4/6, TGFB1, MAPK1.

Figure 5. Up-regulated expression of the top five-decreased was composed of miRNA targets are associated with herpes virus infection. a. The network of miRNAs and their targets which were enriched in herpes virus infection. The yellow circles indicated the genes of which expression were increased more than two-fold after EBV infection. b. The expression alteration of targets associated with herpes virus infection after EBV infection. P-value was calculated using Fisher exact test in R. c. The genes of which expression were increased more than twofold after EBV infection. The gene expression data were calculated using RNA-Seq data with or without EBV infection from SRA.

Figure 6. Upregulated expression of the top five-decreased miRNA targets associated with cell cycle. a. The network of miRNAs and their targets which were enriched in cell cycle. The yellow cycles indicated the genes of which expression were increased more than two-fold after EBV infection. b. The expression alteration of target genes associated with cell cycle after EBV infection. P-value was calculated using Fisher exact test in R. c. The genes of which expression were increased more than twofold after EBV infection. The gene expression data were calculated using RNA-Seq data with or without EBV infection from SRA.

Figure 7. Upregulated expression of the top five-decreased miRNA targets associated with oncogenesis. a. The network of miRNAs and their targets which were enriched in cancer pathway. The yellow cycles indicated the genes of which expression were increased more than two-fold after EBV infection. b. The expression alteration of targets associated cell cycle after EBV infection. P-value was calculated using Fisher exact test in R. C. The genes of which expression were increased more than twofold after EBV infection. The gene expression data were calculated using RNA-Seq data with or without EBV infection from SRA.

The network between the top five-deceased miRNAs through interaction with circrnaome of EBV and small compounds

The expression of some miRNAs was decreased after EBV infection including miRNAs sponged by EBV circRNAome. According to the miRNet database, these miRNAs including has-miR-30c-1-3p, has-miR15b-5p and has-miR-424-5p were associated with 27 small compounds forming a miRNA-compound network (Figure 8). Interestingly, some of the compounds such as etoposide, 5-fluororacil in the network were extensively used to treat cancers including EBV-associated non-Hodgkin’s lymphoma and nasopharyngeal carcinoma [32–35]. These data suggested that these chemotherapeutical compounds might exert their treatment by interfering the miRNAs sponged by EBV circRNAome.

Figure 8. The network between the top five-deceased miRNAs which interacted with circRNAome of EBV and small compounds. The network was generated using miRNet database. Some compounds as medicines were usually used to treat cancers in clinic including tumors associated with EBV infection.

Discussion

Recently, R-resistance RNA sequencing identified the circRNAs encoded by EBV genome, which composed the EBV circRNAome [11,12]. The functions of the circRNAs are not well studied yet. In the current research, we tried to investigate the function of cicrRNAs as a sponge of host miRNAs using miRNA database, miRNA-seq data, and bulk RNA-seq data from EBV-infected cells.

Circrnaome of EBV sponged mirnas of host cells involved in malignancies

CircRNAs could bind to miRNA so that circRNAs competitively inhibit miRNAs to bind to their targets which prevented the function of miRNAs. Through miRNA binding prediction, miRNA- and RNA-seq data analysis, miR-15b-5p, miR-30c-1-3p, miR-30c-2-3p, miR-424-5p, and miR-4286 were identified as top possible candidates sponged by EBV circRNAome. The overexpression of miR-15b promoted apoptosis and arrested cell cycle of glioma by targeting cyclin D1 [36] while the knockout of miR-15b promoted B-cell malignancies by modulating the expression of cyclin D1 and D2 [37]. EBV was associated with B lymphoma such as Burkitt, Hodgkin and Diffuse Large B cell Lymphomas (DLBCL) [3]. The expression of miR-30c was positively correlated to survival of breast cancer patients [38,39]. Loss of miR-424 was associated with aggressiveness of breast cancers in human. Knockout of miR-424 promoted breast tumorigenesis in mice [40]. miR-4286 was also highly expressed in chemo-sensitive tumors [16]. Since these miRNAs had such important functions in the multiple aspects of tumor biology, the inhibition of these miRNAs by EBV circRNAs would promote tumorigenesis, cancer aggressiveness and chemotherapy resistance.

The targets of miRNAs sponged by EBV circrnaome regulated herpes virus infection

CircRNAs sponged miRNAs to increase the expression of miRNAs’ target genes. Most of the targets associated with herpes virus infection were from miR-15 such as FADD, TAB2, SRSF1, EIF2AK2 and SRPK1. Their expressions elevated more than twofold in EBV-infected cells. Kaposi’s-sarcoma-associated human herpesvirus-8 encoded v-FLIP could interact with FADD through death-effector domains and inhibit apoptosis signaling by death receptors [41]. The protein complex TAB2/TAK1/TRAF6 was recruited by N-terminal of Epstein-Barr virus latent membrane protein 1 (LMP1) to promote LMP1-mediated activation of JNK [42]. Herpes simplex virus (HSV) infected cell culture polypeptide 27 (ICP27) interacted with SRSF1 or SRPK1 to regulate host pre-mRNA aberrant splicing, which contributed to virus-induced host shutoff required for efficient viral growth [43,44]. HSV-1 protein US11 physically interacted with EIF2AK2 to negatively regulate the autophagy-constrained virus infection [45]. These evidences supported that the up-regulated targets of miRNAs sponged by circRNAs could regulate a variety of functions associated with virus infection including viral growth, avoidance of limitation by host cells, host cell activation and survival.

The targets of miRNAs sponged by EBV circrnaome regulated cell cycle

Except for miR-30c-1, the targets of the other four miRNAs (Figure 3b) were involved in cell cycle regulation. These targets fall into the following categories: 1) cyclin and cyclin-dependent kinase such as CCNE1/2, CDK4/6, 2) cell division cycle phosphatase such as CDC25A, CDC27, 3) minichromosome maintenance protein complex such as MCM7, 4)YWHA family such as YWHAB/H/Q, 5) the others such E2F3, GSK3B, SMC1A. Cyckin E controlled G1-S phase of cell cycle in complex with CDK2 [46]. CDK4/6 controlled G1 phase of cell cycle in complex with cyclin D [46]. Both CDC25A and CDC27 promoted G1/S-phase transition [47,48]. MCMs primed chromatin for DNA replication through binding origins of DNA replication during the late M to early G1 phase transition of the cell cycle [49]. 14-3-3 proteins (YWHA family) functioned at several key points in G1-S and G2-M transition by binding to regulatory proteins and modulating their functions [50]. E2F3 controlled the expression of genes which rately limited the initiation of DNA replication [51]. GSK3B regulated the expression of cyclin D, which was involved in controlling G1 phase cell cycle in CDK4/6 complex [52]. Knockdown of SMC1A led to arresting at G0/G1 phase [53]. These studies suggested that the up-regulated targets of miRNAs which were sponged by circRNAs promoted the cell cycle progress. In contrast, the negative regulators of cell cycle such as CDKN1A and TGFB1were up-regulated. This might be the negative feedback to balance the over-activation of cell cycle.

The targets of mirnas sponged by EBV circrnaome regulated cancer pathway

Except for miR-30c-2, the targets of other four miRNAs were involved in oncogenesis. Since cell cycle are tightly associated with cancers, a lot of targets in cancer pathway overlapped with the targets in cell cycle including cyclin and cyclin-dependent kinases. In addition to the targets associated with cell cycle, the other targets including ITGA3, PTCH1, CUL2, CRK, FOXO1, were involved in other aspects of cancers. Knockout of ITGA3 prevented skin tumor formation by promoting epidermal turnover and depletion of slow-cycling cells [54]. Overexpression PTCH1 drove skin carcinogenesis [55]. Overexpression of CUL2 drove cervical carcinogenesis which might be associated with miR424 [55]. Overexpression of CRK induced the prominent tyrosine phosphorylations of scaffolding molecules such as p130^Cas and paxillin, ultimately contributing to the increased motility and aggressive behavior of cancer cells [56]. FOXO1 could promote tumor metastasis by up-regulating phosphorylation and MMP expression [57]. Thus, in addition to targets that regulate cell cycle, lots of targets also regulated other aspects of cancers ranging from carcinogenesis to migration, invasion and metastasis.

Applications and perspective

Currently, tons of data from next generation sequencing (NGS) were deposited or collected by various databases such SRA to be used for free by the research community. It is a huge challenge on how to utilize these data to dig out useful information so as to aid the biomedical research. In the current research, integrative data from R-resistance RNA-seq, miRNA-seq, and RNA-seq using EBV-associated biomaterials found the miRNAs sponged by cicrRNAome of EBV as well as the target genes that regulated these miRNAs. This research provided an example of how to use the NGS data stored in the databases to discover the real biomedical project. The network among EBV circRNAome, host miRNAs, and miRNAs’ targets was beneficial for the investigation of malignancies associated with EBV.

Conclusion

R-resistance RNA-seq defined the circRNAome of EBV. The miRNA-seq and RNA-seq data supported that circRNAome of EBV as miRNA sponge regulated miRNA functions to increase gene expression which was involved in carcinogenesis and cancer progression.

Highlights

• EBV cicrRNAome was involved in EBV-associated malignancies.

• EBV circRNAome was a sponge of host cells’ miRNA.

• The targets of EBV circRNAome sponged miRNAs constituted to the comprehensive network.

• The targets of EBV circRNAome were involved in multiple categories including virus infection, cell cycle, and oncogenesis.

• CCNE1/2, CDC25A, E2F3, and TAB2 genes were the top-ranked targets of miRNAs which were sponged by EBV circRNAome.

Availability of data and material

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Author’s contributions

WY contributed to the study design; all authors collected the data and performed the data analysis; all authors prepared the manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.