Long-read sequencing uncovers a complex transcriptome topology in varicella zoster virus

Abstract

Background

Varicella zoster virus (VZV) is a human pathogenic alphaherpesvirus harboring a relatively large DNA molecule. The VZV transcriptome has already been analyzed by microarray and short-read sequencing analyses. However, both approaches have substantial limitations when used for structural characterization of transcript isoforms, even if supplemented with primer extension or other techniques. Among others, they are inefficient in distinguishing between embedded RNA molecules, transcript isoforms, including splice and length variants, as well as between alternative polycistronic transcripts. It has been demonstrated in several studies that long-read sequencing is able to circumvent these problems.

Results

In this work, we report the analysis of the VZV lytic transcriptome using the Oxford Nanopore Technologies sequencing platform. These investigations have led to the identification of 114 novel transcripts, including mRNAs, non-coding RNAs, polycistronic RNAs and complex transcripts, as well as 10 novel spliced transcripts and 25 novel transcription start site isoforms and transcription end site isoforms. A novel class of transcripts, the nroRNAs are described in this study. These transcripts are encoded by the genomic region located in close vicinity to the viral replication origin. We also show that the ORF63 exhibits a complex structural variation encompassing the splice sites of VZV latency transcripts. Additionally, we have detected RNA editing in a novel non-coding RNA molecule.

Conclusions

Our investigations disclosed a composite transcriptomic architecture of VZV, including the discovery of novel RNA molecules and transcript isoforms, as well as a complex meshwork of transcriptional read-throughs and overlaps. The results represent a substantial advance in the annotation of the VZV transcriptome and in understanding the molecular biology of the herpesviruses in general.

Background

Varicella zoster virus (VZV) is one of the nine herpesviruses that infect humans. It is the etiological agent of chickenpox (varicella) caused by primary infection and shingles (zoster) due to reactivation from latency, which is established in the spinal and trigeminal ganglia [1]. VZV belongs to the Varicellovirus genus of the subfamily Alphaherpesvirinae, and is closely related to the pseudorabies virus (PRV) a veterinary Varicellovirus and to the herpes simplex virus type-1 (HSV-1) belonging to the Simplexvirus genus of alphaherpesviruses. The investigation of the VZV pathomechanism is not easy due to its highly cell-associated nature and the lack of animal models, except the SCID mouse model [2, 3]. The VZV virion is ~ 80–120 nm in diameter and is composed of an icosahedral nucleocapsid surrounded by a tegument layer [4]. The outer virion component is an envelope derived from the host cell membrane with incorporated viral glycoproteins [5]. The VZV genome is composed of a ~ 125-kbp double-stranded DNA molecule containing at least 70 annotated open reading frames (ORFs). The viral DNA consists of two unique genomic regions (UL and US), and a pair of inverted repeats (IRs) surrounding the US region. Three genes (ORF62, ORF63 and ORF64) located at the IR regions are therefore in duplicate [6, 7]. VZV encodes six genes (ORFs S/L, 1, 2, 13, 32, and 57), which are not present in HSV [8]. It has been long held that VZV does not express the RNA molecule homologues to the HSV latency-associated transcripts (LAT), however, a recent study reported the identification of VZV latency-associated transcript (VLT) overlapping the ORF61 gene [9]. It has also been thought that VZV lacks the α-TIF protein encoded by the HSV’s vp16 gene that activates the transcription of immediate-early (IE) genes during the initial events of the virus life cycle [10]. VZV ORF10 together with ORF4, 62 and 63 are transcriptional transactivators that are all present in the virus tegument [11, 12]. ORF10 is supposed to substitute the function of VP16 in the transactivation of ORF62 (homologue of HSV ICP4), but unlike VP16, it is unable to control the expression of ORF4, ORF63, and ORF61 [13, 14].

The viral replication is primarily controlled by the regulation of transcription, which is carried out through a sequential activation of viral genes. First, the IE genes are expressed, which is followed by the activation of the early (E) and then the late (L) kinetic classes of genes [15]. The IE protein ORF62 is the major transactivator of the viral genes, which recruits the general transcriptional apparatus of the host cell and thereby controls the expression of other viral genes [16, 17]. The VZV E genes encode the proteins required for the DNA replication, while L genes code for the structural elements of the virus. Polycistronism represents an inbuilt characteristic of the herpesvirus genome [18, 19]. However, the herpesvirus multigenic transcripts differ from the bacterial polycistronic transcripts in that the downstream genes on the viral RNA molecules are untranslated. Exceptions to this rule were found in the Kaposi’s sarcoma-associated herpesvirus [19, 20] and in HSV-1 [21]. The functional significance of this organization of the herpesviral transcripts has not yet been ascertained. Alternative splicing expands the RNA and protein repertoire with respect of the one gene, one RNA/protein situation. In contrast to the beta and gammaherpesviruses, splicing events are rare among the alphaherpesviruses [22, 23]. In VZV only a few genes have been shown to produce spliced mRNAs so far, which include the ORF0 (also referred as ORFS/L) located at the termini of UL region, the UL15 homologue ORF42/45, ORF50 the glycoprotein M (gM), and the newly discovered, multiple spliced VLT [6, 9, 23,24,25,26].

RNA editing includes the adenosine deaminase acting on RNA1 (ADAR1) enzyme that catalyzes the C-6 deamination of adenine (A) to inosine (I). The I is recognized as guanine (G) by the reverse transcriptase (RT) enzyme [27], therefore it can be detected with cDNA sequencing methods. In hyper-edited sites the ADAR1 transforms multiple adjacent As on an RNA strand. A G-enriched neighborhood has been described at the RNA editing sites [28] with an upstream uracil (U) exerting a stabilizing effect on the RNA-ADAR complex [29]. Additionally, RNA hyper-editing has been shown to play a crucial role in the viral life cycle by dodging the host’s immune system [30] and also in the control of DNA replication [31].

The translation of eukaryotic mRNAs occurs according to the scanning model, where the 40S ribosomal subunits scans the RNA strand in 5′ to 3′ direction and initiates the translation at the first AUG they encounter [32]. Mammalian mRNAs contain an essential sequence context for translation initiation, known as the Kozak consensus sequence [33]. A purine at − 3 and G at + 4 position has the strongest binding effect for the translational machinery, while AUGs with a different context tend to be overlooked by the ribosomal subunits (leaky scanning) [34]. Upstream ORFs in the 5′-untranslated regions (UTRs) of the RNAs have been shown to exert a regulatory effect on the protein synthesis through a process called translation reinitiation [32, 35]. In short upstream (u)ORFs translation reinitiation shows a positive correlation with distance between stop codon of uORF and the AUG of downstream ORF [36].

The hybridization-based microarray and the second-generation short-read sequencing (SRS) techniques have revolutionized genome and transcriptome research, including herpesviruses [37,38,39,40]. However, both techniques have limitations compared to long-read sequencing (LRS), including Pacific Biosciences (PacBio) and Oxford Nanopore Technologies (ONT) platforms. Microarray and SRS techniques perform poorly in the detection of multi-intronic transcripts, transcript length variants, and polycistronic RNA molecules, as well as transcriptional overlaps. Additionally, the LRS sequencing platforms are capable of determining the transcript ends with base-pair precision, without the need of any supplementary method [41, 42].

Similar to other LRS techniques, ONT cDNA sequencing is afflicted by sample degradation, resulting in false transcriptional start sites (TSSs). This problem can be mitigated by using cap selection. Another source of non-specificity is the presence of nucleotide sequences complementary to the MinION strand switching oligonucleotide, which can lead to either template switching [43, 44] or false priming [45]. False priming with the oligo(dT) primer can also occur at A-rich regions of the transcripts. These errors contribute to artifactual 3′-ends of the reads.

LRS techniques have already been used in genome and transcriptome studies of viruses [46,47,48,49], including herpesviruses [18, 22, 23, 50,51,52]. These reports have revealed a hidden transcriptome complexity, which included the discovery of long RNA molecules and a large variation of transcript isoforms. ONT and PacBio-based studies have also detected a number of embedded transcripts with in-frame truncated (t)ORFs in several herpesviruses, such as PRV [18, 46], HSV-1 [53], and human cytomegalovirus (HCMV) [22]. Additionally, a large extent of transcriptional read-through and overlaps has also been uncovered by these studies [18, 22, 23]. It has been proposed that the transcriptional overlaps may be the result of a transcriptional interference mechanism playing a role in genome-wide regulation of gene expressions [54].

In this work, we used the ONT MinION sequencing technique to investigate the structural aspects of the polyadenylated VZV transcripts. We report numerous novel transcripts, transcript isoforms, and yet unknown splice events. This study also explores a complex meshwork of transcriptional overlaps. Additionally, we report and analyze a hyper-editing event in a VZV transcript.

Methods

Cells and viral infection

Human primary embryonic lung fibroblast cell line (MRC-5) was obtained from the American Type Culture Collection (ATCC) and grown in DMEM supplemented with antibiotic/antimycotic solution and 10% fetal bovine serum (FBS) at 37 °C in a 5% CO₂ atmosphere. The live attenuated OKA/Merck strain of varicella zoster virus (VZV) was cultured at 37 °C in MRC-5 cell line, the cells were harvested by trypsinization, when the monolayers had displayed specific cytopathic changes. For subsequent propagation of the virus, infected cells were used to inoculate MRC-5 cultures previously grown to full confluence at a ratio of 1:10 infected to uninfected cells. The cultures were then incubated at 37 °C for 5 days, when the cytopathic effect was near 100%.

RNA purification

Total RNA was isolated immediately after collecting the infected cells, using the Nucleospin RNA Kit (Macherey-Nagel) according to the manufacturer’s instruction. Briefly, infected cells were collected by centrifugation and the cell membrane was disrupted using the lysis buffer provided in the kit. Genomic DNA was digested with the RNase-free rDNase solution (supplied with the kit). Samples were eluted in a total volume of 50 μl nuclease-free water. To eliminate residual DNA contamination a subsequent treatment with the TURBO DNA-free Kit (Thermo Fisher Scientific) was conducted. The RNA concentration was measured with a Qubit 2.0 Fluorometer, with the Qubit RNA BR Assay Kit (Thermo Fisher Scientific).

Poly(a) selection

Thirty-five μg of total RNA was pipetted in separate DNA LoBind Eppendorf tubes (Merck). The poly(A)⁺ RNA fraction was isolated from the samples using the Oligotex mRNA Mini Kit (Qiagen). RNA samples were stored at − 80 °C until use.

Oxford Nanopore MinION cDNA sequencing and barcoding

The cDNA library was prepared using the Ligation Sequencing kit (SQK-LSK108; Oxford Nanopore Technologies) according to the modified 1D strand switching cDNA by ligation protocol. In short, the polyA(+) RNA fraction was reverse transcribed using an oligo(d)T-containing primer [(VN)T20 (Bio Basic, Canada)]. The RT reaction was carried out as recommended by the 1D protocol, using SuperScript IV enzyme (Life Technologies); a strand-switching oligo [containing three O-methyl-guanine RNA bases (PCR_Sw_mod_3G; Bio Basic, Canada)] was added to the sample. The cDNA was amplified by using KAPA HiFi DNA Polymerase (Kapa Biosystems) and Ligation Sequencing Kit Primer Mix (part of the 1D Kit) following the ONT 1D Kit’s manual. End repair was carried out on the samples using NEBNext End repair / dA-tailing Module (New England Biolabs), which was followed by “barcoding”: the C11 barcode (ONT PCR Barcoding Kit 96; EXP-PBC096) was ligated to the sample according to the 1D PCR barcoding (96) genomic DNA (SQK-LSK108) protocol, Barcode Adapter ligation step. The barcoded- sample was amplified by PCR using KAPA HiFi DNA Polymerase, as well as the C11 PCR barcode according to the 1D PCR barcoding protocol. The PCR product was end-repaired, then it was followed by adapter ligation using the sequencing adapters supplied in the kit and NEB Blunt/TA Ligase Master Mix (New England Biolabs). The cDNA sample was purified between each step using Agencourt AMPure XP magnetic beads (Beckman Coulter). The concentration of cDNA library was determined using a Qubit 2.0 Fluorometer through use of the Qubit (ds)DNA HS Assay Kit (Thermo Fisher Scientific). Samples were loaded on R9.4 SpotON Flow Cells, while the base calling was performed using Albacore v2.1.10.

Cap selection, cDNA synthesis and sequencing

For the precise determination of TSSs, the ONT’s 1D strand switching cDNA by ligation protocol was combined with a 5′-cap specific protocol. The cDNA sample was prepared from total RNA using the TeloPrime Full-Length cDNA Amplification Kit (Lexogen). Reverse transcription (RT) reaction was performed at 46 °C for 50 min according to the kit’s recommendations by using an oligo(d)T-containing primer (5′ - > 3′: TCTCAGGCGTTTTTTTTTTTTTTTTTT). The RT product was purified by using spin column-based method (silica columns are from the Lexogen kit). A double-strand specific ligase enzyme (Lexogen Kit) was used to ligate an adapter to the 5’C of the cap of the RNA. The reaction was carried out at 25 °C, overnight. The sample was washed by applying the silica-column method. The Second-Strand Mix and the Enzyme Mix (both from the TeloPrime Kit) as well as a Verity Thermal Cycler (Applied Biosystems) were used to produce double-stranded cDNAs according to the Kit’s guide. In order to generate sufficient amount of cDNAs for MinION library preparation, samples were amplified by endpoint PCR following the TeloPrime Kit’s manual. The generation of the sequencing-ready library from this sample is based on the 1D strand switching cDNA by ligation protocol from ONT. The Ligation Sequencing kit (SQK-LSK108, ONT) and NEBNext End repair / dA-tailing Module (New England Biolabs) was used to repair the cDNA ends. This step was followed by the 1D adapter ligation, which was carried out according to the 1D protocol, using the NEB Blunt/TA Ligase Master Mix (New England Biolabs). We used barcoding for the better in silico identification of the transcripts’ ends. We found this method useful because it helped the base-pair precision mapping of TSSs. The library was sequenced on an ONT R9.4 SpotON Flow Cells.

Targeted sequencing of putative ONT transcripts

Validation of the TSSs of three ‘near the transcription origin’ (NTO) transcripts, and the spliced isoforms of the region was carried out using targeted long-read sequencing. The cDNA library was prepared using the ONT 1D protocol, as described in the ‘Oxford Nanopore MinION cDNA sequencing and barcoding’ section, with the following modifications: we used ribodepletion instead of poly(A) selection, and target-specific primers (Integrated DNA Technologies) were applied instead of oligo(dT)-priming for the cDNA preparation. The target-specific primers used in this study are listed in Table 1.

Table 1 Target-specific primers used for the detection of NTO transcripts

Data analysis and alignment

Reads resulting from ONT sequencing were aligned to the reference genome of VZV (GeneBank accession: NC_001348.1) and the host cell genome (Homo sapiens - GRCh37, BioProject number: PRJNA31257) using GMAP v2017-04-24 with the default settings [55]. In order to annotate the 5′ ends of the transcripts, we used our custom script (https://github.com/zsolt-balazs/LoRTIA). In short: the last 16 nt of the MinION 5′ strand switching adapter or the last 16 nt of the cap selection adapter was aligned in a window of -10 nt and + 30 nt from the first mapped position of a read using the Smith-Waterman algorithm with a match cost of + 2 or − 3, a gap opening cost of − 3, and a gap extension cost of − 2. Read ends with a score below 14 were considered false 5′ ends and were discarded. A 5′-end position was considered a TSS if the number of reads starting at this position was significantly higher than at other nucleotides in the region surrounding this start position. For this the Poisson-probability (Poisson [k₀;λ]) of k₀ read starting at a given nucleotide in the -50 nt to + 50 nt window from each local maximum was calculated with $\lambda =\frac{\sum \limits_{i=-50}^{50}{k}_i}{101}$. The 5′-ends of the rare long reads were inspected individually using the Integrative Genome Viewer (IGV) [56]. Poly(A) tails were defined using the same algorithm and parameters, by mapping 15 nt of homopolymer As or Ts to the soft-clipped region on the read’s end. Read ends with a score below 14 or with more than 5 As/Ts directly upstream their positions were considered artefacts of false priming, and were discarded. Transcription end sites (TESs) were defined using the same criteria for the Poisson-probability as in case of TSSs. The computational pipeline is outlined in Additional file 1.

Reads passing both criteria and present in more than three copies were considered as “certain transcript isoforms”. Those with less than three copies or the longest unique reads with questionable 5′ ends were considered “uncertain transcript isoforms”. Reads with sequencing adapters or poly(A) tails on their both ends were discarded, except for complex transcripts, which were individually inspected using IGV. Reads with a larger than 10 nt difference in their 5′ or 3′ ends were considered novel length isoforms (L: longer 5’ UTR, S: shorter 5’ UTR, AT: alternative 3′ termination). Short length isoforms harboring a truncated version of the known open reading frame (ORF) were considered novel putative protein coding transcripts, and designated as ‘.5′. If multiple putative protein coding transcripts were present, then the one with the longest 5’ UTR was designated ‘.5′ and its shorter versions were labelled in an ascending order. Transcripts with TESs located within the ORFs of genes (therefore lacking STOP codons) or with TSSs within the coding regions without in-frame ORFs were both considered non-coding transcripts.

Multigenic transcripts containing at least two genes standing in opposite were named complex transcripts. If a TSS was not obvious at these transcripts, we assumed that they start at the closest upstream annotated TSSs. Splice junctions were accepted if the intron boundary consensus sequences (GT and AG) were present in at least ten sequencing reads and if the frequency of introns was more than 1% at the given region. The relative abundance of the transcript isoforms was calculated by dividing the number of reads with the total number of reads of the dataset; both resulting from the computational pipeline described before. A ± 10 nt-variation was allowed at both the TSS and TES of sequencing read for considering them as certain transcript isoforms. To assess the homology of protein products possibly translated from the ORFs of novel transcript isoforms we used the online BLASTP suite [57], with an expected threshold of 10.

To evaluate the effect of hyper-editing on the secondary structure of RNAs, we used the ‘RNAstructure Software’ suite [58] with the following parameters: temperature = 310,15 K, Max % Energy Difference = 10, Max Number of Structures = 20, Window Size = 0. To simulate the presence of inosine, we changed the edited adenines to guanine. To assess the secondary structure of the NTO3-ORF62 dsRNA hybrid we used the first 467 bases of the ORF62 labeled ORF62–5′ (from genomic position 109,204 to 108,738) as one of our starting sequences, and the whole sequence of NTO3 as the other one. Reads were visualized using the Geneious [59] software suite and IGV. GC-, CCAAT- and TATA-boxes were annotated using the Smith-Waterman algorithm for canonical motifs [60].

Results

Analysis of the VZV transcriptome with third-generation sequencing

In this work, the ONT MinION RNA-Seq technique was used for profiling the genome-wide expression of the lytic VZV transcripts (Fig. 1). Both cap-selected and non-cap-selected protocols were applied for the analyses of the whole transcriptome. Additionally, targeted sequencing was used with the non-cap-selected protocol for validating the TSS and splice junctions of the NTO transcripts. Sequencing of the non-cap-selected non-targeted sample yielded 619,687 reads of which 66,455 mapped to the VZV reference genome (GeneBank Accession: NC_001348.1) with an average mapped read length of 1349 bp, and an average coverage of 625, while 553,232 reads mapped to the host cell genome (Homo sapiens - GRCh37, BioProject number: PRJNA31257).

Fig. 1

The VZV transcriptome. The transcriptome of the varicella zoster virus was analyzed using long-read cDNA sequencing. The reads were annotated using our custom script (https://github.com/zsolt-balazs/LoRTIA). Transcriptional start sites and transcriptional end sites of the annotated transcripts are the 5′ and 3′-ends with significantly higher abundance than any other read end in their ±50 vicinity, according to Poisson-probability, and those that are not a result of strand switching or false priming

Targeted sequencing yielded 1,792,059 reads of which 328,112 mapped to the VZV genome, with an average mapped read length of 651 bp, and an average coverage of 1710, and 1,463,947 reads mapped to the host cell genome. The 2.41% of the targeted reads mapping to the VZV genome included the target-specific primers and mapped to the correct position, whereas, 24.112% of the reads included the primer but mapped off-target, while for 73.477% the target-specific primer was missing. These latter two categories could be a result of nonspecific priming or template switching.

Sequencing of the cap-selected samples yielded 6,918,054 reads, of which 614,192 mapped to the viral genome, with an average coverage of 1200, and 6,303,862 reads mapping to the host cell. This latter technique performed poorly in VZV, which was indicated by the short average aligned read length (294 bp) (Fig. 2 panel a). Intriguingly, we obtained the same poor result with other herpesviruses, including PRV [23] HSV-1 (our unpublished results), however, the cap-selection technique performed very well in the baculovirus Autographa californica multiple nucleopolyhedrosis virus (AcMNPV), and vaccinia virus (our unpublished results). These results together with an analysis of the nucleotide distribution at the vicinity of the 5′-ends in the cap-selected datasets (Fig. 2 panel b) demonstrate that the read length of cap-selected RNA-s is not determined by the GC-content (PRV: 73%, HSV-1: 68%, VZV: 46%), but by a yet unknown factor.

Fig. 2

The read length distribution of cap-selected and non-cap-selected RNAs. a. The frequency of the reads of the cap-selected and non-cap-selected sequencing binned by 200 nucleotides shows that the cap-selected reads are skewed towards shorter read lengths. b. The fraction of each nucleotide in the vicinity of the 5′ ends of the cap-selected reads shows random distribution, suggesting that the short read lengths are not caused by the presence of a specific nucleotide on the RNA but other, yet unknown reason

In total, we detected 1124 5′-ends and 255 3′-ends in the non-cap-selected dataset, and 1428 5′-ends and 279 3′-ends in the cap-selected dataset (Additional file 2: Table S1-S4). The number of 5′-ends qualifying as a putative TSS was 10.86% of the total 5′-end positions, while 32.95% of the total 3′-end positions turned out as TES (Table 2). We excluded 49 5′-ends and 16 3′-ends from the non-cap-selected and 21 5′-ends and 16 3′-ends from the cap-selected datasets because they proved to be the products of false priming.

Table 2 The number of read ends following each step of filtering

Earlier results with primer extension, S1 nuclease assay and Illumina sequencing determined the transcript ends of some of the RNA molecules of VZV. Using the ONT sequencing platform, we confirmed 18 previously known TSSs and nine TESs. Additionally, we annotated 124 new putative TSSs and 71 TESs listed in Additional file 2: Table S5. The sequences upstream of the TSSs and TESs were analyzed in silico to detect putative GC-, CCAAT- and TATA-box motifs, and Poly(A) signals (PASs), by motif searching (Table 3).

Table 3 Putative promoter sequences and Poly(A) signals (PASs) of the VZV transcripts

Twenty-two and five transcript isoforms were identified with the non-cap-selected, and cap-selected dataset, respectively with a relative abundance higher than 1% (Additional file 2: Table S5).

In this work, we detected an enrichment of As and Us in an interval of up to − 10 nt upstream of the putative TES, and an abundant GT-rich region downstream in the + 10 nt interval (Fig. 3.). We annotated altogether 240 novel transcripts of which 143 were confirmed by at least three reads, additionally, we mapped the TSS and TES of 37 previously detected transcripts lacking former mapping of RNA end positions (Additional file 2: Table S5-S6).

Fig. 3

The frequency of nucleotides in the vicinity of transcriptional end sites (TES). The ±10 nt vicinity of N = 61 unique TESs was analyzed using the online WebLogo v3.6.0 service. The Weblogo shows an enrichment of A and U bases upstream, while an enrichment of G and U bases downstream of the TES. This pattern is akin with the sequence surroundings of mammalian TESs

Putative mRNAs

Transcripts embedded into larger RNA molecules are easily detected by the LRS techniques. If such a transcript contains an in-frame ORF shorter [truncated (t)ORF] than that of the host ORF, it can be considered as a putative mRNA potentially coding for an N-terminally abridged polypeptide. Earlier in silico approaches have annotated the VZV ORFs [6, 7] one of which (the orf33.5) is a tORF. In this study, we report the identification of 25 embedded viral transcripts containing tORFs. We could identify promoters for 10 of these transcripts, located at a mean distance of 93.85 ± 32.18 nt (mean ± SD) from their TSSs (Additional file 2: Table S7). Three of these transcripts (ORF13.5, ORF54.5–53 and ORF62.3) contain strong Kozak consensus sequences near their AUGs, while 19 have at least one of the two essential nucleotides present on their − 3 or + 4 positions. Eighteen of the possibly protein coding VZV transcripts contain a PAS at a 19.44 ± 8.91 (mean ± SD) distance upstream their TESs (Additional file 2: Table S7).

Determination of the ends of mRNAs with annotated ORFs

In this work, we determined the TSSs and TESs of 25 mRNAs whose ORFs had earlier been described (Fig. 1, Additional file 2: Table S5). The TSS and TES of these monocistronic mRNAs were not characterized before by any other means. Twelve of these transcripts have a promoter sequence 51.57 ± 33.31 (mean ± SD) upstream their TSSs, and 19 have a canonical PAS upstream their TESs at a mean distance of 14.96 ± 11.59 (mean ± SD).

Non-coding transcripts

The ncRNAs are transcripts without the ability to encode proteins. They can overlap the coding sequences of the genes in the same orientation, or in the opposite orientation [antisense (as)RNAs], or they can be located at the intergenic regions. The 5′-truncated (TR) transcripts have their own promoters but share the PASs and TESs with the ‘host’ transcript, while the 3′-truncated (NC) RNAs are controlled by the same promoters as the canonical transcript in which they are embedded but have no in-frame ORFs. Twenty-three of the novel non-coding transcripts are long non-coding (lnc)RNAs exceeding 200 bps in length per definition, while five of them are small non-coding (snc)RNAs with a size below 200 bps.

In this work, we identified 18 embedded ncRNAs, of which eight were 5′-truncated, while ten were 3′-truncated transcripts (Additional file 2: Table S5). We also detected one intergenic and seven antisense ncRNAs, all of them being controlled by their own promoters. In total, 17 of the novel ncRNAs contain a canonical promoter sequence 60.4 ± 35.19 nt (mean ± SD) upstream their TSSs, while 21 have a PAS at a 16.35 ± 8.87 (mean ± SD) distance of their TESs. ORF42–43-AS and ORF35-AS overlap multiple mRNAs. ORF42–43-AS stands in antisense orientation with respect of ORF42/ORF45-SP-1 and ORF42/ORF45-SP-2, while in tail-to head orientation with ORF43 and ORF43–44. ORF35-AS stands in antisense orientation with the polycistronic transcript ORF35–34–33 and in tail-to-head orientation with ORF36 and ORF36-S and ORF36–37 (Fig. 1). The LAT RNA has been described in every member of the alphaherpesvirus subfamily [61, 62], including VZV [9]. We confirmed the existence of four previously detected lytic isoforms of VLT (VLT_ly) of which two are TSS isoforms, one is a splice variant and one is both a TSS isoform and a splice variant (Fig. 4).

Fig. 4

The genomic region of VLT. The transcriptome of the varicella zoster virus was sequenced using long-read cDNA sequencing. Reads of the VLT_ly transcript isoforms were visualized using IGV. Because of their low abundance, their 5′ ends are feathered on the annotation, resembling uncertain TSSs

5′- and 3’-UTR isoforms

The 5’-UTR isoforms (TSS variants) start upstream or downstream of the TSS of the earlier annotated transcripts, and their expression is regulated by their own promoters, while 3’-UTR isoforms (TES variants) contain distinct PAS and their polyadenylation occurs upstream or downstream of the TES of the main transcript. In this report, we detected 18 novel 5’-UTR length variants, 7 being shorter and 11 being longer than the earlier annotated transcript isoforms. We found canonical promoter sequences in 10 of the 5’-UTR isoforms at a distance of 94.93 ± 38.58 (mean ± SD). Additionally, we detected eight 3’-UTR length isoform, all with a canonical PAS 22 ± 7.65 (mean ± SD) upstream their TESs. An intriguing finding is the putative TSS at the genome position + 4 belonging to the rare transcripts ORF0–1-L-C and ORF0–1-2-L-C at the extreme termini of UL region (TR_L), which suggests the existence of a promoter located on the other terminal repeat (TR_S) of the genome. This putative promoter is supposed to be active in the circular genome.

Splice sites and splice isoforms

Reverse transcription can produce false introns between repetitive sequences of the template RNA due to the phenomenon of template switching. In order to exclude these artefacts, we removed sequencing reads with low abundance (≤ 1%) and with a repeat of more than six nucleotides next to one of the splice junctions. From the initial set of 10,064 unique splice acceptor and donor site candidates 24 matched our criteria resulting in a total of 16 splice isoforms being above the 1% intron depth level of acceptance.

We detected two novel splice variants encoded by the ORF42/45 gene. Furthermore, we identified twelve novel splice sites and confirmed the existence of nine previously described spliced transcripts, all with a consensus GT at the splice donor site and AG at the splice acceptor site.

ORF63, homologue to the HSV-1 and the PRV us1 gene, is one of the main regulator of VZV transcription. In this work, we discovered ten novel splice variants of ORF63. Similarly to the HSV-1 US1 mRNA, NTO1v1 and the NTO1v3 harbors one intron in its 5’ UTR, while the NTO1v4 is spliced twice just as the PRV US1 mRNA [50, 53] (Fig. 5). Ten splice junctions of ORF63 splice variants coincide with those of the VLT_ly splice variants [9], thus they were labeled as VLT-ORF63-C (Fig. 6). These were confirmed using targeted sequencing. It is noted that one of the splice donor sites, present in both NTO1v1 and NTO1v2 is GC, which differ from the canonical splice donor sequence.

Fig. 5

Structurally similar regions of three alphaherpesvirus transcripts neighboring OriS. The ORF63 of VZV, similarly to the HSV-1 US1 has a two-exon-containing splice isoform, while akin with PRV’s US1 the three-exon-containing splice isoforms. Both VZV and PRV express non-overlapping non-coding RNAs in the proximity of OriS. The target-specific primers used to confirm the TSSs of VZV’s NTO2, NTO3, NTO4 are shown as green arrows

Fig. 6

Types of near-replication-origin (nro)RNAs of three alphaherpesviruses. The nroRNAs can be classified in four distinct types according to their position to the replication origin and coding capacity. Type 1 nroRNAs are non-coding RNAs that do not overlap the Ori. Type 2 nroRNAs are ncRNAs that overlap the Ori. Type 3 nroRNAs are mRNAs with alternative TESs and long 3’UTRs and overlapping the Ori, while type 4 nroRNAs are mRNAs with alternative TSSs and long 5’UTRs overlapping the Ori

We detected all of the three ORF50 splice isoforms previously described [24], but ORF50C occurs in relatively low abundance, below our acceptance threshold. We detected another splice variant in low abundance in the ORF50 cluster, and named it ORF50D (Table 4).

Table 4 Splice junction sites of the VZV transcriptome

In four isoforms, splicing affects the sizes of the ORFs, all producing hypothetical proteins truncated near the N-terminal. Roviš and colleagues [63] have identified the majority of VZV proteins using monoclonal antibodies. In some cases, lower-size bands appeared in their Western blots, which they explained to be the result of putative post-translational modification, or proteolytic cleavage. We propose splicing as possible option for these nonspecific protein isoforms, listed in Table 5.

Table 5 The proteins and hypothetical proteins produced by spliced transcripts

Near-replication-origin (nro)RNAs – A novel class of transcripts

In our earlier work, we reported [64] that PRV expresses transcripts located near the replication origins (Oris): the CTO family (including CTO-S, CTO-S-AT, CTO-M and CTO-L) at the OriL and the PTOs (PTO and PTO-US1) [18] at the OriS. Expression of nroRNAs has also been described in HSV-1 (OriS-RNA: [65]) and HCMV (OriLyt: [66]; RNA4.9: [67]). The VZV genome contains exclusively OriS, and lacks the replication origin at the UL segment. We identified nine nroRNAs starting in the proximity of VZV OriS. The NTO1v1, NTO1v3 and NTO1v5 are the spliced long TSS variants of ORF63 while the NTO1v2, NTO1v4 and NTO1v6 are the spliced long TSS variant of ORF64 (Fig. 5 panel a). NTO2 starts at the same TSS as NTO1v1 but terminates 30 bp downstream of its splice donor site. The 5′ end of the reads of NTO3 and NTO4 are positioned downstream of NTO2. These transcripts were present in our cap-selected but not in the non-cap-selected sample presumably because of their low abundance. In order to augment the detection of these transcripts, we used target-specific primers (Fig. 5 panel a.) for sequencing library preparation, followed by MinION cDNA sequencing. The starting positions detected in the cap-selected sample were not confirmed by the targeted sequencing of these transcripts, thus we hypothesize that NTO3 and NTO4 starts at the closest upstream confirmed TSS, that of NTO1v1. The arrangement of the NTO transcripts is similar to that observed in PRV, where the TSS of the PTO-US1 (positionally similar to the NTO1v1) is located at the same genetic locus as PTO (locationally similar to the NTO2, NTO3 and NTO4, but not homologous) (Fig. 5 panels a and c). Our in silico analysis revealed that the poly(A) cleavage site of NTO2 is located within a canonical CA element, with a GU-rich region starting 13 nt downstream from the TES. The NTO3 has a CTTAAA poly(A) signal [68] starting 20 nt upstream from its TES, with a cleavage site in a homopolymer C run and with a consensus GU-rich region starting 25 nt downstream from the TES. The NTO4 has an ATATAA poly(A) signal [69] starting 24 nt upstream the TES, the cleavage site being marked by a consensus CA signal, followed by a GU-rich region 33 nt downstream. No enrichment of adenines was observed in the genomic region of the TES-s in any of the three NTO transcripts, which would be the source of strand switching or false priming. Based on these results, we can distinguish four types of nroRNAs: (1) ncRNAs that do not overlap the Ori (such as CTO-S, CTO-S-AT, PTO, as well as NTO2, NTO3 and NTO4); (2) ncRNAs that do overlap the Ori (such as CTO-M), (3) mRNA isoforms with very long alternative TES (such as CTO-L and CTO-L2); and (4) mRNA isoforms with very long TSS variant [such as PTO-US1, US1-L (PRV), OriS-RNA2 (HSV-1), and the now discovered NTO1 isoforms] (Fig. 6).

Polycistronic and complex transcripts

A major issue of SRS, as wells as microarray and quantitative PCR approaches is that they have severe limitations in distinguishing between mono- and polycistronic transcripts. In contrast, LRS sequencing is suitable for making this distinction, and it is particularly superior in the detection of low abundance multi-genic transcripts. In this work, we identified 33 novel multigenic RNAs, including 29 polycistronic and 4 complex transcripts. Complex transcripts are multigenic RNAs that contain one or more genes in opposite orientations. Antisense sequences on the complex transcripts are unable to encode proteins. We also detected ten complex transcripts in low abundance in the region of VLT, seven of which are co-terminal with ORF63 and three with ORF64. These transcripts overlap with several oppositely oriented coding sequences and are spliced in a similar manner as the VLT_ly isoforms (Fig. 7). In silico analysis detected an in-frame ORF incorporating the coding sequence of ORF63 (it has an upstream AUG). This results in VLT_ly-ORF63-C1, VLT_ly-ORF63-C4, VLT_ly-ORF63-C5, VLT_ly-ORF63-C6 and VLT_ly-ORF63–64-C1 encoding hypothetical proteins whose N-terminal is longer with 88 amino acids (aa) than the one encoded by the orf63 gene, while the VLT_ly-ORF63-C2; VLT_ly-ORF63-C3 and the VLT_ly-ORF63–64-C2 encoding hypothetical proteins with 179 aa longer than those coded by orf63 (Fig. 7). The N-terminal overhangs of the putative proteins show no homology with any other known proteins in online databases.

Fig. 7

Splice isoforms of ORF63 and ORF64 with a similar splicing pattern as the VLT_ly. The ORF63 and ORF64 splice isoforms encompass two spliced ORFs initialized by two upstream AUGs. The vlt-orf63–1 is translated to a hypothetical protein product 88 amino acid longer, while the vlt-orf63–2 is translated to a hypothetical protein product 179 amino acid longer than those encoded by orf63

Transcriptional overlaps

Transcripts can overlap each other in parallel (tandem or head-to-tail), convergent (tail-to-tail) or divergent (head-to-head) manner (Fig. 8). RNAs identified and annotated with a certain TSS and TES in this work form a total of 481 overlaps, of which 15 are head-to-head, 450 are head-to-tail and 16 are tail-to-tail (Additional file 2: Table S8). The overlaps can be full or partial. Full overlaps can be formed between the RNA molecules encoded by polycistronic transcription units (Fig. 8 panel a.), between embedded and host mRNA molecules, between mRNAs and 3′ as well as 5′ truncated transcripts, between the TSS and TES isoforms, etc. Partial overlaps can be formed between every transcript type. An overlap is ‘hard’ if two genes can only produce overlapping transcripts, or ‘soft’ if both overlapping and non-overlapping transcripts are formed. The soft overlap can be the result of alternative promoter usage (TSS isoforms) or transcriptional read-through (TES isoforms). An example for the latter case is the ORF17 and ORF18, which produce non-overlapping transcripts, but the TES isoform of ORF18 (ORF18-AT), with its additional 89 bps in the 3’UTR, overlaps with ORF17 in a tail-to-tail manner (Fig. 8 panel c.).

Fig. 8

Types of overlaps between VZV transcripts. Reads were visualized using IGV. Blue arrow-rectangles represent transcript annotations. The overlaps are framed in orange. a. Parallel overlaps of the ORF5–4-3 cluster. b. Divergent overlaps of the ORF8 and ORF9 and ORF9A-9 transcripts. c. Convergent overlaps of the ORF17 and ORF18-AT transcripts

Upstream ORFs

Using in silico methods, we detected 44 potential (u)ORFs on the 5’-UTRs of 81 VZV transcripts (Additional file 2: Table S9). Five of the longer TSS variants contain uORFs, while their shorter isoform does not. We have previously described this phenomenon in HCMV [52] and PRV [23]. The average size of an uORF was 54.9 ± 44.64 nt (mean ± SD, median = 35 nt), while the average distance of the uORF stop codon from the protein coding ORF’s AUG was 174.84 ± 143.58 nt (mean ± SD, median = 110 nt). This space between the two uORF and the canonical ORF is enough for a potential reinitiation event. We identified Kozak consensus sequences in five uORFs (Additional file 2: Table S9).

RNA editing

The sequencing reads of NTO3 transcript show a very high frequency of A to G substitution, which is not present in the overlapping reads of other transcripts. We found that 58% of all substitutions are A- > G (Fig. 9 panel a and b) in the reads of NTO3, which is significantly higher than the 12.98% in the overlapping transcripts in the same region (p < 0.0001, Fisher’s exact test) (Fig. 9 panel c), making 22.07% of all As of the transcript edited. The targeted sequencing dataset revealed a similar A- > G substitution pattern in the NTO3 overlapping reads as the cap-selected data. This suggest a hyper-editing event in NTO3. Using the MEME software suite [70] we could detect a slight enrichment of the GU dinucleotide preceding the editing site resulting in a GUAG motif (the editing site is underlined) (Fig. 9 panel d).

Fig. 9

A to I hyper-editing of NTO3. a. Reads of NTO3 and the overlapping transcripts mapping to the VZV genome visualized with IGV. The orange dots represent G mismatches, indicating editing events. b. Substitution matrix of the NTO3 reads (n = 703 substitutions / 7749 nt, p < 0.0001, Fisher’s exact test). c. The position and frequency of A- > G substitutions on the genomic sequence corresponding to the NTO2 transcript, showing both NTO3 and the overlapping transcripts. Substitutions with high frequencies indicate A to I editing events, while those with low frequency are sequencing errors. d. The motif surrounding the editing sites. The motif was found using the MEME software suite with an E-value of 0.58 and a log likelihood ratio of 46. The edited adenosine is marked with a *

Using the ‘RNAstructure’ software suite [58], we predicted the secondary structure with the lowest free energy of the NTO3 ssRNA (Additional file 3 a and b) and of the NTO3-ORF62–5’ dsRNA hybrid, using both the unedited and edited forms of the asRNA (Additional file 3 c and d). When forming an intramolecular secondary structure, the unedited form has a higher free energy state than the edited form (− 143.2 kcal/mol compared to − 169.4 kcal/mol), which suggests that hyper-editing confers thermodynamic stability to the secondary structure of the RNA. Twelve of the 17 Is may aid in the formation of stem structures, while in the unedited RNA only four of the As in the same position have a complementing nucleotide. Contrarily, the secondary structure of the dsRNA formed by the edited NTO3 and the ORF62 is slightly less stable than that formed by the unedited form (− 818.7 kcal/mol compared to − 822.9 kcal/mol), but allows the formation of identical secondary structures.

Discussion

Until now, the VZV transcriptome have been analyzed by Northern blot, primer extension, microarray and Illumina sequencing [17, 71,72,73,74,75,76,77,78]. These techniques have generated useful data, but they have limitations to provide a comprehensive list of the VZV transcripts. In this work, we used the ONT MinION LRS technique for the investigation of the poly(A) + fraction of the VZV transcriptome. Our results increased the number of full length transcripts in VZV by one order of magnitude. We identified 108 novel transcripts, including novel mRNA molecules, monocistronic transcripts, UTR isoforms, as well as multigenic transcripts. Novel splice sites and splice variants were also detected. Additionally, we discovered five sncRNAs and 21 lncRNAs and annotated their TSSs and TESs with base pair precision. The precise annotation of the viral transcriptome is essential for understanding the genetic regulation of VZV. Different TSS isoforms can be regulated by alternative promoters, like for the ul44 gene of HCMV, where two transcript isoforms initiate downstream of a canonical TATA promoter and are expressed at early time points of the viral infection, while the third promoter is non-canonical, and the mRNA isoform starting downstream from it is transcribed at a later phase [79]. Furthermore, alternative TESs may influence RNA metabolism by allowing or forbidding the binding of micro (mi)RNAs to the viral RNA [80] or by facilitating or preventing deadenylation [81].

The enrichment of As and Us upstream of PASs in mammalian systems is well-established [82,83,84]. These homopolymer A stretches may also cause false priming and template switching events during reverse transcription, which results in false TESs. In our work, we excluded these artefacts by the use of our analysis pipeline, and demonstrated the verity of the annotated TESs by the presence of the canonical GU-rich region in the + 10 interval downstream of TESs [85].

Similar to other herpesviruses [22, 53], this study also revealed a complex meshwork of transcriptional read-throughs and overlaps. It has been known earlier that the polycistronic transcription units include co-terminal multigenic RNA molecules, which represent a large extent of parallel overlaps along the entire viral genome. We can raise the question as to whether the downstream sequences of multicistronic transcripts have any other functions besides their possible coding potential, or if they are mere random read-through products representing transcriptional noise. We can address the same question for the convergently and divergently overlapping RNA segments, and also for the alternative transcriptional overlaps. We have put forward a hypothesis that explains the potential role of this phenomenon which is based on transcriptional interaction between the RNA polymerase molecules at the overlapping regions. Since essentially every gene produces overlapping transcripts and therefore the pairwise interactions can spread along the genome thereby forming a transcriptional interference network (TIN), which alongside the promoter-transcription factor system determines the global gene expression pattern of the viral DNA [54].

Understanding the factors influencing the expression of orf63 is of particular interest as its transcript is translated predominantly during latent infections of human ganglia [9, 86]. It has been previously shown that orf63 has no trans-regulatory effect on the expression of orf62, which encodes the major immediate-early transactivator protein of VZV [87]. Nevertheless, deletion of orf63 has demonstrated that it is indispensable for VZV replication, and its product has a protein interaction site with the translated orf62 [86]. Other experiments have shown that deletion of orf63 increases the expression of orf62 [88], suggesting that there is a link between the regulation of the two genes. A possible explanation may be the transcriptional interference caused by the head-to-head overlaps between ORF62 and the isoforms of ORF63 (NTO1v1 and NTO1v2) in the late phase of the viral life cycle.

Polycistronic (especially bicistronic) transcripts are also common in eukaryotic organisms, but their generation is explained by trans-splicing and not by transcriptional read-through [89]. However, unprocessed transcripts are difficult to detect because of their short existence, therefore, it cannot be excluded that these transcripts are the result of a transcription read-through mechanism. The predominant occurrence of adjacent genes in the chimeric transcripts suggests that transcriptional read-through followed by cis-splicing may be the case, that is, the existence of a large extent of transcriptional read-throughs may be not restricted to the herpesviruses but they may represent a general phenomenon. Furthermore, the antisense RNAs produced from their own promoters or by transcriptional read-throughs may hybridize with their sense counterparts thereby initiating RNA interference [90]. Very long complex transcripts form a distinct category among multigenic transcripts because of their oppositely oriented ORFs and increased size. Similarly to other herpesviruses [22, 53] they are present in very low abundance in VZV, however their existence is ambiguous, as they can be formed during reverse transcription by template switching [91].

It has been previously thought that the OriS of VZV lacks any overlapping transcripts [6]. Although Davison and colleagues (1986) hypothesized the existence of non-coding RNA in the intergenic region between orf62 and orf63 [6]. We detected several transcripts overlapping OriS, and a small non-coding RNA (the NCO3) between the before-mentioned two genes. This region of the viral transcriptome is structurally similar to the PRV transcriptome. In the vicinity of OriS, we identified two additional novel non-coding transcripts (NTO2 and NTO4) and two 5′ elongated and spliced version of ORF63. These nroRNAs are supposed to play a role in the regulation of the viral replication [92] through the interplay between the transcriptional and replication apparatuses [54]. Specifically, the role of the transcription of these RNA molecules might be to interfere with the replication machinery, in order to force the replication fork to an unidirectional progression [64]. Furthermore, the intergenic region between orf62 and orf63 contains binding sites for T cell proteins including heat shock factor 1 (HSF-1), nuclear factor 1 (NF-1), NF-κB/Rel, and octamer binding proteins (Oct-1), as well as cyclic AMP-responsive elements (CRE) [93]. The presence of NTO and ORF63 isoforms overlapping this region could interfere with the binding of these host cell proteins. Herpesviruses can evade host cell interferon production and T cell response by producing miRNAs derived from lncRNAs, or antisense RNAs [94,95,96]. Short noncoding RNAs have been recently documented in the orf61–62-63 region of VZV [78]. The NTO and ORF63 isoforms detected in our study could be a possible source of these sncRNAs.

We detected a hyper-editing process in the novel NTO3 transcript. This phenomenon was previously observed in members of the herpesvirus family including alphaherpesviruses [97,98,99]. This process plays a crucial role in the cell’s innate immunity [100], while it can be hijacked by some viruses to evade inactivation [30]. Additionally, in hepatitis delta virus hyper-editing is indispensable for replication [31]. A to I editing decreases the affinity of the antisense transcript to the sense RNA, destabilizing their interaction, which may affect the binding of dsRNA enzymes like RNase III homologues [101]. Our in-silico analysis shows that despite elevating the level of free energy of the sense-antisense RNA hybrid, hyper-editing does not result in a change of the secondary structure, but in the formation of I·U wobble pairs which are significantly more resistant to possible Dicer cleavage inhibiting RNA interference [102]. It is also possible that I·U pairs play a role in the cleavage and degradation of ORF62 by a process mediated by the Tudor staphylococcal nuclease (Tudor-SN) [103, 104]. Further investigation is needed to elucidate the significance of hyper-editing on the NTO3. An evaluation of gene expression at different time points could shed light on the presence of the asRNA and on the amount of editing in different stages of the viral life cycle. Additionally, miRNA assays of the viral transcriptome in different time points of the infection could prove the effect of the RNAi on the dsRNA formed by the sense-antisense pair.

Splice events are thought to be rare in alphaherpesviruses, however, they appear to be underestimated in the light of LRS techniques [105]. In this work, we enriched the list of spliced transcripts of the lytic phase of the viral infection, and we identified the novel combinations of splice sites in ORF42/45. The scarcity of the novel splice isoforms implies further validation, using techniques capable of detecting rare transcripts like targeted direct RNA sequencing.

At least six VZV-transcripts have been shown to be expressed in latency [106, 107], however recent target enrichment SRS data suggests that VLT and ORF63 are the only two expressed transcripts, maintaining the dormant phase of VZV. In this work Depledge and coworkers showed that the VLT is spliced differently in the latent and lytic phase of the viral lifecycle, and identified several length and splice isoforms during productive infection [9]. We confirmed the existence of four introns of this transcript in the longer, lytic forms of VLT. Additionally, we showed a splice isoform of ORF63 and ORF64 possessing the same introns as the VLT. This suggests that occasionally VLT, ORF63 and ORF64 can form a single transcriptional unit. We found that a novel splice isoform of ORF63 produces a long transcript, which has non-canonical splice donor site (GC, instead of GU). Another rare but distinct splice variant of ORF63-L has the same, non-canonical splice donor site. The reason for this could be simply the higher GC content of these regions. Another explanation may be that the splice isoforms of ORF63 are recognized by host cell spliceosomes differentially, in order to perform varying splicing patterns during the life cycle of virus. Similar alternative splicing has been described in VLT, assuming distinct functions for it in lytic and latent phase [9, 105].

The uORFs are supposed to play a regulatory role in the translation of eukaryotic mRNAs. Our results suggests, that at least some of the uORFs could play a role in the production of N-terminally truncated proteins, while most of them have a large-enough space between their stop codons and the protein coding ORF’s AUG for a reinitiation event, thus resulting in unaltered protein translation [108]. Further studies implying ribosome profiling and mRNA and protein expression studies are needed to determine the precise function of these uORFs. Nevertheless, the use of alternative promoters for producing TSS variants with or without uORFs may have a role in providing a differential control of translation at distinct stages of the viral lifecycle [22]. TSSs and TESs detected by former S1 nuclease and primer extension experiments agree with our results with base pair precision, demonstrating that our comprehensive LRS transcriptomic survey of VZV complemented with the proposed bioinformatics pipeline can provide a useful background to design other viral transcriptome reannotation analysis by ONT platform. It is also important to know which transcripts are affected by a genetic modification of the virus. The organization principles of the viral genome and the viral replication can only be understood with a detailed knowledge of the transcriptome.

Conclusions

This study substantially redefines the VZV transcriptome by identifying a large number of novel RNA molecules and transcript isoforms, as well as revealing a complex pattern of transcriptional overlaps. The discovery of novel TSS and TES isoforms can lead to a better understanding of the regulation of the viral gene expression through alternative promoters or through RNA turnover, whereas uORFs discovered in this work can modulate translation. The extensive transcriptional overlaps may indicate an interaction between the transcription machineries. Additionally, the discovery of nroRNAs suggests an interference between the replication and transcription apparatuses. Besides the significant advance in transcriptome annotation, these data, especially the data on novel putative mRNAs and ORFs, may also provide useful information to target transcript candidates for controlling this virus.