A Genome-Wide Epstein-Barr Virus Polyadenylation Map and Its Antisense RNA to EBNA

Active EBV expression in JSC-1, Raji, and Akata cells revealed by PA-seq.

To map the genome-wide pA sites and their usage of EBV transcripts, three EBV-positive cell lines, EBV- and KSHV-coinfected JSC-1 (47), EBV nonproducer Raji (48), and EBV producer Akata (49), from latent and lytic infections, were used for the study by PA-seq analysis. The three-EBV-genome alignment in Fig. S1 in the supplemental material shows that the Raji EBV genome has two large deletions, first from nucleotides (nt) 87069 to 90217 (3,148 bp) and then from nt 163986 to 166643 (2,657 bp) (50), but fewer repetitive sequences from nt 170351 to 172550 (22). The Akata EBV genome has fewer repetitive sequences from nt 96351 to 97100 and has two small deletions, first from nt 88144 to 88248 and then from nt 90706 to 90733 (51). The EBV reference genome (GenBank accession no. NC_007605) (19, 22, 23) also has two small deletions, first from nt 88617 to 88772 and then from nt 90706 to 90733, compared with the Raji and Akata genomes (Fig. S1).

EBV lytic reactivation in JSC-1 cells was induced by n-butyrate, that in Raji cells was induced by a combination of both butyrate and 12-O-tetradecanoylphorbol-13-acetate (TPA) (52), and that in Akata cells was induced by surface receptor cross-linking by the addition of goat anti-human serum to the culture medium (49, 53). The most effective method for lytic EBV induction in Akata cells, under multiple tested conditions, was to treat the cells for 24 h with a 1:200 dilution of whole serum, as measured by the expression of the lytic BMRF1/2 gene. The total RNAs from both uninduced (latent) and induced (lytic) cells were harvested for construction of PA-seq libraries (46). The PA-seq libraries from JSC-1 cells were obtained from our previous study and were successfully used to determine KSHV genome-wide pA sites with high levels of accuracy and specificity (25).

To identify EBV pA sites, we first extracted EBV-specific reads by alignment of the PA-seq reads derived from these libraries to the EBV reference genome (GenBank accession no. NC_007605) (23). We identified 2,001,977 sequence reads in JSC-1 cells, 1,701,835 sequence reads in Raji cells, and 4,925,618 sequence reads in Akata cells uniquely mapped to the EBV genome (Table S1). Among these EBV reads, 3,233 reads in JSC-1, 27,875 reads in Raji, and 29,557 reads in Akata cells were identified from uninduced cells with latent EBV infection, and 1,998,744 reads in JSC-1, 1,673,960 reads in Raji, and 4,896,061 reads in Akata cells were from induced cells with lytic EBV infection (NCBI accession no. SRP157138), representing 0.02% and 1.07% to 13.11% of the total PA-seq reads, respectively (Fig. 1A and Table S1). We noticed that the numbers of EBV reads in JSC-1 cells were smaller than the numbers of KSHV reads obtained from the same libraries. Although the total numbers of EBV-specific reads were lower than those mapped to KSHV, the ratio (>500) of lytic EBV reads from the butyrate-treated cells over the latent EBV reads from the untreated cells was much higher than that (>200) of KSHV lytic reads over the latent reads (Fig. 1B). We also noticed that the level of active expression of EBV lytic genes in JSC-1 cells was much higher than that in Raji and Akata cells upon lytic induction, despite all three cell lines exhibiting similar levels (∼0.02%) of latent gene expression.

Genome-wide distribution of EBV polyadenylation sites.

To determine the nucleotide position and usage frequency of individual EBV pA sites, the mapped EBV reads from the cells with latent or lytic infection were first visualized against the EBV reference genome (GenBank accession no. NC_007605) for both viral DNA strands by using IGV software. The EBV reference genome is a hybrid genome sequence composed of the reference EBV B95-8 strain and the Raji strain. As expected, we observed only a few viral reads from each cell line with latent EBV infection but a large number of viral reads from the cell lines with lytic induction, as represented by JSC-1 cells (Fig. 2A). Although the viral reads were scattered along the viral genome, a cluster of the EBV reads can be visualized as a distinct peak on both EBV DNA strands immediately downstream of the majority of annotated open reading frames (ORFs) (Fig. 2). In uninduced JSC-1 cells, a majority of the reads were derived from transcripts of BART, a known EBV latency-associated transcript (54, 55). Analysis of the selected peaks showed that the obtained EBV reads specifically distribute only downstream of polyadenylated Pol II transcripts like LMP-2A/2B. None of the reads were mapped downstream of the nonpolyadenylated Pol III-derived EBER1 and EBER2 RNAs (Fig. 2B, left), two abundant small noncoding EBV RNAs of 167 and 172 nt, respectively (56, 57). Importantly, the clustered reads’ peak could be annotated on the EBV genome in a strand-specific manner and usually downstream of an annotated ORF, such as BBRF3 and BBLF1 (Fig. 2B, middle). In general, the peak heights could be correlated with the expression levels of individual viral genes, as exemplified by BBRF3 by comparing its reads in latent versus lytic infections. Most of the clustered read peaks were mapped to a previously identified pA site (29–34), like the peak downstream of BMLF1-BSLF2 at the position of a previously mapped pA site at nt 70442 (minus strand [−]) (32) (Fig. 2B, right). Thus, this analysis showed that the clustered peaks of PA-seq reads represent the pA sites of EBV RNA transcripts.

Subsequently, we performed a peak-calling analysis using F-Seq (25, 46). Obtained peaks represent the clusters of PA-seq reads from a specific region where heterogeneity of cleavage sites exists, as normally seen in KSHV and mammalian RNA polyadenylation (25, 58). Each identified peak was then redefined by two additional parameters: peak mode and usage. The peak mode corresponds to the nucleotide position with the highest number of reads within a clustered individual peak. This represents the most prevalent cleavage site used for polyadenylation and, therefore, is further referred to as a “pA site.” The pA site usage is a sum of all reads within each clustered peak and correlates with the relative expression levels of the RNA transcripts. After peak-calling analysis, we eliminated all peaks resulting from internal priming to genomic A-rich regions and those peaks having fewer than 50 reads in all samples. Consequently, the total numbers of PA-seq reads mappable to the identified pA sites were 1,387,309 reads (including 1,725 latent EBV reads) for the JSC-1 EBV genome, 1,621,448 reads (including 18,940 latent EBV reads) for the Raji EBV genome, and 4,797,216 (including 28,590 latent EBV reads) for the Akata EBV genome (see Table S1 in the supplemental material). Using these criteria, we identified a total of 62 pA sites from the JSC-1 EBV genome, 60 pA sites from the Raji EBV genome, and 53 pA sites from the Akata EBV genome (Table S1). We used IGV software to visualize the distribution of individual pA sites across the EBV reference genome (GenBank accession no. NC_007605), mapping 27 pA sites to the plus strand (+) and 35 sites to the minus strand of the JSC-1 EBV genome, 26 pA sites to the plus strand and 34 pA sites to the minus strand of the Raji EBV genome, and 26 pA sites to the plus strand and 27 pA sites to the minus strand of the Akata EBV genome (Fig. 3A). Among them, 42 pA sites were consistently found in all three EBV genomes (Fig. 3B). We also determined genome-specific pA site positions after mapping PA-seq libraries to the complete Raji (GenBank accession no. KF717093) and Akata (GenBank accession no. KC207813) EBV genomes (Table S2) and compared them to the reference EBV genome (GenBank accession no. NC_007605). Although the mapped genome positions of the individual pA sites vary from the reference EBV genome to the Raji and Akata genomes, and to the few previously reported pA sites based on the B95-8 genome (29–34), these data indicate that the mapped pA sites represent a comprehensive atlas of EBV pA sites used for the expression of their corresponding genes in these transformed B cells.

Comparing the mapped pA sites in the JSC-1 EBV genome (Fig. 3B), the Raji EBV genome utilizes the 48 pA sites as seen in the JSC-1 EBV genome but misses 6 pA sites from the plus strand at nt 10281, 48435, 79250, 89500, 99451, and 165728 and 8 pA sites from the minus strand at nt 10130, 35713, 66588, 87158, 89217, 89495, 94664, and 164368, whereas the Akata EBV genome shares 49 pA sites as seen in the JSC-1 genome but misses 3 pA sites from the plus strand at nt 37739, 65684, and 99451 and 10 pA sites from the minus strand at nt 35713, 39843, 39961, 66588, 82967, 86049, 87185, 89217, 89495, and 94664 (Table S1). A cluster of the eight mapped pA sites from nt 39691 (−) to nt 40045 (−) in the BHLF1 coding region seen in Raji, but not Akata, cells suggests that they are not essential for lytic EBV induction. A similar implication could also be applied to the pA sites at nt 7126 (+) and nt 7181 (+) seen only in Raji but in neither JSC-1 nor Akata cells. The heterogeneity in the compositions of EBV pA sites among the three cell lines may be due to viral genome polymorphism, variable numbers of internal repeats, and cell line-specific deletions.

Usage of EBV pA sites from viral latent to lytic infection in JSC-1, Raji, and Akata cells.

Based on the number of reads mapped to each peak, referred to as pA site usage, we next assessed the relative expression levels of EBV transcripts polyadenylated at each identified pA site. The usage for each pA site was first normalized to library coverage by dividing the number of raw reads by the total number of reads of each library, in millions, to obtain RPM (reads per million). Normalized RPM were then used to compare pA site usages within or between individual samples from the three cell lines. As shown in Fig. 4A and Table S1 in the supplemental material, all three cell lines exhibited a distinct expression pattern of latent and lytic genes, with a limited number of reads in latent infection compared to a large number of reads in lytic induction. Although Raji cells (type III latency) (59) exhibited more expression of the EBV latent genes than JSC-1 (type II latency) (47) or Akata (type I latency) (35, 49) cells, the pA site usage at nt 160531 (+) was found in all three cell lines for the expression of latent BART transcripts (also called RPMS1/A73/BARF0 transcripts). BART transcripts encode a cluster of EBV-derived miRNAs essential for the establishment and maintenance of latency (60). The top five most-used pA sites quantified as a percentage of the total EBV reads from cells with latent infection are shown in Fig. 4B (left): for JSC-1 cells, 41% at nt 160531 (+) for BART expression, 23% at nt 38009 (+) for BHLF1, 5% at nt 69910 (+) for BMRF2, 3% at nt 77117 (−) for BLLF1/2, and 2% at nt 99518 (+) for BKRF4; in Raji cells, 41% at nt 43251 (+) for BHRF1, 11% at nt 37739 (+) for EBNA-2/EBNA-LP, 9% at nt 97669 (+) for EBNA-1, 7% at nt 7181 (+) for an unannotated gene, and 6% at nt 140693 (−) for LF3; and in Akata cells, which displayed the most bias usage of the three pA sites, 41% at nt 140693 (−) for LF3, 33% at nt 38009 (−) for BHLF1, and 16% at nt 160531 (+) for BART.

Expression during lytic viral infection was more consistent among all three tested cell lines, although JSC-1 cells displayed usage of the most pA sites for expression of the corresponding lytic genes compared to Raji or Akata cells (Fig. 4A and Table S1). Three of the top five pA sites (Fig. 4B, right) were located at nt 38009 (−) for BHLF1, nt 140693 (−) for LF3 and nt 69910 (+) for BMRF2 transcripts, which accounted for 5% to 48% of the total pA-seq reads. When combined, the read counts from the three pA sites were 58% for JSC-1, 73% for Raji, and 95% for Akata cells (Fig. 4B and Table S1). Surprisingly, in Akata cells, the reads from nt 38009 (−) for BHLF1 and from nt 140693 (−) for LF3 over the total pA site reads accounted for 48% and 47%, respectively. In contrast, variations in the usage of the nt 69910 (+) pA site for BMRF2 expression were noticed from Akata (1%) and JSC-1 (5%) to Raji (22%) cells (Fig. 4B). These pA sites, also expressed during latent infection, displayed a several-100-fold increase in the number of PA-seq reads between latent and lytic infections. Thus, the data suggest that they are truly inducible genes during lytic EBV infection. We also observed that lytic induction of Raji cells, but not JSC-1 or Akata cells, blocked the pA site usage at nt 97669 (+) for the expression of latent EBNA-1 and at nt 37739 (+) for EBNA-2 (Fig. 4A and Table S1).

We next compared the fold changes (FCs) in viral gene expression associated with the switch from latent to lytic infection by calculation of the normalized pA site usage from EBV lytic infection over that from latent infection as a ratio of RPM_Lytic/RPM_Latent. As shown in Table S1 in the supplemental material, expression of the pA site at nt 137448 (−) for BILF2 (viral glycoprotein) had the highest increase of 5,183-fold, followed by the second highest increase of 4,698-fold for the pA site at nt 140693 (−) for LF3 transcripts after EBV reactivation in JSC-1 cells. However, the highest FCs from the latent-to-lytic pA site switch were for the pA sites at nt 89839 (−), with a 2,503-fold increase, for BRLF1 (Rta)/BZLF1 (Zta) and at nt 128592 (−), with a 1,570-fold increase, for BXLF1 (thymidine kinase)/BXLF2 (gH protein) in Raji cells and for the pA sites at nt 38009 (−), with a 262.3-fold increase, for BHLF1 (noncoding RNA) and at nt 140693 (−), with a 210-fold increase, for LF3 (noncoding RNA) in Akata cells. The lowest expression level changes for two pA sites clearly distinguishable from the rest in all three cell lines were the pA sites at nt 97669 (+) for BKRF1/EBNA-1 and at nt 160531 (+) for BART transcripts with FCs of only 2.9- and 3.4-fold in JSC-1, 3.1- and 14.6-fold in Raji, and 1.7- and 1.6-fold in Akata cells after virus reactivation. Therefore, EBNA-1 and BART transcripts are two true latent transcripts in all three cell lines. The pA site usages at nt 166466 (−) and at nt 5856 (+) for the expression of two classical EBV latent genes, LMP-1 and LMP-2A/B (61), were increased by 1,331- and 455-fold in JSC-1 but 11- and 7-fold in Raji and 68- and 33-fold in Akata cells, respectively, during virus lytic induction. These differences in the increased expression of LMP-1 and LMP-2A/B in lytic induction could be due to the induced expression of lytic BNLF2 and BNRF1 for their usage of the same pA site. Similarly, the pA site at nt 43251 (+) for the expression of BHRF1 (bcl-2 homolog)/miBH1-1-3 from latent to lytic infection was increased only 3.7-fold in Raji but 930-fold in JSC-1 and 147-fold in Akata cells. The pA site at nt 82957 (+) for EBNA-3A was relatively inducible in all three cell lines, but the pA sites at nt 89500 (+) for EBNA-3B and -3C expression were also relatively inducible only in JSC-1 and Akata but not in Raji cells with a deletion in this EBV region (62).

To verify if pA site usage truly mirrors EBV lytic expression, we performed Northern blot analysis using total RNA from all three cell lines prepared under the same conditions as those used for PA-seq. The transcripts polyadenylated at the most prevalent pA sites at nt 38009 (−), nt 140693 (−), and nt 69910 (+) were detected by antisense oligonucleotide probes located in BHLF1, LF3, and BMRF1/2 loci. As shown in Fig. 4C, the highly induced expressions of BHLF1, LF3, and BMRF1/2 were detected as single species of RNA transcripts of the expected size from all three cell lines with lytic EBV induction, although more host shutoffs by lytic EBV induction or toxicity of chemical inducers were observed for JSC-1 and Raji cells, as indicated by glyceraldehyde-3-phosphate dehydrogenase (GAPDH) RNA levels. Based on the expected size, the detected BMRF transcripts were viewed as bicistronic BMRF1/2 transcripts because the mapped pA site at nt 69489 (+) for the expression of BMRF1 in Raji cells was not used in JSC-1 cells or little used in Akata cells (Table S1). Moreover, BHLF1 RNA in Akata cells showed a larger band size than that from JSC-1 or Raji cells (Fig. 4C) because of the presence of an additional 125-bp GC-rich repeat in this region in the Akata EBV genome (Fig. S1).

RNA cis elements in regulation of EBV expression and polyadenylation.

To determine the cis elements in regulation of EBV polyadenylation, we analyzed the sequence flanking each identified pA site within a range of −50 to +50 nt from each pA site. First, we performed an alignment of all sequences using the sequence logo algorithm to identify any consensus motifs surrounding the mapped pA sites. As shown in Fig. 5A, an A/T-rich sequence motif was found to be highly conserved around −20 to −15 nt upstream of the mapped pA site, and a T/G-rich sequence motif was enriched +10 to +20 nt downstream of the mapped pA site. These motifs coincide with the PAS and DSE commonly found in mammalian transcripts (63, 64). Since the PAS plays a critical role in pA site selection, by scanning 50 nt upstream of each mapped pA site for a PAS signal sequence, we identified the presence of a PAS in 59 of the 62 mapped pA sites in the JSC-1 EBV genome, 48 of the 60 pA sites in the Raji EBV genome, and 52 of the 53 pA sites in the Akata EBV genome (see Table S1 in the supplemental material). Among the mapped pA sites, 49 (79%) are the canonical “AAUAAA” PAS and 10 (16%) are noncanonical PAS variants (AAUACA, AUUAAA, GGUAAA, GAUAAA, AGUACA, CAUAAA, CAUAUA, and UAUAAA) in JSC-1 cells, 41 (68%) are canonical and 7 (11%) are noncanonical in Raji cells, and 44 (83%) are canonical and 8 (15%) are noncanonical in Akata cells. We were unable to identify any PAS or PAS variants in 3 pA sites in JSC-1 cells, 13 (including those in the BHLF1 region) in Raji cells, and only 1 in Akata cells (Fig. 5B and Table S1). The overall PAS prevalence in association with the mapped pA sites is close to those observed in KSHV and eukaryotic pA sites (25, 64, 65). Among the top 10 most-used pA sites identified in lytic EBV infection in all three cell lines, 90% bear a canonical PAS upstream. However, in the least-used bottom 10 pA sites, only 50 to 60% in JSC-1 and Akata cells, and 16% in Raji cells, use a canonical PAS (Fig. 5B). We observed a positive correlation between the pA site expression level and PAS sequence conservation in expression of the EBV genome (Table S1). A similar correlation was observed for KSHV expression (25). Together, these data indicate that EBV polyadenylation sites are defined by common regulatory elements found in host transcripts and that the sequence composition surrounding the mapped pA site directly correlates with the level of viral gene expression.

Verification of PA-seq-identified pA sites.

To confirm mapping accuracy, we compared the identified pA sites with several previously reported pA sites, including pA sites for LMP-2A/2B, BZLF1, EB2 (BMFL1-BSLF2), and LF3 transcripts (29–32). In all cases, we observed a high correlation between the PA-seq-identified pA site and the reported pA site. In addition, we experimentally confirmed the pA sites at nt 5856 (+) for LMP-2A/2B expression and at nt 160531 (+) for BART (RPMS1/A73/BARF0) expression by 3′ RACE on total RNA from three EBV-positive cell lines, JSC-1, Raji, and Akata. In all cases, two pA sites determined by 3′ RACE were matched to those mapped by PA-seq and were identical in all three cell lines with lytic EBV infection (Fig. 6), indicating that the EBV pA site usage in JSC-1 cells also coinfected with KSHV is not different from that of Raji or Akata cells infected with only EBV. Altogether, we have confirmed a high accuracy of pA sites identified by PA-seq.

Taking into account the nucleotide position and strand specificity of the pA sites mapped by PA-seq, we next analyzed the position of each identified pA site in the context of annotated EBV genes. As expected, the majority of pA sites were mapped to the regions immediately downstream of the assigned ORFs on both plus (sense) and minus (antisense) DNA strands. Interestingly, many were mapped downstream of a gene clustering locus containing multiple colinear genes. Their unique usage for RNA polyadenylation in this gene cluster leads to the production of multiple polycistronic RNA transcripts commonly observed in all herpesviruses (25) and papillomaviruses (66, 67). To confirm the production of multiple polycistronic RNA transcripts, we chose JSC-1 cells, which exhibit the highest usage of the most pA sites among all three cell lines, and performed Northern blot analysis of five gene clusters on total RNA extracted from JSC-1 cells with or without lytic virus induction using oligonucleotide probes located upstream of the mapped pA sites (Fig. 7). These five pA sites were chosen because of their high expression levels with thousands of mapped reads in the corresponding gene transcripts during lytic EBV infection. First, we investigated the locus at the head-tail junction region of the episomal EBV genome encoding the latent LMP-1 oncoprotein, BNLF2a, and BNLF2b from its minus DNA strand. While BNLF2a functions as an expressible viral immunomodulator (68), the expression and function of BNLF2b remain hypothetical. Northern blotting using a sense oligonucleotide probe upstream of the mapped pA site at nt 166466 (−) revealed the expression of two major transcripts of ∼0.8 and 2.8 kb from this region in JSC-1 cells with only lytic EBV induction and not with latent infection (Fig. 7A). The data are consistent with the expression level determined by PA-seq reads (26,785 lytic reads versus 18 latent reads) in this pA site (see Table S1 in the supplemental material). Based on RNA size, these transcripts represent LMP-1, with a predicted size of ∼2.4 kb (from the ATG codon to the identified pA site), and BNFL2a/b, with a size of ∼0.6 to 0.8 kb. The second locus selected consisted of latent LMP-2A/2B and colinear BNRF1 (major tegument protein) from its plus DNA strand. Using an antisense oligonucleotide probe located in the terminal exon of the LMP-2A/2B locus, we detected one major transcript, BNRF1, with a size of ∼4.5 kb only in JSC-1 cells with lytic EBV induction (Fig. 7B) for expression of the EBV major tegument protein (69, 70). The expression of LMP-2A/2B transcripts was hardly detectable from this gene cluster region (Fig. 7B). These data indicate that the PA-seq reads detected in JSC-1 cells were mainly derived from BNRF1 transcripts expressed from this latent locus.

Next, we determined the expression of polycistronic transcripts from two lytic loci. One contained BFRF1A (viral DNA packaging protein), BFRF1 (virion egress protein), BFRF2 (uncharacterized protein), and BFRF3 (small capsid protein) and had two pA sites, at nt 48435 (+) and at nt 49799 (+), which were identified by PA-seq. The other locus used encodes BORF1 (minor capsid protein), BORF2 (viral ribonucleoside-diphosphate reductase, large subunit), BaRF1 (viral ribonucleoside-diphosphate reductase, small subunit), BMRF1 (viral DNA polymerase processivity factor), and BMRF2 (glycoprotein) and has three pA sites, at nt 65684 (+), nt 66617 (+), and nt 69910 (+) in JSC-1 cells, which were identified by PA-seq (Table S1). Northern blotting with an oligonucleotide probe upstream of the nt 49799 pA site detected two major transcripts only in the lytic JSC-1 sample, an ∼1-kb BFRF3 transcript (predicted size, ∼0.6 kb) (71) and an ∼3.7-kb BFRF1-3 transcript (72), and one minor ∼3-kb BFRF2/3 transcript (73) (Fig. 7C). Although a mapped pA site with few (only 417) reads at nt 48435 (+) was found downstream of BFRF1, the lack of a canonical poly(A) signal for this pA site usage (Table S1) would make most BFRF1 transcripts skipping the nt 48435 (+) pA site to polyadenylate at the nt 49799 (+) pA site, resulting in the production of a bicistronic BFRF1-3 RNA with a size of ∼3.7 kb, as detected by others with a BFRF1 probe (72) (Fig. 7C). Like the nt 48435 (+) pA site, the mapped nt 65684 (+) pA site also lacks a canonical poly(A) signal and was mapped with few (only 155) reads in the locus (Table S1). Our detection was thus focused only on two highly used pA sites at nt 66617 (+) and nt 69910 (+) in this locus. Northern blotting with a BORF2 probe (oVM376) upstream of the nt 66617 (+) pA site revealed the expression of two major transcripts corresponding to monocistronic BORF2 and bicistronic BORF1/BORF2, with predicted sizes of ∼2.5 and ∼3.7 kb, respectively (Fig. 7D). Two additional bands representing large RNA species (>10 kb in size) were also detected, which were presumably polycistronic BURF1 or BURF1-like RNA transcripts (36) (Fig. 3). The probe oVM377 upstream of the nt 69910 (+) pA site identified 3 major lytic transcripts: the less abundant monocistronic BMRF2 transcript, the most abundant bicistronic BMRF1/BMRF2 transcript, and the less abundant polycistronic BaRF1/BMRF1/BMRF2 transcript, for which the observed sizes were in close agreement with the corresponding predicted sizes of ∼1.1, ∼2.3, and ∼3.3 kb (Fig. 7D).

In conclusion, we have generated a comprehensive atlas of EBV polyadenylation sites and their contribution to the EBV transcriptome in the context of current gene annotation. This study also carefully elucidated that EBV produces multiple bicistronic and polycistronic transcripts from several gene cluster regions by sharing a single pA site for their RNA polyadenylation.

Identification of novel pA sites for polyadenylation of EBV antisense transcripts.

Although most of the mapped pA sites could be assigned downstream of the annotated EBV ORFs, we also identified a subset of pA sites either positioned within the annotated ORFs or not associated with any ORF (Fig. 3; see also Table S1 in the supplemental material). For example, the mapped pA sites in JSC-1 cells assigned to the known ORF regions include the pA site at nt 65684 (+) within the coding region of BORF2 and the pA sites at nt 39843 (−) and nt 39961 (−) in BHLF1, nt 99451 (+) in BKRF4, and nt 140693 (−) in LF3 loci (Fig. 3). In general, these pA sites are characterized by few expression reads, lack a canonical PAS and/or DSE, and are not consistently detected in all three cell lines. Therefore, they represent a group of cryptic pA sites in the coding region that could be activated by noncanonical poly(A) polymerases (74), presumably for regulation of the full-length RNA expression and coding functions. An exception is the pA site at nt 140693 (−), which was mapped to 123 nt upstream of the LF3 stop codon at nt 140570 (−) (Fig. 3). In contrast to other pA sites mentioned above, this pA site was defined by a large number of the mapped reads (Table S1) and has a canonical PAS. In contrast, the pA site at nt 139398 (−) located ∼50 bp downstream of the predicted LF3 protein stop codon (75) was not detectable in either JSC-1 or Akata cells and expressed only 143 reads in lytic Raji cells by using a noncanonical PAS. Therefore, the mapped pA site at nt 140693 (−) is authentic for the production of a noncoding LF3 transcript (30, 38, 51).

The mapped pA sites not associated with any annotated ORF in the EBV genome are pA sites at nt 6793, 7126, 7181, and 79250 from the plus DNA strand and at nt 1764, 10130, 35713, 66588, 82967, 86049, 87185, and 94664 from the minus DNA strand (Fig. 3 and Table S1). A cluster of pA sites at nt 6793 (+), nt 7126 (+), and nt 7181 (+) in proximity of the EBER locus was detected mainly in Raji cells with low-level usage because they lack a PAS or PAS-like signal upstream. The pA site at nt 79250 (+) without a PAS was detected with 100 reads in JSC-1 and with 44 reads in Akata cells with lytic induction (Table S1) and is located 2.1 kb downstream of the highly expressed BLRF2 pA site at nt 77143 (+). Thus, these pA sites are likely all cryptic pA sites. From the minus EBV DNA strand, the pA site at nt 1764 was mapped only in Akata cells. Although not highly used, it contains a canonical PAS for an antisense transcript to LMP-2A/2B and BNRF1. Interestingly, the remaining other six pA sites from the minus DNA strand are positioned in the antisense direction from individual EBNA transcripts. More precisely, the pA site at nt 10130 (−) is upstream of the EBNA Cp promoter; the pA site at nt 35713 (−) is upstream of the EBNA-2 and EBNA-LP terminal exons; the pA site at nt 66588 (−) is in the intron region of EBNA-3 and EBNA-1 and also upstream of the BMRFs and BaRF1 ORFs; the pA sites at nt 82967 (−), nt 86049 (−), and nt 87185 (−) are in the terminal exons of EBNA-3A, -3B, and -3C; and the pA site at nt 94664 (−) is upstream of the EBNA-1 terminal exon and BKRFs (Fig. 3). Their usage was highly varied among individual pA sites, from a low-use pA site at nt 94664 (−) to a high-use pA site at nt 82967 (−) (Table S1). In addition, all reads mapped to these pA sites came almost exclusively from cells with virus lytic induction. Analysis of the sequences surrounding these pA sites showed that a majority contain a canonical “AAUAAA” PAS, except for the pA sites at nt 86049 (−) with no detectable PAS and the pA site at nt 87185 (−) with a noncanonical “GATAAA” PAS. Their expressions were highly variable among three cell lines. While all pA sites were mapped from JSC-1 cells, only two (pA sites at nt 82967 and 86049) were detected from Raji cells, and one (pA site at nt 10130) was detected from Akata cells (Table S1). These data indicate the presence of several potential polyadenylated transcripts antisense to various forms of EBNA transcripts, as suggested by several laboratories (36, 76, 77).

Given the importance of EBNA proteins in EBV biology and cell transformation, we decided to further verify the mapped pA sites for the expression of EBNA-antisense transcripts by 3′ RACE analysis. The pA site at nt 94664 (−) was excluded from the analysis due to a very low number of associated reads detectable only in JSC-1 cells. Total RNAs from all three cell lines with lytic virus infection were analyzed by 3′ RACE with an EBV-specific primer (oVM284, oVM285, oVM286, oVM287, or oVM288) upstream of the tested pA sites (Fig. 8A). We confirmed the expression of all mapped pA sites in JSC-1 cells but partially in Raji and Akata cells (Fig. 8B). In particular, the pA site at nt 10130 (−) upstream of the EBNA Cp promoter (78) was not detectable or was expressed at a very low level in Raji cells with type III latency (59). The amount of 3′ RACE products was comparable with the calculated pA site usage, with the pA sites at nt 82967 (−) (oVM286 product) being the most prominent in JSC-1 and Raji cells (Fig. 8B). The identity of the final RACE products and pA site positions in JSC-1 cells were further confirmed by Sanger sequencing (Fig. 8C). Raji cells exhibited the main usage, with similar levels, as mapped by PA-seq, of two adjacent pA sites at nt 82967 (−) (oVM286 product) and nt 86049 (−) (oVM287 product) (Fig. 8B). Akata cells showed a usage of two pA sites at nt 10130 (−) (oVM284 product) and nt 35713 (−) (oVM285 product), with a weak 3′ RACE product from the pA site at nt 35713 (−) (oVM285 product). The pA site at nt 35713 (−) was not mapped by PA-seq, probably due to its low expression level. Together, the existence of EBNA-antisense pA sites and their usage during EBV lytic infection in a cell type-dependent manner were verified by using 3′ RACE.

Detection of polyadenylated transcripts antisense to EBNA in JSC-1 and Raji cells by Northern blot analyses.

To further understand the nature of EBNA-3-antisense transcripts, we performed Northern blot analysis on poly(A)-selected RNA from three cell lines, with poly(A)⁺ RNA from the BCBL-1 cell line without EBV infection as a negative control. To determine transcript expression kinetics, we compared their expressions in these cells with latent (no chemical induction) or lytic (with chemical induction) infection. Three specific-sense oligonucleotide probes corresponding to the same locations as those of sense oligonucleotides used for 3′ RACE (Fig. 9A; see also Table S3 in the supplemental material) were used for Northern blot analyses for detection of the antisense transcripts to EBNA-3B and -3C, where the associated pA sites at nt 82967 (−), nt 86049 (−), and nt 87185 (−) expressed a large amount of mapped reads (Table S1). As shown in Fig. 9B, the oVM293 probe located upstream of the pA site at nt 82967 (−) provided a positive signal, as expected, only in JSC-1 and Raji cells. However, the antisense RNA transcripts varied in size and between the two cell lines, with two antisense transcripts with sizes of ∼1.4 and ∼5.5 to 6.0 kb in JSC-1 cells but a single antisense transcript of ∼7.5 kb in Raji cells. The oVM294 probe upstream of the pA site at nt 86049 (−) detected two distinct antisense RNA transcripts only in Raji cells with sizes of ∼3.5 and ∼5.8 kb. In contrast, the oVM295 probe specific for the pA site at nt 87185 (−) detected two antisense RNA transcripts only in JSC-1 cells with estimated lengths of ∼2.0 and ∼3.5 to 4.5 kb but not in Raji cells, as expected, due to the deletion in this region (from nt 86838 to 89830) in the Raji EBV genome (62, 79). The observed variation in the expression of the full-length antisense transcripts to EBNA-3B and -3C between JSC-1 and Raji cells was in agreement with our 3′ RACE results (Fig. 8). According to PA-seq read counts, all three species of antisense RNA transcripts were detected only in cells with lytic infection. Due to their expression characteristics, we named these species of seven antisense RNA transcripts lytic anti-EBNA-3B transcript a (∼1.4 kb), transcript b (∼5 to 6 kb), transcript c (7.5 kb), transcript d (∼3.5 kb), and transcript e (∼5.8 kb) and lytic anti-EBNA-3C transcript a (∼2.0 kb) and transcript b (∼3.5 to 4.5 kb). No such transcripts were detected from EBV-negative BCBL-1 cells, indicating their origin from the EBV genome.

Considering that the antisense transcripts to EBNA-2 were also detectable from total RNA extracted from all three cell lines during lytic EBV reactivation (Fig. 8), we further showed that these EBNA-2 antisense transcripts could be enriched from poly(A)-selected RNA and became easily detectable by 3′ RACE from all three EBV-positive cell lines only with lytic EBV infection but not from EBV-negative BCBL-1 cells (Fig. 9C). The lytic induction of EBNA-2 antisense transcripts was much stronger than that of the latent RPMS1/BART miRNA transcripts in all three cell lines (Fig. 9C). By Northern blot analysis we also demonstrated that antisense EBNA-2 transcripts were similar in size, at ~2.7-kb, in both Akata and Raji cells with lytic EBV infection (Fig. 9D, left). Although EBNA-2 was inducible in both JSC-1 and Akata cells (Fig. 9C), the highly expressed EBNA-2 in Raji cells during latent infection was largely reduced in lytic EBV infection, as determined by 3′ RACE (Fig. 9C) and Northern blot analyses (Fig. 9D, right). Interestingly, this reduction of EBNA-2 RNA transcripts in Raji cells with lytic EBV infection was accompanied by the appearance of a smaller product, presumably EBNA-LP (Fig. 9D, right). This inverse expression correlation between EBNA-2 and antisense EBNA-2 transcripts suggests that the expression of antisense EBNA-2 transcripts during EBV lytic infection in Raji cells might contribute to the decreased expression of latent EBNA-2 (Fig. 9C and D).

Cells.

KSHV- and EBV-double-positive JSC-1 (47), EBV-single-positive Raji (92) and Akata (49), and KSHV-single-positive BCBL-1 (93) cell lines were cultivated in HEPES-buffered RPMI 1640 medium (catalog no. 22400-089; Thermo Fisher Scientific) supplemented with 10% HyClone defined fetal bovine serum (catalog no. SH3007.03; GE Healthcare) and 1× penicillin-streptomycin-glutamine (catalog no. 10378-016; Thermo Fisher Scientific). To induce viral lytic replication, all cell lines were adjusted to 5 × 10⁵ cells/ml. To induce lytic replication, JSC-1 cells were treated for 48 h with a final concentration of 3 mM sodium salt of n-butyric acid (sodium butyrate; Bu) (catalog no. B5887; Sigma-Aldrich), Raji cells were treated for 48 h with a combination of 3 mM sodium butyrate and 20 ng/ml of phorbol 12-myristate 13-acetate (TPA) (catalog no. P8139; Sigma-Aldrich), BCBL-1 cells were treated for 48 h with 1 mM sodium salt of valproic acid (VA) for BCBL-1 (catalog no. P4543; Sigma-Aldrich), and Akata cells were treated for 24 h with goat anti-human serum (catalog no. 55087; MP Biologicals) at a 1:200 dilution. After induction, the cells were harvested by centrifugation and lysed in TriPure reagent (catalog no. 11 667 165 001; Roche) for total RNA isolation. The untreated cells were harvested in parallel for analysis of latent gene expression.

PA-seq and data analysis.

The polyadenylation sites were determined by using a PA-seq method described in detail previously (25, 46). The obtained sequence libraries were mapped to human (version hg19; UCSC), KSHV (GenBank accession no. U75698), and EBV (GenBank accession no. NC_007605) reference genomes by the Burrows-Wheeler alignment (BWA) tool (94) and further processed by SAMtools (95). The reads from Raji cells were also mapped to the Raji-specific genome sequence (GenBank accession no. KF717093), and reads from Akata cells were mapped to the Akata-specific genome (GenBank accession no. KC207813). All identified pA cleavage sites were extracted and converted to wiggle (.wig) files for visualization of the nucleotide positions and frequencies using IGV software (http://software.broadinstitute.org/software/igv/). The uniquely mapped reads were then analyzed for viral polyadenylation sites by peak calling using the F-Seq program on combined libraries (25, 46, 96). To eliminate false peaks, we removed all peaks containing a stretch of “A’s” resulting from internal priming and low-abundance peaks with fewer than 50 reads per peak. The peak mode, representing the nucleotide position with the highest number of reads, was then defined as the pA site. To obtain the relative expression levels during latent and lytic infection, the reads were mapped back to the original libraries. The heat map expressing read abundances in individual samples was generated by conditional formatting in Excel software (Microsoft) after normalization to library size and expressed as reads per million (RPM) with maximum scales of 50 RPM for latent and 1,000 RPM for lytic samples. The pA site usage in individual samples was also calculated as a percentage of reads mapped to each pA site from total reads mapped in that particular sample. To identify cis-regulatory sequence elements surrounding viral pA sites, we extracted DNA sequences ±50 nt from each pA site. Sequence conservation was determined by sequence logos (see http://weblogo.berkeley.edu/logo.cgi) (97). The frequencies of PASs were determined by manual inspection of the upstream 50 nt of each pA site.

3′ RACE.

3′ RACE was performed by using the SMARTer RACE 5′/3′ kit (catalog no. 634858; Clontech/TaKaRa), as recommended by the manufacturer (see Table S3 in the supplemental material for primer sequences used). After reverse transcription, 3′ RACE products were amplified for 35 cycles by PCR using a universal primer in combination with a gene-specific primer. The obtained PCR products were analyzed on an agarose gel, gel purified, and directly sequenced by the Sanger sequencing method. Their sequences were mapped back to the EBV genome, and a nucleotide immediately upstream of a pA tail was considered a pA site (cleavage site).

Northern blot analysis.

Northern blotting was carried out as described previously (98), using 5 μg of total RNA or 1 μg of the polyadenylated RNA fraction isolated by using an Illustra QuickPrep mRNA purification kit (catalog no. 27-9254-01; GE Healthcare). After separation on a 1% formaldehyde-containing agarose gel, the RNA was transferred to a nylon membrane and hybridized with an antisense DNA oligonucleotide probe generated by [γ-³²P]ATP end labeling (see Table S3 in the supplemental material). The specific signals were collected and processed by using the Typhoon Trio imaging system (GE Healthcare). The riboprobes for detection of EBNA-2 transcripts (probe B) and antisense transcripts to EBNA-2 (probe A) by Northern blotting were synthetized by in vitro transcription using Promega’s T7 riboprobe system (catalog no. P1440) on PCR templates in the presence of [α-³²P]GTP and gel purified, and 5 × 10⁶ counts per minute (CMP) were used for hybridization. PCR products were generated by amplification of part of the Raji EBV genome using oVM434 and oVM285 for riboprobe A and oVM432 and oVM433 for riboprobe B.