Single-Cell Transcriptomics Reveals a Heterogeneous Cellular Response to BK Virus Infection

Abstract

BK virus (BKV) is a human polyomavirus that is generally harmless but can cause devastating disease in immunosuppressed individuals. BKV infection of renal cells is a common problem for kidney transplant patients undergoing immunosuppressive therapy. In cultured primary human renal proximal tubule epithelial (RPTE) cells, BKV undergoes a productive infection. The BKV-encoded large T antigen (LT) induces cell cycle entry, resulting in the upregulation of numerous genes associated with cell proliferation. Consistently, microarray and transcriptome sequencing (RNA-seq) experiments performed on bulk infected cell populations identified several proliferation-related pathways that are upregulated by BKV. These studies revealed few genes that are downregulated. In this study, we analyzed viral and cellular transcripts in single mock- or BKV-infected cells. We found that the levels of viral mRNAs vary widely among infected cells, resulting in different levels of LT and viral capsid protein expression. Cells expressing the highest levels of viral transcripts account for approximately 20% of the culture and have a gene expression pattern that is distinct from that of cells expressing lower levels of viral mRNAs. Surprisingly, cells expressing low levels of viral mRNA do not progress with time to high expression, suggesting that the two cellular responses are determined prior to or shortly following infection. Finally, comparison of cellular gene expression patterns of cells expressing high levels of viral mRNA with those of mock-infected cells or cells expressing low levels of viral mRNA revealed previously unidentified pathways that are downregulated by BKV. Among these are pathways associated with drug metabolism and detoxification, tumor necrosis factor (TNF) signaling, energy metabolism, and translation.IMPORTANCE The outcome of viral infection is determined by the ability of the virus to redirect cellular systems toward progeny production countered by the ability of the cell to block these viral actions. Thus, an infected culture consists of thousands of cells, each fighting its own individual battle. Bulk measurements, such as PCR or RNA-seq, measure the average of these individual responses to infection. Single-cell transcriptomics provides a window to the one-on-one battle between BKV and each cell. Our studies reveal that only a minority of infected cells are overwhelmed by the virus and produce large amounts of BKV mRNAs and proteins, while the infection appears to be restricted in the remaining cells. Correlation of viral transcript levels with cellular gene expression patterns reveals pathways manipulated by BKV that may play a role in limiting infection.

Keywords: BKV; polyomavirus; single-cell transcriptomics.

FIG 1

Viral gene expression varied greatly among BKV-infected RPTE cells at both 2 and 5 dpi. (A) Using results from the RPTE bulk RNA-seq experiment at 2 dpi, cellular genes in the BKV sample were placed into different groups based on their transcripts per kilobase million (TPM) values. (B) The mean and median percentages of cells expressing genes in each TPM group for the 2-dpi SCT were calculated and plotted (M mean, mock mean percent expression; M median, mock median percent expression; B mean, BKV mean percent expression; B median, BKV median percent expression). Note that the results between mock and BKV are similar, and percent expression remains low (<20%) for genes in <100 TPM groups. (C) Variations in viral early region (ER) and late region (LR) gene expression across 2,199 filtered BKV-inoculated cells at 2 dpi. The numbers of BKV-inoculated cells (y axes) were plotted against the percentages of BKV ER (left) or LR (right) gene counts relative to the total UMI counts in histograms. The areas of the histograms in the red boxes were enlarged to better visualize the distribution of cells with 0 to 0.5% ER gene counts (left, indicated with a block arrow) and cells with 0 to 25% LR gene counts (right, indicated with a block arrow). (D) Variations in viral ER and LR gene expression across 2,205 filtered BKV-inoculated cells at 5 dpi. Histograms of the numbers of cells with specific percentages of BKV ER (left) or LR (right) counts relative to the total UMI counts are shown. The enlarged area of the ER histogram shows the distribution of cells with 0 to 2% ER gene counts (left, indicated with a block arrow), and the enlarged area of the LR histogram shows the number of cells with 0 to 25% LR gene counts (right, indicated with a block arrow). (E) Percentage of total BKV-inoculated cells in each group at 2 and 5 dpi. BKV-inoculated cells were divided into low-, medium-, and high-BKV groups based on the proportions of UMIs that mapped to the viral late transcript relative to those that mapped to cellular transcripts (low-BKV cells, ≤1%; medium-BKV cells, >1% and ≤10%; high-BKV cells, >10%).

FIG 2

BKV-infected RPTE cells differed widely in their expression of viral protein products. RPTE cells were inoculated with BKV at an MOI of 1, and the presence of early (LT) and late (VP1) viral products was monitored by immunofluorescence with specific antibodies. (A) Overlapping images of DAPI (4′,6-diamidino-2-phenylindole) (blue) and viral products (magenta) reveal many different levels of intensity in cells from the same culture and viral inoculation, at 2 and 5 dpi. (B) Box plot of staining intensities of individual cells according to the specific viral product and time postinfection. The table indicates the average, median, and standard deviation values for each of the four experimental conditions shown. (C) Distribution of staining intensity according to the number of cells showing similar protein levels.

FIG 3

Distribution of cells into specific stages of the cell cycle according to viral load. RPTE1 cells from three different experimental SCT sets (2 dpi high confluence [HC], 2 dpi, and 5 dpi) were grouped according to their viral load (mock [none] or low, medium, or high expression of late viral transcripts). Within each group, the percentages of cells at different stages of the cell cycle (G1, S, or G₂/M) are depicted.

FIG 4

BKV infection reprogramed cellular gene expression in RPTE cells. (A, left) The principal-component analysis of mock- and BKV-inoculated cells at 2 dpi was visualized using a t-distributed stochastic neighbor embedding (tSNE) plot. Each dot represents a cell (BKV-inoculated cells and mock-inoculated). The tSNE plot for mock- and BKV-inoculated cells shows a mixture of cells from the two populations. (Middle and right) Expression levels of the viral early region (ER) (middle) and late region (LR) (right) are also depicted using tSNE plots. Expression levels of ER and LR transcripts in BKV-inoculated cells are colored based on LogFC values. BKV cells with low levels of LR transcripts largely overlapped mock cells. (B) tSNE plots of mock- and BKV-inoculated cells at 5 dpi were also generated and show that even the low-BKV-expressing cells were completely separated from mock cells. (C) Mock- and BKV-inoculated cells at 5 dpi were ordered into a total of 5 distinct states based on their transcriptional profiles (S1 to S5) using the Monocle pseudotime function. (D) The distributions of mock cells and BKV cells in the low, medium, and high groups were visualized separately by their corresponding pseudotime states (top). The cell numbers from the mock and BKV groups in each state are listed in the summary table.

FIG 5

Cellular genes whose levels covaried with BKV late gene expression levels had functions in diverse biological processes. A total of 80 common covarying genes were identified across 4 independent experiments. The changes in RNA levels for these genes are summarized in a heat map generated using LogFC values of these covarying genes as group markers (low-, medium [Med]-, or high-BKV group versus mock) and/or markers by cell cycle (CC) phase. The markers by cell cycle phase were identified by comparing cells from each of the 4 groups of BKV-inoculated cells against mock cells in matching cell cycle phases. Due to the high similarity between the Low* group cells and cells with no detectable viral UMIs, only the Low* group results are presented for markers by phase. Each row corresponds to a gene. The official gene symbols and corresponding cellular processes are listed on the right. Specific BKV groups with or without cell cycle phase specifications and the time points of infection are indicated on the top. The black box highlights the results of bulk RNA-seq (LogFC of BKV/mock) at 2 dpi.

FIG 6

Confirmation of gene expression trends observed by SCT through analysis of specific protein markers. (A) Simplified tSNE plots indicating the expression of late viral products or the indicated markers in RPTE cells 2 dpi after inoculation with BKV. The expression of CCND1 is reduced in cells showing high levels of viral products, while RRM2 expression directly correlates with the expression of viral products. (B) Immunofluorescence analysis of CCND1 protein expression in comparison to VP1 protein in cells inoculated with BKV at 2 dpi. Combined overlapping images of DAPI, CCND1, and VP1 are shown on the left, while overlaps between CCND1 and VP1 are shown on the right. (C) Immunofluorescence analysis of RRM2 protein expression in comparison to VP1 protein in cells inoculated with BKV at 2 dpi. Combined overlapping images of DAPI, RRM2, and VP1 are shown on the left, while overlaps between RRM2 and VP1 are shown on the right. (D) Quantification of CCND1 and VP1 costaining. (E) Quantification of RRM2 and VP1 costaining. Note that CCND1 was a statistically significant covarying gene in all 4 independent SCT experiments, although its adjusted P value was 0.025 in 1 experiment.

FIG 7

Broad and strong upregulation of ERGs by BKV infection. (A) More extensive and more robust upregulations of ERGs occurred at the later time points. The overlapped box-and-dot plot shows all ERGs with a LogFC of  ≥0.5 in each BKV group (group 1, undetected or 0%, cells from which no viral UMIs were detected; group 2, Low*, cells with viral UMIs that comprised >0% and ≤1% of total cellular UMIs; group 3, Med, cells with viral UMIs comprising >1% and ≤10% of total cellular UMIs; and group 4, High, cells with viral UMIs comprising >10% of total cellular UMIs). (B) Markers by cell cycle phase show greater upregulation of many ERGs in BKV-infected cells than those in the normal cell cycle (columns M_G2M and M_S) at 5 dpi. Due to the high similarity between the Low* group cells and cells with no detectable viral UMIs, only the Low* group results are presented for markers by phase. Note that the markers by phase for the BKV group cells were determined by comparing BKV cells in each group with mock cells in the matching cell cycle. Therefore, the positive LogFC values of BKV markers visualized by the heat map show further upregulations on top of the increased levels of the corresponding genes during the normal cell cycle.

FIG 8

Dynamic changes of cytoplasmic ribosomal protein (RP) genes were associated with different levels of BKV late gene expression. Genes involved in cellular energy metabolism were extensively regulated during BKV infection. (A) A heat map was generated using LogFC values of group markers at 2 and 5 dpi without scaling. The time points and BKV groups are indicated on the top. RP genes (RPL [ribosomal large subunit] and RPS [ribosomal small subunit]) are labeled on the left. The black box in the middle lane highlights the results of bulk RNA-seq (LogFC of BKV/mock) at 2 dpi. (B) The general trend of changes in cytoplasmic RP genes represented by their average LogFCs is summarized in a chart. (C) Heat maps were generated using LogFC values of group markers involved in energy metabolism at 2 and 5 dpi. The time points and BKV groups are indicated on the top. Gene symbols are listed on the right side of the heat maps.

FIG 9

BKV infection resulted in the downregulation of detoxification genes. (A) Expression of detoxification genes in BKV-inoculated RPTE cells according to viral expression levels. A heat map was generated from LogFC values without scaling (red is upregulation, and blue is downregulation). Regulation of detoxification-related transcripts at the bulk RNA level at 2 dpi is included for comparison. (B) Average changes in expression levels (LogFC) of detoxification genes in different cells according to viral load. (C) Simplified tSNE plots indicating cell identity (mock or BKV) and the corresponding expression of late viral products or selected detoxification genes. (D) Relative expression of selected detoxification genes in control and BKV-inoculated cells after 2, 5, and 8 dpi. Values were determined by real-time PCR and normalized against the levels of an endogenous control (glyceraldehyde-3-phosphate dehydrogenase [GAPDH]).

FIG 10

BKV infection orchestrated the regulation of multiple immune genes. The heat map depicts variations in immune genes, generated using LogFC values of group markers at 2 and 5 dpi. The time points and BKV groups are indicated on the top. Gene symbols corresponding to individual rows are listed on the right. All genes encoding acute-phase proteins are enclosed in the left bracket. Genes highlighted in red are the top regulated genes (genes with LogFC values greater than or equal to 0.75 or less than or equal to −0.75 in at least one group). The black box in the middle lane highlights the results of bulk RNA-seq (LogFC of BKV/mock) at 2 dpi.

FIG 11

Downregulation of members of the TNF pathway in BKV-inoculated RPTE cells. (A) Regulation of TNF-related transcript levels according to viral expression in two SCT experimental sets. The heat map was generated from LogFC values. Red indicates upregulation, and blue reflects downregulation. Regulation of TNF-related transcripts at the bulk RNA level at 2 dpi is included for comparison. (B) Average expression levels (LogFC) of TNF-related genes in different cells according to viral load. (C) Expression of different TNF-related factors in supernatants from mock- or BKV-inoculated cells at 2 dpi, as determined by a multiplex enzyme-linked immunosorbent assay (ELISA). For each factor, values are expressed as picograms per milliliter of culture supernatant.