Characterizing resistant cellular states in nasopharyngeal carcinoma during EBV lytic induction

Abstract

The pervasive occurrence of nasopharyngeal carcinoma (NPC) is intricately linked to Epstein-Barr virus (EBV) infection, making EBV and its associated pathways promising therapeutic targets for NPC and other EBV-related cancers. Lytic induction therapy, an emerging virus-targeted therapeutic strategy, capitalizes on the presence of EBV in tumor cells to specifically induce cytotoxicity against EBV-associated malignancies. Despite the expanding repertoire of compounds developed to induce EBV lytic reactivation, achieving universal induction across all infected cells remains elusive. The inherent heterogeneity of tumor cells likely contributes to this variability. In this study, we used the NPC43 cell line, an EBV-positive NPC in vitro model, and single-cell transcriptomics to characterize the diverse cellular responses to EBV lytic induction. Our longitudinal monitoring revealed a distinctive lytic induction non-responsive cellular state characterized by elevated expression of SOX2 and NTRK2. Cells in this state exhibit phenotypic similarities to cancer stem cells (CSCs), and we verified the roles of SOX2 and NTRK2 in manifesting these phenotypes. Our findings reveal a significant challenge for lytic induction therapy, as not all tumor cells are equally susceptible. These insights highlight the importance of combining lytic induction with therapies targeting CSC-like properties to enhance treatment efficacy for NPC and other EBV-associated cancers.

Fig. 1. Heterogeneous response of NPC43 cells to lytic induction treatment.

A UMAP representation of NPC43 cells at three distinct treatment time points, highlighting distinct clustering patterns. B UMAP plots illustrating the cellular states at different treatment time points. C Heatmap depicting the correlation between pseudo-bulks of different cellular states at various treatment time points, based on the top 1000 variable genes. Heatmap was hierarchically ordered. D Violin plots showing the activation levels of the EGFR, MAPK, PI3K, and NF-κB pathways in NPC43 cells. Pathway activity scores were inferred using PROGENy. Wilcoxon signed-rank test was performed. E Dot plot illustrating the expression levels of markers associated with different cellular states, demonstrating their conservation across treatment time points. F Bar plots showing the percentile distribution of cellular states at different treatment time points. Composition analysis was performed with scCODA, setting FDR to 0.05, and using cycling cells as the reference cell type. The proportions of NR cells, ECM-related cells (2), and Differentiated cells showed credible changes. G Scatter plots depicting the activity scores of Gene Ontologies (GOs) differentially activated in cellular states. H Scatter plot highlighting the most specific Differentially Expressed genes (DEs) of NR cells at T48. I Dot plot showcasing the expression levels of surface markers of NR cells in different cellular states, indicating conservation across different treatment time points.

Fig. 2. Prevalence of SOX2&NTRK2 in tumor samples from NPC patients.

A UMAP representation of different cell types in a published NPC dataset. B UMAPs displaying the expression levels of SOX2 and other surface markers in this published NPC dataset. C Dot plot illustrating the expression levels of SOX2 and other surface markers in malignant cells from various published NPC datasets. D Immunostaining images demonstrating the detection of SOX2 and TrkB in patient slides. E RNA-scope images revealing different abundances of SOX2 and NTRK2 in patients with long survival (over 150 weeks) or short survival (less than 50 weeks). SOX2 is labeled in red, NTRK2 is labeled in green, and nuclei are labeled in blue. F Box plots presenting the quantitative results of RNA-scope; the normalized expression levels of SOX2 and NTRK2 in long and short survival patient groups. G Scatter plot displaying the correlation between SOX2 and NTRK2 using a bulk RNA-seq dataset. Patients were separated into two groups based on SOX2/NTRK2 co-expression. H Kaplan–Meier progression-free survival curves for the two groups of patients with high SOX2/NTRK2 co-expression or low co-expression.

Fig. 3. Experimental validation of single-cell RNA-seq results.

A UMAP representation of the top 10% of cells expressing NTRK2, with a bar plot quantifying the proportional distribution of cellular states. B FACS detection of TrkB signaling events under different conditions, including an unstained control (unstained), UT, T24, and T48 NPC43 cells. C NTRK2 and SOX2 expressions determined by qRT-PCR. Cells from different conditions were sorted into NTRK2-high and NTRK2-low groups to extract RNA for qRT-PCR. Significant differences (p < 0.05) were observed between T24-high vs. T24-low and T48-high vs. T48-low groups. D Western blot showing SOX2 and TrkB abundance in TrkB-high and TrkB-low cells in UT and T48 conditions. GAPDH was used as a loading control. E qRT-PCR analysis of NTRK2, SOX2, and BZLF1 expression in FACS-sorted TrkB-high and TrkB-low populations (T48 condition). Significant differences (p < 0.05) were detected between T48-high and T48-low groups. F Western blot showing TrkB and Zta abundance in TrkB-high and TrkB-low cells in T48 condition. GAPDH was used as a loading control. G Relative copy numbers of BZLF1, EBER1 (2 EBV genes), and SOX2 determined by qPCR. Cells from different conditions were sorted into TrkB-high and TrkB-low groups to extract DNA for qPCR. H Expressions of NTRK2 and BZLF1 determined by qRT-PCR. Cells were sorted into TrkB-high and TrkB-low groups first (UT-high and UT-low). Then, sorted cells underwent lytic induction treatment for 48 h (UT-high-T48 and UT-low-T48). Significant differences (p < 0.05) were observed between UT-high vs. UT-low and UT-high-T48 vs. UT-low-T48 groups.

Fig. 4. Characterization of NR Cells and the function of SOX2 and NTRK2.

A GSEA results from DEs between NR cells and Keratinized cells in the UT dataset. B Violin plots showing the expression or activity score of genes (SOX2 and CD44) and gene sets (Stemness, p-EMT, EMV-IV, and Hypoxia). Wilcoxon signed-rank test was performed. C Images showing the size of tumorspheres. Cells with different TrkB abundance were seeded for this tumorsphere-forming assay. D Bar plots showing the numbers of tumorspheres (size over 10 µm). Cells with different TrkB abundance were seeded for this tumorsphere-forming assay. E Fluorescent immunostaining of tumorspheres showing the existence of the SOX2 and TrkB. F Images showing the size of tumorspheres. TrkB-high NPC43 or C17 cells undergo SOX2 or NTRK2 knockdown before seeding for the tumorsphere-forming assay. G Bar plots showing the numbers of tumorspheres (size over 10 µm). TrkB-high NPC43 or C17 cells undergo SOX2 or NTRK2 knockdown before seeding for tumorsphere-forming assay. H Heatmap showing the correlation of SOX2, CD44, and NTRK2 with EMT-related gene sets. Three gene sets defined from three different studies consistently show higher correlation with NTRK2. I Images showing the wound healing assay with different cells. NPC43 cells were knocked down with a shRNA vector. NTRK2-KD-sh shows slower migration compared to the control. J Violin plots showing ELF3 expression as a marker of differentiated cells in NPC43. K qPCR analysis showing increased expression of ELF3 and BZLF1 following NTRK2 knockdown and lytic induction treatment. Compared to NC-sh of UT, paired t-test p-values were less than 0.0383 for ELF3 and 0.00942 for BZLF1. L FACS analysis of Zta+ cells indicating an increased proportion of cells entering the lytic activation state after NTRK2 knockdown. A two-proportion z-test showed a p-value smaller than 10⁻⁵.

Fig. 5. Characterizing the regulation of SOX2 and NTRK2.

A Western blot showing TrkB, SOX2, E-cad (EMT-related marker), pERK (MAPK-related marker), and pS6 (PI3K-related marker) abundances in NPC43 cells undergoing SOX2-KD or NTRK2-KD. B Heatmap showing the activity of predicted TFs in different conditions. Predicted TFs were inferred from scRNA-seq data with the SENIC algorithm. C Rank plots showing the most specified TFs in NR cells at UT, T24, and T48. SOX2 or NTRK2-related TFs were highlighted. D Network showing how SOX2 is connected with NTRK2. TFs that are activated in NR cells were visualized in the network. E FACS plots showing the detection of SOX2 and TrkB in NP460 cells with SOX2 overexpression vector. F Western blot showing TrkB and SOX2 in NP460 cells with SOX2 overexpression.

Fig. 6. NTRK2 and SOX2 serve as robust features of NR cells across different NPC models.

A Heatmap displaying the top features correlated with NTRK2 expression, including genes and pathways. Only features recurrent across multiple datasets are shown in the heatmap. B UMAP of integrated scRNA-seq data from C666-1, with cells colored by treatment conditions and cellular states. Expression levels of NTRK2 and SOX2 are also visualized on the UMAP. C UMAP illustrating representative pathway terms associated with different cellular states. Stemness terms are enriched in NR-like cells; cytolysis and response to UV-B are enriched in cytolysis cells; unfolded protein response and inhibitory synaptic assembly are enriched in Prelytic cells. D Violin plot showing scores for EMT, Adhesion-related, hypoxia, and stress pathways. Only UT and T48 samples are included to compare NR-like cells. E UMAP shows the density of EBV gene detection and the expression of the lytic activation marker BZLF1.