Multiomics analysis reveals the genetic and epigenetic features of high-risk NK cell-type chronic active EBV infection

Abstract

Chronic active Epstein-Barr virus (EBV) infection (CAEBV) is an orphan disease characterized by the proliferation and infiltration of EBV-infected T/natural killer (NK) cells into multiple organs. Although CAEBV is a heterogeneous disease with diverse clinical courses, its pathogenesis remains poorly understood. In this study, we explored the molecular mechanisms underlying CAEBV by performing a comprehensive multiomics analysis, including genome, transcriptome, epigenome, and single-cell transcriptome and surface proteome analyses, of 65 patients with CAEBV. Methylation analysis identified 2 distinct subtypes of NK cell-type CAEBV based on the CpG island methylator phenotype (CIMP). In CIMP-positive CAEBV, regions associated with enhancer of zeste homolog 2 binding sites and histone H3 lysine 27 trimethylation exhibited increased DNA hypermethylation, resulting in downregulation of tumor suppressor and antiherpesvirus genes. CIMP-positive CAEBV had a particularly poor prognosis and displayed a "neoplastic" phenotype with a DNA methylation pattern similar to that of extranodal NK/T-cell lymphoma, a higher tumor mutation burden, and frequent copy number alterations. In addition, both in vitro and in vivo functional assays demonstrated that 5-azacytidine, a hypomethylating agent, was a potentially effective agent for high-risk CIMP-positive CAEBV. Finally, we established a method to effectively detect EBV-infected cells in single-cell analysis, suggesting that EBV-infected NK cells have tissue-resident properties and that innate and adaptive immunity to EBV is compromised in patients with CAEBV. The present findings provide insight into the complex molecular features of CAEBV and suggest potential molecular therapies.

Graphical abstract

Figure 1.

Classification of NK cell–type CAEBV based on the DNA methylation status. (A) Unsupervised hierarchical clustering of promoter-associated CpG island methylation profiles in NK cells from in-house samples (NK cell–type CAEBV, n = 32; IM, n = 4; non-EBV–IM-like, n = 2; healthy volunteers, n = 4) and deposited data (ENKTL, n = 15; normal NK, n = 21). Case/control, disease, and cluster are indicated in different colors (top). DNA methylation heat map of the 3000 most variable probes associated with the promoter (bottom). (B) Principal component analysis of the samples described in panel A using the 3000 most variable probes associated with the promoter. (C) Kaplan-Meier OS curves of 12 CIMP-positive and 20 CIMP-negative CAEBVs. The table below the plot indicates the number of patients at risk at each time point. (D) Forest plots of multivariate analyses of OS performed using the Cox proportional hazard model (hazard ratio ± 95% confidence interval [CI]). CIMP status and known clinical prognostic factors were included in the model. AIC, Akaike information criterion; PC1, principal component 1; PC2, principal component 2.

Figure 2.

Integrative analysis of the epigenome and transcriptome. (A) Integrative scatterplots contrasting expression differences with DNA methylation differences in CIMP-positive (left) and CIMP-negative (right) CAEBV compared with the normal control. Genes and probes associated with the promoter showing significant differences in gene expression (|log₂FC| ≥1.5 and adjusted [adj] P < .05) and DNA methylation (|ΔM| ≥2 and adj P < .05) are highlighted. Tumor suppressor genes and antiherpesvirus genes are highlighted in green and blue, respectively. (B) Pathway analysis of DNA hypermethylated silenced genes in CIMP-positive CAEBV (highlighted in the lower right corner in panel A) performed using Metascape software. (C) LOLA enrichment analysis of 100 differentially hypermethylated regions in promoters. Significantly enriched categories from CISTROME or ENCODE entries of the LOLA Core databases are shown. (D) Expression of EZH2 in NK cells from CIMP-positive, CIMP-negative, EBV-IM, and control cases (healthy volunteers and non-EBV–IM-like). Box plots show median (lines), interquartile ranges (IQRs) (boxes), and ±1.5 × IQR (whiskers). Significance was assessed using the Dunnett T3 test for multiple comparisons against controls. (E) CUT&RUN profile differences in CIMP-positive (S_C_10 NK_LCL) (top) and CIMP-negative (T_H_01 NK_LCL) (bottom) samples compared with the normal control (Control 1 NK cell) for H3K27me3 (left) and H3K4me3 (right) over the average gene body for all genes or hypermethylated and silenced genes in CIMP-positive CAEBV. GO, gene ontology; Log₂FC, log₂ fold change; TES, transcription end site; TPM, transcripts per million; TSS, transcription start site.

Figure 3.

Mutational landscape of CAEBV according to CIMP status. (A) Gene mutations and CNAs in NK cell–type CAEBV with different CIMP status. (B) Pairwise associations among mutated genes found in at least 3 patients with NK cell–type or T-cell–type CAEBV and PD-related death outcome. Only significant correlations (q < .05) are shown together with their odds ratios. ARID1A and KMT2D mutations and CD274 3′-UTR SVs are defined as risk genes and highlighted in red. (C) Kaplan-Meier survival curves for NK cell–type or T-cell–type CAEBV according to risk gene mutations. The table below the plot indicates the number of patients at risk at each time point. (D) Forest plots of multivariate analyses of OS in NK cell–type or T-cell–type CAEBV constructed using the Cox proportional hazard model (hazard ratio ± 95% CI). Risk genes and known clinical prognostic factors were included in the model. AIC, Akaike information criterion; WES, whole-exome sequencing; WGS, whole-genome sequencing.

Figure 4.

Enrichment of tissue-resident signatures in EBV-infected NK cells. (A) UMAP embedding of scRNA-seq data for 49 878 cells from PBMCs of 3 patients with CIMP-positive CAEBV, 3 patients with CIMP-negative CAEBV, 2 patients with IM, and 2 healthy volunteers. Twenty-two cell types were defined according to the RNA expression of marker genes (supplemental Figure 8A). (B) UMAP plot colored according to EBER1 expression level. (C) Violin plot showing the expression of EBER1 in each cell type. (D) The expression of NK-lineage-defining surface markers is displayed using violin plots in 3 NK clusters (NK, EBV-infected NK, and CD56^bright NK). (E) UMAP embedding of NK cells colored by CD49⁺ peripheral blood NK (CD49⁺ pbNK) up (upper) and down (lower) gene score. The score was calculated using a published gene list (supplemental Table 12). (F) Representative images of dual immunostaining for EBER (brown) and CD49a (red) in a human stomach biopsy sample (K_C_01, left), a skin lesion biopsy sample (T_H_01, middle), and a liver sample from a PDX mouse engrafted with NK-LCLs (S_C_10, right). cDC, conventional dendritic cell; dnT, double negative T; gdT, γδ T; HSPC, hematopoietic stem and progenitor cell; MAIT, mucosal associated invariant T; pDC, plasmacytoid dendritic cell; PD-L1, programmed cell death 1 ligand 1; PDX, patient-derived xenograft; Treg, regulatory T.

Figure 5.

Innate and EBV-adaptive immune compromise in CAEBV. (A-B) (Top) UMAP embedding of PBMCs colored by IFN-α response score (A) and IFN-β response score range (B). The score was calculated using gene sets termed “GOBP_ GOBP_RESPONSE_TO_INTERFERON_ALPHA (GO:0035455)” and “GOBP_RESPONSE_TO_INTERFERON_BETA (GO:0035456).” (Bottom) Heat maps depicting the average scores in each of 22 cell types. The module score of “EBV-infected NK” was compared with the “NK” score using a 2-sided Welch t test. (C) UMAP embedding of T cells colored by clonal expansion size. Clonal expansion divided into 3 categories (left) and clonal expansion sizes ranging from 0 to 250 (right) are shown. (D) Distribution of the clone status of T cells suspected to be specific to EBV based on the CDR3 amino acid sequence in CAEBV (n = 6), IM (n = 2), or healthy volunteers (n = 2). (E) Rate of EBV-specific T cells with the TCR specific for LMP2 in CAEBV (n = 6), IM (n = 2), or healthy volunteers (n = 2). The difference in the rate of TCR specific to LMP2 between CAEBV and IM was evaluated using the Fisher exact test. cDC, conventional dendritic cell; dnT, double negative T; gdT, γδ T; HSPC, hematopoietic stem and progenitor cell; MAIT, mucosal associated invariant T; pDC, plasmacytoid dendritic cell; Treg, regulatory T.

Figure 6.

Efficacy of 5-Aza in CIMP-positive NK-LCLs. (A) Number of live cells in CIMP-positive NK-LCLs (n = 2), CIMP-negative NK-LCLs (n = 1), and normal NK cells from healthy volunteers (n = 2) treated in vitro with 4 different concentrations of 5-Aza ranging from 8.0 × 10^–3 (×1/125) to 1 μg/mL (×1) for 7 days relative to those treated with vehicle (DMSO). Three wells per concentration were used in each experiment (mean ± standard deviation [SD]), and the experiment was repeated twice for validation. The overall P value was calculated using the Kruskal-Wallis test. The P values for pairwise comparisons were calculated using the Mann-Whitney U test. (B) M values were calculated for CIMP-positive NK-LCLs (S_C_10) and NK cells from a healthy volunteer (Control 4) at probes where silencing by DNA methylation was observed in CIMP-positive CAEBV. The ΔM values, representing the difference between the average M values from 3 wells of 5-Aza (1×) treatment and 3 wells of the vehicle control (DMSO), were plotted for each probe (mean ± 95% CI; n = 855 probes). Significance was assessed using the Mann-Whitney U test. (C) Hypomethylated (ΔM < –0.5 and adj P < .05) and upregulated genes in CIMP-positive NK-LCLs (S_C_10) after 7 days of in vitro 5-Aza treatment vs vehicle. Significantly upregulated genes (log₂ fold change ≥ 1 and adj P < .05) are highlighted in red. (D) Pathway analysis of DNA hypomethylated and upregulated genes in CIMP-positive NK-LCLs after in vitro 5-Aza treatment was performed using Metascape software. (E) Schematic showing the protocol for the in vivo study of 5-Aza using NOG-IL2 mice injected with CIMP-positive NK-LCLs (S_C_10). (F) Percentage of hCD45⁺ cells in the peripheral blood of 5-Aza vs vehicle-treated NOG-IL2 mice at 42 days after transplantation of the CIMP-positive NK-LCLs (S_C_10) (mean ± SD; n = 5 mice per group). Differences were evaluated using 2-way analysis of variance. (G-H) Spleen weight (G) and percentage of hCD45⁺ cells in the spleen (H) on day 42 in mice treated with 5-Aza or vehicle (mean ± SD; n = 5 mice per group). The difference was evaluated using the Mann-Whitney U test. DMSO, dimethyl sulfoxide; dsRNA, double-stranded RNA; GO, gene ontology; R-HSA, reactome-Homo sapiens.