Single-Cell Transcriptomic Analysis of Kaposi Sarcoma

Abstract

Kaposi Sarcoma (KS) is a complex tumor caused by KS-associated herpesvirus 8 (KSHV). Histological analysis reveals a mixture of "spindle cells", vascular-like spaces, extravasated erythrocytes, and immune cells. In order to elucidate the infected and uninfected cell types in KS tumors, we examined twenty-five skin and blood samples from sixteen subjects by single cell RNA sequence analyses. Two populations of KSHV-infected cells were identified, one of which represented a CD34-negative proliferative fraction of endothelial cells, and the second representing CD34-positive cells expressing endothelial genes found in a variety of cell types including high endothelial venules, fenestrated capillaries, and endothelial tip cells. Although both infected clusters contained cells expressing lytic and latent KSHV genes, the CD34+ cells expressed more K5 and less K12. Novel cellular biomarkers were identified in the KSHV infected cells, including the sodium channel SCN9A. The number of KSHV positive cells was found to be less than 10% of total tumor cells in all samples and correlated inversely with tumor-infiltrating immune cells. T-cell receptor clones were expanded in KS tumors and blood, although in differing magnitudes. Changes in cellular composition in KS tumors after treatment with antiretroviral therapy alone, or immunotherapy were noted. These studies demonstrate the feasibility of single cell analyses to identify prognostic and predictive biomarkers.

Fig 1. Landscape of Primary KS.

Single cell suspensions from viably frozen primary KS blood and enzymatically separated skin tumor samples processed for scRNAseq resulted in A) an unsupervised UMAP clustering projection of the merged data set containing 248,741 cells from 23 samples representing a landscape of 25 clusters. Between 15 and 25% of cells from each sample were removed when filtering out dying cells (>20% mitochondrial genes), low quality cells (proportion of UMI > 93^rd percentile), and doublets (identified by DoubletFinder package). Seurat’s default settings were used for normalization and scaling. Principal component (PC) analysis was based on 22 high variable genes with clustering of cells set at a resolution of 0.6. The box around clusters 2, 13, and 19 indicates they are endothelial cells in skin tumors with 95% of KSHV-infected cells mapping to cluster 13. B) UMAP clustering projection color coded based on individual sample distribution across the merged dataset. Samples were submitted in three batches indicated by the color of the sample name in the legend (batch 1, red; batch 2, blue; batch 3, black) and the percent of KSHV+ cells contributed by each sample and the percent of KSHV+ cells present in each sample is shown in S1C Fig. C) UMAP clustering projection color-coded based on tissue source of each sample demonstrated that cells obtained from skin (n=175,472) and peripheral blood mononuclear cell (PBMC) (n=73,269) form clearly delineated clusters. D) UMAP clustering projection color-coded based on reads corresponded to the KSHV genome in which 3,118 KSHV+ cells (with 2 or more reads of viral transcripts) in the merged object corresponding to an average of 1.8% of total tumor cells. 2,970 (95%) of KSHV+ cells were in cluster 13. KSHV+ cells were not detected in PBMC samples. E) UMAP clustering projection color-coded were based on cell annotation using SingleR and the Human Primary Cell Atlas (HPCA) reference dataset. The inset focuses on vascular (red) vs. lymphatic endothelial cells (green) within clusters 2, 13, and 19 identifying cluster 13 as predominantly lymphatic endothelial cells.

Fig 2. Latent and Lytic KSHV in Primary KS.

A “latent” cell is defined as any cell expressing at least two reads from any of the following KSHV latency transcripts: LANA (ORF73, gp81), Kaposin (K12, gp79), vFLIP (ORF71, gp80), K15 (ORF75, gp85), vOX-2 (K14, gp83), or vIRF-2 (gp65) AND expressing no other KSHV gene. A “lytic” cell is defined as any cell expressing at least one read from a latency transcript AND at least one read from any other viral gene. When all samples are combined in a single data set as in Fig 1, KSHV infected cells populate in a single cluster (cluster 13). A) Lytic vs. Latent: Within the UMAP clustering projection described in Fig 1, 2,970 KSHV+ cells are in cluster 13; of which the virus is in latency in 1,257 cells (42%) and lytic replication in 1,713 cells (58%). However, when samples are evaluated individually, unsupervised clustering divides KSHV infected cells into two distinct clusters. B) Two clusters: 10X Cell Ranger, unsupervised, graph-based t-SNE plots of KS6B, KS8, and KS9 skin tumors with KSHV+ (red) clusters (A and B) highlighted. Remaining cluster colors are not coordinated between samples. C) An exploded diagram of the t-SNE plot of sample KS6B in which clusters A (brown) and B (purple) are magnified and color coded based on Log2 KSHV gene expression and lytic vs latent cell identity demonstrate that viral gene expression and lytic vs latent replication are not the defining features that distinguish the cluster A from cluster B. Samples shown in panels B and C were those samples with >2% KSHV-positive cells (S1C Fig).

Fig 3. Two Populations of KSHV Infected Cells in Primary KS.

To characterize infected tumor cells, a new unsupervised analysis of the 2,970 KSHV+ cells within cluster 13 was performed (with number of principle components = 10 and resolution = 0.3) in which cells were annotated using SingleR and Human Primary Cell Atlas (HPCA) dataset resulting in A) a UMAP projection of 10 KSHV+ clusters (numbered 13-0 through 13-9). Although cells in cluster 13 were identified as lymphatic endothelial, small subclusters of KSHV+ cells also express keratinocyte markers (cluster 13-8, n=19 of 135), fibroblast markers (subcluster 13-5, n=23 of 204), vascular endothelial markers (subcluster 13-6, n=40 of 182), or markers of mitosis, including MKI-67, RRM2, and DIAPH3 (subcluster 13-9, n=76). B) The remaining subclusters can be divided into two groups with subclusters 13-1, 13-2, and 13-3 representing cluster B from Fig 2 and subclusters 13-0, 13-4, 13-7 representing cluster A. A UMAP projection color-coded based on sample reveals that subclusters 13-3 and 13-4 are exclusively from sample KS6B and subclusters 13-2 and 13-7 are exclusively from KS8. C) A FeaturePlot of Log2 expression values representing normalized read counts for GAPDH, CD34, EP300 and CREBBP identify cluster A (subclusters 13-0,4,7) as CD34⁺GAPDH^HI and cluster B (subclusters 1,2,3) as CD34^-GAPDH^LO.

Fig 4. Differential Expression Analysis of KSHV

+ cells. A) Differential expression analysis of genes expressed in KSHV+ cells included 1,118 genes with a p value less than 1e-300. When graphed as the average Log2 fold change of gene expression in KSHV+ vs. KSHV negative cells against the ratio of percent of KSHV+ cells expressing that gene (Pct.1) to the percent of KSHV- cells expressing that gene (Pct. 2) the KSHV latency cluster downstream of the LTd promoter (designated KSHV LAT, including viral genes v-FLIP, v-CYC, and K12) along with the lytic transcript K5 were the top viral biomarkers. The voltage-gated sodium channel SCN9A (Nav1.7) emerged as the top, non-KSHV biomarker among a number of previously described factors including FLT4, PROX1, STC1 and CD36.

Fig 5. Immunity in Primary KS.

A) t-SNE plots of primary KS tumor samples in two groups, three samples with greater than (>) 2% KSHV positive cells (KS6B, KS8, KS9) and three samples with less than (<) 2% KSHV positive cells (KS7A, KS10A, KS12) color- coded to highlight Macrophages (Black dots; defined as all cells with elevated Log2 sum expression >2 for CD14, CD80, CD86, TLR2, MARCO, IL10, IL1B, IL1A, and ITGAX (CD11C) also see S9 Fig), CD8+T cells (Blue dots; defined as all cells with elevated Log2 sum expression >2 for CD2, CD3E, CD8A, NKG7, LAG3, GZMH, PTPRCAP, and negative for CD4), and KSHV+ tumor cells (Red dots, defined as all cells expressing KSHV genes). All samples with >2% KSHV-positive cells are shown, and samples shown with <2% KSHV-positive cells are representative of the entire group of samples with that characteristic. B) Graph representing the ratio of Macrophages to KSHV+ cells (black bars), and CD8+T cells to KSHV positive cells (blue bars) in samples are shown in panel A. C) Percent Macrophages, Percent CD8+ T cells and Percent KSHV+ cells in samples from HIV+ vs. HIV- donors with p values shown for comparisons of the results from each group. D) Dot plot of the ratio (Fold) of Macrophages/ Total; CD8+ T cells/ Total; and KSHV+ cells/ Total, for each skin sample. Correlation coefficient for Macrophage vs. KSHV = -0.273 (1 tailed p value = 0.14); CD8+T cell vs. KSHV = -0.356 (one tailed p value = 0.074); Macrophage vs. CD8+ T cell = 0.69 (one tailed p value = 0.00007).

Fig 6. Low CD4:CD8 Ratio in KS Blood.

A) UMAP plots for 4 primary peripheral blood mononuclear cell (PBMC) samples in which multiomic scRNAseq and T-cell receptor (TCR) sequencing was performed (for batch 1 samples including KS6A PBMC, TCRseq was not performed). The combination of TCR reads with gene expression for CD4 and CD8A clearly defined clusters representing the following cell types: TCR+CD4+ (Blue; CD4+ T cells), TCR-CD4+ (Orange; monocytes) also see S9 Fig, TCR+CD8- (Green; CD8+ T cells), TCR-CD8+ (Red; NK), TCR+DP (Purple; double-positive T cells), TCR+DN (Brown; double-negative T cells). The ratio of CD4+ T cells to CD8+ T cells is indicated, as well as the percentage of CD4+ and CD8+ T cells out of total TCR+ cells. B) Quantitation of the percent CD4+ TCR+ and percent CD8+ TCR+ cells in KS peripheral blood mononuclear cell (PBMC) samples (three AIDS-KS and one HIV-negative iatrogenic KS sample are compared to publically available PBMC data from a normal reference data set from ScaleBio [80]. C) Data reported in S4 Fig of Ravishankar et al. [33] in which the CD4:CD8 ratio in AIDS-KS is similar to the samples shown in panel B.

Fig 7. Characterization of T cell diversity in KS.

A) Diversity of T cell clonal expansion in peripheral blood mononuclear cell (PBMC) samples from 4 patients, characterized by proportion of space occupied by clones. Clonal expansion was ascribed to one of 4 relative abundance categories, based on a clone size score calculated by scRepertoire. Clonal size cut points are listed in the panel legend [81]. B) Number of samples expressing the most common complementarity-determining region 3 (CDR3) motifs found in GLIPH2 clustering (Grouping Lymphocyte Interactions by Paratope Hotspots). Among 289 motifs, only one was found among PBMC samples from all 4 subjects. C) The graph depicts the frequency of T-cell receptor β chain joining segment (TRBJ) usage in the 150 most abundant T-cell receptor (TCR) clones in primary KS samples with matched blood and skin samples (KS6B, KS10, KS11, KS12). D) Three sequence logos; top) 9 CDR3 sequences in abundant TRBJ1-1 clones in KS PBMC samples from 4 patients described in panel E; middle) 19 CDR3 sequences from TRB1-1 clones found by GLIPH2 clustering described in panel B; compared to bottom) the anti-open reading frame (ORF) 25 P88 sequence reported by Roshan et al. in [40]. The y-axis BITS indicates the height method used by the sequence logo generation package. E) Alignment of 9 CD3 TRBJ1-1 sequences from abundant T cell clones from 4 KS subjects to the sequence reported in [40]. ‘Max score’ refers to the highest bit score, or sequence similarity score, for the alignment. ‘Perc. Ident’ refers to the percent of total sequence amino acids that are identical to the Roshan et al. sequence [40]. ‘E value’ refers to the significance score. Similar to p values, E values closer to zero represent fewer chances of a result being due to chance, and hence higher significance. The sequence color code refers to the E values, with gray marking an alignment score of <40, blue marking a score of 40-50, and red marking identical residues (score >200) [82]. Samples shown in panels C-E are all samples for which TCR analysis was performed on both skin and PBMCs.