Expression of the Antisense-to-Latency Transcript Long Noncoding RNA in Kaposi's Sarcoma-Associated Herpesvirus

An updated map of known and potential KSHV mRNAs, miRNAs, and lncRNAs.

Arias et al. published KSHV 2.0, which mapped poly(A)-containing viral RNAs and assessed their translational status using ribosomal profiling (13), while Dresang et al. combined a tiling microarray with liquid chromatography and tandem mass spectrometry to explore the protein-coding capacity of the transcriptome (20). A high-throughput transcriptional mapping approach was taken by Majerciak et al. (21), whereas Chandriani et al. focused on tiling microarray supported by Northern blotting or 3′ rapid amplification of cDNA ends (RACE) (19). From these studies at least 16 potential lncRNAs have been reported (13, 19–21), although there are discrepancies between the different approaches (Table 1). Seven potential lncRNAs found by tiling microarray (19, 20) were not detected by subsequent RNA-seq studies (13, 21, 59, 60). We sought to integrate these and prior studies (1, 13, 19–21, 50, 59–61) into an updated, comprehensive map of KSHV transcription annotations (Fig. 1).

At a global level, our analysis of existing data sets found evidence of extensive overlapping and antisense transcription across the entire KSHV genome which was not previously appreciated, as very few of the initially defined open reading frames (ORFs) overlap in KSHV. Nearly the entire genome is transcribed during lytic replication. Several of the KSHV lncRNAs overlap not only mRNAs but also each other, either on the same strand or on opposite strands. For instance, lncRNAs T1.5 and PAN RNA completely overlap T6.1 on the same strand, while PAN RNA is antisense to the anti-PAN lncRNA (Fig. 1). A large number of KSHV lncRNAs are reported to have ill-defined boundaries due to limitations in the methods used for their initial discovery. lncRNAs whose 5′ and/or 3′ ends have not yet been determined were classified as ill-defined in Table 1 and colored gray in Fig. 1, while lncRNAs with well-defined boundaries across multiple independent studies were annotated as such (Table 1; colored purple in Fig. 1). Two lncRNAs, as-vIL6-K4.2 and as-K9-ORF62, were annotated solely based on a verbal description in the results section of a single report (13), while three others, as-ORF65-69, as-vGPCR, and as-ORF75, were mapped by estimating their coordinates based on a published figure (19).

KSHV lncRNA ALT is an early lytic transcript.

ALT is too large and not sufficiently abundant to be reliably assessed by Northern blotting. Conventional reverse transcription-PCR (RT-PCR) methods of detecting RNA are difficult for the ALT lncRNA due to the highly abundant latency transcripts in an overlapping region on the opposite strand (see Fig. 3). To improve sensitivity and specificity, we designed pools of biotinylated antisense oligonucleotides (ASOs) to capture strand-specific RNA by using streptavidin-coupled magnetic beads (62, 63). To assess this enrichment, we used RT-PCR of control RNAs and regions unique to ALT. This strategy is applicable to lncRNAs whose ends are not well known, since capture oligonucleotides can be designed antisense to any region with evidence of transcription.

To validate this approach, we used the well-characterized human lncRNA MALAT1 (metastasis-associated in lung adenocarcinoma transcript 1; also called NEAT2, for nuclear-enriched abundant transcript 2). First, we tested our ability to reverse transcribe full-length MALAT1, as lncRNAs are often difficult templates to extend due to their extensive length and secondary structure. We isolated total RNA from 2 million cells of the PEL cell line TREx BCBL1-RTA (64). Using oligo(dT) to initiate reverse transcription, we verified that three commercial reverse transcriptase (RT) enzymes could generate the full ∼6.6 kb of MALAT1 cDNA from its 3′-poly(A) tail to its 5′ end (as ascertained by primer MALAT#1). Importantly, of the three we tested, only ThermoScript RT also avoided false priming (65), which is critical to this study, as false priming obscures the strand origin for cDNA sequences and can lead to artifacts in stranded RNA-seq approaches. All oligonucleotide primer pairs for RT-PCR experiments amplified the appropriately sized products from genomic DNA (see Fig. 4B).

For strand-specific enrichment of MALAT1, we designed a pool of ASOs complementary to the 3′ end of MALAT1 transcripts starting near the transcription termination sequence (TTS) and spaced roughly every 200 nucleotides (Fig. 2A, dotted black lines, with numbering as shown in Table 2). Total RNA was incubated with the MALAT1-specific pool of biotinylated ASOs, and then ASO-bound RNA was isolated on streptavidin-coupled magnetic beads. RNA-ASO-bead complexes were washed on a magnetic stand to remove nonspecific interactors, and then the RNA of interest was eluted from ASO-bead complexes by high-temperature incubation. To assess the enrichment of MALAT1 transcripts, we performed RT-PCR on input or captured RNA using three MALAT1 amplicons (Fig. 2A, gray shading) and a negative-control amplicon for the highly abundant GAPDH mRNA. All RNA was incubated with or without RT, and no amplification was observed in the non-RT control samples (Fig. 2B). The RT-PCR products were visualized by a Caliper LabChip GX microfluidic instrument, which allows for high-accuracy sizing and quantification. All three MALAT1 regions were amplified from captured RNA at levels similar to input RNA, while GAPDH mRNA was significantly reduced in the captured RNA compared to the input (Fig. 2B). Here and in subsequent experiments, we added random hexamers to the RT reaction to increase consistency across captured lncRNAs of unknown size (Fig. 2B). These results establish that target RNAs, up to an estimated length of 6,600 nucleotides (nt), can be selectively enriched and converted to full-length cDNA using this method.

As the name suggests, ALT lies on the opposite strand of the latency (Lat) transcripts (11, 30, 66), completely overlapping most of them (Fig. 3). In addition, ALT is on the same strand and is coterminal with a bicistronic lytic transcript containing ORF K14 (v-OX2, a homolog of cellular surface receptor OX2) and ORF74 (vGPCR, for viral G protein-coupled receptor) (19, 67). With such extensive overlapping transcription, careful placement of oligonucleotide primers is paramount when trying to detect, capture, and distinguish individual transcript isoforms in this region. To distinguish between the ALT lncRNA and the KSHV latency locus mRNAs, we used primers that span or lie within introns of the latency locus mRNAs and only permit amplification of transcripts from one strand or another (Fig. 3, blue shading). Two such strand-specific primer pairs (ALT specific or LT2, as listed in Table 2) were designed to only generate a single amplified product from one of the two opposing strands. The first, the ALT-specific primer pair (Fig. 3), lies completely within the introns of various latent LANA (latency-associated nuclear antigen) transcripts or upstream of a lytic transcript, labeled LANA LytT, on the bottom strand. It is also upstream of the K14-ORF74 bicistronic transcript on the top strand. Therefore, the ALT-specific primer pair only amplifies ALT cDNA. The second primer pair, the LT2 primer pair, has the potential to produce three amplicons, two from spliced latency transcripts and one from unspliced ALT lncRNA (Fig. 3). PCR amplification of cDNA from the Lat transcript isoform LT2 would result in a 392-bp amplicon, while relatively long amplicons would be produced from either a Lat isoform with a short intron (LT1; 3,812 bp) or the ALT lncRNA (4,311 bp). To ensure that the LT2 primer pair amplifies only the 392-bp product spanning the intron of isoform LT2, we varied the extension times during PCR on PEL cell DNA to determine conditions under which the two long amplicons are not amplified. Thus, we chose an extension time of 45 s, which underrepresents long products and yields a single LT2-specific product from cDNA (Fig. 4A and B).

Additionally, two other primer pairs (LT1 or K14-ORF74, as listed in Table 2 and Fig. 3) were designed to amplify multiple products of different sizes in order to differentiate between spliced or unspliced transcripts from either the same strand or opposite strands. The first, the LT1 primer pair (Fig. 3), spans an intron from Lat isoform LT1 and amplifies a product of 156 bp from that strand but a full-length amplicon of 655 bp from the ALT lncRNA on the opposite strand. The second, the K14-ORF74 primer pair (Fig. 3), spans an intron that is present in a subpopulation of the K14-ORF74 and ALT lncRNA transcripts (19), and it amplifies either a 233-bp product from unspliced RNA or an 85-bp spliced product. Finally, the ORF72 primer pair amplifies a single non-strand-specific product that can originate from either the Lat transcripts or the ALT lncRNA. All of these primer pairs gave the expected product using BCBL1 genomic DNA (Fig. 4B).

KSHV transcription expression is altered during lytic replication. Therefore, to determine the best conditions for maximal expression of ALT, we interrogated latent and lytic conditions. TREx BCBL1-RTA cells carry the potent inducible RTA/ORF50 expression plasmid, and upon addition of doxycycline, the cells synchronously initiate lytic transcription (64). Our laboratory has extensively characterized transcription from the latency locus in this cell line, and therefore we chose the time points of 12 and 24 h postinduction to assess lytic changes in transcription.

To investigate the expression levels of transcripts on both strands of the latency locus, we used oligo(dT) and random hexamers to reverse transcribe mRNA from either latent or lytic cells and then assessed transcripts using the strand- or isoform-specific primer pairs (Fig. 4A). Poly(A) RNA was quantified using the RNA HS assay on a Qubit 3.0 fluorometer, and the same mass of RNA was added to each RT reaction. Consistent with our prior studies (68), two LANA transcripts, LT1 and LT2, were expressed at steady levels during the latent and lytic stages, as seen with the Lat-specific LT2 primer pair (Fig. 4A, lanes 10, 12, and 14) or with the 156-bp LT1-specific amplicon from the LT1 primer pair (Fig. 4A, lanes 18, 20, and 22). Given the massive changes in cellular transcription during reactivation and the reproducibility of LT2 expression in this 24-h period, we used LT2 as a control for transcript expression. The ALT-specific primer pair (Fig. 4A, lane 4) and the LT1 primer pair (i.e., the 655-bp ALT-specific amplicon in Fig. 4A, lane 20) each amplified ALT, and we observed the highest expression of ALT relative to LT2 at 12 h postlytic reactivation. We were able to detect both the spliced and unspliced isoforms from the K14-ORF74 primer pair at all stages of the viral life cycle (Fig. 4A, lanes 26, 28, and 30), which supports the results from the ALT-specific primer pair. No product was observed in the non-RT negative-control samples (Fig. 4A, lanes 3, 5, 7, 11, 13, 15, 19, 21, 23, 27, 29, and 31) or the nontemplate PCR control (NTC) samples (Fig. 4A, lanes 8, 16, 24, and 32). The bands at or below 50 nucleotides represent single-stranded primers and double-strand primer-dimers, which are only visible due to the use of a high-sensitivity microfluidic device. It verifies that primer was added in excess in all reactions.

Using the LT1 primer pair, we noticed an unexpected amplicon of ∼400 bp that is absent from latent cells, appears at 12 h postreactivation, and is expressed at the highest level at 24 h postreactivation (Fig. 4A, lanes 18, 20, and 22). We found the same ∼400-bp amplicon in JSC1 and BC1 cells, confirming this is not a cell line-specific phenomenon. Further, we determined that this amplicon corresponds to a new splice isoform of a LANA-containing transcript. This isoform has a canonical splice acceptor site at genomic position 127888, which is upstream of the known acceptor site in latency transcript LT1 at position 127640. This yields a corresponding intron that is shorter than LT1, as shown in graphical form as LANA new (Fig. 3) and in a multiple-sequence alignment (Fig. 4C) using the online tool MUSCLE, version 3.8 (69). The intron of this new splice isoform lies entirely within the 5′ untranslated region (UTR), as with the intron of LT1, and therefore does not affect the LANA coding DNA sequence (CDS).

To improve the detection of the ALT lncRNA, we next sought to capture ALT prior to reverse transcription. Because ALT was expressed at its highest level at 12 h postreactivation (Fig. 4A, lanes 4 and 20), this time point was used for all subsequent experiments. To ensure the most complete lytic reactivation from doxycycline-inducible TREx BCBL1-RTA cells, we induced with doxycycline and added the histone deacetylase inhibitor valproic acid (VPA). Unlike other inducers, such as n-butyrate, 12-O-tetradecanoylphorbol-13-acetate (TPA), or vorinostat/suberanilohydroxamic acid (SAHA), VPA exhibits no toxicity in PEL even after prolonged incubation. We also added foscarnet/phosphonoformate (PFA) to abolish KSHV DNA replication, nuclear reorganization due to the formation of replication compartments, and premature cell lysis. This increases the chances of detecting authentic, transcriptionally regulated lytic lncRNAs, rather than fragmented, abortive transcription across replicating episomes. To establish strand specificity, we sought to enrich ALT and the Lat transcripts in separate ASO pools. We designed these pools of biotinylated ASOs clustered near the 3′ end of both sets of transcripts, avoiding repetitive regions (Fig. 3, dotted black lines, with numbering as shown in Table 2). For ALT, ASOs were distributed evenly throughout the predicted length of the lncRNA, while ASOs for Lat transcripts were grouped within a common region at the 3′ end to maximize the capture of all Lat transcripts.

To assess the enrichment of ALT and Lat transcripts, we performed RT-PCR on input or captured RNA using the strand- or isoform-specific amplicons within the latency locus (Fig. 3, blue shading) and a negative-control amplicon within the glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene. No amplification was observed in the negative-control samples incubated without RT (Fig. 5). Capture ASOs from both strands were effective, as the ORF72 primer pair, whose amplicon is not strand specific, equally amplified a single product in ALT- and Lat-captured RNA at levels similar to those of input RNA (Fig. 5A). The topography of the KSHV genome is highly repetitive and highly GC rich, leaving few targets for optimal design of ASOs. This, coupled with the low temperature at which ASO binding occurs, means that the expected outcome for this assay is enrichment, not isolation, of our target. Thus, we observe that some GAPDH is amplified in the captured RNA from each ASO pool (Fig. 5B), but at a much smaller amount than the intended target RNAs (e.g., ORF72) (Fig. 5A).

ALT was enriched in ALT-captured RNA compared to Lat-captured RNA (Fig. 5C, E, and F, left versus right side), as seen with either the ALT-specific primer pair (Fig. 5C), the 655-bp ALT-specific amplicon from the LT1 primer pair (Fig. 5E), or both K14-ORF74 amplicons (Fig. 5F). Conversely, the Lat transcripts were only amplified efficiently in Lat-captured but not ALT-captured RNA relative to input RNA (Fig. 5D and E, right versus left side), as seen with the Lat-specific LT2 primer pair (Fig. 5D) and the 156-bp Lat-specific amplicon from the LT1 primer pair (Fig. 5E). These data establish the presence of KSHV ALT by a direct, strand-selective approach independent of potential interpretation biases inherent to short-read sequencing and current heuristic, short-read mapping algorithms.

Identification of ALT 5′ and 3′ ends using RNA-seq.

Having established the specificity of our enrichment method, we sought to define the ends of captured lncRNA by RNA-seq analysis. To this end, we created strand-specific RNA-seq libraries from TREx BCBL1-RTA cells 12 h postreactivation with or without ASO enrichment. We next aligned all sequence reads from ALT-enriched RNA to the KSHV strain JSC1 genome (GenBank accession number GQ994935.1). To identify 3′ ends and transcription termination sequences (TTS), we searched for genomic regions where an abundance of sequence reads contain a poly(A) tract that does not match the reference genome and identified the initial nonmatching adenine. We then determined the percentage of all reads in that genome position that contain adenine or any other nucleotide on either strand. As expected (11, 70), all poly(A)-containing reads in the latency locus data set aligned to known TTS of the Lat transcript, one at JSC1 genomic position 122510 (30% of all reads at that position), which is 3′ to ORF71, and one at position 117247 (82% of reads), which is 3′ to Kaposin/K12 (Fig. 6B).

We identified a TTS downstream of the common ALT and K14-ORF74 polyadenylation (pA) signal. For this TTS shared by the coterminal ALT lncRNA and K14-ORF74 transcripts, 94% of all reads contain ≥5 adenine residues starting at position 130873 (Fig. 6A), while 6% of reads span this position in either direction. These results were obtained from RNA-seq analysis of ALT-enriched RNA, but similar percentages were observed at the same position in libraries from non-ASO-enriched RNA (all raw reads are available in the NCBI Sequence Read Archive under accession code SRP082699). As the two genome-encoded adenine residues upstream of position 130873 interfere with determining the exact start of this poly(A) tract, this result is consistent with the published TTS at positions corresponding to 130870 to 130872 in the JSC1 genome (13, 19, 21, 59, 60).

To examine the TSS of the ALT lncRNA, we mapped a continuous stretch of reads in the same strand orientation as the 5′ region of the ALT lncRNA (Fig. 6B). These reads covered the first ∼1,500 nucleotides of ALT but were not present in the region upstream of the published ALT TSS. Reads were clustered near the two reported 5′ ends of the ALT lncRNA, including the long predominant isoform starting at a position corresponding to 120739 in the JSC1 genome and the short isoform at 121012 (19). We also examined other properties of the ALT lncRNA. To explore the potential splicing of ALT, we used two different tools (71, 72) that predicted >200 potential splice donor and splice acceptor sites within the >10 kb of this lncRNA. Of note, each of these tools confirmed the known splice donor/acceptor pair at JSC1 genomic positions 129545 and 129694, respectively, which represents the only known intron shared by the K14-ORF74 and ALT lncRNA transcripts (19); however, we found no evidence in our RNA-seq data of intra-ALT splice sites.

Although the ends of ALT are clearly defined in our data, we note that there are a few interior areas that have poor read coverage, particularly in repetitive regions. Lack of coverage in the repetitive region within ALT is due to an inherent inability of the BBMerge algorithm to extend reads through repeats. Lack of coverage in other areas is likely due to reduced efficiency of cDNA initiation during reverse transcription or sequence-dependent reduced efficacy of sequencing chemistry used in the Illumina machine. In sum, our experiments confirmed the existence of the ALT lncRNA in KSHV by direct, strand-specific, sequencing-independent means (RT-PCR), uncovered a novel splice site in the 5′ UTR of LANA, and verified the termination and start site of the ALT lncRNA.

Cell culture and reagents.

The KSHV-infected primary effusion lymphoma (PEL) cell line TREx BCBL1-RTA (64) was cultured in RPMI 1640 medium (Gibco), incubated at 37°C in 5% CO₂, and supplemented with 10% tetracycline-free fetal bovine serum (Omega Scientific), l-glutamine (2 mM; Gibco), 100 g/ml streptomycin sulfate, and 100 U/ml penicillin G (Corning). Lytic reactivation of KSHV was induced with doxycycline (1 μg/ml) and valproic acid (1 mM) in the presence of foscarnet (phosphonoformate; 100 μM).

RNA isolation.

RNA was isolated from TREx BCBL1-RTA, JSC1, and BC1 cells. Doxycycline, valproic acid, and foscarnet were added when cultures reached the desired population of ≥10 million cells for total RNA or ≥100 million cells for poly(A) RNA. The culture was either harvested immediately or growth was continued for 12 or 24 h postinduction. Cells were pelleted by centrifugation at 300 × g for 5 min, and supernatants were removed. Cell pellets were resuspended in TRI Reagent (Molecular Research Center) and lysed according to the manufacturer's protocol, except an additional extraction was performed using acidic phenol-chloroform-isoamyl alcohol (125:24:1; Fisher). RNA pellets were dissolved in nuclease-free water and treated with 0.12 U/μl of TURBO DNase using a TURBO DNA-free kit (Life Technologies) according to the manufacturer's protocol. Removal of DNA was assessed by RT-PCR (described below), and RNA quality was determined using an RNA 6000 NanoChip on a Bioanalyzer instrument (Agilent). RNA was quantified using the RNA HS assay on a Qubit 3.0 fluorometer (Invitrogen). Poly(A) RNA was selected using Oligotex mRNA kits (Qiagen) according to the manufacturer's protocol, except that bead-bound RNA was first eluted twice with the same aliquot of nuclease-free water heated to 70°C and then eluted with an additional four aliquots of heated water. RNA in each eluted sample was quantified using Qubit 3.0. Poly(A) RNA was used to assess the expression of ALT and the Lat transcripts by RT-PCR (Fig. 4A) and to generate RNA-seq libraries from one set of biological replicates, while total RNA was used in all other RT-PCR or RNA-seq experiments.

Enrichment of specific RNAs by biotinylated ASOs.

Oligonucleotides antisense to the RNAs of interest were designed using the Stellaris online tool ChIRP (chromatin isolation by RNA purification) Probe Designer, based on work from the Chang laboratory (62, 63). ASOs were designed to bind to nonrepetitive regions of RNA, and biotin was added to the 3′ end along with a triethyleneglycol spacer to allow for more efficient binding to magnetic beads (Eurofins). ASO-mediated enrichment of poly(A) RNA was used to generate a single RNA-seq library (ALT03) (Fig. 6), while ASO-enriched total RNA was used in other RT-PCR experiments (Fig. 2B and 5). For enrichment, ≥2.0 μg poly(A) RNA from ≥70 million cells or ≥20.0 μg total RNA from ≥2 million cells was incubated with 100 pmol of pooled biotinylated ASOs in ChIRP hybridization buffer (750 mM NaCl, 1% SDS, 50 mM Tris-HCl, pH 7.0, 1 mM EDTA, 15% formamide) (63) and RiboLock RNase inhibitor (Thermo Fisher Scientific) at 37°C for 2 h. ASO-bound RNA was then isolated by adding 100 pmol streptavidin-coupled magnetic beads in binding solution from a Dynabeads kilobaseBINDER kit (Invitrogen) and incubating on a rotator at room temperature for 2 h. Beads were washed and RNA was eluted in 10 mM Tris-HCl, pH 7.0, heated at 85°C for 2 min. The concentration of enriched RNA was determined using a Qubit 3.0 fluorometer (Invitrogen) as mentioned in “RNA isolation” above.

Assessment of RNA expression or ASO enrichment by RT-PCR.

Poly(A) RNA was used to assess the expression of ALT and Lat transcripts by RT-PCR (Fig. 4A), and ASO-enriched poly(A) RNA was used to generate the RNA-seq library ALT03 (Fig. 6). ASO-enriched or non-ASO-enriched (input) total RNA was used in RT-PCR experiments (Fig. 2B and 5). Total RNA, poly(A) RNA, or ASO-enriched RNA was reverse transcribed using ThermoScript reverse transcriptase (Invitrogen) with the following slight modifications to the supplier's protocol. The annealing step was performed with 90 to 140 ng of RNA and either with oligo(dT)₂₀ or a mix of oligo(dT)₂₀ and random hexamer primers at 75°C for 5 min. Reverse transcription was performed at 25°C for 10 min and 55°C for 15 min, and then heat inactivation was performed at 85°C for 5 min. The resulting cDNA was amplified using GoTaq Green master mix (Promega) and the pairs of oligonucleotide primers listed in Table 2. PCR primers were first tested on genomic DNA isolated from BCBL1 cells (Fig. 4B). For RT-PCR, cycling conditions included an initial denaturation step of 95°C for 30 s, amplification for 35 cycles (denaturation at 95°C for 20 s, annealing at 56°C for 20 s, and extension at 72°C for 30 s), and a final extension at 72°C for 1 min. Cycling conditions for genomic DNA were the same, except the extension time was 45 s. The reactions were run by capillary electrophoresis on a Caliper LabChip GX microfluidic instrument (PerkinElmer).

Preparation of RNA for high-throughput sequencing.

Strand-specific RNA-seq libraries were made using total RNA, poly(A) RNA, or ASO-enriched RNA (22 to 85 ng) from PEL cells 12 or 24 h postreactivation and the Ovation Universal RNA-seq system (NuGEN) with the following modifications to the supplier's protocol. cDNA was fragmented to an approximate size range of 200 to 500 bp using a TruSeq protocol (i.e., parameters included 5% acoustic factor, 175-W peak incident power, and 200 cycles per burst) with a modified treatment time of 40 s on a Covaris E220 focused ultrasonicator instrument. Human rRNAs were depleted using the provided NuGEN insert-dependent adapter cleavage primers. Libraries were amplified for 20 to 21 cycles. Final cDNA libraries were quantified using the dsDNA HS assay on a Qubit 3.0 fluorometer (Invitrogen), while the size distribution was determined using a high-sensitivity DNA chip on a Bioanalyzer instrument (Agilent). All libraries were paired-end sequenced using an Illumina HiSeq 2500 instrument for 50-bp or 100-bp read lengths.

RNA-seq data analysis.

Sequence reads were trimmed using bbduk.sh (ktrim = r, k = 23, hdist = 1) from the BBMap suite v.36.20 (https://sourceforge.net/projects/bbmap/) and the adapters.fa file included in the BBMap software that contains sequences for commonly used adapters from next-generation sequencing library preparation. Since the 50-bp reads used in this study do not overlap, we used a feature in the BBMerge program (included in the BBMap suite) that makes use of a tadpole (a read assembling program, also included in the BBMap suite). BBMerge tries to extend the 3′ end of each read via assembly using kmer frequencies. BBMerge was run on input read files (R1/R2) with these parameters (k = 31, extend2 = 200) to allow extension of reads up to 250 bp (Brian Bushnell, personal communication, and unpublished BBMerge manuscript). Reads (adapter-trimmed original or merged using BBMerge) were mapped with BBMap (with maxindel = 4000, intronlen = 10, ambig = random, qin = 33, sam = 1.3; the rest of the parameters were default) to the KSHV strain JSC1 genome (GQ994935). BBMap is able to generate alignment results directly as BAM files when SAMtools (v.1.3.1) software is available on the server. We obtained between 3 and 13 million sequence reads that aligned to the JSC1 genome.

We identified the 3′ ends of RNA transcripts in the latency locus of KSHV by using Integrative Genomics Viewer (IGV) v.2.3.79 (85). We first found genomic regions for which the majority of sequence reads contain ≥2 adenine residues that do not match the reference genome and form a poly(A) tract of ≥5 nucleotides, and we visually inspected each genomic position in that region to identify the initial nonmatching adenine. For that particular genomic position and its surrounding positions, we then used bam-readcount v.0.7.4 (https://github.com/genome/bam-readcount) to report the number and percentage of reads that contain each nucleotide base on either strand. We visualized 5′ ends by using Genomics Workbench software v.9.0.1 (CLC bio) to identify collections of sequence reads all in the same strand orientation without any upstream reads.

Accession number(s).

The sequence data have been deposited in the NCBI Sequence Read Archive under accession code SRP082699, and the annotated .gff file is available at https://www.med.unc.edu/orfeome/downloads/annotated-target-genomes.