Single-cell RNA-sequencing of Herpes simplex virus 1-infected cells identifies NRF2 activation as an antiviral program

Herpesvirus infection initiates a range of perturbations in the host cell, which remain poorly understood at the level of individual cells. Here, we quantified the transcrips of single human primary fibroblasts during the first hours of lytic infection with HSV-1. By applying a generalizable analysis scheme, we defined a precise temporal order of early viral gene expression and found unexpected bifurcations and bottlenecks. We identified individual host cell genes and pathways relevant in early infection by combining three different computational approaches: gene and pathway overdispersion analysis, prediction of cell-state transition probabilities as well as future cell states. One transcriptional program, which was turned on in infected cells and correlated with increased resistance to infection, implicated the transcription factor NRF2. Consequently, Bardoxolone methyl, a known NRF2 agonist, impaired virus production, suggesting that NRF2 activation restricts the progression of viral infection. Our study provides novel insights into early stages of HSV-1 infection and serves as a general blueprint for the investigation of heterogenous cell states in virus infection.

activation restricts the progression of viral infection. Our study provides novel insights into early stages of HSV-1 infection and serves as a general blueprint for the investigation of heterogenous cell states in virus infection.

Introduction
Herpes simplex virus 1 (HSV-1) is one of nine known Herpes viruses that affect humans. Although an estimated 80% of the worldwide population is infected in a quiescent, latent form, the virus may cause a variety of diseases during lytic replication and reactivation [1]. A hallmark of HSV-1 is the way it alters the cellular RNA metabolism and RNA content on many levels. On one hand, viral mechanisms activate transcription of viral genes and interfere with splicing regulation, RNA polymerase II transcription, and RNA stability to inhibit synthesis of cellular proteins [2][3][4][5][6][7][8]. On the other hand, virus entry affects a range of cellular pathways, which may in turn lead to transcriptional activation or repression of downstream target genes [9]. Viral infection is a dynamic process driven by the interplay of anti-viral cellular pathways and viral mechanisms which evolved to suppress them. Incoming HSV-1 virions bind to receptors on the cell surface, including NECTIN1/NECTIN2, the TNF receptor superfamily member 14 (TNFRSF14) [10] and the members of the integrin family [11], resulting in activation of cellular pathways such as NF-κB, or of the transcription factors IRF3 and IRF7 [11]. After entry into the host cell, a range of pattern recognition receptors sense viral DNA or RNA, such as the cyclic GMP-AMP synthase MB21D1 (also known as cGAS) [12,13]. Eventually, these pathways can lead to the induction of inflammatory cytokines and interferons [10].
The characterization of cellular heterogeneity due to activation of different host pathways and the progression of viral infection is of great interest. Recent singlecell RNA-sequencing (scRNA-seq) efforts provided an unbiased characterization of virus-host interactions in individual cells which are masked at the population level [14][15][16][17][18][19][20][21][22]. However, deeper insights into unique molecular signatures and discovery of specific cell subsets could be obtained by increasing sequencing depth and the application of advanced analytical approaches to study viral infections.
Here, we profiled deep transcriptomes of tens of thousands of individual cells harvested before and at several times post HSV-1 infection. Our results related the progression of infection to cell cycle phases, and defined a precise temporal order of viral gene expression. The depth of data and using unspliced mRNA as a predictor for future cell states allowed us to connect the course of infection to the activity of specific host cell genes and pathways. Particularly, we investigated the relationship of HSV-1 infection and the transcription factor NRF2, which is activated during infection, and demonstrated that the NRF2 agonist Bardoxolone methyl impairs a productive viral replication. Overall, our study provides novel insights into early stages of HSV-1 infection, and provides an analytical framework to study viral infections using scRNA-seq.

Single-cell RNA-sequencing of HSV-1 infected primary human fibroblasts
To investigate the heterogeneity of molecular phenotypes in the first hours of viral infection, we infected primary human fibroblasts (NHDF) with HSV-1 at an MOI of 10 ( Fig. 1a,b) and profiled the transcriptomes of uninfected cells as well as cells harvested at 1h, 3h and 5h post infection using the droplet-based singlecell sequencing (Drop-seq) [23,24]. For further analysis, only cells with more than 2000 detected genes were used, a threshold that has been previously shown to reduce technical variability [25]. An overview of the dataset (Table S1), number of characterized cells (Table S2), distribution of unique molecular identifiers (nUMI) and number of detected genes (nGene) ( Supplementary Fig.   1a), as well as correlation between scRNA-seq and bulk RNA-seq (Fig. S1B) are provided in Supplemental data. Low-reproducibility genes (Table S3) were subsequently omitted or flagged.
The analyzed cells clustered based on harvesting time point, cell cycle markers, and the amount of viral mRNA, suggesting that the strongest contributors to cellular variability were cell cycle state and the progression of infection (Fig. 1c).
However, cells did not group by biological replicates, indicating that replicates provided comparable and reproducible data (Fig 1d).
The distribution of the viral gene expression per single cell at the different harvesting time points indicated the progression of infection over time. (Fig. 1d).
Separating cells based on their cell cycle state (G1 vs. non-G1) showed that, for a given harvesting time point, non-G1 cells generally contain more viral transcripts (Fig. 1e), suggesting that S, G2 and M phase cells are more susceptible to viral infection, and/or that the infection progresses faster in these cells.
Consequently, at 5 hpi we observed that cells bearing high levels of HSV-1 mRNA (8-30%) showed a lower nUMI count when related to the number of detected genes ( Supplementary Fig. 1c), indicating less complex transcriptomes due to large number of viral transcripts and/or reduction of fibroblast encoded mRNAs likely as a consequence of the host cell shut-off [2].

Precise definition of the temporal order of viral gene expression
Within lytic infection, viral genes have been classified as immediate early (alpha), early (beta), and late (gamma1, gamma2) [26,27]. Since the temporal order of genes expressed from the virus genome is intrinsically encoded in their single-cell expression profiles, we developed an approach to derive a precise model of the early viral gene expression cascade.
As a proxy of the sequential appearance of virus-encoded transcripts, we counted, for each gene, the percentage of infected cells in which it was detected ( Fig. 2a). In order to focus on early viral transcriptional events, only cells harvested at 1hpi and 3hpi were used (cells harvested at 5hpi are depicted in Supplementary Fig. 2b). Interestingly, viral gene transcription appears to start mostly around the internal repeats, and around gene UL23 (Fig. 2a).
For the first eight viral genes as defined by this ordering, normalized expression values were represented in Fig. 2b. Cells were first sorted according to appearance of the eight genes, and then by the abundance of US1 (coding for ICP22), or, if absent, UL54 (ICP27). Under the assumption that the progression of infection, at least in an early phase, is correlating with the accumulation of viral transcripts, this represents a pseudo-time course of the lytic infection. In agreement with previous findings [26,27], US1 and UL54 were the first genes to appear.
Early in infection, there are likely three distinct sets of cells, with either US1 or UL54, or both genes expressed. A comparably high number of cells in these three first states suggested that viral transcription progresses slowly until a number of US1 and UL54 transcripts were present, indicating that the two earliest genes US1 and UL54 were transcribed before the UL23/UL50 transcripts are expressed.
As shown above (Fig. 2a), we observed a bimodal distribution of the relative densities of HSV-1 transcript amounts per single cell. We connected the bimodal distribution of viral transcripts per single cell (Fig. 1e) to the sequential expression of viral genes in two ways. First, cells in the first peak bear mainly transcripts from one or two viral transcripts (Fig. 2b, Supplementary Fig. 2b).
Second we generated smoothened two-dimensional densities of cells with plotting the number of detected genes vs. the percentage of viral transcript per cell ( Supplementary Fig. 2c). Together, this indicates that the bimodality in Fig. 1e arises from the transition between US1/UL54 expression and later stages of infection.
Calculating the likelihood of transitions between single cells based on diffusion maps [28,29] indicated "bottlenecks" in the accumulation of HSV-1 transcripts, as represented by clusters of high transition probabilities ( Supplementary Fig.   2d). Cells in the two-gene state (US1 and UL54) showed a higher likelihood to transition into the four-gene state, deriving a early gene expression cascade shown in Fig. 2c.

The small GTPases RASD1 and RRAD are expressed in infected cells and reduce virion production
To study changes in viral gene expression upon infection, we first examined differential gene expression by bulk mRNA sequencing ( Supplementary Fig. 3a).
To identify genes differentially expressed only in infected cells, we plotted the correlation of host gene expression to viral transcripts in single cells against the maximal fold change in the population mRNA-seq data (Fig. 3a).
Among host genes, activated only in infected cells, we identified hemoglobin alpha genes, HBA1 and HBA2 (Fig. 3a), which were previously shown to be induced by the viral transcription factor ICP4 [7]. Additionally, two genes encoding Ras-related small GTPases, RASD1 (Fig. 3b) and RRAD (Fig. 3c), showed strong positive correlations with viral transcripts. Similarly, we found that RASD1 and RRAD were up-regulated in previously published microarray, RNA-seq, and ribosome profiling datasets [5,6,30,31], indicating that the induction of these two genes is independent of the cell type and virus strain. To investigate the influence of RASD1 and RRAD on viral production, we depleted RASD1 and RRAD by siRNA knockdown prior to infection (Fig. 3d,e), and observed increased viral DNA in the supernatant 24 hpi, suggesting that both genes provide some antiviral activity.

Distinct subpopulation of cells marked by mutually exclusive NQO1 and SULF1 gene expression
To identify distinct subpopulations of cells, defined by a particularly high or low expression of specific marker genes, we used an overdispersion analysis [32].
Not surprisingly, cell cycle, viral genes and host cell genes correlating with viral infection appeared as the strongest markers ( Fig. 4a-e). In addition, two other gene sets, not directly related to the progression of infection and the cell cycle, outlined distinct subpopulations. Importantly, the genes discussed here do not necessarily by themselves influence the infection, but rather are indicators of a specific cellular state that could favor or impair the infection.
The first subpopulation was marked by high mRNA levels of the sulfatase SULF1 Whereas SULF1 was previously shown to be induced by TNFα in MRC-5 fibroblasts [33], we found no clear activators for OXTR, but indications that proinflammatory cytokines such as IL-1β, IL-6, and TNFα might induce its transcription [34]. For the activation of NQO1 and the correlating genes FTH1, FTL, and SQSTM1, several reports pointed to the transcription factor NFE2L2 (also known as NRF2) [35][36][37].
Pairwise gene expression analysis revealed a strong anticorrelation of SULF1 and NQO1 expression in single cells, whereas FTH1 and FTL correlated well with NQO1 (Fig. 4h). On the other hand, OXTR expression levels correlated with SULF1 but not with NQO1, indicating the existence of subpopulations of cells with distinct marker gene expression.
Next we related these subpopulations to the progression of infection. To this end, we used graph abstraction (PAGA) [38], which calculates transition probabilities between groups of cells (Fig. 4i). Cells with high NQO1 levels, and therefore high preceding NRF2 activity, had a relatively low probability to progress into the infection, compared to cells from groups D and F.

RNA velocity shows interruption of the cell cycle by HSV-1 infection
To infer expression dynamics in individual cells, we used RNA velocity [39]. RNA velocity uses the ratio of spliced to unspliced mRNAs to derive transcriptional rates. RNA velocity values are shown for the genes mentioned above ( Fig. 4j-n).

The transcription factor NRF2 is activated upon viral infection
We used the inferred RNA velocity (transcriptional rates) as a precise read out of activation of upstream signaling pathways. To visualize pathway activity, we mapped the cells on a two-dimensional embedding not based on the mature mRNA expression levels as in Fig. 4, but based on the RNA velocity values (Fig.   5a,b). Cells with high SULF1 transcription (Fig. 5c, lower panel), which could reflect TNFα-responsive cells, clustered together in the lower left part. These cells are mostly harvested at later time points but do not show high expression of HSV-1 transcripts, reflecting the established antiviral activity of TNFα [41,42].
Cells with larger amounts of viral transcripts showed relatively high NRF2 activity as deduced from high NQO1 RNA velocity (Fig. 5d, lower panel) suggesting that NRF2 is activated as a part of the cellular defense against the progressing infection.
We again applied the PAGA algorithm to cells clustered based on transcriptional activity (Fig. 5e). The highest probability to proceed into infection are cells in the S/G2/M phases (groups F, N, Ja, supporting again that these phases of the cell cycles favor infection. On the other side, cells with high SULF1 transcription are in cluster C, which shows a comparably lower probability to transition into infection, confirming the lower susceptibility of cells with an activated TNFα response to HSV-1. To summarize, we made two observations regarding NRF2 activity and HSV-1 infection. The analysis of mature mRNA distribution showed that cells with a high level of transcripts of NRF2-driven genes, and therefore high preceding NRF2 activity, have a low transition probability into later stages of infection ( Fig.   4). Looking at RNA transcription however showed that cells at later stages of infection appear to respond by increasing transcription of NRF2 target genes ( Fig. 5), which could reflect a cellular defense mechanism against HSV-1 infection.
In order to strengthen our interpretation of the data, we also analyzed an independent experiment with again two biological replicates, where a procedure

The NRF2 agonist Bardoxolone methyl restricts HSV-1 infection
Our analysis suggested an activation of NRF2 in a subset of infected cells. Under physiological conditions, NRF2 is repressed by KEAP1, which sequesters NRF2 and facilitates its ubiquitination and degradation [44]. Upon disruption of this interaction in response to oxidative stress, NRF2 translocates into the nucleus, and induces transcription of a number of target genes, including NQO1. Using the NRF2 agonist Bardoxolone methyl [45] at final concentrations of 0.1-0.4 µM, we observed an increase in NQO1 mRNA expression in primary fibroblasts (Fig. 6a) and in HEK 293 cells ( Supplementary Fig. 6a), as expected for NRF2 activation [46,47]. At the concentrations applied here, Bardoxolone methyl did not have any apparent effect on cell growth (data not shown). Since, as described above, the progression of infection depends on the cell cycle, we also probed mRNA levels for the cell cycle markers GMNN and TOP2A introduced in the previous sections ( Supplementary Fig. 6d). In the primary fibroblasts but not HEK 293 cells, the levels of these mRNAs were somewhat reduced, indicating that Bardoxolone methyl could also dampen cell cycle progression or promote cell cycle exit.  Fig. 6c).

Discussion
We performed single-cell RNA-sequencing of primary human fibroblasts at early stages of infection with HSV-1. Using the resolution of scRNA-seq, we showed that cells in the S/G2/M cell cycle phases bear more viral transcripts. The relationship between the cell cycle and the progression of infection has been studied for decades [48][49][50][51], with most reports pointing to a G1 or G1/S arrest upon HSV-1 infection. Single-cell RNA sequencing now allows for perturbationfree analysis of the relationship between the cell cycle and progression of infection, reducing the risk of experimental artifacts. We have seen that cells in S/G2/M phases on average bear more viral transcripts (Fig. 1f). This is corroborated by recent findings that the activity of CCNE1 (also known as cyclin E) and CDK2, which promote the G1/S transition, correlates with productive infection [52]. Similarly, high levels of the G1/S transition marker GMNN (also known as geminin) favored infection in microscopy-based experiments [53]. We Virus infection leads to changes in host cell transcription. We have identified a number of deregulated genes in bulk-RNA sequencing, but only using the scRNAseq data we could distinguish on how the deregulation is related progression of infection in the same cell. Two small GTPases, RASD1 and RRAD, were directly induced by the infection. We have observed that both possesses antiviral properties (Fig. 3). Interestingly, these two genes are usually not expressed in epithelial cells. They are however putative targets of the primate specific gene DUX4 [54]. This germline transcription factor was recently shown to be induced in HSV-1 infected cells [55]. Correlating host gene expression with viral gene expression generally emerges as a focus in scRNA-seq studies of virus infected cells [15,17], and comparative studies might reveal common topics and modulators of viral infection.
The emergence of pseudo-time inference algorithms for scRNA-seq data [56], enabled us to order cells along defined trajectories and made it possible to deduce relationships between biological mechanisms, without artefact-prone perturbances such as gene knockdowns/knockouts. We first used overdispersion analysis [32] to find genes with inhomogeneous expression across the entire dataset, which might therefore play a role in the infection (Fig. 4). The relationship of these genes to the progression of infection was then analyzed using RNA velocity [39] and graph abstraction [38]. We focused on two groups of genes, one represented by NQO1 (target of the transcription factors NRF2), and the other by SULF1 (associated with the TNFα response). Importantly, in bulk RNA-seq data, these genes do not appear as differentially expressed and would likely be omitted from subsequent studies.
NRF2 has previously been linked virus infections [57], however without a clearly defined role. Particularly, in a recent Influenza virus scRNA-seq study, NRF2 activity has been associated with high amounts of viral transcripts [15]. For Rotavirus, a dsRNA virus, NRF2 activation was recently shown to reduce virus production [58]. Our analysis suggests that NRF2 activity in a subset of cells leads to an "escape state", diverging from a progressing lytic infection (Fig. 4), and that NRF2 becomes active upon infection (Fig. 5). For Herpes infections, reports on the role of NRF2 so far are mixed. Treatment of mice with the NRF2 agonist tert-butylhydroquinone protected them from MCMV infection [59]. In contrast, HSV-1 infection in A549 was drastically reduced upon a two-day RNA interference of NRF2 [60]. This treatment however induced the antiviral genes STING and IFI16, which could confound the direct effect of the NRF2 depletion.
We followed up on our observations by using the NRF2 agonist Bardoxolone methyl, which is currently in phase 3 clinical trials for treating chronic kidney disease [62,63]. Interestingly, Bardoxolone methyl impaired the infection. At which stage of the infection viral production is impaired, and whether by mere presence of nuclear Nrf2 or via its target genes, as well as the relevance of reported NRF2 interaction with the NF-kappaB pathway, remains to be investigated.
In summary, our study provides a detailed analysis of the events at the beginning of HSV-1 infection in primary human fibroblasts. We could relate the activity of single genes to the progression of infection. Using overdispersion and cellular trajectory analysis allowed the identification of biologically relevant pathways, demonstrating how scRNA-seq can provide an unprecedented insight in the heterogeneity in viral infections.

Acknowledgments
The authors wish to thank Julia Madela, Sebastian Voigt, Lars Dölken, Benedikt

Competing financial interests
The authors declare no competing financial interests.

Data availability
Raw sequencing reads as well as raw read counts for both the bulk RNA-seq and the single-cell RNA-seq, along with normalized counts and tSNE coordinates and cell cycle information, are available in the NCBI GEO repository, accession number GSE123782.

Code availability
Custom R, Python, Perl, awk and bash scripts are readily available from the authors upon request.

Supplementary Tables
Supplementary Table 1: sample overview   Supplementary Table 2: variable genes   Supplementary Table 3:             Note that only cells with more than 2000 detected genes were used, thus the lower apparent limit for the green dots. Light and dark gray bars designate low SULF1 cells and high SULF1 cells according to the bimodality of the plot. The horizontal gray line denotes the cutoff expression between high and low. Only high SULF1 cells were used to calculate the correlations in Fig. 4h.
Wild-type HSV-1 was derived from a patient isolate and was obtained from Department of Virology, Saarland University Medical School, Homburg, Germany.
The virus was propagated in NHDF cells. At 72 hpi, viral supernatant was collected and sterile-filtered through a 0.45 μm pore size filter and stored at −80°C. Viral titers were determined by plaque assay as described previously [5].

Infections
Cells were seeded in 10 cm dishes (for scRNA-seq) and in 6 well plates (for bulk RNA samples and slides for immunofluorescence). For the non-synchronized infection at 37 °C, half of the medium (conditioned medium) was removed, and then cells were incubated for 30 min with HSV-1 (MOI 10) or were mock infected. The supernatant was removed and cells were washed with PBS, followed by re-applying the conditioned medium. We harvested uninfected cells and at early time points (1hpi, 3hpi, and 5hpi) in two biological replicates.

DropSeq single-cell RNA-sequencing
Cells were fixed and used for Drop-seq single-cell sequencing as previously described [24]. Monodisperse droplets of about 1 nl in size were prepared on a self-built Drop-seq set up following closely the instrument set up and library generation procedure as described [23].

Bulk RNA-sequencing and analysis
RNA was extracted from cells using Trizol (Thermo Fisher) and the RNA clean and concentrator kit (Zymo). Sequencing libraries were prepared using the TruSeq Stranded mRNA kit (Illumina) according to the manufacturer's instruction and sequenced on a HiSeq4000 device using 2x76nt paired-end sequencing. Sequencing reads were aligned to the hg38 version of the human genome using tophat2 [64]. Standard differential expression analysis was performed using quasR [65] for counting reads (using the Gencode gene annotation version 28) and edgeR [66].

Single-cell RNA-seq data processing
The scRNA-seq data was processed using the PiGx pipeline [67]. In short, polyA sequences are removed from reads. The reads are mapped to the hg38 version of the human genome using STAR [68] with gene models from Gencode version 28.
The number of cells, for each sample, is determined using dropbead [24]. Finally, a combined digital expression matrix is constructed, containing all sequenced experiments. This table is available at the GEO entry. Cells containing less than tSNE coordinates were calculated using the Seurat package [69]. Cell cycle was assigned to each cell using the CellCycleScoring function. Cell cycle gene sets were taken from [70]. Normalized total HSV-1 transcription was calculated by summing up raw counts of all detected viral genes, divided by the total number of raw counts, multiplied by a scaling factor of 10000 and log(2) transformed.

Diffusion maps for Herpes
Transition probabilities between cells were estimated based on the diffusion pseudo-time distances [29], implemented in the destiny Bioconductor package [28]. Diffusion pseudo-time distance between two cells represents the cumulative probability of travelling from one cell to the other on a probabilistic graph. The graph is constructed by embedding cells using gaussian kernels.
Obtained distances were converted to transition probabilities using an exponential probability distribution.

Binning and correlations
Correlations of host cell gene expression with HSV-1 gene expression (Fig. 4), SULF1 or NQO1 expression (Fig. 5) were calculated using binned cells to reduce noise. For this, cells with "high" expression (see Fig. S3a and S4c) were first sorted according to normalized values of HSV-1 transcripts, SULF1 or NQO1.
Then, cells were binned with 20 cells per bin. Each bin was then treated as a "metacell". Normalized gene expression values for the metacella were calculated by summing up, per gene, all raw values in the bin, dividing by the total raw count in the entire metacell, multiplied by a scaling factor of 10000 and log (2) transformed. Then, the linear Pearson correlation coefficient and slope for each gene or antisense transcript group with HSV-1 transcripts, SULF1 or NQO1, respectively, were calculated.

RNA interference
NHDF cells were plated in 6-well plates and incubated at 37°C. 24

RT-qPCR
RNA was extracted from cells using Trizol (Thermo Fisher) and the RNA clean and concentrator kit (Zymo). DNase treatment using DNase I amplification grade (Thermo Fisher) and RT using SuperScript III (ThermoFisher) was performed according to the manufacturer's protocol. For all experiments, four measurements were done from two biological replicates. Power SYBR Green PCR Master Mix on a StepOnePlus system (both Thermo Fisher Scientific) was used for qPCR. Data were normalized to the indicated timepoint/sample and the GAPDH signal using the ∆∆Ct method [73]. Primers for qPCR were designed using the Universal ProbeLibrary (Roche Life Sciences), tested for efficiency and a single amplicon and listed in supplementary table S4.

UMAP
For visualization purposes, we employed in Fig. 4 and 5 the UMAP dimensionality reduction algorithm [74] instead of tSNE as in the previous figures. The reason was that the directionality, which the RNA velocity algorithm calculates based on all detectable genes, and then projects as arrow clouds on the two-dimensional map (Fig. 4C), could not be visualized when using tSNE.
Variable genes were define using the FindVariableGenes function with the default parameters. First twenty principle components were calculated using the defined variable genes. Dimensionality reduction was performed with UMAP [74]. The cells were clustered using Louvain algorithm with a set of resolution parameters ranging from 0.5 -2.

Velocity and graph abstraction
The data was pre-processed using the Velocyto CLI pipeline. Only cells with more than 500 intronic and more than 500 exonic reads were kept for further analysis.
Intronic signal represented about 3% of all sequencing reads.
The following parameters were used within the Velocyto analysis:  (40,40), n_neighbors=100) For the nascent RNA clustering, the Velocyto estimated transcriptional rates were processed using scanpy [75]. The rates were used to calculate the principal components, UMAP embeddings, cluster the cells using the Louvain algorithm, and compute PAGA [38].    UL1  UL3  UL5  UL6  UL9  UL10  UL14  UL15  UL17  UL21  UL22  UL23  UL24  UL25  UL28  UL29  UL30  UL32  UL33  UL36  UL37  UL38  UL39  UL41  UL42  UL43  UL44  UL47  UL48  UL49A  UL50  UL51  UL52  UL54  UL55  UL56  US1  US2  US4  US6  US8   Color scales are shown to the right of these panels. areas containing cells with relatively high expression levels were marked with a light green ellipse. The dark green and blue ellipse denote S and G2/M phase cells, respectively. h, correlation of gene expression with NQO1 expression (horizontal axis) and SULF1 expression (vertical axis) in cells harvested at 3hpi. Every dot represents a gene, colored by the slope of the linear correlation with NQO1. Selected genes were labeled by name. i, Groups of cells according to PAGA were labeled A to Q. The line thickness of the connections (graph edges) between the cell clusters (graph nodes) indicate the likelihood that cells can move from one cluster to the other (right). Areas of interest are marked with light green ellipses. j-n, Cells colored by RNA velocity values of the indicated genes. Cells for which no value could be calculated were colored in light gray. o, RNA velocity arrows projected on the UMAP. Cells containing at least one detected viral transcript were colored in light orange, all others in light gray.   a, Cells were projected on a two-dimensional map using UMAP according to RNA velocity values, with cells colored by harvesting time point. Main distinguishing features were marked in grey. b, Coloring by amount of HSV-1 transcripts. c and d, Cells colored by expression values (top panels) and RNA velocity values (bottom panels) of the indicated genes. Color scales are shown to the right of these panels. For expression values, areas containing cells with relatively high expression levels were marked with a light green ellipse. For RNA velocity, cells for which no value could be calculated were colored in gray. e, Groups of cells according to PAGA were labeled A to R. Transition probabilities between clusters are proportional to the line thickness between the clusters (right).