Participant and KS sample characteristics

Thirty Ugandan participants contributed 41 cutaneous KS tumors for RNA-Seq analysis, with 11 participants contributing 2 samples (S1 Table). The majority of participants were men (21/26; 80%) with a median age of 34 years (range, 23–56 years) (Table 1). At the time of KS diagnosis, all had advanced T1 tumor stage per AIDS Clinical Trials Group staging criteria, which includes those with extensive oral KS, visceral KS or tumor-associated edema [27]. All patients were naïve to antiretroviral treatment. Median CD4 T-cell count was 183 cells/mm³ (IQR, 58,331 cells/mm³) and median HIV plasma RNA was 5.5 log10 copies/mL (IQR, 5.1, 5.6 log10 copies/mL). Tumor samples represented a range of morphotypes, including 24 macular (58%), 13 nodular (32%), and 4 fungating (10%) lesions. Total nucleic acids were extracted from the 41 KS biopsies and cDNA libraries were prepared from poly-A-selected RNA and subjected to RNA deep-sequencing on the Illumina platform (S1 Table). 37 independent KS samples were sequenced for 50 bp from paired-end non-stranded libraries. Four additional KS samples were sequenced from stranded libraries to distinguish sense and anti-sense RNA transcripts. Total reads ranged from 81–124 million for the paired-end libraries and 35–52 million for the stranded libraries, which were analyzed at a lower depth of sequencing. RNA reads mapping to the human genome HG19 were subtracted from the libraries and the remaining reads were mapped onto the KSHV NCBI reference sequence NC_009333, strain GK18, with mapped KSHV reads ranging from 13 to 158,924, with a median of 10,232 (Fig 1A). Five of the KS biopsies had greater than 100,000 KSHV-mapped reads. While seven of the KS tumors had less than 1,000 total KSHV-mapped reads, the level of reads mapping to the human genome in these samples was comparable to the other KS samples (S1 Table). A comparison of KSHV mRNA expression in the different tumor morphotypes showed no significant differences in total KSHV-mapped reads (Fig 1B). The KSHV genome copy number per cell was determined for four of the KS samples (001_C, 006_B, 026_B, and 029_B). A comparison of the Ct values obtained from qPCR assays targeting KSHV and the cellular gene oncostatin M, essentially as described in [28], revealed similar KSHV genome copy numbers ranging from 0.5–0.9 KSHV genomes/cell (mean = 0.6). A comparison of the number of mapped KSHV reads/KSHV genome copies showed a variance of 8.6% across the four KS biopsies.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Comparison of reads mapping to the KSHV genome in the RNAseq libraries from different Ugandan KS tumor biopsies.

A) RNA reads from the RNAseq libraries prepared from 41 KS biopsies (S1 Table) were initially mapped to the human genome (hg19) using TopHat2. The set of reads not matching with the human genome were further mapped to the KSHV reference sequence NC_009333 for the KSHV GK18 strain and the reads from the different KS biopsies were compared. B) The RNA read data were grouped according to the tumor morphology classification.

https://doi.org/10.1371/journal.ppat.1007441.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Participant Demographic and Clinical Characteristics.

https://doi.org/10.1371/journal.ppat.1007441.t001

Visualization of KSHV gene expression in KS tumors

The RNA-seq reads mapping to the complete KSHV genome were visualized with the Integrated Genome Viewer (IGV), using a linear scale to provide an overview of the highly expressed regions of the genome. Representative data from KS tumors with 2,432 to 158,924 total KSHV-mapped reads are shown (Fig 2), while data from the seven KS tumors with total KSHV-mapped reads less than 1,000 are provided in S1 Fig The highest level of RNA reads mapped to the region of the T0.7 RNA transcript within the latency locus at the right end of the KSHV genome (Fig 2). T0.7 encodes ORF K12 Kaposin A, a small transmembrane protein implicated in cell adhesion and transformation [29, 30]. High levels of reads mapping to the T0.7/K12 region were consistently detected in all the KS tumor samples, regardless of the level of total KSHV-mapped reads in each sample. High levels of RNA reads also mapped to a small region of the KSHV genome located immediately to the right of ORF72, observable as a single vertical line of mapped reads (Fig 2; indicated with an asterisk in the bottom graphic). A lower but consistent level of reads mapped to the adjacent region containing ORFs 71 and 72. ORFs 71 and 72 encode vFLIP, which functions to promote cell survival and inhibit KSHV lytic replication, and vCYC, which regulates cell-cycle progression, respectively [31]. Downstream of the latency locus at the right end of the genome, a moderate level of reads mapped consistently to the region encoding ORF75, a large tegument protein essential for viral lytic replication [32], and K15, a multiple-pass membrane protein that modulates cellular signaling pathways associated with KSHV-induced angiogenesis [33–35](Fig 2). At the left end of the genome, a moderate level of reads mapped inconsistently to the multifunctional regulatory polyadenylated nuclear (PAN; T1.1) RNA transcript, with higher levels in the tumors with lower total KSHV-mapped reads shown in the upper part of Fig 2. One fungating tumor 023_B showed high levels of both PAN and total KSHV-mapped reads. Some tumor samples contained moderate levels of reads mapping to ORFK2, the viral interleukin-6 (vIL-6) homolog, or ORFK5, the ubiquitin ligase modulator of immune response (MIR2). One tumor, 008_B, showed high level expression of a sharply delineated region of the KSHV genome extending from ORFK3 to ORF19, with very high level of reads mapping to ORFK5 (Fig 2, bracketed). This transcript pattern was unique to this tumor and was not seen in the paired tumor from the same individual (008_C) or any other KS tumor. This sample is further described below.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. KSHV transcriptome analysis in Ugandan KS tumor biopsies—Linear scale.

RNA-seq analysis was performed on 41 KS tumor biopsies and the RNA reads were aligned to the complete reference genomic sequence of KSHV, strain GK18, (NC_009333) and visualized using IGV from the left to the right end of the genome. The vertical axis represents the number of reads aligned to each nucleotide position of the KSHV reference genome (autoscaled) and the visible axis range is shown on the left. Data from 13 representative biopsies from 12 patients are shown in ascending order based on total KSHV mapped read levels. The positions of various KSHV genes across the KSHV genome are indicated at the top and highly expressed genes are indicated at the bottom. An asterisk in the bottom panel indicates the position of a highly expressed region of the KSHV genome flanking the ORF72 coding sequence (P3-exon1). Two different KS biopsies from the same individual (008) are shown. KS biopsy 008_C shows a transcript pattern that is similar to the other KS biopsies, while biopsy 008_B shows high level expression of a region of the left end of the KSHV genome flanking the region encoding PAN (brackets). As discussed in the text, the KSHV strain in KS biopsy 030_C contains the M-allele of K15 and no reads mapped to the P-allele of K15 present in the GK18 reference sequence used for mapping (B: 030_C lane, arrow).

https://doi.org/10.1371/journal.ppat.1007441.g002

To analyze the complete range of read depths across the KSHV genome, the mapped reads were also visualized using a log-based scale (Fig 3). This analysis revealed the presence of RNA reads mapping across the complete KSHV genome with concentrations in regions encoding specific genes associated with lytic replication (ORFs 6, 59, 60), gene regulation (ORFs K8, 57), virion structure (ORFs 11, 17.5, 27, 33, 38, K8.1, 65, 66) and immune modulation (ORFs 4, K2, K3, K4, K5, 45, K10) (Fig 3). Similar patterns of transcription were visually observed in the different KS tumors regardless of the total number of KSHV-mapped reads in each sample (see for example, tumors 11_D and 013_C, with 5,928 and 158,924 total KSHV-mapped reads, respectively). This similarity was also observed in the 7 KS tumors containing less than 1,000 total KSHV-mapped reads (S1 Fig), suggesting that overall the majority of KSHV-infected cells in the KS lesions expressed the same basic transcriptome pattern. Of note, occasional KS samples, such as 030_B, had minimal levels of reads mapping to ORFK15, even though high levels of reads mapped to the adjacent ORF75 (Fig 3, arrow). Different alleles of ORFK15 have large sequence variation extending into the ORF75 sequence [36, 37], which could have affected the ability of ORFK15 reads in the Ugandan tumors to map to the NC_009333 reference sequence.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. KSHV transcriptome analysis in Ugandan KS tumor biopsies—Log scale.

RNA-seq analysis using a linear scale described in Fig 2 was also performed using a log scale to visualize the complete range of mapped reads. See legend to Fig 2.

https://doi.org/10.1371/journal.ppat.1007441.g003

Phylogenetic analysis of ORF75 reveals a wide diversity of KSHV strains in Ugandan KS tumors

The RNA-seq mapping studies were performed by aligning the RNA reads from the different Ugandan KS tumors to the GK18 NCBI Reference sequence (NC_009333) present in a patient with a case of the classic (HIV-negative) KS. The IGV view of the aligned reads from the Ugandan KS tumors revealed numerous mismatches indicating the presence of different KSHV strains in the KS tumor samples. To perform an initial phylogenetic analysis of these strains, we compared the RNA sequences of the large open reading frame encoding ORF75, which was highly expressed in nearly every KS tumor sample. Although most previous studies of KSHV phylogeny have used the widely divergent K1 and K15 sequences [36, 38], these genes, especially K1, were too minimally and inconsistently expressed in the KS tumors to obtain sufficient sequence for phylogenetic analysis. The complete coding sequences of ORF75 (3891 bp) were assembled from the 50 bp reads for 23 of the Ugandan KS samples. These sequences were aligned with published OR75 sequences from 16 unique KSHV genomes from Zambia [39] and ORF75 sequences in the NCBI database from KSHV strains in several KS biopsies and PEL cell lines obtained in Western countries. Maximum likelihood analysis revealed two major (A, C) and four minor ORF75 clusters (B, D-F) (Fig 4). A BLAST alignment with ORF75 partial sequences that have been previously typed [40] indicated that Cluster A corresponded to subtype [B], Clusters B and C corresponded to subtype A/[B], Cluster D corresponded to subtype R/A, Cluster E corresponded with subtype R/M and Cluster F corresponded with subtype N. The vast majority of Zambian sequences were distributed in the major clusters A and C. Previous studies have shown that ORF75 sequences are linked to the adjacent K15 subtypes, of which 3 distinct alleles P, M and N have been detected [36]. The majority of the ORF75 sequences, including those in Clusters A-D were linked to the K15 P-allele (Fig 4), which was present in the GK18 reference sequence. In contrast, the U030 and ZM123 sequences, like the BC1 sequence, were linked to the M-allele, while U012, ZM095 and ZM128 [39], were linked to the K15 N-allele. The sequence differences between the K15 P-allele of GK18 and the more distantly related M and N alleles resulted in essentially no RNA reads from the U012 and U030 tumor samples aligning with the K15 region of GK18 (see Fig 3). Although our phylogenetic analysis was based on a single gene, it is clear that the KSHV strains infecting the different Ugandan KS tumors show high variability, with strong similarity to KSHV strains identified previously in the Zambian KS samples [39]. Overall, the RNAseq analysis indicated the presence of a single KSHV strain in each tumor. In all cases, but subject 008 (discussed below), independent tumors from the same patient contained the same KSHV strain.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Phylogenetic analysis of the KSHV ORF75 sequences in 23 Ugandan KS tumor biopsies.

The complete 3891 bp sequences of the KSHV ORF75 from the KSHV strains in 23 Ugandan KS tumor samples were assembled from overlapping 50 bp RNA reads and compared to the published ORF75 sequences of KSHV strains from 14 Zambian KS biopsies: ZM004 (KT271453), ZM027 (KT271454), ZM091 (KT271455), ZM095 (KT271456), ZM102 (KT271457), ZM106 (KT271458), ZM108 (KT271459), ZM114 (KT271460), ZM116 (KT271461), ZM117 (KT271462), ZM118 (KT271463), ZM121 (KT271464), ZM123 (KT271465), ZM124 (KT271466), ZM128 (KT271467), and ZM130 (KT271468). In addition, ORF75 sequences from the KSHV reference sequence strain GK18 from an HIV-negative individual with classical KS (NC_009333), the strain JSC1 from an HIV infected male with a lymphomatous peritoneal effusion (GQ994935), the strain DG1 from HIV infected Caucasian male with KS (JQ619843) and the strain BC-1 from a primary effusion lymphoma (U75698) are compared. The sequences were aligned using ClustalX and the phylogenetic relationship was determined by maximum likelihood. Separate clusters of sequences (A-F) were identified. Sequences associated with the M-, N-, and P-alleles of the adjacent K15 subtypes are shown. KS tumors 008_B and 008_C were obtained from the same patient. The scale of divergence is indicated.

https://doi.org/10.1371/journal.ppat.1007441.g004

Transcript mapping in the highly expressed KSHV latency locus

To accurately correlate the RNA reads to specific mRNA transcripts in the latency region, we previously developed a map of transcripts that had been identified in the literature [25]. Since these transcripts were initially characterized in different KSHV strains with different sizes, their positions were mapped onto the sequence of the NCBI reference sequence for KSHV, strain GK18 (NC_009333), and the sizes of the GK18 transcripts were predicted from the transcription start and polyadenylation (poly-A) termination sites (S2 Table). A map of the spliced and unspliced latency region transcripts is shown in Fig 5B, with the transcripts grouped according to their use of the poly-A termination sites at bp 117,553 (Group A), bp 122,342 (Group B), and bp 123,015 (Group C). Due to heterogeneity in the DR5, DR6 and DR7 repeat regions, the sizes of the transcripts from the GK18 strain vary from the published transcript sizes determined for other KSHV strains, such as BCBL-1, as indicated in the corresponding references (summarized in S2 Table).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Quantitation of RNA reads mapping to the KSHV latency region.

RNA reads from 34 Ugandan KS tumor biopsies were mapped to the KSHV genomic region from ORF69 to K15 using TopHat2. A) Representative spliced transcript patterns (Sashimi plots) for four KS tumors are shown for genes in the latency locus. The vertical axis represents the number of reads aligned to each nucleotide position of the KSHV reference genome (autoscaled for the height of reads mapping to the 5’exon downstream of promoter P3(LTd), described in the text; the visible axis range is shown on the left). The total number of reads mapping to the KSHV genome is indicated at the right. The number and position of split reads detected by TopHat2 indicating donor and acceptor sites for spliced transcripts are shown and the corresponding excised intron region is labeled. The positions of poly-A sites of transcript termination are shown: ORF69 (bp 117,509); T0.7 and related Group A transcripts (bp 117,553); ORF71 and related Group B transcripts (bp 122,342); ORF72 monocistronic Group C transcript (bp 123,015), as described in the text. B) The proposed transcription pattern for the latency locus is shown (black = sense strand; red = antisense strand). The proposed non-coding antisense ALT RNA and the positions of the direct repeats DR5, 6 and 7, the long inverted indirect repeat (LIR2) and the microRNA region (miRs) are indicated. The positions of the promoters for the latency transcripts are shown: P1 (LTc; bp 128,089), P2 (LTi; bp 127,789), P3 (LTd; bp 124,131), P4 (~bp 119,035), and P5 (~bp 118,234) promoters and the known latency transcripts are graphically represented. Coding sequences are boxed, adjacent 5’ and 3’ regions are thick lines and introns (a-f) are indicated as thin dotted lines (intron a: bp 123,842…119,048; intron b: bp 127,961…124025; intron c: bp 127,961…119,048; intron d: bp 123,842…123,108; intron e: bp 127,961…127,463; intron f: bp 127,961…127,619). The transcripts terminating at the poly-A site downstream of the Kaposin coding sequence (panel A; bp 117553) are grouped together (Group A), while those terminating at the poly-A site downstream of ORF71 (panel A; bp 122342) are grouped together (Group B). The monocistronic ORF72 Group C transcript terminates at the poly-A site bp 123,015. The major T1.7A spliced transcript in Group A encoding the K12 A/B/C complex is indicated with a double asterisk, while the major T1.7B bicistronic transcript in Group B encoding ORF71/72 is shown with a single asterisk. The first exon (P3-Exon1) in the transcripts derived from promoter P3(LTd) is indicated. The sizes of the transcripts are derived from the GK18 strain (NC_009333 reference sequence) (see S2 Table). The DR5, DR6 and DR7 repeat regions contain different numbers of repeats in different KSHV strains, yielding variation in lengths of the Group A and B transcripts [44]. C) The non-overlapping UCDS features for quantitating transcripts are shown. UCDS features for ORFs 69, 71, 72 and K14 correspond to the coding sequences for these ORFs. The K12A UCDS feature corresponds to the T0.7 transcript, while the K12Aa UCDS feature corresponds to the P3-Exon1 downstream of the P3 (LTd) promoter (bp 124,054‥123,843). The miR UCDS feature identifies reads from unspliced mRNAs, such as the 6.4Kb transcript, which map to the genomic region containing the right long inverted repeat (LIR2) encoding the microRNAs K1-9. The miR UCDS feature does not detect processed miRNAs, as these RNAs are not polyadenylated and thus are not present in the sequenced RNA-seq library. The ORF73a and ORF73b UCDS features target the non-repetitive regions of the ORF73 transcript, and reads are combined for quantitation. The DR6 UCDS feature targets an extremely GC-rich repetitive region. D) mRNA reads mapping to the latency region extending from the flanking ORF69 to the terminal K15 were quantitated for 34 KS tumors using the UCDS features and normalized (TPM), and the median and interquartile range was plotted using an R-based BoxPlot algorithm, see Materials and methods. The reads mapping to the K12Aa, DR6, DR5 and K12A UCDS features are boxed in red to indicate their relationship within the T1.7A spliced transcript, while the reads mapping to the ORF71/ORF72 and ORF75/K15 are boxed in blue and purple, respectively, to indicate their relationship within bicistronic transcripts T1.7B (single asterisk, Panel B) and T6.3, respectively.

https://doi.org/10.1371/journal.ppat.1007441.g005

To analyze the read depth and splicing events in the latency region, Sashimi plots were produced by IGV from the TopHat2 analysis of the RNA reads for the KS tumor samples. Sashimi plots graphically present the read depth across a selected region and show the location and quantity of split RNA reads that define the presence of a splicing event within an mRNA transcript. Sashimi plots of four representative KS tumor samples (Fig 5A) were aligned with the different spliced and unspliced transcripts generated from the latency region (ORFs 69-K14) (Fig 5B). This alignment clearly showed a high level of RNA reads mapping to the T0.7 RNA (herein designated as T0.7A) from the P5 latency promoter, which encodes Kaposin A. In addition, a high level of reads mapped to the adjacent DR5 and DR6 repeat regions, even though the TopHat analysis limited reads to a single map position. A number of spliced and unspliced transcripts that terminate at a common poly-A site (bp 117,553) downstream ORF K12 contain the repeat regions. These transcripts are derived from different latency transcript (LT) promoters and include the spliced T1.6A and T1.8A transcripts from the P1/LTc (constitutive) promoter, the spliced and unspliced T1.7A and T6.5A transcripts from the P3/LTd, (downstream) promoter, and the unspliced T1.5A transcript from the P4 promoter (Fig 5B).

Analysis of the Sashimi plots across the latency region revealed a high level of split reads mapping between the highly expressed DR6 repeat region and the highly expressed genomic region downstream of the P3 (LTd) promoter, herein designated as P3-exon1 (Fig 5A and 5B), which was indicated in Fig 2 with an asterisk. Examination of the corresponding sequences in the KSHV genome revealed classical splice donor (AG|gt) and acceptor (ttacgcccccttcgcag|G) sites at bp 123,843 and bp 119,047, respectively. These sites define the presence of a 4,796 bp intron (Fig 5B; labeled “a”), which is spliced from the pre-mRNA for transcripts T1.7A and T1.8A, and encompasses the sequences encoding ORFs 72 and 71, the right origin of replication and flanking long inverted repeat (LIR2) and the microRNAs miR K1-9 and 11. A high level of reads split across intron “a” were detected in the vast majority of KS tumors (Table 2). The depth of the split reads across intron “a” (for ex. see Fig 5A, ranging from 100 reads (030_C) to 1,307 reads (013_C); Table 2) was similar to the depth of the reads in the flanking exons suggesting that the majority of the reads mapping to the K12/DR5/DR6 region were derived from the T1.7A spliced transcript (Fig 5B, double asterisk). This conclusion was confirmed by RT-PCR amplification of the spliced T1.7A transcript from four KS biopsies using primers derived from the K12 region and the exon junction spanning the spliced “a” intron (S2 Fig). While there was evidence for the splicing of intron “b” in some tumors (see for example Fig 5A: 013_C- “b” = 32 reads; Table 2), compatible with the processing of transcript T1.8A from the upstream P1 (LTc) latency promoter, there was no corresponding accumulation of RNA reads mapping to the 5’ exon of this transcript (Figs 3 and 5A). Thus, the splicing data indicates that the high level of split reads across the K12/DR5/DR6 region and P3-Exon1 are derived from the spliced T1.7A transcript initiating from the P3 (LTd) promoter, which is highly expressed in all the KS tumors.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Quantitation of split reads defining spliced latency transcripts.

https://doi.org/10.1371/journal.ppat.1007441.t002

Quantitation of KSHV transcripts

Quantitation of transcripts from complex genomes, such as KSHV, has been difficult due to the compact nature of the genome and the presence of numerous overlapping transcripts that are differentially expressed. We developed a novel approach to quantitate RNAseq reads mapping to specific genes and gene features in the highly complex KSHV genome using unique coding sequence (UCDS) features specific for all known KSHV ORFs and transcriptional regions [25]. Published information and transcript data from RNAseq analysis of primary latent KSHV infections of several endothelial cell types and long-term latent infections in PEL cells were used to identify transcription start and termination sites and globally map mRNA transcripts. The UCDS features were devised to be non-overlapping, separated by the length of a read (50 bp) so that algorithms, such as HTSEQ, using the “intersection_nonempty” setting could specifically identify and distinguish reads mapping to the different targeted gene features on both DNA strands. RNA read depth was determined by quantitating the reads aligning to the UCDS features in a simplified gene feature file (KSHV NC_009333 UCDS ver 020116.gff; S1 File) based on the KSHV reference genome NC_009333 [25]. The read count was normalized to the size of the feature and number of total KSHV-mapped reads in each tumor sample, yielding a relative transcript level, in transcripts per million KSHV reads (TPM), as described in Materials and Methods. Read data and normalization for each KS sample is provided in S3 Table This normalization allowed the expression of each KSHV gene to be compared as a proportion of the total number of KSHV-mapped reads in each sample. In this and subsequent quantitative analyses of read depths across the KSHV genome, the analyses were limited to 34 of the 41 KS tumor samples, which had sufficient levels of RNA reads mapping to the KSHV genome (>1000) for objective comparisons (see Fig 1). The median normalized level of transcripts targeted by our set of UCDS features across the entire KSHV genome for the 34 samples is provided in S4 Table. To more accurately reflect the protein coding potential for genes in regions with overlapping polycistronic transcripts, we determined the level of the primary transcripts for each ORF derived from its associated promoter, in which the ORF would be the first in the transcript to be translated through 5’ CAP-dependent initiation. The contribution of overlapping transcripts derived from distal promoters was subtracted, yielding an estimation of the primary transcript levels for each ORF (S4 Table), as discussed previously Bruce et al [25].

In the latency and flanking regions, UCDS features were identified targeting the coding sequences for the ORFs K12A (T0.7), 72, 71, 73, K14, 74, 75 and K15, as well as additional areas of interest, including the direct repeats DR5 and DR6, the region encoding microRNAs miR-K1-9, 11 (miR-region), and the short P3-Exon1 downstream of the ORF72 P3(LTd) promoter (K12Aa) (Fig 5C). Quantitation of the reads mapping to the latency region from the 34 KS tumor samples revealed very high transcript levels for the K12/T0.7 region (UCDS K12A) and flanking DR5 (UCDS DR5) and DR6 (UCDS DR6) repeat regions (median = 162,467, 310,389 and 65,018 TPM, respectively) (Fig 5D, boxed in red; S4 Table). Although the known transcripts in this region span both the DR5 and DR6 region, transcribed right to left (Fig 5B), significantly less reads mapped to the DR6 UCDS feature, which would detect the 5’ region of such transcripts. The DR6 repeat and upstream flanking regions are known to have significant sequence heterogeneity between different KSHV strains [41], suggesting that the relatively lower level of reads mapping to the DR6 UCDS feature could be due to mismatches between the reads from the Ugandan KSHV strains and the NC_009333 KSHV reference sequence used for the mapping. Very high levels of transcripts (median = 123,915 TPM) were detected using the K12Aa UCDS feature (Fig 5D, boxed in red), which targets the small P3-Exon1 downstream of the P3(LTd) promoter (Fig 5B and 5C). The similarity in the transcript levels detected across the K12A, DR5, DR6 and K12Aa UCDS features (Fig 5D, boxed in red) correlates with the presence of the common spliced transcript (T1.7A) derived from the P3(LTd) promoter, which encodes the Kaposin A/B/C complex (Fig 5A, double asterisk), as described above. This was confirmed by RT-PCR of RNA from 4 KS tumor biopsies (S2 Fig).

Moderate levels of transcripts (median = 15,794 and 14,430 TPM) were detected using UCDS features targeting ORF71 and ORF72, respectively (Fig 5D; S4 Table). Previous studies have determined that the major transcripts encoding ORFs 71 and 72 are either bicistronic (5’ ORF72/ORF71 3’) initiating from the P3(LTd) promoter or tricistronic (5’ ORF73, ORF72, ORF71 3’) initiating from either the P1(LTc) or P2(LTi) promoters (Fig 5B; S2 Table). The similarity in the read depths detected using the ORF71 and ORF72 UCDS features suggests that the majority of these reads are derived from the major unspliced bicistronic transcript (T1.7B) from the P3(LTd) promoter (Fig 5A, single asterisk). While these reads could also be derived from the common spliced bicistronic transcript (T1.8B) from the P1(LTc) promoter (see Fig 5B), only a small number of reads mapped to the 5’ exon of this transcript near the P1(LTc) promoter (see Fig 5A). Only a low level of transcripts (3,300 TPM) were detected with the miR UCDS feature, which targets transcripts containing the miR region flanking the LIR2 (Fig 5C), suggesting that the majority of the transcripts encoding ORF71 terminated at the poly-A site (bp 122,342) downstream of ORF71 (Fig 5B). Quantitation of the reads mapping to the non-repetitive regions of ORF73 using the ORF73A and ORF73B UCDS features revealed low levels of the tricistronic transcripts T5.2B, T5.4B, T5.5B and T5.7B encoding ORFs 73, 72 and 71 (Fig 5D; S4 Table). This supports the conclusion that the majority of the ORF71 and ORF72 reads map to the unspliced bicistronic T1.7B transcript from the P3(LTd) promoter and not from transcripts from the P1(LTc) promoter.

To further analyze transcription across the ORF72/71 locus, we determined the ratio of transcripts detected using the ORF72 and ORF71 UCDS features. The ORF71/ORF72 ratios in the different tumors ranged from 0.7 to 2.1, suggesting differential expression of the two ORFs that was not strictly compatible with the single bicistronic T1.7B transcript (S3 Fig). Seven tumors showed higher levels of transcripts encoding ORF72 (ORF71/72 ratios from .7 to 0.9) (S3A Fig) indicating the presence of a ORF72 monocistronic transcript T1.0C (Fig 5B), identified previously by Sarid et al [42]. Twenty-one tumors showed higher levels of transcripts encoding ORF71 (ORF71/72 ratios from 1.1 to 2.1) (S3A Fig), indicating the presence of a ORF71 monocistronic transcript T0.9B (Fig 5B), identified previously by Grundhoff and Ganem [43].

Conserved nucleotide changes in the Ugandan KSHV strains alter the protein coding potential of the highly expressed T1.7A transcript

Kaposin A is part of a complex translational program that generates multiple novel proteins from the K12 locus. The Kaposin A sequence is downstream of the DR5 and DR6 repeat sequences, which encode variant translation products initiating at alternate CUG initiation codons that have been detected in cultured PEL cells [41, 44]. The CUG initiation codons occur in the sequence region upstream of the DR6 repeat region, which is highly variable in different KSHV strains [41]. Previous studies identified a non-spliced transcript from a promoter, herein referred to as P4, corresponding to T1.5A (see Fig 5B), whose major translation product in the KSHV strains present in BCBL-1 and JSC-1 cells was a CUG-initiated open reading frame encoding a highly repetitive protein sequence, with an “LAH” N-terminal sequence (Kaposin B)[44](Fig 6B). A minor translation product, Kaposin C, initiated from an alternate downstream CUG initiation codon with a “LQY” N-terminal sequence. Kaposin C contained a similar repetitive protein sequence as Kaposin B but was fused to Kaposin A [44] (Fig 6B). Functional studies have shown that the BCBL-1 Kaposin B can modulate mRNA turnover by stabilizing cytokine mRNAs [45]. However, the functions of Kaposins A and C are as yet unclear, as functions attributed to Kaposin A, such as transforming potential [30], have been further attributed to the microRNA miR-K10, which is embedded within the Kaposin A coding sequence [46].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Sequence variants of the latency region transcripts.

The structure of the KSHV genome across the latency region is shown in the direction of transcription, with the positions of the encoded major ORFs, small upstream ORF (uORF), microRNAs (miRNA), and direct repeats (DR) indicated. The positions of the known promoters and the regions of sequence variability in the repeat regions are also shown. The proposed spliced T1.7A transcripts from the P3 (LTd) promoter and non-spliced transcripts T1.5A and T0.7A from the P4 and P5 promoters, respectively, are shown. The AUG initiating codon for the upstream “MGA” encoding uORF sequence and the CUG initiating codons for the two different downstream repeat open reading frames “LAH” and “LQY” are indicated, as is ORFK12 encoding Kaposin A. A) Putative translation products from the GK18 reference strain are shown: Kaposin D containing the repeat sequence initiating with “LQY” that lacks the Kaposin A sequence and Kaposin E containing the repeat sequence initiating with “LAH” that is fused in tandem with Kaposin A. B) The previously described translation products from the BCBL-1 cell line are shown: Kaposin A alone, Kaposin B containing the repeat sequence initiating with “LAH” (major gene product) and Kaposin C containing the repeat sequence initiating with “LQY” fused in tandem with Kaposin A (minor gene product)[44]. The BCBL-1 genome contains a 2 bp deletion between the DR5 and DR6 repeat regions compared to GK18 (red dotted line). C, D) Putative translation products from the Zambian KSHV strains, ZM114 and ZM004, are shown, which contain 1 bp deletions upstream of the “CUG” initiation region and 1 bp deletions between the DR5 repeat and the ORFK12 coding sequence (red dotted lines). D) Putative translation products from the Zambian KSHV strain, ZM004 are shown, which contains additional nucleotide changes eliminating the “LAH” and “LQY” ORFs, as described in the text. The Ugandan KS tumors with the same sequence changes as ZM114 and ZM004 are indicated.

https://doi.org/10.1371/journal.ppat.1007441.g006

Comparison of the K12 locus of the BCBL-1/JSC KSHV strain with the GK18 reference strain revealed a 2 bp deletion in the BCBL-1/JSC-1 sequence between the DR6 and DR5 repeat regions (Fig 6A and 6B). This insertion alters the open reading frames such that the “LAH” initiated open reading frame from the upstream CUG codon in GK18 is now fused with Kaposin A, suggesting that the major translational product of the GK18 K12 locus would be a novel protein, herein termed “Kaposin D” (Fig 6A). This indicates that the minor translation product of the GK18 K12 locus from the downstream CUG initiation codon would be a novel protein “Kaposin E”, with the N-terminal sequence “LQY”. Unlike Kaposin C in the BCBL-1 strain, Kaposin E would lack the downstream Kaposin A fusion. The putative DNA and encoded protein sequences in this region are provided in S4 Fig.

The RNA sequence data in our study revealed several nucleotide differences between the Ugandan KSHV strains and both the GK18 and BCBL-1/JSC-1 strains. A single nucleotide deletion was detected downstream of the P4 promoter before the CUG initiation codons in all the Ugandan KS samples (Fig 6C and S4 Fig). All Ugandan KSHV strains also had an additional nucleotide deletion between the DR5 repeat and the Kaposin A reading frame (“A” bp 118,228, NC_009333). The deletion of this base, which is also detected in all the published Zambian KSHV strains [39], changes the open reading frame through the repeat region that is contiguous with the Kaposin A sequence downstream. In a subset of Ugandan KSHV strains, from KS tumors 001, 005, 006, 015, 018, 020, 023, 026, 028, 029 and 034, this deletion results in neither open reading frame from the two CUG initiation codons creating a fusion with Kaposin A, as exemplified by the Zambian KSHV strain, ZM114 (Fig 6C). Thus, these KSHV strains would encode Kaposin B as the major CUG-initiated translation product and Kaposin E, as the minor-CUG initiated product, with Kaposin A as a downstream AUG-initiated product. The remaining Ugandan KSHV strains, from KS tumors 003, 004, 007, 008, 009, 010, 012, 013, 022, 024, 030, 037, 088, 099, and 101 had additional sequence differences that altered the possible protein products expressed from this locus. In these strains, the major CUG-initiated ORF, Kaposin B, was eliminated by a change in the putative “CUG” initiation codon to “CGC” (Fig 6D and S4 Fig). In addition, the minor CUG-initiated ORF, Kaposin C, was eliminated by a change in the codon flanking the “CUG” initiator from “CGA” to the stop codon “TGA”, thus immediately terminating translation. A similar situation was seen in the Zambian KSHV strain, ZM004 (Fig 6D), and was confirmed by PCR amplification and sequencing.

Small 5’ upstream ORFs (uORFs) are known to regulate expression of downstream ORFs by capturing or slowing the progression of ribosomes scanning the transcript for favorable initiation codons [47]. Previous studies have shown that a small 24 aa uORF is encoded upstream of ORF72 in the initial P3-exon1 of the T1.7B bicistronic ORF72/71 transcript from the P3(LTd) promoter (Fig 7A). This uORF attenuates the expression of the ORF72 10–20 fold in the T1.7B bicistronic transcript [43], and would be expected to attenuate ORF72 expression in the T1.0C monocistronic transcript and ORF71 expression in the T0.9B monocistronic transcript (Fig 7A). This uORF is also encoded in the initial 5’ exon of the highly expressed T1.7A spliced transcript from the same promoter and is positioned upstream of the ORFs encoding the Kaposin A/B/C complex (Fig 7A). Thus, this uORF would be expected to also attenuate translation of the Kaposin B/C complex proteins, especially since these proteins initiate from CUG codons.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Analysis of highly expressed transcripts from the P3(LTd) latency promoter.

A) A diagram of the coding potential of the spliced and unspliced transcripts derived from the P3(LTd) latency promoter is shown. The presence of a regulatory uORF, introns and coding ORFs are shown. The regulatory uORF has been shown to downregulate ORF72 expression in the ORF71/72 bicistronic transcript, as discussed in the text, and is predicted to downregulate KapB/C in T1.7A and vFLIP in T0.9B. The variable region upstream of the Kaposin B/C (KapB/C) coding sequence and the position of RNA editing affecting the coding sequence of the Kaposin A (KapA) and miR-K10a seed sequence are indicated, as are the median expression levels for the different transcripts in the KS tumor samples. B) Alignment of the Kaposin A variants. Published sequences: (UG-107: ACZ92060; ZM123: KT271465; GK18: NC_009333; BAC-36: HQ404500; Retroperitoneal fibromatosis herpesvirus Macaca nemestrina (RFHVMn): AGY30757); Sequences this study: (Uganda A group: 003, 007, 009, 010, 013, 022, 024, 030, 032, 037, 099, 101; Uganda B group: 001, 004, 006, 008, 011, 012, 014, 015, 018, 020, 023, 026, 028, 029, 034, 088; BCBL-1). Additional unique variations from the Uganda group A and B sequences included: Q14-P14 (011), W20-G20 (023), P40-A40 (034) and unknown sequence variation (no mapped reads) after R54 (006). C) The percent of reads containing the edited base at position 118,096 of NC_009333, resulting in the change from miR-K10a to miR-K10b (panel A) and the change from Ser39 to Gly in Kaposin A (panel B) is shown for all 34 KS tumor samples, with median and interquartile range. Edited reads: median = 6.5, IQR = 0–89; Total reads: median = 79, IQR = 5–161.

https://doi.org/10.1371/journal.ppat.1007441.g007

The highly expressed T1.7A transcript undergoes RNA editing and splicing to liberate an intron encoding KSHV microRNAs

The processing of the pre-mRNA for the T1.7A spliced latency transcript generates an RNA intron that is believed to be the source of the major KSHV microRNA species, miRs-K1-9 and 11 [48–50] (Figs 5B and 7A). Therefore, the high levels of the T1.7A spliced transcript could generate elevated levels of these microRNAs in the KS tumors. Since the other KSHV microRNAs, miR-K10 and miR-K12, are located within the retained exon in the T1.7A spliced transcript, such transcripts could also be processed to produce these microRNAs. Processed microRNAs are not detected in our RNAseq protocol, since they do not contain the 3’ poly-A region used for RNA purification. Thus, purification of small RNA species in the KS tumor samples would be needed to confirm the presence of KSHV microRNAs in the KS tumors.

RNA editing of bp 118,096 (NC_009333) (Fig 7A) has been shown to convert the transformation-associated miR-K10a containing an adenine in its seed sequence into the non-transforming miR-K10b containing an inosine [46]. We examined the RNA reads from the KS tumors for the presence of an edited base in the RNA reads. The fraction of ORFK12 transcripts containing an edited base ranged from 0 to 59% across the 34 KS tumors with a median of 13% (Fig 7C), suggesting that the transforming miR-K10a would be the major microRNA form produced in most of the KS tumors. Since miR-K10 is embedded within the mRNA transcripts encoding Kaposin A, RNA editing also affects the encoded Kaposin A protein sequence, changing a serine to glycine [51](Fig 7A and 7B). Our analysis revealed that only a fraction of the ORFK12 transcripts in the majority of KS tumors would encode the altered Kaposin A.

Transcripts encoding ORF73 LANA are minimally expressed in KS tumors

The spliced and non-spliced Group B latency transcripts transcribed from the P1(LTc) promoter and P2(LTi) promoter terminate at the poly-A site downstream of ORF71 (Fig 5B, S2 Table). Although previous in vitro studies indicated that the P1 and P2 transcripts were the major latency-associated transcripts in PEL cells, there was evidence for only minimal expression in the KS tumors (Figs 3 and 5A). Only two tumors had evidence for splicing of intron “e” within the major latency-associated spliced ORF73/72/71 tricistronic transcript (T5.2B), with only 3 reads detected (Table 2). No evidence was detected for splicing of intron “f” (T5.4B), which is similar to intron “e” but from an adjacent splice acceptor site (“f” acceptor: bp 127,626 [6] compared to “e” acceptor: bp 127,462 [8]). Transcripts initiating from the P1(LTc) promoter upstream of ORF73 LANA were quantitated using the ORF73A and ORF73B UCDS features (Fig 5C). These features flank the large DR7 repeat region within the ORF73 coding sequence, which was not used for quantitation due to its repetitive nature. Consistent low levels of LANA reads (median = 2,195 TPM) were detected in the KS tumors (Fig 5D; S4 Table), indicating low levels of the T5.2B, T5.4B, T5.5B, and T5.7B ORF73/72/71 tricistronic spliced and unspliced transcripts derived from the P1(LTc) or P2(LTi) promoter. The low level of transcription from the P1(LTc) promoter correlates with the very low levels of split reads detected for introns “e” and “f” in the Group B transcripts (Fig 5A and Table 2). RT-PCR analysis of RNA from four of the KS tumor samples with the highest levels of KSHV-mapped reads failed to detect any of the tricistronic transcripts from the P1(LTc) or P2(LTi) promoters, confirming the RNAseq data (S2 Fig).

ORF75 and K15 are consistently expressed at moderately high levels in KS tumors

Additional genes adjacent to the latency locus showed consistent and high-level expression in the KS tumors. These included ORF75 and K15 (Figs 2 and 3), which are expressed in the same orientation as the latency genes described above. While the vast majority of the KS samples showed comparable levels of reads mapping to both ORF75 and K15 (Fig 3), tumors from two individuals, 012 and 030, had high levels of reads mapping to ORF75 with no reads mapping to K15 (Fig 3, arrow). The different alleles of K15 have significant sequence variation extending into the ORF75 sequence. The phylogenetic analysis of the ORF75 sequences indicated that the KSHV strains in tumors 012 and 030 contain the K15 alleles N and M, respectively, thus reads from these strains do not map to the P-allele containing GK18 reference sequence. UCDS features were developed targeting the large ORF75 open reading frame (3,891 bp) (ORF75 UCDS) and the largest exon coding the C-terminal domain of K15 (464 bp) (K15a UCDS) from the heavily spliced K15 gene, which appears to be present in all K15 transcripts. Quantitation of the reads mapping to the ORF75 and K15a UCDS features revealed moderately high transcript levels across 32 KS samples (median = 21,490 and 15,866 TPM respectively) (Fig 5D; S4 Table). The KS samples from 012 and 030 were excluded, as their true read counts were not captured in this analysis. A single poly-A transcription termination signal has been identified for these transcripts downstream of ORF75 [52], suggesting that all of the transcripts detected by the K15a UCDS feature would be bicistronic, encoding K15 at the 5’ end and ORF75 at the 3’ end. Thirty of thirty-two tumor samples showed higher levels of ORF75 transcripts compared to K15 (Fig 5D; S4 Table; S3B Fig), indicating the presence of an ORF75 monocistronic transcript initiating upstream of the AUG initiator of ORF75. Quantitative analysis revealed that approximately two-thirds of the transcripts encoding ORF75 would be bicistronic, with K15 as the primary CAP-dependent translation product. The remaining transcripts would be monocistronic, with ORF75 as the primary translation product. It is not known whether ORF75 translation could occur through IRES-mediated initiation in the bicistronic transcripts, as was observed for ORF71 [43].

Genomic translocation into the latency region induces high level expression of lytic genes

One KS tumor, 008_B, exhibited an unusual high level of reads mapping to the genomic region from ORFK3 to ORF19, as indicated above (Figs 2 and 3; bracketed in black) and in more detail (Fig 8A). The RNA reads at the borders of this region abruptly stopped in the middle of ORFK3 at one end and ORF19 at the other end, with no correspondence to possible RNA termination sites. High level expression of this region was not detected in KS tumor 008_C (Figs 2, 3 and 8A), which was isolated from a different location in the same patient, suggesting the presence of a genomic rearrangement in the KSHV strain in the 008_B tumor. By decreasing the alignment threshold in the initial Bowtie2 alignments, we identified reads containing partial sequences aligning to the ORFK3 and ORF19 regions flanking this highly expressed region. Sequence analysis of these reads revealed that a 14,813 bp ORFK3-ORF19 genomic region had been translocated from its original position at the left end of the genome (bp 19,168 to 33,980) to the right end of the genome within the long-inverted repeat (LIR2) (bp 119,504, numbering from the NC_009333 sequence) (Fig 8). A model genomic sequence of NC_009333 with this translocation was created and used to map the 008_B reads (as indicated in Fig 8). This analysis revealed high levels of spliced transcripts derived from the highly expressed latency promoter P3(LTd) upstream containing the 5’ P3-Exon1 (Fig 8D and 8E). Novel transcripts were detected in which the P3-exon1 was spliced to various acceptor sites in downstream genes within the translocated genome fragment. A spliced transcript encoding ORFK5 was expressed at the highest level. The majority of other genes contained in the translocated region were also highly expressed with non-spliced transcripts driven by their original promoters. Since high level expression of these genes was not observed in the paired 008_C tumor, it appears that the translocation of this genomic region to the highly expressed latency region was responsible for the increased gene expression.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. High level expression of a genomic fragment flanking LIR1, which was translocated into the latency region within LIR2 in KS tumor 008_B.

A) Structure of the KSHV genome in KS tumor 008_C showing the position of a genomic fragment (ORFK3 to ORF19; green) that was translocated into the long-inverted repeat (LIR2) within the latency region in the KS tumor 008_B isolated from the same individual. B) Proposed structure of the KSHV genome in KS tumor 008_B with the ORFK3-ORF19 translocation. While the exact junction of the ORFK3-ORF19 translocated fragment was determined from the sequence of RNA reads mapping across the junction in LIR2, the junction of the genomic sequence across the excision point between ORFs K2 and 20 was not determined due to the absence of RNA reads in this region. C) Sashimi plot showing RNA read levels of genes between ORFs K12 and 72 mapping to a model genome of the KSHV reference sequence NC_009333 containing the proposed translocation. (Vertical axis: Range of RNA reads 0–1000). The position and quantity of split reads detected by TopHat2, indicative of transcript splicing, as discussed in Fig 5, is shown. D) Proposed position and orientation of genes translocated into LIR2 in the 008_B KSHV genome is indicated. The 3’ end of ORF-K3 (K3tr) and the 5’ end of ORF19 (19tr) are truncated at the LIR2 insertion point. The positions of the highly expressed P3(LTd) latency promoter and the initial 5’ exon (P3exon1), described in Fig 5, are indicated. E) Proposed transcripts derived from the P3 promoter with novel splicing patterns are shown. The highly expressed T1.7A and T1.7B transcripts detected in all the KS tumors, as indicated in Fig 5, are boxed. The numbers of reads split across each intron are shown in parentheses, and the transcript levels quantitated using specific UCDS features are indicated as transcripts per million (TPM) KSHV reads.

https://doi.org/10.1371/journal.ppat.1007441.g008

Variable gene expression is detected outside of the latency region

At the left end of the KSHV genome, high levels of RNA reads mapped to the PAN (T1.1) mRNA transcript in some KS tumors (Figs 2 and 3; ex. KS tumors 011_D and 023_B). A UCDS feature targeting the PAN RNA was used to quantitate the RNA reads mapping to PAN. While transcript levels reached nearly 300,000 TPM in several KS samples, with a median of 37,413 TPM, significant variation was observed across the KS tumors (Fig 9A and 9D; S4 Table). The PAN transcript overlaps the longer T6.1 transcript encoding ORFK7 (Fig 9B). A UCDS feature was developed to target the ORFK7 transcript upstream of PAN (Fig 9C) to avoid overlap issues in which high levels of ORFK7 expression were observed in many previous studies that were actually due to PAN. While minimal expression of ORFK7 was detected, with a median of 267 TPM, a strong correlation was observed between the expression levels of ORFK7 and PAN in the KS tumors (R = 0.6401, P<0.0001). Using UCDS features, moderate levels of reads mapping to ORFK2 and ORFK5 were consistently detected in the vast majority of the KS tumors (median = 10,138 and 4,373 TPM, respectively) (Fig 9D; bracketed in red). Quantitative analysis of RNA reads mapping to the other KSHV genes in this region revealed 10 to 100-fold lower levels of transcripts with high variability across the 34 KS tumors (Fig 9D and S4 Table).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Analysis of the expression of PAN and flanking genes at the left end of the KSHV genome.

A) The IGV view of reads mapping to spliced and non-spliced transcripts across the region encoding PAN and flanking genes is shown for 5 representative KS tumor samples. Tumors 008_B and 008_C were obtained from the same individual. The position and depth of mapped reads (linear scale) are shown in line (1) with the range indicated at the far right. The position of split reads mapping to splice junctions is indicated in line (2), where the height and thickness of the arc are proportional to the read depth up to 50 reads. The labeled introns are indicated on the transcript map compiled from published literature and predicted transcripts in (B). ORFs on the sense strand are shown in red, while ORFs on the complementary strand are shown in black. C) Position of UCDS features used for read quantitation. D) Total reads mapping to the UCDS features from 33 KS tumors (without 008_B) were quantitated and normalized (TPM). ORFs with high consistent levels of expression are indicated and boxed in red.

https://doi.org/10.1371/journal.ppat.1007441.g009

Hierarchical clustering analysis reveals three distinct groups of KS tumors with common patterns of gene expression

The transcript levels corresponding to each UCDS feature across the KSHV genome (described previously [25]) were compared for 34 of the KS tumor samples using a hierarchical clustering algorithm implemented in CIMminer [53]. This analysis was limited to the KS tumor samples with more than 1,000 total KSHV-mapped reads. The initial analysis showed expression of the genes in their order within the KSHV genome using the “equal width” binning method to map the TPM expression. This method divides the weightrange of data values into equal width intervals and each interval is mapped to one color for display in the clustered image map (Fig 10; left end genome-top, right end genome-bottom). Obvious high levels of transcripts were detected at the right end of the genome in the latency region, mapping to UCDS features for ORFK12 (UCDS-K12A), the adjacent direct repeat regions DR5 (UCDS-DR5) and DR6 (UCDS-DR6), the P3-Exon1 flanking ORF72 (UCDS-K12Aa), and the terminal ORF75 (UCDS-75) and ORFK15 (UCDS-K15a) (Fig 10). Elevated levels of PAN and ORFK2 transcripts were detected in some of the tumors at the left end of the genome. The hierarchical clustering analysis identified three clusters of KS tumors. Cluster I tumors displayed high level expression of transcripts containing the DR6 repeat region and minimal expression of PAN, while Cluster III tumors displayed high level expression of PAN and lower levels of expression of DR6 (Fig 10). Cluster II tumors were intermediate.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Hierarchical clustering of KSHV transcripts—“Equal width” binning.

Hierarchical clustering using CIMminer was applied to KSHV gene expression data obtained by mapping RNA reads to UCDS gene features from 34 KS tumors with total KSHV-mapped reads greater than 1,000. Reads mapping to each UCDS feature were quantitated, normalized (TPM) and the expression values were color-coded for display in the cluster maps using the “equal width” binning method, as described in the text. The unsupervised clustering is shown at the top and the different KS tumor samples are indicated at the bottom of the map. The UCDS features are listed at the right side in the linear order of the KSHV genome from left (top) to right (bottom).

https://doi.org/10.1371/journal.ppat.1007441.g010

A second analysis was performed using the “quantile” binning method which divides the weightrange of expression data values into intervals each with approximately the same number of data points, spreading out the color differences in the image map. The KS tumors in Cluster I displayed low levels of gene expression throughout the KSHV genome outside of the latency region (Fig 11). In contrast, tumors in Cluster III showed moderate levels of lytic gene expression throughout the KSHV genome. Abnormally high expression of the genomic region ORFK3 to ORF19 was clearly observed in the KS tumor 008_B in this cluster. This is the region that was translocated downstream of the P3(LTd) latency promoter in this tumor. An obvious block of expressed genes was observed in the Cluster III tumors between the long-inverted repeat (LIR1) and PAN at the left end of the genome (Fig 11, top). The Cluster II tumors were intermediate with elevated gene expression across the majority of the KSHV genome. Both of the KS tumors from patient 030 showed unusually high expression of the genomic region from ORF39 through the right end of the genome. This patient was the only one to carry a KSHV strain with a K15 M-allele (Fig 4), as shown by the lack of reads mapping to the K15 UCDS.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Hierarchical clustering of KSHV transcripts—“Quantile” binning.

The hierarchical clustering of KSHV gene expression using CIMminer shown in Fig 10 was also performed using the “quantile” binning method, as described in the text. The position of poly-A transcription termination sites for polycistronic ORFS and the direction of transcription is indicated. The pairs of tumor biopsies from the same individual are shown at the bottom (red dot—macular tumor; blue dot—nodular tumor). The 008_B tumor sample with the unusual genomic rearrangement is indicated with a (*). The gradient of latent to lytic gene expression observed in these tumors is shown at the bottom.

https://doi.org/10.1371/journal.ppat.1007441.g011

Quantitative analysis reveals variable gene expression across the KSHV genome with higher levels of reads mapping to ORFs adjacent to poly-A sites in polycistronic loci

The hierarchical clustering analysis with quantile binning showed an obvious pattern of transcription including expression of transcripts detected by UCDS features for ORFs K11, K2, K4, LIR1, K5, PAN, 38, 45, K8, K8.1, 57, 58, 59, and 65, in addition to latency transcripts K12A, DR5, DR6, miR, 71, 72 and K12Aa (Fig 11, labeled with red text). Other genes, such as ORFs 17.5, 27, 33 and 69 were expressed in a subset of tumors. Higher levels of reads mapped to genes that flanked poly-A transcription termination sites (Fig 11, indicated at right). Most of these genes are present in loci with bicistronic or polycistronic transcripts that terminate at the same site. Thus, the high level of reads mapping to the ORFs adjacent to the poly-A site, in most cases, is due to overlapping polycistronic transcripts that contain the sequences encoding the poly-A-flanking ORF, as we have shown previously [25]. Thus, the high levels of reads mapping to these poly-A-flanking ORFs would not correlate directly with their CAP-dependent protein expression.

No association is observed between gene expression levels and KS tumor morphotypes

The cutaneous KS tumor biopsies had been collected to represent different morphologies, sometimes from the same patient, including macular tumors with a flat appearance, nodular tumors with a distinct 3-dimensional structure and fungating tumors showing ulcerations and necrosis [54]. No correlations were observed between the morphological state and the level of total KSHV-mapped reads (Fig 1B), and the hierarchical clustering analysis revealed no obvious correlations between the morphological state and the patterns of gene expression (Fig 11). Cluster I was composed of 7 different KS tumors (5 macular; 1 nodular; 1 fungating) from 4 different individuals with a median read count of 38,605 (range 15,978–134,856). Cluster II was composed of 18 KS tumors (10 macular; 8 nodular) from 15 different individuals with a median read count of 25,538 (range 2,303–158,924). Cluster III was composed of 9 KS tumors (5 macular; 3 nodular; 1 fungating) from 8 different individuals with a median read count of 3,136 (range 2,548–146,773). Surprisingly, we observed a strong correlation between the gene expression patterns of KS tumor samples from the same patient, regardless of morphotype (Fig 11, bottom; red dot = macular, blue dot = nodular). Seven of 8 paired samples from the same individuals showed clustered gene expression profiles. The only paired sample that did not cluster was from individual 008. The 008_B sample in this pair (Fig 11, indicated with an asterisk at the bottom) exhibited the unusual high-level expression of the region between ORFK3 and ORF19, which was due to the translocation not present in the paired 008_C sample. Gene expression in the translocated region was mainly driven by the P3(LTd) promoter in the latency region, rather than by the normal gene specific promoters.

Principal component analysis (PCA) was used to reduce the complexity of the expression data from the 34 KS tumors. Principal components 1 (PC1) and 2 (PC2) captured 88.4% and 6.2% of the variation in the gene expression data, respectively. This analysis confirmed the grouping of gene expression detected by hierarchical clustering, where the tumor samples are indicated in blue text (Cluster I), black text (Cluster II) and orange text (Cluster III) (Fig 12A). As seen in the hierarchical clustering analysis, no similarities were detected between the patterns of gene expression and tumor morphotype in the PCA analysis (Fig 12A: morphotype of each tumor is color-coded). The PCA analysis revealed a gradient of KSHV gene expression from the more latent phenotype in Cluster I (upper right- hand quadrant) to the more lytic phenotype in Cluster III (lower left-hand quadrant). The similarity in the pattern of gene expression between the paired tumor samples from the same individual was confirmed, with the exception of sample 008_B (labeled with *), as discussed above (Fig 12A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. PCA analysis of normalized KSHV gene expression.

A) the normalized KSHV transcript levels from the 34 KS tumors analyzed in Fig 11 were analyzed by principal component analysis (PCA). The individual KS tumor samples are color coded for morphological type. The gene expression clusters identified in Fig 11 are shown (Group I-blue text and outline; Group II-black text and outline; Group III-orange text and outline). Two tumors 10B and 88C fall outside the PCA clusters, as indicated. The paired KS tumors from the same individual are shown with a black dashed line. B) the RNA reads mapping to the PAN and K12A UCDS features targeting transcripts containing PAN and the K12 Kaposin latency region transcripts, respectively, were normalized to the length of the UCDS feature, expressed as reads per kilobase (RPK), and the ratio of PAN to K12A was plotted for the 34 KS tumors with greater than 1,000 total KSHV-mapped reads (circular dots) and 7 KS tumors with KSHV-mapped reads less than 1,000 (triangles). C) the total KSHV-mapped reads detected in the RNAseq libraries were normalized (KSHV-mapped reads/100 million total reads) and compared between the KS tumors in the gene expression clusters identified in Fig 11: (34 KS tumors, reads > 1,000; circular dots) and the 7 KS tumors with less than 1,000 mapped reads (triangles), which were grouped within the lytic Cluster III due to the high PAN/K12A ratio identified in Panel B. The morphological types are color coded: red-fungating, green-macular, blue-nodular. The median and interquartile range are shown and the gradient from latent to lytic gene expression determined from Fig 11 is indicated. Significance was determined using the unpaired T-test. The 023_B KS sample contained a very high level of KSHV-mapped reads, and due to the unusual fungating morphology of this tumor, the data was plotted (in parenthesis) but not included in the statistical calculations.

https://doi.org/10.1371/journal.ppat.1007441.g012

Since PAN expression is considered a marker of lytic activation, we compared PAN expression with the expression of the major T1.7A latency transcript. The RNA reads mapping to the PAN UCDS feature (size = 1.126 Kb) were compared to those mapping to the K12A UCDS feature (size = 0.183 Kb), which detects the Kaposin region in the major T1.7A latency transcript and overlapping transcripts. The reads were normalized to compare transcript levels using reads per kilobase (RPK), plotted as circular dots using the color code for morphotypes (Fig 12B). Because read data for these two features were present in all 41 of the KS tumor samples, we also compared the PAN/K12A ratio of the 7 KS tumor samples with less than 1,000 total KSHV mapped reads that were excluded in earlier expression analysis (plotted as triangles). Although the total number of KSHV-mapped reads was low in these samples, they showed an expression pattern that mirrored those seen in the lytic Cluster III (S1 Fig). The tumor sample 008_B with the altered genome was excluded. The three expression clusters showed statistically significant differences in the PAN/K12A ratio with median ratios of 0.046 (Cluster I; IQR = 0.009:0.084), 0.1545 (Cluster II; IQR = 0.0138:0.297) and 1.635 (Cluster III; IQR = 1.026:2.092) (Fig 12B), indicating that the PAN/K12A ratio could be used to distinguish the different tumor clusters. Notably, three of the four fungating samples showed high PAN/K12A ratios, compatible with a lytic phenotype.

The number of KSHV-mapped RNA reads present in each tumor sample was compared between tumors in the different expression clusters. KSHV-mapped reads were normalized to the total read depth of each tumor library, indicated as KSHV-mapped reads per 100 million library reads. The 7 KS tumors with less than 1,000 total KSHV-mapped reds were included in lytic Cluster III (triangles) due to the similarities in the PAN/K12A RPK ratios with the Cluster III tumors and the presence of lytic gene expression (S1 Fig). Both the latency Cluster I and intermediate Cluster II tumors showed high levels of KSHV-mapped reads in the tumor libraries (Cluster I median = 38,995, IQR = 18,393:139,027; Cluster II median = 25,059, IQR = 5,307:92,593). The lytic Cluster III tumors showed 23 and 10-fold lower levels of KSHV-mapped reads, respectively (Cluster III median = 1,694, IQR = 99:3,490) (Fig 12C). While the extremely high level of KSHV-mapped reads for the fungating tumor 023B (red dot) is shown in parentheses, it was excluded from the column analysis. The presence of high levels of KSHV-mapped reads in the KS tumors with latent phenotypes and the absence of KSHV-mapped reads in the KS tumors with lytic phenotypes is summarized graphically in Fig 12A.

Correlation of gene-gene expression

A set of genes that were highly expressed across the set of 34 KS tumors was identified. Median transcript levels above 5,000 TPM were observed for 11 UCDS features targeting the latency region T1.7A (Kaposin A/B/C) spliced transcript and the ORF72/71 and ORF-K15/75 bicistronic transcripts, as well as the lytic region transcripts for PAN, K2 and K5 (Fig 13A). Median transcript levels of 1–5,000 TPM were observed for 14 other UCDS features targeting ORFs 11, OLAP, K4, 38, 45, 50, K8, K81, 57, 58, 59, 63, miR and 73 (Fig 13A). To identify biologically relevant KSHV gene regulation modules, the Pearson correlation was determined for each pair of highly expressed genes. The analysis was limited to the highly expressed genes as the read level across most of the other genes was low and variable. The 008_B tumor, which showed the unusual genomic translocation, was not included. Hierarchical clustering revealed groups of co-expressed genes that showed similar expression correlation profiles across the 33 KS tumor, indicated by cohesive purple squares along the diagonal (Fig 13B). These clusters of genes represent modules undergoing similar regulation of gene expression in the tumors. The orange areas showed clusters of genes that were negatively correlated. A strong correlated expression of the different transcripts derived from the latency P3(LTd) promoter was observed in Cluster A, including the spliced Kaposin A/B/C T1.7A transcript (detected with the K12A and K12Aa UCDS features), the unspliced bicistronic ORF71/72 T1.7B transcript (detected with the ORF 71 and 72 UCDS features), and the unspliced T6.5A transcript (detected with the miR UCDS feature (see Fig 5 for transcript details). The large Cluster C contains a number of co-regulated genes implicated in the lytic replication cycle, including genes involved in regulation (ORFs 45, 50, K8, 57, and PAN), immune modulation (ORFs K4, K5), replication (ORF59) and virion structure (ORFs 38, K8.1, 65). Other genes in Cluster C with unknown function include ORF58 and OLAP [55]. The expression of these genes is positively correlated with the expression of other genes in Cluster C, but negatively correlated with the expression of the other highly expressed genes. Cluster B contained a number of genes showing variable co-regulation, including K2 (vIL-6), ORF11, K15, ORF73 and ORF75. The expression of these genes negatively correlated with the expression of the P3 (LTd) transcripts in Cluster A and the lytic replication-associated transcripts in Cluster C.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 13. Analysis of the highly expressed genes in the KS tumors.

A) The median expression levels of the most highly expressed KSHV genes were compared, with each dot representing the normalized transcript level (TPM) detected by each KSHV UCDS feature. Genes with transcript levels above 5,000 TPM (blue line) are specifically indicated. Genes with transcript levels between 1,000 (red line) and 5,000 TPM are listed together. Genes with transcript levels below 1,000 are grouped as “all other”. Clusters of UCDS features detecting common transcripts are indicated: Red—major T1.7A transcript in the latency region; Blue—K15/ORF75 bicistronic transcript; Green—ORF71/72 bicistronic transcript. An asterisk marks the genes proximal to poly-A sites. B) hierarchical clustering of gene-gene expression correlations for those genes with median transcript levels above 1,000 TPM. Pearson correlation coefficients (R) were calculated for pairs of genes using normalized transcript levels. The negative correlation (R = -1.0) is shown in orange and the positive correlation (R = 1.0) in purple, with the absence of correlation (R = 0) in white. The unsupervised hierarchical clustering is indicated at the top and the major clusters of genes (A-C) are indicated.

https://doi.org/10.1371/journal.ppat.1007441.g013

Clinical samples

HIV infected adults with KS were recruited from the Uganda Cancer Institute (UCI)/Hutchinson Center Cancer Alliance in Kampala, Uganda. Eligible participants had to be at least 18 years of age and have late stage (T1) KS by AIDS Clinical Trials Group staging criteria and be naïve for antiretroviral therapy (ART). All KS samples were obtained with written informed consent.

RNA-sequencing

Total nucleic acids were extracted using RLT buffer (Qiagen) and RNA was isolated using the RNeasy mini kit with DNAse treatment step. Total RNA integrity was checked using an Agilent 2200 TapeStation (Agilent Technologies, Inc., Santa Clara, CA) and quantified using a Trinean DropSense96 spectrophotometer (Caliper Life Sciences, Hopkinton, MA). Unstranded RNA-seq libraries were prepared from 300 ng of total RNA using the TruSeq RNA Sample Prep Kit v2 (Illumina, Inc., San Diego, CA, USA). Four KS tumors libraries were prepared using the TruSeq Stranded mRNA Library Kit (Illumina) from 100 ng of total RNA. Library size distributions were validated using an Agilent 2200 TapeStation (Agilent Technologies, Santa Clara, CA, USA). Additional library QC, pooling of indexed libraries, and cluster optimization was performed using Life Technologies’ Invitrogen Qubit 2.0 Fluorometer (Life Technologies-Invitrogen, Carlsbad, CA, USA). The unstranded RNA-seq libraries were pooled (5-plex) and the stranded libraries were pooled (4-plex) and each pool was clustered onto a flow cell lane. Sequencing was performed using an Illumina HiSeq 2500 in “High Output” mode with a paired-end, 50 base reads (PE50) sequencing strategy for the unstranded libraries and non-paired end for the stranded libraries. Image analysis and base calling was performed using Illumina’s Real Time Analysis v1.18 software, followed by ‘demultiplexing’ of indexed reads and generation of FASTQ files, using Illumina’s bcl2fastq Conversion Software v1.8.2.

RT-PCR analysis

RNA from four KS tumors was analyzed by RT-PCR to detect different spliced and unspliced transcripts from the latency region of KSHV (S2 Fig). Forward PCR primers from the 3’ end of the ORF73 sequence (P1/2F, _bp124121 5’ CCCTGCCATTAACCCAGCCAG 3’ _bp124101) and the P3-exon1 sequence downstream of the P3 promoter (P3F, _bp123950 5’ ACCCATCTACCTCAACTGAAC 3’ _bp123930) were develop to use in conjunction with a reverse primer from the 3’ end of the ORF72 sequence (72R, _bp123757 5’ CGATCCTCACATAGCGTGGGA 3’ _bp123777) to amplify a 406 bp fragment of the tricistronic ORF73/72/71 (T5.2–5.7B) transcript and a 194 bp fragment of the bicistronic ORF72/71 (T1.7B) transcript, respectively. Reverse PCR primers from the 3’ end of the K12 sequence (KapR1, _bp118037 5’ CGTTGCAACTCGTGTCCTGAATG 3’ _bp118059) and the sequence spanning the exon junction of transcript T1.7A (P4R, _bp119036 5’ TTATAGCGTTTC 3’ _bp119047: _bp123843 5’ CTGTAGAGCCTG 3’ _bp123854), were developed to amplify 1,042 and 120 bp fragments of the T1.7A spliced transcript encoding the K12 Kaposin A/B/C region, respectively, using the P3F forward primer. A forward PCR primer from the sequence downstream of the P4 promoter flanking the DR6 repeat region (P4F, _{bp 118973} 5’ TTGGATTTACACGTATCGAGG 3’ _{bp 118953}) was developed to amplify a 186 bp fragment of K12 using the KapR1 reverse primer. PCR reaction conditions were developed using BCBL-1 DNA as template. cDNA was prepared from KS tumor RNA using Tetro reverse transcriptase (Bioline) and random decamers at 45 °C in the presence of RiboSafe (Bioline), and PCR amplification was performed using MyTaq HS polymerase (Bioline) and forward and reverse primers (5.0 μM), after denaturing at 95° 2.5 minutes with 45 cycles of 95° 15 sec/60° 30 sec/ 72° 60 sec). The PCR products were visualized by electrophoresis using ethidium bromide. Reactions without reverse transcriptase were performed to control for DNA contamination.

Data analysis

RNA reads mapping to the human genome (hg19) were removed using the Bowtie2 program (version 2.2.6) [65]. The remaining RNA reads were aligned to the KSHV reference sequence NC_009333 for the KSHV GK18 strain using TopHat2 (version 2.0.14) [66] in a local instance of Galaxy [67]. Mapped reads were visualized using the Integrated Genome Viewer (IGV; version 2.3.75) [68]. For quantitation purposes, the reads from paired-end libraries (non-stranded) were analyzed as unpaired (single-end data) to allow each read of the pair to map unambiguously to a single gene feature. The reads from both strands of the stranded libraries (non-paired end) were either concatenated and analyzed together (librarytype = unstranded) for comparison to the non-stranded paired end libraries or were analyzed separately (librarytype = FR) to show strand specificity. TopHat2 was used to detect splicing events ab initio. The default presets were used except that the maximum intron length was decreased to 10,000 and the maximum number of alignments allowed was decreased from 20 to 1, to avoid overcounting reads to repetitive regions. HTSeq (version 0.6)[69] was used to quantitate the reads mapping to the unique set of UCDS gene features within the novel revised gene feature file “KSHV NC_009333 UCDS ver 020116.GFF“(S1 File) [25]. The “intersection (non-empty)” setting in HTSeq was used to count all reads mapping completely or partially to a UCDS feature to maximize read count (featuretype = UCDS; IDattribute = gene). No reads were eliminated by ambiguity since the UCDS features were 50 bp apart, the length of a read. The read count was expressed as transcripts per million (TPM) by first normalizing the read count to reads per kilobase (RPK) by dividing the read counts by the length of the UCDS gene feature, in kilobases. The “per million” scaling factor was determined by summing all of the RPK values in a sample and dividing by 1,000,000. Each RPK value was then divided by the “per million” scaling factor to give TPM of mapped KSHV reads. Hierarchical clustering of TPM normalized expression levels was performed using the algorithm implemented in CIMminer [53]. Hierarchical clustering of the gene correlation matrix was performed by calculating the Pearson correlation between the normalized transcript levels (TPM) associated with each pair of UCDS gene features, using a script in R to create and output the correlation matrix. Shiny web applications were developed for R-based principal component analysis in Fig 12A (available at https://efg-ds.shinyapps.io/pcaApp/) and the boxplot analysis of gene expression levels in Fig 5C (available at https://efg-ds.shinyapps.io/boxplotApp/).

Phylogenetic analysis

Phylogenetic analysis of the KSHV strains in the KS tumors was performed on the complete coding sequences of ORF75 (3,891 bp), which were assembled from the 50 bp RNA reads for 23 of the KS tumors with IGV, using maximum likelihood analysis. The sequences have been deposited in Genbank.

Ethics statement

The biological samples in the study were obtained specifically for the study. All participants provided written informed consent. The protocol was approved by the institutional review boards at the Fred Hutchinson Cancer Research Center, the Makerere University School of Medicine, and the Uganda National Council for Science and Technology.