Epstein-Barr virus reactivation induces divergent abortive, reprogrammed, and host shutoff states by lytic progression

Viral infection leads to heterogeneous cellular outcomes ranging from refractory to abortive and fully productive states. Single cell transcriptomics enables a high resolution view of these distinct post-infection states. Here, we have interrogated the host-pathogen dynamics following reactivation of Epstein-Barr virus (EBV). While benign in most people, EBV is responsible for infectious mononucleosis, up to 2% of human cancers, and is a trigger for the development of multiple sclerosis. Following latency establishment in B cells, EBV reactivates and is shed in saliva to enable infection of new hosts. Beyond its importance for transmission, the lytic cycle is also implicated in EBV-associated oncogenesis. Conversely, induction of lytic reactivation in latent EBV-positive tumors presents a novel therapeutic opportunity. Therefore, defining the dynamics and heterogeneity of EBV lytic reactivation is a high priority to better understand pathogenesis and therapeutic potential. In this study, we applied single-cell techniques to analyze diverse fate trajectories during lytic reactivation in two B cell models. Consistent with prior work, we find that cell cycle and MYC expression correlate with cells refractory to lytic reactivation. We further found that lytic induction yields a continuum from abortive to complete reactivation. Abortive lytic cells upregulate NFκB and IRF3 pathway target genes, while cells that proceed through the full lytic cycle exhibit unexpected expression of genes associated with cellular reprogramming. Distinct subpopulations of lytic cells further displayed variable profiles for transcripts known to escape virus-mediated host shutoff. These data reveal previously unknown and promiscuous outcomes of lytic reactivation with broad implications for viral replication and EBV-associated oncogenesis.


INTRODUCTION
Viral infections lead to heterogeneous cell fate outcomes including resistance, abortive infection, latency, or full virion amplification often leading to cell death.Cells that resist viral infection often display elevated pre-existing anti-viral responses [1][2][3][4] .Likewise, cell responses that enable survival following virus replication can prime for further anti-viral responses 5,6 .
Herpesviruses are large double-stranded DNA viruses that provide a unique and complex infection paradigm to model the heterogeneity of viral infection as they reactivate from a latent state in response to diverse stimuli.
Epstein-Barr virus (EBV) was the first oncogenic human virus to be discovered 7 .Since its isolation from endemic Burkitt Lymphoma (BL) cells in 1964, EBV infection has been linked to an expansive set of human cancers and, more recently, autoimmune diseases [8][9][10][11] .EBV infection in immunosuppressed individuals can lead to post-transplant lymphoproliferative disease (PTLD) 12 and HIV-related diffuse large B cell lymphomas (DLBCL) 13 as well as up to 40% of Hodgkin Lymphoma (HL) 14 .and rare individuals with chronic active EBV (CAEBV) can develop T and NK cell lymphomas 15,16 .Beyond these hematologic malignancies, EBV infection is associated with epithelial cancers such as nasopharyngeal carcinoma (NPC) 17 and gastric carcinomas 18 .
Collectively, EBV causes is or otherwise associated with nearly 2% of all cancers diagnosed annually 8 .This prevalence in malignant disease vastly underrepresents the success of EBV as a human pathogen.Globally, it is estimated that over 95% of adults are infected with EBV 19 .EBV is transmitted via saliva, which enables the virus to traverse oral epithelial tissues and infect B lymphocytes within the tonsils 20 .EBV infects B cells via the surface receptor CD21 (CR2) 21,22 and rapidly induces B cell adaptive immune programs to mimic germinal center (GC)-like dynamics 23- 26 .Successful evasion of antiviral defenses, immune tolerance checkpoints, and growth-induced damage [27][28][29] allows memory B cells latently infected with EBV to exit from this virus-manipulated GC reaction.Viral latency establishment within the memory B cell compartment yields lifelong persistence 30,31 .Lytic reactivation from this latent state triggers the production of new virions and is essential to the replicative cycle and transmission between hosts.The lytic gene program is transcriptionally orchestrated by two immediate early (IE) lytic genes: BZLF1 (encodes for the transcription factor Zta / Z / ZEBRA) and BRLF1 (encodes for the transcription factor Rta / R) [32][33][34] .
While Zta and Rta both play essential roles in lytic reactivation, Zta is the master lytic transactivator in B cells.BZLF1 expression is induced upon cell differentiation and stress 35,36 , a prototypical example being post-GC B cell differentiation into plasmablasts 37 .Host cell transcriptional regulators of plasma cell generation including XBP1 and BLIMP1 (PRDM1) induce .CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.It is made The copyright holder for this preprint this version posted June 14, 2024.; https://doi.org/10.1101/2024.06.14.598975 doi: bioRxiv preprint EBV lytic reactivation via direct transactivation of the BZLF1 promoter [38][39][40] .Zta then transactivates subsequent expression of early and late lytic genes by binding at Z-responsive elements (ZREs) throughout the viral genome 41 .As an AP-1 family homolog 33 , Zta also binds loci throughout the host genome 42 and has characteristics of a 'pioneer' transcription factor.Consistent with this, BZLF1 expression and the early stages of EBV reactivation cause considerable alterations to the host cell epigenome and resulting gene expression 43,44 .
Prior work suggests that lytic gene expression is functionally important for tumorigenesis.
Notably, viral strains that carry the NFATc1-responsive Z promoter variant Zp-V3 exhibit increased lytic replication and are enriched in EBV-associated cancers relative to strains with prototypical Zp 45 .In SCID and NSG mouse models with reconstituted human immune systems, significantly fewer animals developed EBV + lymphomas after infection with BZLF1 knockout virus versus a wild-type (WT) control strain 46 .Further, infection with a Zta-overexpressing strain that failed to complete reactivation (i.e., abortive lytic) promoted tumor growth in mice similar to WT EBV 47 .Recent experiments in immunocompromised mice confirmed the tumorigenic role of abortive lytic infection by using EBV lacking the BALF5 gene, which encodes a viral DNA polymerase subunit essential for lytic replication 48 .These studies demonstrated that expression of BZLF1 (and possibly other early lytic genes) contributes to tumorigenesis in vivo regardless of the potential for horizontal infection of bystander cells by new virions.While detailed insights regarding the oncogenic effects of successful or abortive lytic replication are limited, tumor microenvironment inflammatory conditioning by cytokines secreted from reactivating cells has been proposed [49][50][51][52][53] .
Another complication in the relation between viral reactivation and oncogenicity stems from observations that a significant proportion of EBV-infected tumor cells are resistant or otherwise refractory to lytic reactivation.In Burkitt Lymphoma-derived P3HR1 and Akata cells, high expression of the oncoprotein c-Myc promotes viral latency maintenance and suppresses lytic reactivation via direct interaction with the origin of lytic replication (oriLyt) and inhibition of chromatin looping to activate BZLF1 expression 54 .Accordingly, MYC suppression facilitates BZLF1 expression and the subsequent induction of viral lytic genes.It is noteworthy that constitutive oncogene expression favors viral genome propagation through proliferation of latently infected host cells whereas lytic replication becomes a more advantageous strategy in its absence.Similarly, BL-derived cells refractory to lytic reactivation have also been found to express high levels of STAT3 [55][56][57] , which functions as an oncogene in B cells and inhibits apoptosis via induction of BCL2 expression.Beyond simply being expressed by refractory cells, STAT3 antagonizes lytic reactivation of EBV + cells through the functions of its transcriptional targets 56 .In fact, LCLs derived from patients with autosomal dominant hyper-IgE syndrome (AD-HIES), a disease that leads to non-functional STAT3 activity, went lytic at a higher rate than LCLs derived from healthy donors 58 .Given the therapeutic potential of drug-induced lytic reactivation followed by viral DNA synthesis inhibition to treat EBV-latent cancers, investigators are actively exploring means to make refractory cells more sensitive to lytic induction [59][60][61] .However, such efforts should be weighed against the known associations between the EBV lytic cycle and oncogenesis, which remain to be fully elucidated.
Many EBV gene products contribute to virus-driven malignancies by mediating functions associated with cancer hallmarks including uncontrolled proliferation, tumor suppressor inhibition, epigenetic reprogramming, genome instability, apoptotic resistance, and immune evasion 62 .EBV + cells with cancer stem cell (CSC) features have also been reported in NPC and gastric carcinoma 63,64 , suggesting the potential for cellular self-renewal associated with infection.In the CSC model, a small subset of tumor cells retain the capacity for self-renewal and proliferation through activation of signaling pathways (e.g., Wnt, Notch), transactivators of the epithelial-tomesenchymal (EMT) transition, and critical regulators of pluripotency (e.g., SOX2, OCT4).CSCs may serve as progenitors for other tumor cells, especially in lymphoid malignancies that are derived from cells of origin that intrinsically retain self-renewal properties to support immunologic memory [65][66][67] .Aberrant expression of self-renewal genes and other CSC biomarkers 68 may originate from significant (epi)genomic reprogramming and result in cellular phenotypic plasticity.
Lytic replication of EBV (and DNA viruses from several other families 69 ) clearly constitutes a major reprogramming event for the host cell.Nuclear chromatin is globally disrupted by IE gene expression, the formation of viral replication compartments, and the accumulation of viral DNA 43,70 .Moreover, preferential binding of BZLF1 to methylated promoters can reverse epigenetic silencing of both EBV and cellular genes through nucleosome eviction, resulting in heterochromatin-to-euchromatin conversion 44,[71][72][73] .While evidence for stem-like reprogramming and CSC gene expression during the EBV lytic cycle has not been reported to our knowledge, it is noteworthy that reactivation of HSV-1 (another herpesvirus) induces embryonic development programs including Wnt/-catenin activity that licenses late viral gene expression 74 .
These previous studies demonstrate that EBV reactivation from latency is a complex process that culminates in heterogeneous host cell responses germane to the progression of virusassociated cancers.Single-cell sequencing techniques are particularly well suited to dissect the inherent complexity of host-virus interactions and their effects on cell fate [74][75][76][77] .In recent studies of early EBV infection 25,26 and established latency 78,79 , we have used single-cell sequencing to successfully resolve and study diverse phenotypes arising from complex host-pathogen dynamics.We reasoned that a similar high-resolution experimental and informatic approach would clarify distinct courses of lytic reactivation, provide essential data for future studies of viral pathogenesis, and inform potential therapeutic strategies to address EBV-driven oncogenesis.To this end, we performed time-resolved single-cell RNA sequencing (scRNA-seq), flow cytometry, and RNA Flow-FISH (fluorescence in situ hybridization) in P3HR1-ZHT cells to define initial cell state diversity, differential fate trajectories, and previously unknown lytic response phenotypes within this widely used EBV + Burkitt Lymphoma model.Cellular transcriptomic responses to lytic reactivation were investigated with respect to IE and early versus late viral gene programs and subsequently validated in the B958-ZHT LCL.

Heterogeneous responses to EBV lytic reactivation in individual cells
P3HR1-ZHT cells are an inducible model of EBV lytic reactivation (Fig. 1A).This model system constitutively expresses the EBV immediate early lytic transactivator Zta (encoded by the BZLF1 gene) fused with a modified murine estrogen receptor hormone binding domain.While the encoded fusion protein is normally rapidly degraded, addition of 4-hydroxytamoxifen (4HT) stabilizes it and promotes its nuclear translocation, whereupon the Zta domain binds and transactivates Zta-responsive elements (ZREs) in both host and viral genomes.Because Zta has positive regulatory control of its own promoter via ZRE binding 80 , 4HT treatment also leads to expression of endogenous BZLF1, thus initiating viral lytic reactivation.Although all cells in the P3HR1-ZHT line express the inducible construct, it has been observed that complete EBV lytic reactivation occurs only in a subset of 4HT treated cells 81,82 .
We confirmed inducible yet non-uniform viral reactivation of P3HR1-ZHT cells in response to 4HT treatment using FACS staining for the viral glycoprotein gp350, which was expressed in cells that reached the late stage of lytic reactivation.Unstimulated P3HR1-ZHT cells expressed minimal gp350 (1.1%), but treatment with 100 nM 4HT for 24 hours resulted in gp350 expression in 19.2% of cells.When we simultaneously treated cells with 4HT and PAA, an inhibitor of viral DNA replication, we observed a significant reduction in gp350 expression by 24 hours (Fig. 1B, Fig. S1).These results indicated that cells exhibited heterogenous responses to viral lytic reactivation and that completion of the full lytic cycle was dependent upon successful viral DNA replication, which has been previously described in herpesviruses [83][84][85][86][87] .We expanded upon these gp350 FACS results using RNA Flow-FISH assays to detect viral RNAs from genes expressed at different stages of the lytic cycle: the immediate early lytic gene BZLF1, the early lytic gene BGLF4, and the late lytic gene BLLF1.After 24 hours of 4HT treatment, we observed a significant increase in expression of all three lytic transcripts compared to mock treated cells.However, there was a stepwise decrease in expression level between early and late lytic genes (Fig. 1C, Fig. S2).These results confirmed that a significant proportion of Z-HT induced P3HR1 cells were refractory to full lytic reactivation.
Since we observed heterogeneous responses upon lytic reactivation, we applied timeresolved single-cell RNA sequencing (scRNA-seq) to study the concurrent cellular responses in the P3HR1-ZHT system after 24, 48, and 72 hours of 4HT treatment compared to untreated cells (Fig. 1D).UMAP projection of samples by timepoint demonstrated that substantial transcriptomic changes occurred after 4HT stimulation (Fig. 1E).Cells expressing high levels of viral reads clustered together, however there was a distinction between cells expressing immediate early, early, and late viral transcripts (Fig. 1F).Analysis of all EBV transcripts identified genes with high, moderate, and low expression; however, all 4HT-treated samples expressed more viral transcripts compared to untreated cells (Fig. 1G).These results confirmed heterogeneous responses to lytic reactivation observed by flow cytometry and enabled subsequent genome-wide analyses.

Identification of distinct EBV reactivation response clusters
Cells from integrated timecourse scRNA-seq libraries were hierarchically clustered by host and viral transcriptome similarity, which led to the identification of five main clusters (Fig. 2A).
Unstimulated cells were mostly present in clusters A and B, while clusters C, D, and E primarily comprised 4HT-treated cells across the experimental time course (Fig. 2B) and displayed elevated viral gene expression compared to clusters A and B (Fig. S3).Further examination of these clusters revealed differences in the number of total and unique RNAs, the percentage of viral RNAs, and the percentage of mitochondrial RNAs (Fig. 2C).These differences in unique and total RNA features suggested major phenotypic differences both in unstimulated and reactivated cells.Therefore, we scored the clusters based on cell cycle state and found that there was a decrease in G2/M specific gene expression and an increase in G1 gene expression after 24 hours of 4HT treatment, consistent with EBV lytic reactivation occurring in a pseudo-S phase 88,89 (Fig. S4A).We confirmed this finding using BrdU/7-AAD staining of untreated versus 4HT-treated cells (Fig. S4B).Consistent with induced cell cycle arrest, lytic reactivation upon 4HT treatment led to a reduction of S phase cells (43.2% vs. 54.3%)and modest increase in G0/G1 cells.Because pulsed BrdU staining does not discriminate cellular and viral DNA synthesis, a portion of S phase 4HT-treated cells were likely undergoing viral but not cellular DNA synthesis.This was further evidenced by a significant fraction of gp350 + cells within the gated S phase population (Fig. S4B).
We also assayed MitoTracker signal stratified by gp350 expression and found that gp350 + cells had lower mitochondrial content (Fig. S4C).

Cells traverse heterogeneous biological response trajectories during lytic reactivation
Next, we analyzed differentially expressed genes by cluster and grouped them by ontology using a combined approach with software-based annotation tools 90 and primary literature searches (Fig. 2D).Unstimulated cells were almost exclusively present in clusters A and B, which were distinguished from each other by total transcripts and unique features per cell (Fig. 2B-C).

Unstimulated cells with high RNA and feature counts (cluster A) exhibited a germinal center (GC)
B cell profile including MME (CD10) 91 , BCL6 92,93 , BCL11A 94 , POU2F2 (OCT2) 95 , and AICDA (AID) 96,97 .Along with high MYC expression, this phenotype is consistent with the profile of endemic BL from which P3HR1-ZHT is derived.In contrast, unstimulated cells with low RNA and features counts (cluster B) exhibited a cell stress expression signature that included slight enrichment of genes for ribosomal subunits (RPL34, RPS27), nuclear-encoded components of mitochondrial respiratory complexes (COX7C), and the apoptotic resistance genes PTMA 98 and GSTP1, the latter of which also mediates oxidative stress 99 .Cluster C, which was comprised of 4HT-treated samples, displayed antiviral restriction (APOBEC3G, PPP1R15A, TRIM14, FURIN), inflammatory (CCL4L2, CCL3L1, NKG7), and NF-B signaling (NFKBIA, ICAM1, CD83, BCL2, BCL2A1) signatures.Cluster D had a similar gene expression pattern to cluster B with the addition of lytic transcripts and several long noncoding RNAs from R-loop "hot spots" (C1orf56, AC092069.1,AC005921.2,AC106707.1)associated with genomic instability related to unscheduled gene expression or DNA synthesis (in contexts including herpesviral reactivation) 100- 104 .Finally, cluster E primarily contained cells that had entered the lytic cycle after 4HT treatment.
We next focused on individual genes that are differentially expressed between the clusters.
We specifically chose STAT3 and MYC because they have been established as key regulators of EBV lytic reactivation 54,56,57,59 (Fig. 2E).In line with these published results, MYC expression was Given the observed heterogeneity of phenotypic states before and after lytic induction, we aimed to better understand the distinct response trajectories of EBV-infected cells using pseudotemporal ordering (Fig. 2G).Pseudotime analyses 111

Abortive lytic cells are characterized by high NF-B pathway gene expression
Abortive lytic replication, or the initiation of the lytic cycle without expression of late lytic genes / proteins, has been identified in various systems 47,112 .We sought to characterize this replication sub-state further through analysis of the abortive lytic cells in the cluster C phenotype.Using markers identified in Fig. 2D we were able to clearly distinguish unstimulated, abortive lytic, and lytic cells using CD38, BCL2A1, and BLLF1 expression, respectively (Fig. 3A).STAT3 + cells in the BZLF1 + abortive lytic state (cluster C) notably co-expressed BCL2A1 and other NF-B pathway target genes (Fig. 3B).RNA Flow-FISH for CD38, BCL2A1, and BLLF1 in cells treated with DMSO (control), 4HT (lytic), and 4HT + PAA (an abortive lytic model due to inhibited viral DNA synthesis) confirmed these distinct response states (Fig. 3C, Fig. S6).This experiment confirmed that CD38 RNA was primarily expressed in unstimulated cells and decreased upon 4HT treatment.BLLF1 (gp350) RNA was almost exclusively expressed in 4HT treated cells, and its expression was blocked upon PAA treatment as expected.BCL2A1 RNA was significantly elevated in 4HT + PAA-treated cells, especially by 48 hours post-treatment (Fig. 3D).Thus, these markers reliably delineated latent, abortive, and lytic phenotypes identified from scRNA-seq as clusters A, C, and E.
Because EBV LMP-1 partially mimics the activated CD40 receptor that induces NF-B signaling, we reasoned that LMP-1 might be associated with the abortive lytic phenotype.
However, LMP-1 expression was largely restricted to cluster E (Fig. 3E), consistent with its transcription during the lytic cycle 113,114 .This observation suggested that the abortive lytic phenotype and associated NF-B signaling was not dependent upon LMP-1 expression.We confirmed this finding through FACS detection of gp350 (lytic cells) and ICAM1, a surfaceexpressed proxy for NF-B pathway transcriptional activation (and in this context, abortive reactivation).Untreated P3HR1-ZHT cells did not express gp350 or ICAM1 (Fig. 3F, Fig. S7).
Treatment with 4HT induced expression of both gp350 and ICAM1; notably, expression of these proteins was observed in distinct cell subpopulations, supporting our finding that NF-B signaling was primarily active in cells that had not entered the full lytic cycle.Accordingly, co-treatment with 4HT + PAA to induce an abortive lytic state by blocking viral DNA synthesis led to increased ICAM1 + cell frequency consistent with the BCL2A1 upregulation observed in Fig. 3D.Conversely, co-treatment with 4HT and an inhibitor of IKK (a key component of NF-B signaling) eliminated ICAM1 expression, but did not increase gp350 expression.These results demonstrated that NF-B signaling is a feature of abortive lytic cells that is independent of LMP-1 activity, but does not restrict late viral gene expression.

Lytic subpopulations are reprogrammed to stem-like plasticity during EBV reactivation
We next focused on the lytic fate by analyzing cells in cluster E. Paradoxically, lytic cells in cluster E collectively expressed the most unique genes (i.e., transcript diversity) of any cluster despite having low mRNA density per cell consistent with host shutoff (Fig. 4A).In addition to differences in early and late lytic gene expression across cluster E (Fig. 1F), this observation was consistent with enhanced cell-to-cell variability in gene expression.We therefore subclustered cells at higher resolution to examine heterogeneity among lytic subpopulations (Fig. 4B).This yielded three subclusters of BZLF1 + cellsone with high late gene expression corresponding to complete reactivation (cluster E1) and two with comparatively lower late gene expression (clusters E2 and E3).Differential expression analysis by subcluster revealed remarkably broad cellular plasticity and developmental pluripotency signatures in E2 and E3 (Fig. 4C).Although MYC was downregulated, the master pluripotency regulators POU5F1 (OCT4), SOX2, KLF4, NANOG, and LIN28A were expressed in E2 and E3 [115][116][117][118][119][120] .Intriguingly, many essential transcriptional regulators of pluripotency exit and germ layer specification were also co-expressed with BZLF1 + in E2 and E3 lytic subpopulations.Expression of ALDH1A1, ALPL, ITGA6, CD44, PROM1 (CD133), LGR5, and YAP1 upregulated in the E2 and E3 phenotypes was consistent with cancer hallmarks including cell plasticity, self-renewal, and drug-tolerant persistence 68,[121][122][123][124][125][126] .Related to YAP1 expression, we identified distinct Hedgehog 127 , Notch 128 , and Wnt 129,130 signaling pathway signatures in E2 and E3 lytic phenotypes as well as Hippo-independent YAP pathway 131 components reported in cancer.E3 cells also expressed genes encoding several PIWI-like family proteins, which protect germline cell genomes from transposable element insertion, maintain stemness, and are upregulated in some cancers [132][133][134][135] .
In total, 6,900 of 26,728 cells (25.8%) across all sampled timepoints expressed BZLF1 transcripts.Co-expression of genes including ALDH1A1 and SOX2 in a subset of BZLF1 + cells demonstrated an association between cellular plasticity and EBV lytic reactivation (Fig. 4D).GRN analysis further supported a role for SOX2 transcriptional activity in a fraction of lytic cells (Fig. S9).RNA Flow-FISH validated ALDH1A1 and SOX2 expression in BZLF1 + cells in 4HT treated P3HR1-ZHT cultures (Fig. 4E, Fig. S8).We also used flow cytometry to validate increased expression of the CSC biomarkers CD44, CD133 (PROM1), and CD166 (ALCAM) at the protein level in gp350 + cells (late lytic) relative to gp350 -subsets across treatment conditions (Fig. 4F, Fig. S10).
We next examined whether lytic cycle initiation was sufficient to induce CSC-associated pluripotency expression or if successful viral DNA synthesis was required.To do so, we used RNA Flow-FISH to detect BZLF1, ALDH1A1, and BLLF1.ALDH1A1 was expressed in BZLF1 + BLLF1 + cells following 4HT treatment, consistent with its expression in late stages of lytic reactivation (Fig. 4G, left and middle panels).Consistent with a role for viral DNA replication in CSC gene induction, co-treatment with PAA and 4HT diminished BZLF1 + BLLF1 + cell frequency and ablated ALDH1A1 expression (Fig. 4G, right panel).Collectively, these data support a unique program of cellular plasticity induced in the late phase of EBV lytic reactivation.

Host shutoff escapees in lytic subclusters exhibit distinct ontologies
Because lytic subclusters identified at high resolution displayed distinct cellular transcriptomes, we asked whether host shutoff responses differed among lytic cells.RNA for BGLF5, an early EBV lytic gene that mediates host shutoff 107 , was detected at variable levels across lytic cells and inversely correlated with per cell mRNA feature density as expected (Fig. 5A).Moreover, transcripts for genes previously found to escape host shutoff 109,110 were identified in each lytic subcluster (E1, E2, and E3) (Fig. 5B-C).Host shutoff escapee expression could be broadly categorized by two patternssome escapees (e.g., C19orf66, CDKN1B) were expressed in unstimulated P3HR1-ZHT cells and retained across abortive and lytic cells, whereas other escapees (e.g., IL6, SERPINB2, LHX1, JAG1) were exclusively expressed in lytic cells (Fig. 5C).
Intriguingly, lytic subclusters exhibited different host shutoff escapee profiles.Anecdotally, we also noted that several escapees in clusters E2 and E3 were related to inflammatory responses and overlapped with CSC and developmental pluripotency signatures (Fig. 5D).We applied gene ontology (GO) analyses to differentially expressed genes among lytic subclusters to further investigate potential biological differences.Cells in E2 displayed significant enrichment of GO terms related to mRNA splicing and post-transcriptional regulation and epigenetic regulation versus cells in E3 (Fig. 5E, top panel).RNA processing GO terms were also upregulated in E2 when compared jointly against clusters E3 and A to filter out differences related to transcripts basally expressed in unstimulated cells (Fig. 5E, bottom panel).Conversely, the top enriched GO terms in cluster E3 versus E2 were related to cell-cell adhesion, morphogenesis, and diverse tissue-specific developmental programs (Fig. 5F).Relatively few cellular GO terms were enriched in fully lytic cells (E1), consistent with extensive host shutoff and predominantly viral gene expression (Fig. 5G).

Phenotype validation across viral strain and host background
Finally, we confirmed key findings through additional independent scRNA-seq experiments capturing responses of B958-ZHT cell lines to 4HT treatment (Fig. 6) and technical replication in P3HR1-ZHT (Fig. S13).Unstimulated and 24 h post-4HT B958-ZHT cell libraries were generated and analyzed as in previous experiments (Fig. 6A).High-resolution cluster annotations from P3HR1-ZHT scRNA-libraries were mapped to B958-ZHT cells by anchor feature identification and transfer to evaluate the preservation of biological phenotypes across cell systems (Fig. 6B).Cells corresponding to each high-resolution cluster were identified in the B958-ZHT dataset.Viral IE, early, and late gene expression modules were also scored across B958-ZHT cells and compared against scores for the three lytic subclusters (Fig. 6C).As in the P3HR1-ZHT system, E1 cells exhibited high late gene expression consistent with complete reactivation while E2 and E3 cells displayed reduced late gene scores.In B958-ZHT, the E3 cluster most closely associated with plasticity and self-renewal signatures had the lowest IE, early, and late expression relative to other cells in lytic clusters.Prior findings of viral gene anticorrelation with MYC, STAT3, and BCL2A1 (Fig. 6D) and lytic cell upregulation of cancer-associated stem-like pluripotency and host shutoff escapees (Fig. 6E, Figs. S11-S12) were conserved in B958-ZHT.Thus, our findings in the P3HR1-ZHT system are applicable across EBV strains and host cell genetic backgrounds.

DISCUSSION
The single-cell data presented herein substantially expand and refine transcriptome-wide contours of host-virus dynamics during the EBV lytic cycle.For example, prior studies discovered that EBV-infected BL cells are prone versus resistant to reactivation dependent on STAT3 expression, activity, and functions of its downstream transcriptional targets 56,57,59 .A population of STAT3 -/lo cells in unstimulated P3HR1-ZHT revealed by scRNA-seq (cluster B), which exhibits globally reduced mRNA levels consistent with cellular quiescence, may be more permissive to successful reactivation than cells with basally elevated STAT3 (cluster A).Additionally, cells that undergo abortive replication retain STAT3 expression (and predicted transcriptional activity) after stimulation, while STAT3 and host transcript loads are drastically reduced in fully lytic cells, consistent with host shutoff functions exhibited by diverse viruses [136][137][138][139][140] .Single-cell data are also consistent with the functional importance of c-MYC in regulating EBV latency versus lytic reactivation 54 .MYC expression exhibits cluster-level patterns similar to STAT3, with the notable exception that MYC is more strongly expressed in cluster B cellslikely due to constitutive expression resulting from the chr8:chr14 (Ig-MYC) translocation in BL.Single-cell sequencing and RNA-FiSH results further identify unique upregulation of NF-B and IRF3 pathway transcriptional targets in abortive lytic cells.Paired with STAT3 and MYC activity, we speculate that this concerted response might sustain viability and reinforce latency in cells that fail to meet the lytic switch threshold.
Acquisition of cellular plasticity within lytic cell subsets in multiple EBV + B cell models is particularly striking.Several aspects of the lytic cycle could conceivably contribute to host cell plasticity through reversing epigenetic repression of lineage-ectopic genes.As observed across several DNA virus families, EBV genome replication within intranuclear compartments induces dramatic reorganization of host chromatin 69,70,73,141 .Along with this alteration to nuclear architecture, Zta binding at accessible AP-1 recognition sequences 33 (particularly methylated sites 71,72 ) may reverse epigenetic silencing through supporting nucleosome eviction, enhancement of chromatin accessibility, and recruitment of transactivators to facilitate aberrant gene expression 43,44 .ChIP-seq for Zta has revealed many such potential sites throughout the host genome, including POU5F1 (Oct-4) 42,43 .From the viral perspective, Zta binding across the cellular genome may function as a "sink" that supports bimodal control of the switch between latency (Zta absence or noise-level expression) and lytic reactivation (high Zta) 43 .From the host perspective, our findings suggest that these BZLF1 interactions with cellular DNA and nuclear chromatin remodeling during later stages the lytic cycle have substantialand potentially pathogeniccollateral effects on biological reprogramming.Along these lines, developmental reprogramming associated with Wnt/-catenin signaling has been observed in single-cell study of HSV-1 lytic infection 74 .
Additionally, DNA damage, antiviral nucleic acid sensing, cytoskeletal rearrangements, and other major mechanobiological changes that manifest during reactivation may activate intrinsic responses to cellular injury leading to NF-B and IRF3 signaling [142][143][144] .Paired with lytic-mediated growth arrest 145,146 , we speculate that this process may engage cellular senescence and injury responses that promote autocrine and paracrine cellular reprogramming.An essential feature of damage-associated induction of cellular pluripotency is upregulation of pro-inflammatory cytokines such as IL-6 147 .In both P3HR1-ZHT and B958-ZHT scRNA-seq datasets, IL6 was exclusive to fully lytic cell subsets.However, IL6R was expressed in abortive cells in P3HR1-ZHT and most latently infected cells in B958-ZHT.Expression of JAK1/2 and STAT3 in latently infected cells from both lines was suggestive of an IL-6 response axis (IL6(R)/JAK/STAT3) known to be activated in hematologic malignancies 148 .This raises the intriguing possibility that cells from one reactivation trajectory and viral replication mode (fully lytic cells) might reinforce the survival and proliferation of tumor cells resulting from an alternative response (abortive, latently infected) through paracrine mechanisms.In addition to its escape from host shutoff 110 , IL-6 autocrine support for latent EBV + B cell proliferation and its depletion in BZLF1-and BRLF1-deificient tumors in murine models of EBV-driven lymphoproliferative disease are especially noteworthy 53,149,150 .A similar effect has been observed during infection with KSHV, which encodes a viral IL-6 homolog.Thus, the developmental pluripotency profiles and responses of lytic cell subsets may be associated with cellular DNA damage responses that have inadvertent pathogenic effects in EBV + tumors.Notably, cytokine production by EBV-infected tumor cells (including abortive lytic cells) has also been proposed to support oncogenesis through microenvironment conditioning, polarization of tumor infiltrating lymphocytes, and evasion of Tcell surveillance 49,50 .
In summary, our findings support a model of differential response trajectories to EBV lytic induction.The first determinant in this model is initial cell state, where ground-state STAT3 and MYC expression and activity predict a 'high-resistance', low-probability path to full reactivation.
Conversely, cells with globally reduced transcription and reduced expression of STAT3 (and MYC) at the time of lytic reactivation traverse a 'low-resistance' path with high probability of complete reactivation.These data have potentially important clinical implications, as they suggest that quiescent EBV + tumor cells may be more sensitive to lytic induction therapies.However, a critical second fate determinant that manifests in lytic cells may complicate this pursuit.To this point, our scRNA-seq and RNA Flow-FISH results are consistent with the previously identified role of lytic cycle induction in tumorigenesis 46,47,53 .Most cells that undergo full reactivation and new virion release are likely to die.However, some lytic cells undergo profound reprogramming to plastic CSC-like states that may promote malignancy through multiple mechanisms, even independent of their own survival.For example, we found transcript-level evidence that lytic cells could reinforce viral latency and survival of abortive or refractory cells via IL6/JAK2/STAT3 signaling.Additional studies are necessary to explore, dissect, and therapeutically perturb the IL-6/JAK/STAT3 pathway in EBV + lymphomas.Given these findings, subsequent examinations of the epigenetic consequences of early EBV reactivation at high resolution should be prioritized, and the possibility of double-edged consequences of oncolytic therapies should be specifically examined in detail.Future single-cell approaches should interrogate the frequency of viable abortive lytic cells 151 and the particular changes in chromatin accessibility as well as other epigenetic features of this phenotype.Similar experimental approaches should be applied to study clinical EBV + tumor specimens to understand oncogenic correlates of lytic reactivation in situ.

MATERIALS AND METHODS
Cell lines, culture, and treatments P3HR1-ZHT cells (derived from the Type 2 EBV + [P3 strain] Jijoye eBL line) and B958-ZHT (a marmoset lymphoblastoid cell line transformed with Type 1 EBV [B95-8 strain]) were used in this study.Each cell line was cultured at 37°C with 5% CO2 in RPMI + 10% FBS (R10) media (Gibco RPMI 1640, ThermoFisher).To induce lytic gene expression, 4x10 5 cells/mL for a given cell line in log-phase growth were treated with 25 nM, 50 nM, or 100 nM 4-hydroxytamoxifen (4HT) in methanol (4HT, Millipore Sigma).Phosphonoacetic acid (PAA, 1 M) was included in parallel with lytic induction treatments to inhibit viral DNA synthesis and prevent complete reactivation in separate experimental groups (i.e., abortive lytic replication).Control groups were prepared via treatment with 0.1% DMSO (and DMSO + PAA).All treatments for flow cytometry and RNA Flow-FISH experiments described below were performed in triplicate (technical replicates) in 6-, 12-, or 24-well culture plates.
Cells were then incubated with target probes for 2h in a 40°C water bath.Cells were washed and stored overnight at 4°C and then incubated with a Pre-amplification buffer for 1.5h in a 40°C water bath followed by a 1.5h incubation in amplification buffer.Cells were then incubated in label probes for 1 hour in a 40°C water bath, washed with FACS buffer and subsequently analyzed on a Cytek Aurora.Spectral flow unmixing was performed with SpectroFlo software and uniformly applied to all samples.Further analysis and gating was completed in FlowJo.
Single-cell sample and library preparation P3HR1-ZHT cells were plated at 4 x 10 5 cells/ mL in 5 mL R10 then treated with methanol (mock-0 h) or with 25 nM 4HT (4-hydroxytamoxifen).The cells incubated in 4HT for 72, 48, and 24 hours then all cells were harvested for library preparation at the same time.The viabilities of the 0, 24, 48, and 72 h samples at time of collection were approximately 90%, 80%, 75%, and 75%, respectively.Harvested cells were resuspended at the recommended concentration to collect approximately 10,000 cells per sample during GEM generation.Single-cell transcriptomes from all four samples were captured and reverse transcribed into cDNA libraries using the 10x Genomics Chromium Next GEM Single Cell 3' gene expression kit with v3.1 chemistry and Chromium microfluidic controller according to recommended protocols (10x Genomics, Pleasanton, CA).All cDNA gene expression libraries were pooled for sequencing.
Sequencing, read alignment, and QC Pooled single-cell libraries were sequenced across two lanes of an S2 flow cell on a NovaSeq6000 (Illumina, San Diego, CA) with 50 bp paired-end reads at a target sequencing depth of 50,000 reads per cell.Output base calls (.bcl) were assembled into sampledemultiplexed reads (.fastq) using cellranger mkfastq with default settings (10x Genomics, Pleasanton, CA).Reads were mapped to a concatenated reference genome package (hg38 + NC_009334 [type 2 EBV]; prepared via cellranger mkref) to generate single-cell expression matrices by running cellranger count (10x Genomics, Pleasanton, CA).Cellranger output files (genes.tsv,barcodes.tsv,matrix.mtx)were used to create Seurat data objects in R [152][153][154] , which were subsequently pre-processed using QC filters.Cells and features were included if they met the following criteria: feature (gene) expression in a minimum of three cells; mitochondrial genes accounting for < 25% of all transcripts; a minimum of 200 unique expressed genes; < 65,000 total transcripts to exclude non-singlets.The elevated mitochondrial transcript and total transcript cutoffs relative to those used for resting PBMC samples 155 were chosen because of the highly proliferative nature of the P3HR1 cell line, the expectation of apoptosis as one outcome to lytic reactivation, and the implementation of viability enrichment prior to library preparation described above.A total of 26,728 cells across the timecourse passed all QC filters (nuntreated = 10,196; n24h = 7,905; n48h = 5,841; n72h = 3146).Data pre-processing, dropout imputation, analysis, and visualization A complete list of loaded packages and versions (RStudio sessionInfo() output) is provided as a supplementary file.Single-cell expression data were analyzed and visualized with R (v4.0.5) / RStudio (v2022.07.1+554) using Seurat v4.1.0.Data from each timepoint were analyzed separated and merged into a single object to support time-resolved analysis.Raw count data were normalized and scaled prior to feature identification (NormalizeData and ScaleData functions).Cell cycle scores and phases were assigned based on annotated gene sets provided in the Seurat package (CellCycleScoring function).Expression data were dimensionally reduced using principal component analysis of identified variable features (RunPCA), and the first 30 principal components were used for subsequent UMAP dimensional reduction (FindNeighbors, RunUMAP).Cell clusters were identified at multiple resolutions for phenotype identification and comparative analysis (FindClusters).
Biological zero-preserving imputation was applied to correct technical read dropout using adaptive low-rank approximation (ALRA) of the RNA count matrix 156 .Data presented throughout this study was generated from imputed read data.Differential gene expression analysis of the merged timecourse RNA and imputed (ALRA) assays was performed at multiple clustering resolutions.Outputs from this analysis are provided as supplementary files.Single-cell gene expression, co-expression, and cluster-averaged expression were visualized with Seurat functions (e.g., DimPlot, FeaturePlot, FeatureScatter, VlnPlot, DotPlot, DoHeatmap).Additional visualization of multi-gene co-expression was generated with the UpSetR package 157 .

Pseudotime analysis
Pseudotime trajectories were calculated for day 0 and merged timecourse datasets using Monocle3 111,158 .Briefly, Seurat objects were adapted as cell dataset objects and used to learn and order cells along pseudotime graphs anchored at manually determined root cells.Calculated pseudotime values were added as a feature to original Seurat objects and used for subsequent gene expression analyses.Pseudotime-gene correlation was plotted and fit via smoothing splines to visualize expression dynamics across clusters (cell phenotypes).

Gene ontology and gene regulatory network analyses
Low and high-resolution cluster gene ontology (GO) enrichment analysis for biological processes was performed using the enrichGO function in clusterProfiler 90 .Statistically significant enrichment results were visualized using the barplot, pairwise_termism, and emapplot functions.
Cluster-level gene regulatory network (GRN) inference of transcription factor activities was conducted using CollecTRI in the R package decoupleR 159,160 .

Statistical analyses
Raw and adjusted p values (Bonferroni correction) were calculated and provided for all identified differentially expressed genes from scRNA-seq data (see supplementary tables 1-4).
For conventional flow cytometry and RNA flow-FISH experiments, statistically significant differences between treatment groups were determined via two-tailed Welch's t test (n = 3 replicates per condition).
Human cell lines used in this study were not accompanied with HIPAA identifiers or PHI.All experiments were thus categorized as non-human subjects research and approved by a Duke University IRB (eIRB #Pro00006262).(F) Biological process GO analysis for genes upregulated in lytic subcluster E3 versus E2.

MAIN FIGURE LEGENDS
(G) Biological process GO analysis for genes upregulated in lytic subcluster E1 versus E2 and E3.
strongly anti-correlated with BZLF1 induction (Fig 2E, bottom left panel).STAT3 expression, which has been previously shown to be upregulated in cells refractory to lytic reactivation 59 , was likewise anti-correlated with expression of BZLF1 (Fig. 2E, bottom middle panel).STAT3 and MYC expression were positively correlated and highest in unstimulated (cluster A) and abortive (cluster C) cells (Fig. 2E, bottom right panel).Prediction of transcription factor activities based on gene regulatory network (GRN) enrichment likewise identified enhanced STAT3 (and NF-B) target expression in cluster C (Fig. S5).RNA Flow-FISH detection of BZLF1 and MYC validated scRNA-seq data and provided additional insight with respect to partial versus complete reactivation indicated by expression of the late lytic gene BLLF1 (Fig. 2F).Specifically, 4HT treatment induced significant increases in BZLF1 + and BLLF1 + cells and a concomitant decrease in MYC + cells relative to DMSO-treated controls (Fig. 2F, top and middle panels).Moreover, the majority of BLLF1 + cells were BZLF1 + /MYC -(Fig.2F, bottom panel).
are preferable over purely chronologic sampling for studying biological state transitions due to initial state variability and asynchronous responses to infection among individual cells 26 .Root cells (pseudotime=0) for the reactivation trajectory graph were chosen within clusters A and B since both of these phenotypes were represented by unstimulated cells (Fig. 2G, top panel).As shown by per cell viral fractions of captured mRNA transcripts, reactivation generally progresses in pseudotime, with limited viral expression in abortive cells at intermediate coordinates and high viral expression in fully lytic cells in late pseudotime (Fig. 2G, bottom panel).Notably, trajectories from both clusters A and B pass through incomplete reactivation states (C and D, respectively) before convening within the lytic phenotype (cluster E) at late pseudotime (Fig. 2G).Collectively, cluster-resolved expression, MYC and STAT3 profiles, and pseudotime trajectory analysis enabled us to construct a state model for lytic reactivation in the P3HR1-ZHT system (Fig. 2H).Unstimulated cells express elevated MYC and STAT3 and may undergo abortive reactivation in response to 4HT in which BZLF1 expression is minimal while MYC and STAT3 levels are largely maintained.Alternatively, cells may proceed to lytic reactivation, during which both MYC and STAT3 expression are severely diminished.Although clusters C and E were connected by a bridge of cells in the UMAP embedding, we cannot make definitive conclusions from these data alone regarding possible interconversion between abortive and lytic states.While global mRNA levels decrease along the transition from A to E consistent with host shutoff, the trajectory from cluster B (unstimulated) through D (intermediate) toward E (lytic) was characterized by relative increases in total and unique host and viral mRNA content.However, reduced MYC expression was also observed along the B to E trajectory.Overall, these results indicated that heterogeneity in unstimulated cells and differential responses to BZLF1 induction each contributed to the generation of distinct cell states during lytic reactivation.Analysis of gene expression along state-specific pseudotime trajectories captured these distinct biological response coordinates (Fig. 2I).For example, trajectories starting from clusters A and B both exhibited upregulated BZLF1 and net MYC reduction.However, STAT3 expression was consistently low across B, D, and E while STAT3 increased from A to C and decreased from A to E. Likewise, dynamic expression of GC B cell and NF-B signature genes along A→(C)→E were not observed along B→D→E.

Figure 1 .
Figure 1.EBV lytic reactivation in the P3HR1-ZHT Burkitt Lymphoma line at single-cell resolution.(A) Schematic of 4HT-inducible BZFL1 (Zta) expression initiating lytic reactivation in the Burkitt Lymphoma-derived P3HR1-ZHT cell line.(B) Flow cytometry validation 24 h after 4HT-induced lytic reactivation and inhibition of complete reactivation by phosphonoacetic acid (PAA) in P3HR1-ZHT.Cellular expression of the viral glycoprotein gp350 (encoded by the late lytic gene BLLF1) serves as a proxy for successful

Figure 2 .
Figure 2. P3HR1-ZHT phenotypic heterogeneity and response trajectories during lytic induction.(A) P3HR1-ZHT cell clusters identified in merged timecourse data via unsupervised methods.(B) Cluster composition of cells from individual timepoints.Cluster colors are coded as in 2A.(C) QC feature distributions by cluster.The total number of mapped reads per cell is given by nCount_RNA.The number of unique RNA features (i.e., genes, lncRNAs) per cell is given by nFeature_RNA.The viral fraction of mapped reads per cell (viral.pct)and mitochondrial transcript fractions (percent.mt)were calculated using the PercentageFeatureSet() function in Seurat 152,154 .(D) Differential RNA expression by cluster.Sequences are annotated by their known biological roles and functions derived from gene ontology (GO) analysis and primary literature.Dot size represents the percentage of cells in each cluster that express a given gene and color encodes average expression across the cluster.(E) UMAP expression profiles (top row) and pairwise correlation plots (bottom row, Pearson R) for BZLF1, MYC, and STAT3.Correlation plots depict individual cells colored by cluster.(F) RNA Flow-FISH validation of reduced MYC expression in BZLF1 + BLLF1 + cells (top panel).Asterisks in the middle panel bar plot denote significantly reduced frequency of MYC + P3HR1-ZHT cells and increased frequencies of BZLF1 + and BLLF1 + cells after 4HT treatment (n=3 per

Figure 3 .
Figure 3. Validation of an abortive response with elevated NF-B activity distinct from full lytic reactivation.(A) Identification of CD38, BCL2A1, and BLLF1 as respective biomarkers for unstimulated, abortive, and lytic P3HR1-ZHT cells.(B) Co-detection of BZLF1 and NF-B pathway transcriptional targets in abortive cells (co-positive cells in red) by timepoint.(C) RNA Flow-FISH validation of full (BLLF1 + ) and abortive (BCL2A1 + ) reactivation as orthogonal responses.DMSO control-treated cells are predominantly CD38 + and exhibit minimal spontaneously lytic (full or abortive) cells (top panel).4HT treatment induces distinct full lytic and abortive subsets (middle panels).Inhibition of viral DNA synthesis with PAA blocks full lytic reactivation and increases the frequency of BCL2A1 + abortive cells (bottom panels).Colored circles denote predicted corresponding model states defined from scRNA-seq.(D) Frequencies of CD38 + and BCL2A1 + cells presented in 3C by treatment condition at 24 h and 48 h.Asterisks denote significantly decreased frequencies of CD38 + cells and increased frequencies of BCL2A1 + cells upon 4HT and 4HT+PAA treatment versus respective control treatments (n=3 per condition; two-tailed Welch's t-test; **p<0.01).(E) EBV LMP-1, which encodes a potent activator of NF-B signaling, is expressed in late lytic cells (left panel) but not associated with abortive cells that exhibit upregulated NF-B transcriptomic signature including BCL2A1 (right panel, Pearson R=-0.06).(F)Flow cytometry analysis of protein biomarkers of full lytic reactivation (gp350) and NF-B activity (ICAM1).Consistent with mRNA measurements, separate gp350 + and ICAM1 + populations are induced following 4HT treatment.Co-treatment with PAA reduces gp350 + cell frequency and increases ICAM1 + fractions.IKK inhibitor co-treatment reduces ICAM1 + cell frequency but does not substantially affect gp350 + cell frequency.

Figure 4 .
Figure 4. Cancer-associated cellular plasticity and self-renewal signature identification in EBV lytic cell subsets.(A) Total mapped RNA reads per cell (top panel) versus total unique genes expressed across each cluster (bottom panel).(B) Unsupervised identification of high-resolution subclusters across P3HR1-ZHT timecourse scRNA-seq data.(C) Differentially expressed genes upregulated in lytic subclusters (E1, E2, and E3).Genes were identified by comparing each subclusters versus all others, summarized by gene ontology methods, cross-referenced against primary literature, and curated by biological annotation.(D) Co-expression of BZLF1 and genes associated with cellular pluripotency and cancer stemness (SOX2, ALDH1A1) in single cells (co-positive cells in red) by timepoint.

Figure 6 .
Figure 6.Lytic subset reprogramming and host shutoff escape signatures are conserved in B958-ZHT lymphoblastoid cells.(A) UMAP representation of scRNA-seq data from the inducible lytic marmoset lymphoblastoid cell line B958-ZHT before (Unstim) and 24 h after 4HT treatment.(B) Mapping of cell subclusters defined from P3HR1-ZHT analyses to B958-ZHT scRNA-seq data via transfer anchor integration (left panel).Subcluster composition is presented for unstimulated and 4HT-treated cell libraries (right panel).(C)Viral expression module (IE, early, late) and mapped lytic subcluster scores in timecourse merged B958-ZHT data.Of note, the assigned subclusters in 6B represent qualitative classifications based on maximum annotation signature scores for each cell.Accordingly, a given cell may score highly for more than one related signature while being assigned to a single classification.The underlying quantitative signature scores for E1, E2, and E3 presented here thus reflect a lytic phenotypic continuum rather than purely discrete states.

Figure S13 .
Figure S13.Independent scRNA-seq replicate validation of key heterogeneous responses in P3HR1-ZHT cells.(A) Overview of P3HR1-ZHT replicate experiment treatments (methanol control and 4HT) and identified clusters.(B) UMAP visualization of global QC metrics (top row), differential abortive and lytic responses correlated with STAT3 and MYC levels (2 nd and 3 rd rows), and upregulated pluripotency signature in lytic cell subsets (4 th and 5 th rows).