Virus-human chromatin interactions reorganise 3D genome and hijack KDM5B for promoting metastasis in nasopharyngeal carcinoma

Abstract

Epstein-Barr virus, the first identified human DNA tumour virus, is detectable in more than 90% of nasopharyngeal carcinoma patients in endemic regions. The 3D chromosome conformation analysis reveals that virus‒host chromatin interactions induce the spatial reorganisation of loops and compartments, resulting in "enhancer infestation" and switch of "H3K27 bivalency" at EBV-interacting regions. Through this mechanism, EBV hijacks KDM5B, a gatekeeper of genome stability, and its binding targets, leading to aberrant activation of an EBVIR-enhancer-KDM5B signature. Cancer cells with this signature present increased MYC activation, DNA damage responses, and epigenetic plasticity of epithelial-immune dual features with metastatic potential. Our multicentre multiomics study confirms that this signature is the prerequisite for chromosome instability and can be utilised as a risk factor for distant metastasis. This study highlights a mechanism in which latent viral episomes can alter the host high-order epigenotype, promoting transcriptional rewiring and metastasis in NPC.

Fig. 1. EBV-associated 3D genome alterations.

A EBV-associated compartment features in NPC cells. EBV+ NPC cell lines, EBV+ PDXs and an EBV+ clinical tissue sample were clustered separately from NPE and EBV− NPC cell lines. Differential compartment analysis based on principal component 1 (PC1) was carried out with the EBV+ NPC versus EBV− NPC cell lines. A total of 279 regions switched from compartment A to B, and 465 regions switched from compartment B to A. B Unsupervised hierarchical clustering of B-to-A and A-to-B compartment switching in NPC cells. C Hierarchical clustering of EBV-associated compartment switching from B to A and A to B in NPE, EBV+ NPC and EBV− NPC cells visualised on the basis of PC1 value (a 50 kbp window size) through the IGV browser. A positive (negative) PC1 value indicates compartment A (compartment B). D Redistribution of AB compartments in two paired EBV+ and EBV− NPC cell lines at chromosome 4. A Pearson correlation matrix was normalised by the Knight‒Ruiz (KR) method. E Visualisation of the 3D genome of chromosome 4 in matched NPC cell line pairs. A selective genomic region showed that compartments A and B were clearly separated from each other without EBV infection, while EBV-infected cells have reorganised compartments with B-to-A transitions. The red colour indicates A (active) compartments, and the blue colour indicates B (inactive) compartments. F Reduction in chromatin loops in EBV-infected NPC cells. Aggregate peak analysis (APA) plots based on the loop interactions called from Moustache. A total of 66,100 regions were identified as loops in NPE, EBV+ NPC and EBV− NPC cells and are plotted in 20 × 20 square pixels. APA plots of paired EBV+ and EBV− NPC cells. P2LL denotes the peak-to-lower-left of the interaction matrices. The higher P2LL score indicates stronger chromatin interactions. G Reduction in chromatin loops in EBV+ cells. P2LL denotes peak-to-lower-left. Two-tailed Student’s t-test was performed EBV+ (n = 7) versus EBV− (n = 3) cells (p-value = 0.018), *p-value ≤ 0.05. The error bar indicates mean ± standard deviation (SD) of the data. NPE (n = 2).

Fig. 2. Identification of three EBV-interacting regions (EBVIRs) associated with enhancer activation and inactivation.

A EBV-interacting events across human chromosomes in EBV+ NPC cells, PDX, and EBVaGC (YCCEL1) cells. The frequency of EBV tethering to the human genome was calculated with a 50 kb window size and visualised as a heatmap on a Circos plot. B Consistent B-to-A compartment switching was found at EBVIRs. The odds ratio (OR) was used to calculate the association of EBVIRs with six types of compartment switching (A-to-B, B-to-A, B–B/A–A enhance, B–B/A–A reduce) in EBV+ NPC cells, PDXs and EBV-infected NPC cells. The NPE cell control was used as a reference. One-tailed Fisher’s Exact Test has been performed to test if EBVIRs are associated with subsets of compartment switches, * p-value ≤ 0.05. The total number of events (C17 = 110658, C666 = 110648, NPC43 = 110660, HK1EBV = 110520 and Xeno32 = 110490) has been examined for the significance. The error bar indicates a 95% confidence interval (CI) of the data. Compartment A-to-B (B-to-A) switching is indicated by a change in the PC1 value from positive (negative) to negative (positive). The enhanced and reduced compartments indicate the magnitude of the change in the PC1 value within the same compartment. The arrows indicate that more EBV was consistently situated at B-to-A compartments than at A-to-B compartments across all samples. C EBVIRs are associated with enhancer activation and enhancer inactivation. Euclidean clustering resulted in two clusters of EBVIRs, which are named enhancer activation and enhancer inactivation EBVIRs. The CUT&RUN-seq signals were normalised in Reads Per Kilobase per Million mapped reads (RPKM). D Discovery of the three most frequent contact regions in the EBV genome. The EBV genome is annotated in the left panel. A two-tailed Z-score was calculated with data from individual samples, Z-score > 1.645, 90% confidence level (p-value ≤ 0.1); number of regions used for statistical test n_C17 = 1232, n_C666 = 18,405, n_NPC43 = 7063, n_HK1EBV = 25,796, n_C15 = 22,348, n_xeno23 = 1,11,888, n_xeno32 = 6464, and n_NPC113 = 31,400. Three regions were significantly enriched across C17, C666, NPC43, Xeno32 and HK1EBV cells. The RPKM normalised chromatin accessibility (ATAC), CTCF binding, and H3K4me3 histone modification (CUT&RUN-seq) with a 10 bp window size and summarised signal tracks along the EBV genome in EBV+ NPC cells are shown.

Fig. 3. Multiomics integrated data analysis identified KDM5B and its targets related to EBV infection.

A The expression levels of LMP1, RPMS1, EBV latency genes and KDM5B in epithelial cells were determined by scRNA-seq. LMP1 and RPMS1 were used to stratify the cells into EBV-high and EBV-low clusters. The normalised KDM5B expression is examined and linked to high expression of EBV latent genes. B Venn diagram following the multiomics integrated data analysis. Upregulated genes in EBV-high clusters were identified from differential gene expression analysis by comparing the EBV-high (n = 5335) and EBV-low (n = 1876) clusters in scRNA-seq data. A two-tailed Mann–Whitney U-test has been performed, with an adjusted p-value ≤ 0.05. C Putative transcription factors involved in regulating upregulated genes in EBVIRs through enhancer activation. TFs were ranked on the basis of their combined score (log₁₀ (p-value) from Fisher’s exact test × the the Z score for deviation from the expected rank). D Integrated IGV browser signal tracks at the KDM5B locus. Histone modifications, including H3K27ac and H3K27me3, CTCF binding signals at the same scale, together with Omni-C chromatin interactions for C666 and NP69 cells, and EBVIRs are included. The orange frame highlights EBVIRs at the KDM5B promoter region, with chromatin contacts forming CTCF chromatin loops in nearby regions. E neoTAD formation at the KDM5B gene locus. Omni-C data were used to visualise TAD formation, and data were compared between EBV+ NPC and NPE and EBV- NPC cells using Juicebox. The plot was further normalised by sequencing coverage visualising at 100 kbp resolution. If the observed deduce control data is greater than zero (red), it indicates more chromatin contacts in C666, and vice versa. F EBVIR–enhancer–KDM5B target expression determined by scRNA-seq. The sum of normalised expression levels in EBVIR–enhancer–KDM5B signature has been calculated in EBV-high and EBV-low cells (n_cells = 7211). The expression levels in the EBV-high (n = 5335) and EBV-low (n = 1876) clusters were statistically evaluated by a two-tailed unpaired t-test, ****p-value ≤ 0.001. G Spatial transcriptomes of KDM5B and EBVIR–enhancer–KDM5B target expression levels in the EBV-high and EBV-low regions in a clinical sample. The average expression levels of genes in EBVIR–enhancer–KDM5B signature were calculated for each circle (n_spots = 264), visualised by image and quantified by the violin plots. The expression levels in the EBV-high (n = 63) and EBV-low (n = 201) clusters were statistically evaluated by two-tailed unpaired t-test (p-value = 0.04), *p-value ≤ 0.05.

Fig. 4. Single-cell spatial transcriptome analysis revealed subclusters of epithelial NPC cells with aberrant upregulation of the EBVIR–enhancer–KDM5B signature.

A Single-cell spatial transcriptome of a metastatic NPC patient sample (n_FOV = 2883) illustrated by the AtoMx platform. Two FFPE slides from one MET NPC patient and one non-MET NPC patient have been used to conduct single-cell spatial transcriptomes. PanCK, CD45 and DNA staining were performed. The epithelial cell clusters (Epi-1, Epi-2 and Epi-3) are shown, and select genes, including PankCK, CD45, RPMS1, CIITA and EPCAM, were visualised (highest expression: yellow, moderate expression: green, and lowest expression: blue). Neighbourhood analysis using Jaccard distance stratified cells into three niches in this selected region of interest. B Pathway activity of EBVIR–enhancer–KDM5B targets in different cell types. The GSVA score for EBVIR–enhancer–KDM5B pathway activity was calculated for each cluster of cells. Red bars show the most significant pathway activity in epithelial clusters 1–3. C Expression levels for selected EBV and human genes in different human cell types. RPMS1 and LMP1 were used as indicators for EBV-high cancer cells, as discovered by scRNA-seq data. Expression levels of the genes upregulated in epithelial clusters 1–3 were selected in the plot. Red indicates a higher expression level, and blue indicates a lower expression level. D Enrichment scores of NPC-related cancer hallmarks from the Human Molecular Signatures Database (MSigDB) for different cell types. The selective hallmarks with statistical enrichment in epithelial clusters 1–3 were visualised in the plot. Red indicates a positive enrichment score, and blue indicates a negative enrichment score. E Ectopic expression of NPC-related genes in RPMS1-expressing EBV+ NPC cells determined via scRNA-seq. EBV-high clusters were identified via scRNA-seq were subclustered (n_cells = 1696) via UMAP into RPMS1-high (n = 1112) and LMP1-high (or RPMS1-low) (n = 584) groups on the basis of expression of the EBV latency genes RPMS1 and LMP1, respectively. Genes commonly found via both scRNA-seq and single-cell spatial transcriptomics were selected. Cells with high expression are coloured blue, and those with low expression are coloured yellow. F Epithelial cell proportions in MET and non-MET patients. The proportion of cells in each epithelial cell cluster was examined for MET (n_cells = 89,743) and non-MET (n_cells = 14,104) patients. Epi-1/−2/−3 clusters account for about 25% of the total epithelial cells in MET patients. N denotes normal epithelial cells, and T denotes tumour epithelial cells.

Fig. 5. KDM5B is a risk factor for MET, and its targets are involved in chromosomal instability and DNA damage response.

A GSEA of EBV tethering-mediated enhancer activation of KDM5B targets in MET patients. Genes were ranked based on differential expression analysis using the bulk RNA-seq data from MET patients and non-MET patients in the HK and GZ cohorts to calculate the enrichment score of EBVIR–enhancer–KDM5B signature. B Human whole-transcriptome atlas (WTA) spatial analysis of 26 clinical specimens (MET: n = 16, nonMET n = 8). The yellow frame indicates the region of interest (ROI). A heatmap shows differentially expressed genes (DEGs) between the MET and non-MET samples. A two-tailed Mann–Whitney U-test with log2FC > 0.6 or < −0.6 and adjusted p-values ≤ 0.05 estimated by permutation test identified DEGs. The DEGs in both cohorts were evaluated. C Putative transcription factors associated with metastasis-related genes. The putative transcription factors among metastasis-associated genes were evaluated in the HK (n_genes = 236) and GZ (n_genes  = 291) cohorts. A scatter plot shows the odds ratio against −log₁₀ (p-value). The combined score was calculated as log₁₀ (p-value) from one-tailed Fisher’s exact test × the Z score for deviation from the expected rank. D Functional annotation of genes expressed in the MET patients. A metastasis-associated gene set (n = 82) common to both cohorts was identified. The hypergeometric test for Reactome signalling pathways and GO terms for biological processes with Benjamini–Hochberg correction has been used to examine the significance. All significant terms had an adjusted p-value ≤ 0.05, and the font size is based on the odds ratio. E Chromosomal instability (CIN) scores for EBV tethering-mediated enhancer activation of KDM5B targets in both the HK and GZ cohorts. A cut-off of 0.1 was assigned and used to identify CIN-high (HK: n = 21 and GZ: n = 30) and CIN-low patients (HK: n = 13 and GZ: n = 29). Student’s t-test was used to test the significance of the difference between the two groups. A two-tailed unpaired t-test was performed to examine the significance between CIN-high and CIN-low patients in both cohorts. The bar plot with error bars indicates mean ± standard deviation (SD) of the data. CIN-high and CIN-low patients have shown a significant difference in HK and GZ cohorts (p-value = 0.04 and p-value < 0.001, respectively).

Fig. 6. Univariate and multivariate logistic regression analyses revealed that KDM5B pathway activity is an independent risk predictor for NPC metastasis.

A The HK (n = 77) and GZ (n = 100) cohorts were included in the analysis. Univariate logistic regression analyses were performed to assess the association between individual variables and the binary outcome. KDM5B pathway activity; OR = 3.43, 95% CI (2.54, 4.63) with p-value = 0.0003, CIN Score; OR = 2.26, 95% CI (1.68, 3.05) with p-value = 0.0076, C-reactive protein (CRP); OR = 3.54, 95% CI (2.31, 5.44) with p-value < 0.0001, Neutrophil–lymphocyte ratio (NLR); OR = 3.57, 95% CI (2.40, 5.31) with p-value = 0.3173, Fibrinogen (FBG); OR = 5.11, 95% CI (3.10, 8.42) with p-value < 0.0001, Baseline EBV DNA copy; OR = 2.37, 95% CI (1.50, 3.75) with p-value < 0.0001, EBV DNA copy after treatment; OR = 7.96, 95% CI (4.75, 13.34) with p-value < 0.0001, Stage III; OR = 1.42, 95% CI (0.87, 2.32) with p-value < 0.0001, Stage IV; OR = 4.5, 95% CI (2.69, 7.54) with p-value < 0.0001, Age; OR = 0.8, 95% CI (0.59, 1.09) with p-value < 0.0001, Gender; OR = 1.72, 95% CI (1.21, 2.45) with p-value < 0.0001. Multivariable logistic regression analysis was conducted to evaluate the independent associations between selected clinical and molecular variables and the binary outcome. Group 1 is the reference group for analysis. KDM5B pathway activity; OR = 4.65, 95% CI (2.98, 7.25) with p-value < 0.0001, EBV DNA copy baseline; OR = 2.62, 95% CI (1.47, 4.67) with p-value = 0.001, EBV DNA copy after treatment; OR = 7.97, 95% CI (4.47, 14.21) with p-value < 0.0001, Stage III; OR = 2.26, 95% CI (1.11, 4.62) with p-value < 0.0001, Stage IV; OR = 6.12, 95% CI (2.83, 13.23) with p-value < 0.0001. The plots show odds ratios (ORs) in the centre with 95% confidence interval (CI) of the data. The size of the red square denotes the sample size.

Fig. 7. EBV hijacks host factor—KDM5B to provide host cancer cell survival.

A Generating KDM5B wild-type (WT) (n = 3) and knockdown (KD) (n = 4) cells in C666. A western blot result for KDM5B at ~170 kDa and α-tubulin control at ~55 kDa demonstrate the success of establishing KDM5B–KD clones. B Quantification of cells in KDM5B–WT and –KD cells. An equal number of cells were cultured for 7 days. A cell count chamber was used to calculate the number of cells in WT and KD cells. Biological replicates have been performed for WT clones (n = 3) and KD clones (n = 4). The number of cells was normalised to WT. An unpaired t-test was used to test the significance, showing KD cells have approximately 60% cell growth reduction with p-value = 0.047. * p-value ≤ 0.05. The bar plot with error bars indicates mean ± standard deviation (SD) of the data. C MTT assay results for 7-day treatment of KDM5B inhibitor (GSK-467) in C666, NPC43 and YCCEL1 cells. The experiments were performed with technical replicates (n = 5) for both DMSO control and KDM5B treatment groups. The KDM5B inhibitor results were normalised to their matched DMSO control. A two-way ANOVA test was performed with Bonferroni correction, ****adjusted p-value ≤ 0.0001. A concentration of 75 µM of KDM5B inhibitor treatment showed significant suppression of cancer cell growth in EBV+ NPC and EBVaGC cells. The bar plot with error bars indicates mean ± standard deviation (SD) of the data. D CRISPR screening revealed that the EBVIR–enhancer–KDM5B signature suppresses tumour cell growth. A genome-wide loss-of-function (LOF) analysis employing the CRISPR screening approach was conducted. CRISPR libraries were transduced into the cells, and cells were subsequently selected by puromycin. The sgRNA sequences were amplified and sequenced by the Illumina platform. The LOF of a gene conferring a growth advantage, characterised by a higher abundance of sgRNAs, was quantified based on the sequencing. Conversely, the LOF of a gene lacking a growth advantage was determined by sequencing with a significant reduction of the sequencing signals (Created in BioRender. Dai, W. (2025) https://BioRender.com/7e447e1). The GSEA result showed significant negative enrichment of EBVIR–enhancer–KDM5B signature (n = 64), NES = −1.36, FDR q-value = 0.039, suggesting LOF of EBVIR–enhancer–KDM5B signature strongly inhibited cancer cell growth.

Fig. 8. EBV hijacks the host factor—KDM5B to maintain its genome copies.

A KDM5B regulates genes in response to the EBV life cycle and cell death. Bulk RNA-seq has been conducted for KDM5B–WT (n = 3) and KDM5B–KD (n = 4) with biological replicate clones. Pre-ranked genes based on log₂ fold change of differentially expressed genes were subjected to GSEA. The downregulated genes in the EBVIR–enhancer–KDM5B signature were extracted for GO–BP analysis. A scatter plot shows the odds ratio against −log₁₀ (p-value). The combined score (log₁₀ (p-value) from one-tailed Fisher’s exact test × the Z score for deviation from the expected rank) is indicated by colour. The terms with p-value ≤ 0.05 are highlighted. B Profile of EBNA1, KDM5B binding, and H3K27ac on EBV genome visualised in IGV browser. The RPKM normalised EBNA1 and KDM5B binding signals show the involvement of EBNA1 and KDM5B at EBV–oriP–EBVIR with H3K27ac signals. C Quantification of EBV copy per cell in KDM5B–WT and –KD cells. qPCR was performed in triplicate for each clone (n_WT = 9 and n_KD = 12). A two-tailed unpaired t-test was used to test the significance, showing a significant reduction of EBV copies in KD cells with p-value = 0.004. ** p-value ≤ 0.01. The error bar indicates mean ± standard deviation (SD) of the data. D Pearson correlation plots between the KDM5B pathway and cellular EBV copies in the patient biopsies. The Pearson correlation tests showed a strong positive correlation of KDM5B pathway activity to EBV copies (n_patients = 22) with r = 0.633, y = 0.08536x − 0.6128, and p-value = 0.0016. Correlation is significant at the 0.05 level (two-tailed). The error bar indicates mean ± standard deviation (SD) of the data. E MTT assay results for 72 h of treatment of BRD inhibitor (JQ1) in C666, NPC43 and YCCEL1 cells. The experiments were performed with technical replicates (n = 5) for both drug treatment and DMSO control. A two-way ANNOVA test was performed with Bonferroni correction, *adjusted p-value < 0.05, **adjusted p-value < 0.01, ***adjusted p-value < 0.001 and ****adjusted p-value < 0.0001. The bar plot with error bars indicates mean ± standard deviation (SD) of the data.