Genome-wide analyses of Epstein-Barr virus reveal conserved RNA structures and a novel stable intronic sequence RNA

Background Epstein-Barr virus (EBV) is a human herpesvirus implicated in cancer and autoimmune disorders. Little is known concerning the roles of RNA structure in this important human pathogen. This study provides the first comprehensive genome-wide survey of RNA and RNA structure in EBV. Results Novel EBV RNAs and RNA structures were identified by computational modeling and RNA-Seq analyses of EBV. Scans of the genomic sequences of four EBV strains (EBV-1, EBV-2, GD1, and GD2) and of the closely related Macacine herpesvirus 4 using the RNAz program discovered 265 regions with high probability of forming conserved RNA structures. Secondary structure models are proposed for these regions based on a combination of free energy minimization and comparative sequence analysis. The analysis of RNA-Seq data uncovered the first observation of a stable intronic sequence RNA (sisRNA) in EBV. The abundance of this sisRNA rivals that of the well-known and highly expressed EBV-encoded non-coding RNAs (EBERs). Conclusion This work identifies regions of the EBV genome likely to generate functional RNAs and RNA structures, provides structural models for these regions, and discusses potential functions suggested by the modeled structures. Enhanced understanding of the EBV transcriptome will guide future experimental analyses of the discovered RNAs and RNA structures.

Despite intense investigation for more than 50 years, the precise roles played by the virus in these diseases remains to be elucidated.
The~170000 bp genome of EBV is a linear doublestranded (ds)DNA that circularizes to form an episome in the host cell nucleus. Infection occurs by entry of the EBV virion into human host epithelial cells and initially proceeds via an aggressive lytic phase. The virus migrates to B cells where it causes persistent lifelong infections marked by extended periods of latency with interspersed lytic reactivation [14]. EBV latency proceeds via three distinct programs, each expressing a different set of coding and non-coding viral gene products. Viral latent gene products rewire B cells to evade the host immune system and propagate the virus [15]. In a manner not yet fully understood, this rewiring increases the tumorigenic potential of EBV-infected cells, including making infected cells resistant to apoptotic pathways that would otherwise kill cancerous cells [16].
The EBV transcriptome is complex, consisting of multiple pre-mRNA transcripts that undergo extensive alternative splicing events and can, at various states of infection, yield different sets of gene products [17]. Two well-studied EBV-encoded ncRNAs (EBER1 and EBER2) are known to be expressed throughout infection [18], as are subsets of as many as 50 different EBV-encoded microRNAs (miRNAs) [19,20]. Although the precise functions of the EBERs have remained obscure, they are the most highly expressed RNAs in EBV-infected cells (~10 7 copies per cell [18]). Roles in cancer have been proposed [21] and the EBERs are capable of inducing tumors in immunesuppressed mice [22]. Both EBERs form ribonucleoprotein complexes [23][24][25][26], including host proteins, and adopt well-defined secondary structures [27]. Except for the EBERs, RNA structure is a relatively understudied aspect of EBV virology; an analysis of RNA secondary structure in EBV could enable significant advances in understanding this important human pathogen.
RNA secondary structure plays important roles in many viruses, mediating functions as diverse as translation initiation [28], catalysis [29], viral genome packaging [30], alternative splicing [31,32], interactions with the host innate immune system [33], RNA lifetime/stability [34][35][36], and regulation of gene expression [37]. The discovery of RNA structure is an active area of research [38][39][40] and is aided by the availability of sophisticated non-coding (nc)RNA-discovery algorithms [41][42][43]. A successful strategy for discovering conserved RNA structures is implemented in the program RNAz [44]. In this approach, homologous sequences are aligned and divided into fragments. For each fragment, calculated parameters that include measures of the thermodynamic stability (z-score [45]) and conservation of secondary structure are used to make predictions based on a support vector machine (SVM) trained on known ncRNAs, which generates a probability classifier (p-class, scale of 0.0 to 1.0) for the fragment containing conserved secondary structure. The RNAz approach has been used to discover RNA secondary structures in a number of viruses [38,[46][47][48].
Across the entire EBV genome the z-score average is -0.19 and the average p-class is 0.13. RNAz discovered 265 regions that map to known transcripts from the EBV genome (named locus 1 to 265 in Additional files 1 and 2) and have strong evidence of locally conserved and stable RNA structure. These regions (average length~200 nt) are comprised of overlapping windows with high predicted probability of generating conserved RNA structure. The average p-class for the set of predictions was 0.91 with minimum and maximum values of 0.61 and 1.00, respectively (Additional file 1). The thermodynamic z-score is an important component in determining the p-class. This parameter reflects the thermodynamic stability of the native sequence predicted structure versus that of random sequence [45] (e.g. measures greater than expected folding stability). For predicted structured regions in EBV RNA, the average z-score was -2.38 with minimum and maximum values of -1.18 and -9.71, respectively. This analysis,  Figure 1 Phylogenic relationships of analyzed viruses. (Left) Phylogenetic tree showing the three herpesviridae sub-families (alpha, beta, and gamma), where the gamma sub-family is further divided into two genera: gamma-1 (Rhadinovirus) and gamma-2 (Lymphocryptovirus). The tree was generated from data in [51]. (Right) Breakdown of the gamma herpesvirus branch shown on the left. The consensus neighbor-joining tree (100 replicates) was based on the conserved herpesvirus DNA polymerase catalytic subunit (DPOL). Kaposi's sarcoma-associated herpesvirus (KSHV), a gamma-1 herpesvirus, is used as an outgroup to illustrate the phylogenetic relationships between four EBV strains (EBV-1, EBV-2, GD1, and GD2) and two other lymphocryptoviruses used in this study [Macacine herpesvirus 4 (MHV4) and Callitrichine herpesvirus 3 (CaHV3)]. The tree was generated using the Geneious program. however, does not account for long-range RNA interactions that may occur outside the overlapped prediction windows.

Identification of known EBV ncRNAs
The genomic region specifying EBER1 occurs within locus 12 (Additional file 1). The structure model derived for EBER1 (Figure 2A) is almost identical to the model in the Rfam database (ID# RF01789). A single base pair in the Rfam model, between nt 6760 and 6770, is shifted to form a pair between 6646 and 6770, and the predicted closing AU pair of this hairpin is open in the Rfam model ( Figure 2A). Additionally, C6649 is in a bulgeloop in the Rfam model, but the insertion of a G residue across from this nt allows a GC pair to form. EBER2, surprisingly, does not occur in an RNAz-predicted locus. The region containing EBER2 was not calculated to have a favorable z-score; however, this is most likely due to the inability of the primary sequence alignment to capture the conservation of structure for this RNA. EBER2 is less conserved in sequence than EBER1 (50.3% sequence conservation vs. 78.3% for EBER1) making primary sequence alignment difficult in this region. Common helices in EBER2 were not properly aligned, thus reducing the predicted free energy of the consensus structure, which resulted in an unfavorable z-score.
All but two of the EBV miRNAs annotated in the miRBase database [53] occur in RNAz-predicted structured regions. Each of the mature miRNA sequences, when mapped to predicted structure models, falls within a hairpin RNA structure that resembles a canonical pre-miRNA (Additional files 1 and 3). For example, Figure 2B shows three hairpins in locus 209, which include the BART3, BART4, and BART1 miRNAs. Each of these hairpins places the 5P and 3P mature sequences in canonical structural contexts: the mature miRNAs are offset and imperfectly base pair with each other. This is also true for the other 36 mature miRNA sequences that fall within predicted structured regions (Additional file 3).
Finally, a reported EBV-encoded viral small nucleolar RNA (v-snoRNA; [52]) also falls within an RNAz-predicted structured region (locus 235 in Additional file 1). The v-snoRNA corresponds to nt 154352 to 154416 (Additional file 3) and, in isolation, has a predicted folding free energy and z-score of -17.8 kcal/mol and -0.81, respectively ( Figure 2C).
In addition to the regions corresponding to EBER1, the v-snoRNA, and 14 loci specifying EBV pre-miRNAs, there are 249 additional EBV regions with predicted conserved RNA structure (Additional file 1).

Structure within introns
Sixty predicted structured regions overlap introns. Introns are a fertile source of structured ncRNAs [54,55] that play roles in a number of important biological processes: snoRNAs [56], miRNAs [57], and piRNAs [58] are expressed from introns. RNA structure within introns can alter accessibility or distance between functional elements and thus regulate splicing [59].
The most extensive regions of predicted structure in introns occur within a repetitive region of the EBV genome (the W repeats). The number of W repeats can vary, and the optimal number appears to be five to eight [60]. The W repeats are transcribed at the 5′ end of the long (~100000 nt) primary transcript encoding the Epstein-Barr nuclear antigen (EBNA) latent proteins ( Figure 3). EBNA transcripts are particularly interesting because they have high nuclear-to-cytoplasmic abundance, suggesting noncoding nuclear functions [61].
Six EBNA proteins are expressed from mRNAs produced by extensive alternative splicing. Transcription of the EBNA RNA can occur from an upstream promoter, the C promoter (Cp), or a downstream W promoter (Wp) [62]. Early in latency III Wp is utilized, where splicing of W0 to W1 generates the start codon for EBNA-LP production [63]; later in latency III there is a switch to Cp where the non-coding C2 exon is joined to W1 to create the EBNA-LP start codon ( Figure 3). Within the W repeats two short coding exons, W1 and W2, which are 66 and 132 nt long, respectively, are joined to form the EBNA-LP protein open reading frame (ORF) by removal of a long (2791 nt) and short (81 nt) intron. RNAz predicts widespread conserved and stable RNA structure covering 49% of the long W repeat intron. There are five predicted regions, which are identical in sequence in the W repeats ( Figure 3 Figure 3). Cp occurs just upstream of locus 17 and the C1 non-coding exon is within locus 17 ( Figure 3).The first copy of Wp overlaps locus 21, which also contains W0 ( Figure 3).
Within loci 18, 23, 28, 33, 38, 43, 48, and 53 is a highly unusual modeled RNA structure, a very long (586 nt) tetraloop hairpin ( Figure 4A). The INFERNAL package [64,65] Figure 3 Cartoon of the primary transcript encoding the EBNA proteins. Exons are indicated with black or grey rectangles and the genomic location of the end of each protein-coding region is indicated above. The identical W1 and W2 coding exons that compose most of EBNA-LP are shown in grey. The Cp and upstream Wp promoter sites are indicated by bent arrows with the genome coordinates of the transcription start sites given above. The boxed region zooms in to show the 5′ end of the pre-mRNA, including the locations of the predicted long hairpin (HP) and ebv-sisRNA-1, colored yellow and green, respectively, with their genome coordinates given in parentheses. RNAz-predicted structured regions are shown in blue and labeled L17 to L23 for loci 17 to 23 (genome coordinates for these loci are in Additional file 1).
Additional file 4). The EBV, MHV4, and CeHV12 hairpins are each capped with a tetraloop closed by a GC pair ( Figure 4A, B, and C) and only the first nt in the loop varies (UUGG, CUGG, and GUGG, respectively). The more distant CaHV3 hairpin, in contrast, is capped with a seven-nt loop ( Figure 4D). Perhaps the identities of the capping nucleotides are of lesser importance to the function of this long hairpin than the overall conserved structure, which is marked by extensively base paired stem regions interspersed with internal (or bulged) loops ( Figure 4). Internal loops may adopt interesting 2D and 3D structures; for example, the three-by-three internal loop that occurs 20 bps from the hairpin loop in EBV [ Figure 4A, also conserved in MHV4 ( Figure 4B)], may fold into a "3 RRs" motif [67] where the opposing internal loop nucleotides form consecutive sheared purinepurine pairs. Additionally, global thermodynamic z-scores were calculated for each hairpin ( Figure 4): all four have very negative z-scores (ranging from -6.3 to -5.9). Thus, the structures generated by these sequences are more than six standard deviations more thermodynamically stable than random sequences with the same dinucleotide content. This indicates that evolution is acting to preserve base pairing within each sequence, as well as the striking structural homology between these viral hairpins.
Very long RNA hairpins that are structurally similar to the EBV W repeat hairpin have been described in humans [68] and in Caenorhabditis elegans [69]. These structures are excellent substrates for adenosine deaminase acting on RNA (ADAR) editing enzymes, which convert adenosine to inosine (A-to-I editing). Interestingly, a previous study of A-to-I editing in EBV [66] found evidence for an editing site that maps to a region of the genome that we predict to be part of the long hairpin ( Figure 4A). ADARs act on structured RNAs and have preferences for adenosines in particular sequence contexts [70], allowing the prediction of potential A-to-I editing sites and of the extent of editing. The reported editing site is predicted to be moderately edited (12.9% by ADAR2, Additional file 4) and is modeled to be in a CA internal loop ( Figure 4A). There are additional sites predicted to be strongly edited in the EBV W repeat hairpin and in the homologous structures found in MHV4, CeHV12, and CaHV3 (Additional file 4). Few of these strongly predicted sites (seven), however, have apparent conservation (Additional file 4), which is unsurprising, as selection of ADAR sites may be partially stochastic [71].
RNA editing in introns, typically within Alu inserts [71], has been proposed to regulate RNA splicing [72]. Inosine can form wobble pairs with G and C residues and can contribute to RNA-protein interactions: thus A-to-I editing can alter important RNA-RNA and RNA-protein interactions [73]. Editing can alter splicing regulatory elements or create/destroy splice sites [72]. Perhaps editing of this long intronic hairpin also plays a role in splicing or in other functions post-splicing.

Discovery of an EBV-encoded sisRNA
A small RNA library (30-200 nt size range) was constructed using nuclear RNA isolated from cultured human B lymphocytes stably infected with EBV (BJAB-B1 cells [74]). These cells exhibit a latency III program of gene expression, the most transcriptionally active program, which also produces latency I and II transcripts. Nuclear RNAs include not only EBER1 and EBER2 [75] but also intronic sequences that are excised during the splicing process. Unsurprisingly, the abundant EBER1 and EBER2 are present in the library and had the highest number of aligned RNA-Seq reads ( Figure 5 and Table 1). A significant peak (~7.8% EBER1 abundance) corresponding to the v-snoRNA ( Figure 5 and Table 1) was also observed. There was a surprisingly large peak originating from the W repeat region ( Figure 5). This peak corresponds to reads that cover the entirety of the small intron (81 nt) that separates the W1 and W2 ORFs ( Figure 3). Next to the EBERs, this was the largest peak observed. Comparing length-normalized reads between the intron and EBER1, the intron is estimated to be present at 21% the level of EBER1 or roughly on a par with EBER2, which is estimated to be 25% as abundant as EBER1 ( Table 1). The presence of the intron in total nuclear RNA samples was confirmed by Northern blot ( Figure 6A) and RT-PCR ( Figure 6B). RT-qPCR then revealed that the intron is~9-fold enriched in nuclear versus cytoplasmic RNA ( Figure 6C).
Stable intronic sequence (sis)RNAs have been described in Xenopus oocytes [78] and, interestingly, in herpes simplex virus-1 (HSV-1, an α-herpesvirus) [79] and in cytomegalovirus (a βherpesvirus) [80]. In these  Figure 5 Small RNA library RNA-Seq reads aligned with the EBV genome. Reads are colored blue; the EBER1, EBER2, and sisRNA peaks are much larger than the represented area (see Table 1). A cartoon of the EBV genome is shown in black (exons as boxes and introns as thin lines). The locations of latently expressed genes [EBNA-LP, LMPs, EBERs (EBER1 and EBER2), BHRF1 miRNAs, and BART miRNAs (clusters I and II)] are indicated at the top of the figure, including the location of the Cp and first Wp promoter sites for producing mRNA for the EBNA-LP and five other EBNA proteins (additional details are in Figure 3). Peaks for the ebv-sisRNA-1 and v-snoRNA are indicated. Images were generated using the integrated genome viewer (IGV, [76,77]) with data from Additional file 5. other herpesviruses a stable intron, known as the latencyassociated transcript (LAT) in HSV-1, is implicated in maintenance of latency [81]. The EBV sisRNA (ebv-sisRNA-1) is the first of this class of ncRNAs to be described in EBV. The ebv-sisRNA-1 differs from the LAT in two major respects: the LAT is much larger (>2 kb), and the functional form of the LAT is believed to be the lariatintron splicing intermediate, whereas the ebv-sisRNA-1 is likely to be a linear molecule. For both ends of ebv-sisRNA-1 to have been free for ligation in the small RNA library construction step, the 5′ end of the released intron could not have been sequestered in a lariat structure. This suggests that ebv-sisRNA-1 has been acted upon by the debranching enzyme, which also linearizes intron-derived snoRNAs during the maturation of this important class of ncRNAs [82]. Introns homologous to ebv-sisRNA-1 were found in other lymphocryptoviruses ( Figure 6D). The ebv-sisRNA-1 sequence and these RNAs have an APSI of 83.3%, with the CaHV3 sequence being most distant (57.6% APSI). In all analyzed EBV strains (EBV-1, EBV-2, GD1, and GD2), the 81-nt ebv-sisRNA-1 sequence is 100% conserved, whereas the APSI of the genome sequences of these strains is 95.6%. Within the CA-rich region an 11-nt sequence, from positions 49 to 59, is 100% conserved in the analyzed lymphocryptoviruses ( Figure 6D). The distribution of dinucleotide frequencies in the ebv-sisRNA-1 sequence is skewed, with CA being the most frequent, followed by AC and CC. The presence of CA-rich regions in the sisRNA is interesting as such sequences are able to bind hnRNP L protein and modulate splicing [83].
Although ebv-sisRNA-1 did not overlap an RNAzpredicted region, the secondary structure for this RNA was modeled on the basis of the eight sequences shown in Figure 6E. Two small hairpins are predicted ( Figure 6E): an upstream hairpin that places the U-rich motif in a loop and a weak (3 base pair) hairpin in the CA-rich region. A compensatory mutation (double point-mutation that preserves base pairing) is observed in the more stable upstream hairpin ( Figure 6E, a conserved GC pair is converted to a UA pair in CaHV3), which suggests that evolution is acting to preserve this structure. It is interesting to find evidence of a conserved hairpin that places the U-rich motif in a loop, as U-rich sequences often bind proteins [84] and hairpins can pre-organize RNA sequence motifs for RNA-protein interactions [85]. The CA-rich region is likely unstructured and free to associate with other interactors.

Additional putative sisRNAs in EBV
Other than the W repeat sisRNAs, EBV possesses 18 additional short introns ( Table 2) interspersed throughout its genome. These short introns range in conservation from 82.9% to 100.0% APSI (average of 97.2%) and are predicted to fold into structures with a range of thermodynamic stabilities (from -0.4 to -59.6 kcal/mol). Four sequences have significantly negative z-scores [< -1.0, from nt 11466 to 11610, 80362 to 80446, 87014 to 87087, and 105929 to 106056 (Table 2)] and one intron has a moderately favorable z-score (-0.54) from nt 1496 to 1573. These short introns also overlap with RNAz-predicted structured regions (Table 2).
Particularly interesting is the sequence from nt 105929 to 106056, which has the most favorable predicted folding free energy and z-score of any EBV intron ( Table 2). This intron occurs within the lytically-expressed BBLF2/3 gene, which generates part of a helix-primase complex essential for replication of the EBV genome [86]. The modeled secondary structure of the intron is shown in Figure 7A; it forms a tetraloop hairpin that includes all of the intron sequence. This intron structure is contained within RNAz-predicted locus 152, where the basal stem of the hairpin extends into coding regions by 25 nt (Additional file 3, nt 66692 to 66844 in the reverse genome sequence) and the splice donor and acceptor are base paired within helices. A similar hairpin structure, identified by INFERNAL, occurs in the MHV4 genome ( Figure 7B). Rather than being joined together via splicing, MHV4 BBLF2/3 is translated from two ORFs that produce BBLF2 and BBLF3 as separate proteins. This same intronic region in EBV partially codes for BBLF2 in MHV4 (the translation stop codon occurs within the UUAA tetraloop) and is also present in the 5′ UTR of BBLF3 ( Figure 7B). It is remarkable that the hairpin structure has been maintained evolutionarily between EBV and MHV4 despite different apparent roles in protein expression: in EBV it brings splice sites closer together in space and may be processed out to form an independent structured RNA, whereas in MHV4 it presents a stop codon in a loop and is not processed out of the mRNA. Four EBV short introns have significant numbers of mapped reads in the small RNA library (Table 1); however, they are lower in abundance than ebv-sisRNA-1. In the EBV latent membrane protein (LMP-2B) gene, for example, there are three short introns with significant numbers of aligned read peaks (nt 459 to 539, nt 789 to 870, and nt 952 to 1025 in Table 1). Indeed, the intron encompassing nt 952 to 1025 has an abundance of~5.7% relative to EBER1, which is similar to the value observed for the v-snoRNA (7.8%, Table 2). The short introns in EBV, including those that are expressed during latency III and those with few or no mapped reads, may represent a family of sisRNAs. It will be interesting to learn if they are differentially expressed in other viral latency programs or during the lytic phase. Short introns are an attractive source of small ncRNAs, as they are released as byproducts of pre-mRNA splicing.

Structure overlapping splice sites
The importance of RNA structure in the regulation of splicing has been highlighted in two excellent reviews [59,87]. Structure can be inhibitory, by sequestering splice sites [88] or regulatory elements [89], or it can enhance splicing by presenting sites in an accessible conformation [90] or by bringing splice sites into closer proximity with each other [91]. Conformational switching between accessible and inaccessible structures may also regulate splicing [92].
In EBV, there are seven predicted structured regions that include splice sites (Additional file 1). Of these, three include a splice donor or acceptor, while the remaining regions include both sites. For example, in the LMP-2B gene, locus 1 overlaps a splice donor site (after nt 788, Additional file 3), which is buried within the stem of a tetraloop hairpin; towards the 3′ end of this structure (nt 827 to 833) is a putative intronic splicing enhancer ( Figure 8A). Locus 2 spans an entire intron and includes both the donor and acceptor sites (after nt 1495 and 1573, respectively, in Additional file 3). Additionally, there is a putative exonic splicing enhancer at the beginning of locus 2 (nt 1475 to 1482). The modeled structure for locus 2 ( Figure 8B) organizes this region into a complex multibranched structure where the splice donor and acceptor sites are each in helices: the donor at the end of the basal stem and the acceptor within a hairpin loop. Interestingly, each splice site occurs directly after or near a predicted UU/UU twoby-two internal loop motif, in which UU mismatches maintain stable [93,94] helical 3D conformations [95]. Internal loops and non-canonical pairs can present complex patterns of chemical functional groups in proteinaccessible regions of a folded RNA [96], which may influence the binding of splicing factors. Additionally, the 3D folding of this structure may bring splice sites close together .

Structure overlapping coding regions
Many (150) predictions overlap EBV coding regions, in some cases including a significant portion of an ORF. In the lytically-expressed EBV FGAM-synthase gene, for example, predictions cover 40% the coding nt (locus 3 to locus 9; Additional file 1). More dramatically, the BGLF3 gene overlaps four predicted structured regions (locus 163 to locus 166; Additional file 1), which span 84% of the coding region including part of the start codon. Interestingly, many viral genomes appear to be enriched for RNA secondary structure [97,98], including within coding regions [38,99,100]. RNA structure in coding regions can have a variety of functions [101], such as altering reading frames via frame-shifting pseudoknots [37] or by inducing ribosomal pausing [102], which can impact numerous co-translational processes [103,104]. The identification of conserved and stable RNA secondary structures within EBV coding regions suggests that they may play important roles in this virus. Twenty-one predicted structured regions include translation start or stop codons. Locus 263, for example, includes the start codon for the latent membrane protein 2A (LMP-2A), which is involved in the regulation of host B-cell signaling and is implicated in EBV-associated tumorigenesis [105]. The first 63 nt of locus 263 fall within the 5′ UTR of the LMP-2A mRNA ( Figure 8C). The start codon is predicted to occur within the stem helix of a multibranch loop. RNA structure can play a role in regulating translation through modulating the accessibility of codons. Strong RNA structure at stop codons, for example, can stimulate stop codon readthrough [106], whereas modulation of start codon accessibility can affect translational efficiency [107]. Perhaps the stable structure (z-score -1.48, Figure 8C) that encompasses the LMP-2A start codon modulates the expression of this important gene.

Conclusion
Bioinformatics analyses of the EBV genome revealed 265 regions with putative conserved RNA structure. The computational screen successfully identified most of the known EBV ncRNAs: EBER1, a virally-encoded snoRNA, and 42 miRNAs. Other predicted structured regions of the EBV genome occur at or near sites that suggest potential functions for RNA structures. A very long hairpin structure in the W repeat region that is conserved in lymphocryptoviruses contains at least one A-to-I editing site, with possible implications for splicing. Additional regions span splice sites or start and stop codons. These structure predictions may prove useful since RNA and RNA structures are attractive targets for chemotherapy [108][109][110][111]. The modeled RNA structures are rich in internal loops, e.g. the UU/UU internal loops near two splice sites in LMP-2B, which are potentially drugable structural motifs [112]. RNA splicing, in particular, is an attractive target for RNA therapeutics [113,114].
RNA-Seq data identified several putative sisRNAs, including one in the W repeat region (ebv-sisRNA-1) that is abundant and localizes in the nucleus. These results lay the groundwork for future experiments to further characterize the proposed structures and to understand the functional significance of the discovered RNAs and RNA structures in EBV. All structural data are available in the Additional files. In addition, data for the W repeat long hairpin and ebv-sisRNA-1 will be deposited in the Rfam database.

Methods
Structure prediction and modeling EBV RefSeq genomes from four different viral strains, EBV-1 (NC_007605.1), EBV-2 (NC_009334.1), GD1 (AY961628.3), and GD2 (HQ020558.1), plus one closely-related genome from Macacine herpesvirus 4 (NC_006146.1), were obtained from the NCBI nt database. Genome sequences were aligned using the MAFFT program [49,50] implementing the FFT-NS-1 alignment strategy. The alignment was processed into 120-nt windows with a 10-nt step-size using the Perl script rnazWindow.pl provided with the RNAz program. Predictions of windows with likely conserved and stable RNA structure were made for both the plus and minus strands with the RNAz 2.1 [42,44] program using a dinucleotide randomization model. The RNAz output was processed to combine windows with higher probability of structure (p-class > 0.5) into loci. These loci were further filtered to include only those that had at least one window with p-class > 0.9. When predictions had ambiguous strand orientations, e.g. occurred in a genome location that could be transcribed from either strand, the program RNAstrand [115] was used to determine if structure was more likely to be in the forward or reverse strand.
From each predicted locus the EBV-2 sequence was extracted and used in a BLAST [116] search against the NCBI nt database. The returned homologous sequences were compiled and each locus was aligned using the MAFFT program implementing the Q-INS-i method. Initial models for these alignments were made using RNAalifold [117]. Alignments and predicted consensus structures were converted to Stockholm format alignments and used to build covariance models with the INFERNAL package [64]. These covariance models were used to search for potentially homologous sequences/structures in a database of all available herpesvirus RefSeq genomes. Returned sequences were added to the alignments and structural predictions were made using RNAalifold as well as the program TurboFold [118] followed by manual model refinement. Resulting models were used to perform additional iterative INFERNAL searches and structure revisions until no further model improvement could be made (finding the structure that was best conserved in the alignment).
2D renderings of modeled RNAs were generated using the PseudoViewer program [119] and manually processed using the open-source graphics program Inkscape and ImageJ [120,121]. The Gibbs free energy of folding at 37°C [122,123] for models was calculated using the efn2 method implemented in the RNAstructure package [124]. The z-scores for sequences were calculated by generating sets of 200 dinucleotide randomizations using the SIMMONICS program [125] and then calculating energies with RNAfold [126]. The difference in native free energy versus the randomized set was normalized by the standard deviation to give the z-score.
Small RNA-Seq analysis EBV-positive BJAB-B1 cells were maintained at 37°C and 5% CO 2 in RPMI media supplemented with FBS (10%) and penicillin [20]. Cells that had been expanded from 1 ml frozen stock to 200 ml culture volume were grown for an additional 48 hrs and were then pelleted by centrifugation at 1000 × g for 5 min at 4°C. Pellets were washed 3× with ice-cold PBS buffer (plus DTT). The pellet, on ice, was resuspended in cell disruption buffer [127] (10 mM KCl, 1.5 mM MgCl2, 20 mM Tris-Cl, 1 mM DTT, 0.1% Triton-X) and moved to a chilled Dounce homogenizer (Kontes, type B). Homogenization proceeded on ice until mostly free nuclei were observed under microscope. The lysate was transferred to centrifuge bottles and spun at 1500 × g for 5 min at 4°C. The supernatant, containing the cytoplasmic fraction, was removed to a fresh tube and a 3× volume of Trizol reagent was added. The nuclear fraction containing pellet was dissolved in 3 ml of Trizol reagent and both fractions were flash frozen in a dry-ice/ethanol bath and stored at -80°C overnight. Each sample was chloroform extracted and then precipitated at room temperature for 10 min using isopropanol. RNA was pelleted by centrifugation at 10000 × g for 10 min at 4°C. Pellets were resuspended in water and treated with RNase-free DNase I according to manufacturer's protocol (Promega). RNA was further purified/recovered by RNeasy RNA purification kit (Qiagen). Samples were treated again using DNase I to remove any traces of genomic DNA contamination and recovered using the RNeasy kit. RNA integrity was checked using an Agilent BioAnalyzer: the RNA integrity number for the undigested samples was equal to 10, indicating that sample annealing did not degrade the RNA. RNA samples were used to build libraries using a small RNA library kit (Illumina), where specific adapters were ligated to the 5′ and 3′ ends, before barcoding and cDNA synthesis. The cDNA was fragmented and size-selected to contain fragments in the range of 30 to 200 bp. RNA-Seq data, 75 bp reads, were acquired on an Illumina HiSeq instrument run in paired-end mode. Reads were aligned to the EBV-2 RefSeq genome (NC_009334.1) using the Bowtie 2 program [128] (data given in Additional file 5).
Analysis of ebv-sisRNA-1 expression 20 μg of total RNA from EBV-positive BJAB-B1 and EBV-negative (but otherwise isogenic) BJAB cells were fractionated by denaturing 10% PAGE. RNA was transferred to an Amersham Hybond-N + positively charged nylon membrane then cross-linked to the membrane using UV radiation. The membrane was incubated at 37°C for 1 hr in pre-hybridization buffer (GE) and washed 3× in double-distilled water (ddH 2 O). The membrane was then incubated overnight at 37°C with 5′-radiolabeled DNA probe complementary to nts 37 to 66 of the 106 nt human U6 snRNA (NR_046494.1), probe sequence: 5′-GCAGG GGCCATGCTAATCTTCTCTGTATCG-3′. The blot was washed 3× in wash buffer (0.005% SDS and 3× SSC), wrapped in plastic film, exposed to a phosphor screen for 16 hrs and imaged on a Storm phosphorimager. The blot was rehybridized with 5′-radiolabeled DNA probe complementary to nt 19 to 46 of ebv-sisRNA-1, probe sequence: 5′-TGTGGTGGAGTGTTGGGCTTA GCAGAAA -3′.
For RT-qPCR a master mix was made of FastStart Universal SYBR Green Master Mix (Roche), ddH 2 O and each primer (final concentration of 300 nM). 45 μl aliquots were added to a 96-well plate before adding 5 μl 100× dilute cDNA (three biological replicates, each run with three technical replicates, including no template and reverse transcriptase controls). In addition to the ebv-sisRNA-1 primers, amplification of nt 8 to 104 of the nuclear-enriched U2 snRNA (NR_002716.3) used primers: FWD 5′-CTCGGCCTTTTGGCTAAGAT-3′ REV 5′-TA TTCCATCTCCCTGCTCCA-3′; and nt 994 to 1071 of the cytoplasmically-enriched 18S rRNA (NR_003286.2) used primers: FWD 5′-CGAAAGCATTTGCCAAGA AT-3′ REV 5′-GCATCGTTTATGGTCGGAAC-3′. The same thermocycler protocol used for PCR was used for the qPCR run, and data was collected on an Applied Biosystems StepOnePlus Real-Time PCR system. Data were exported to Microsoft Excel and analyzed using the Delta Ct method to find the ratio of nuclear to cytoplasmic abundance.