EBV strain interacts with host HLA to drive nasopharyngeal carcinoma risk

2026-05-21Chen,Yanhong,Liang,Jingtong,Zhang,Wanlin,Zhang,Xinyu,Zhang,Xinyi,Wang,Ching-Yuan,Yao,Fei,Zhang,Shanshan,Zhou,Xiang,Ye,Weimin,Feng,Ruimei,Cai,Yonglin,Zhang,Zhe,Ji,Mingfang,Cui,Qian,Lin,Xihong,Li,Jiesen,Xu,Jialei,Zhang,Q

Nature 2026年8115期

关键词：公式

Abstract

Epstein–Barr virus (EBV) infects more than 95% of adults worldwide but is associated with endemic nasopharyngeal carcinoma (NPC) specifically in southern China 1 , 2 , 3 , 4 . Here, through a stepwise host–EBV genome interaction analysis, we identify a genetic interaction between HLA-A*11:01 and the high-risk EBV variant 85841G as a key determinant of NPC risk. Individuals carrying a susceptible HLA-A background (HLA-A*11:01 − or HLA-A*02:07 + ) and infected with the high-risk 85841G EBV form a dual-risk subgroup with substantially elevated, interaction-driven NPC risk, far exceeding the effects of host or virus alone. This dual-risk subgroup comprises 20.5% of the population and accounts for approximately 47% of NPC cases. We show that EBV 85841G encodes an EBNA3B peptide that binds to HLA-A*11:01 and elicits specific T cell responses capable of lysing EBV + B cells transformed by 85841G-carrying strains, and is associated with reduced salivary viral load and lower NPC risk among A*11:01 carriers. Evolutionary analysis reveals that 85841G arose via ancient recombination events between northern and southern EBV and subsequently underwent clonal expansion in southern China, leading to co-enrichment of interacting host and viral risk factors that, in turn, contribute to NPC endemicity. These findings reveal a markedly stratified, interaction-driven risk architecture in NPC and highlight opportunities for precision prevention.

Main

EBV infects most adults worldwide, yet its associated diseases, including NPC, Burkitt lymphoma and multiple sclerosis, are distributed unevenly across populations 1 , 2 . For example, nearly one-half of all NPC cases arise from approximately 2% of the global population residing in southern China, particularly in Guangdong and Guangxi 3 , 4 . Elucidating the causes and mechanisms underlying this epidemiological enigma is essential for effective prevention. NPC is aetiologically linked to the complex interplay between EBV infection, host genetic susceptibility and environmental exposures 4 , 5 . EBV genomic studies have identified high-risk viral lineages enriched in southern China and associated with NPC endemicity 6 , 7 , and human genome-wide association studies (GWAS) have consistently implicated the class I human leukocyte antigen (HLA) genes as the strongest host risk loci 8 , 9 , 10 , 11 , 12 . A candidate–variant study has further suggested that HLA class I variation can modify the risk conferred by high-risk EBV infection 13 . Given the central role of HLA in antigen presentation and antiviral immunity, these studies point to a compelling hypothesis that NPC susceptibility may be shaped by host–EBV interactions rather than independent effects of viral strain and host genotype. However, the causal determinants of such interactions on the host and viral sides remain unresolved.

Motivated by genome-to-genome studies that have mapped critical host–pathogen interactions in infectious diseases 14 , 15 , 16 , 17 , 18 , 19 , 20 , here we applied a stepwise genome-to-genome scan to resolve the genetic basis of the host–EBV interaction in NPC (Extended Data Fig. 1). We further integrated functional and viral population genetic analyses to establish the molecular basis of this interaction and to reconstruct the evolutionary history of the interacting EBV variant. Together, our results show that NPC risk is shaped not by host and viral variants independently, but by their synergistic interactions, providing a foundation for precision prevention strategies targeting susceptible populations.

HLA–EBV interaction drives GWAS signal

To identify the host and viral determinants underlying the host–EBV interaction in NPC susceptibility, we adopted a two-step genome-to-genome design optimized for statistical power (Extended Data Fig. 1 and Supplementary Notes 1 and 2). We screened the human genome for variants whose genetic effects on NPC risk differed by high-risk EBV infection status, using two independent datasets (step 1; discovery: 747 cases and 1,251 controls; validation: 644 cases and 880 controls), and then fine-mapped the EBV genome to pinpoint viral variants driving this interaction in paired host–viral genomic data (step 2; 324 cases and 410 controls). Power simulations indicated more than 80% power for both steps under the observed effect sizes (Methods; Supplementary Note 1).

We defined NPC-associated high-risk EBV subtype by three coding variants in BALF2 (single-nucleotide polymorphisms (SNPs) 162215A>C, 162476T>C and 163364C>T; hereafter BALF2-CCT) 6 or the alternative marker in the non-coding EBER2 gene (SNP 7048A>C; hereafter 7048C) 7 . BALF2-CCT and 7048C are in strong linkage disequilibrium (r 2 = 0.70) and largely tag the same high-risk EBV lineage, present in approximately 78% of NPC cases in southern China and conferring approximately ninefold increased risk (Supplementary Fig. 1 and Supplementary Note 3). BALF2-CCT-defined high-risk EBV infection was used as the interacting factor in step 1.

For step 1, we performed a host genome-wide interaction scan in a discovery dataset and evaluated the lead signals in an independent validation dataset. NPC cases and ancestry-matched controls in step 1 were genotyped using an Illumina genotyping array. BALF2-CCT status was determined from saliva DNA using SNP genotyping 6 (Methods; Extended Data Fig. 1, Supplementary Tables 1–3 and Supplementary Fig. 2). Population structure analyses confirmed that all participants were of southern Chinese ancestry (Supplementary Fig. 3).

We applied a genome-wide gene–environment interaction scan approach, which estimates interaction as the difference between the overall genetic effect and the main genetic effect after excluding the interaction effect 21 . In the discovery phase, we screened for host genetic variants whose effects on NPC risk differed by high-risk EBV status (BALF2-CCT versus other strains), adjusting for age, sex and the top four human multidimensional scaling (MDS) components. Quantile–quantile plots confirmed no inflation of test statistics (Extended Data Fig. 2a–c). This genome-wide scan identified the HLA region as the sole locus with genome-wide significant genotype × EBV interaction signals, with the top interaction signal at rs417162 (P interaction = 4.45 × 10 −10 ; Fig. 1c). rs417162 was in almost complete linkage disequilibrium (r 2 = 0.99) with rs2860580, a top risk variant from previous GWAS 9 , 22 , and also showed strong association with NPC in the current dataset (P overall = 5.73 × 10 −12 ; Fig. 1a). Including the genotype × EBV interaction term abolished association signals across the genome, including those at the HLA locus (Fig. 1b), suggesting that the well-established GWAS signals at the HLA locus may be largely driven by host–EBV interaction.

**Fig. 1: Genome-wide scan of host–EBV genetic interactions in NPC.**

A*11:01 underlies HLA–EBV interaction

To fine-map HLA variants and functional alleles underlying the HLA–EBV interaction in NPC susceptibility, we imputed a common set of 4,161 variants, including classical HLA alleles, amino acids in MHC proteins and SNPs, in the MHC region using SNP2HLA with the Pan-Asian reference panel 23 , 24 , 25 (Supplementary Table 4). In the discovery phase, the strongest interaction signal was observed at HLA-A amino acid position 62, where the presence of glutamine (A62Q) showed a highly significant interaction with high-risk EBV strains (P interaction = 5.88 × 10 −12 ). This residue, located in the peptide-binding groove, has previously been identified as one of the top NPC susceptibility variants 8 , 11 . At the allelic level, HLA-A*11:01, which is tagged by A62Q (r 2 = 0.70), exhibited the most statistically significant interaction (P interaction = 6.55 × 10 −10 ; Extended Data Fig. 2d). Conditioning on either HLA-A*11:01 or A62Q abolished all residual interaction signals within the HLA locus, suggesting that this allele or its tagging residue probably represents the driver of the observed HLA–EBV interaction (Extended Data Fig. 2e,f).

To validate our findings, we analysed the interaction between common HLA-A alleles and the high-risk EBV subtype BALF2-CCT in an independent dataset (validation phase). We confirmed the HLA-A*11:01 × BALF2-CCT interaction on NPC risk (P interaction = 8.48 × 10 −5 ), and meta-analysis of the discovery and validation datasets provided stronger evidence (P interaction = 4.89 × 10 −13 ; Supplementary Table 5). HLA-A*11:01 carries several amino acid residues in the peptide-binding groove that show significant interaction signals with high-risk EBV strains (A62Q, A70Q and A114R; Supplementary Table 6), supporting a potential functional interaction. Of note, among non-high-risk EBV infections, A*11:01 showed only a weak protective trend in meta-analysis (P = 7.18 × 10 −4 ), whereas among high-risk EBV infections, it markedly modified NPC risk: non-carriers of A*11:01 had significantly higher risk than carriers (odds ratio (OR) = 10.60 versus 5.78, P = 1.05 × 10 −9 in meta-analysis; Fig. 1d, Extended Data Fig. 2g and Supplementary Table 7). Moreover, the interaction with high-risk EBV explained the majority (approximately 85%) of the overall protective effect of A*11:01 on NPC risk (Extended Data Fig. 2h).

In addition to A*11:01, we observed a second suggestive host × high-risk EBV interaction for HLA-A*02:07: A*02:07 showed minimal risk effect in non-high-risk EBV carriers, but among high-risk EBV carriers, A*02:07 carriers had substantially elevated risk compared with non-carriers (OR = 10.38 versus 6.20; Fig. 1e and Supplementary Table 7). Nevertheless, after conditioning on the main and interaction effects of A*11:01, evidence for an independent A*02:07 interaction effect in the meta-analysis did not meet genome-wide significance (P interaction = 4.86 × 10 −4 ; Supplementary Table 5). We also detected borderline interaction signals involving HLA-B alleles in the meta-analysis, but these were also much weaker than the HLA-A interaction signals and did not remain significant after multiple-testing correction (Supplementary Notes 4 and 5). Although interesting, these additional suggestive interaction effects will require future studies for validation. Together, these findings indicate that HLA-A*11:01 is the primary genetic modifier of NPC risk associated with high-risk EBV infection.

EBV 85841G underlies A*11:01 interaction

The discovery of host–EBV interactions between HLA-A*11:01 and high-risk EBV strains raised a key question regarding which viral variants drive this interaction. To address this, we conducted an EBV genome-wide interaction scan to pinpoint variants interacting with HLA-A*11:01 in NPC risk. Using capture-based targeted sequencing 6 , we obtained EBV whole-genome sequences from 734 individuals (324 cases and 410 controls with no strain-based selection; Methods; Extended Data Fig. 1, Supplementary Tables 8 and 9 and Supplementary Figs. 4–6). We applied the same interaction approach as in the host genome-wide analysis to scan the EBV genome for variants interacting with HLA-A*11:01. To control viral population structure and relatedness 26 , we fit a generalized linear mixed model with a viral genetic relatedness matrix random effect and fixed effects for age, sex, the top four viral principal components and the top four host MDS components (Methods; Extended Data Fig. 3a,b and Supplementary Note 6). This model adjusted for broad viral clade structure (viral principal components) and residual viral relatedness, as well as host population structure. The quantile–quantile plot showed that this model provided good control of viral population structure with modest inflation (Extended Data Fig. 3c). The most significant interaction was observed at EBV SNP 85841A>G (P interaction = 5.57 × 10 −4 ), surpassing the empirical EBV genome-wide threshold (P threshold = 5.92 × 10 −4 ; Fig. 2a,b and Supplementary Note 7). This significance threshold was derived from 10,000 permutations using a minimum-P procedure to control the family-wise error rate (α = 0.05) across the EBV genome. Permutations were implemented in BRASS, which accounts for sample clonality and relatedness (Methods) 27 . Conditioning on the viral SNP 85841 abolished the interaction between HLA-A*11:01 and EBV variants across the genome (Extended Data Fig. 3e,f). In particular, when we further excluded the top variant (85841) and the variants within its ±1-kb flanking region, the quantile–quantile plot became deflated (Extended Data Fig. 3d), confirming that residual inflation in the overall association results was minimal. Moreover, the interaction effect at SNP 85841 was more significant than that observed for the high-risk EBV subtype BALF2-CCT or 7048C (Supplementary Note 8). Together, these results indicate that SNP 85841 is the primary viral variant driving the observed interactions between HLA-A*11:01 and high-risk EBV strains.

**Fig. 2: EBV genome-wide scan for the interaction with HLA-A*11:01.**

Supporting this finding, stratified analyses revealed that the HLA-A*11:01 interaction was confined to the high-risk EBV strains carrying the 85841G variant (hereafter high-risk 85841G). We quantified the interaction effect using the relative excess risk due to interaction (RERI), a classical epidemiological measure that captures the excess joint effect of carrying both risk factors beyond the sum of their individual effects, thereby reflecting the contribution of the host–virus interaction to the overall disease risk 28 , 29 (Supplementary Note 9). Individuals lacking HLA-A*11:01 and infected with high-risk 85841G strains had RERI 5.44 compared with A*11:01 carriers infected with non-high-risk EBV strains (Fig. 2c). By contrast, high-risk strains carrying the 85841A variant (high-risk 85841A) conferred elevated NPC risk regardless of HLA-A*11:01 status, showing no evidence of interaction (Fig. 2c).

We validated the interaction between HLA-A*11:01 and the high-risk 85841G strains in two independent datasets from southern China and Singapore (Extended Data Fig. 1, Supplementary Table 10 and Supplementary Fig. 4). EBV genotypes for SNP 85841 and the NPC-associated markers BALF2-CCT and 7048C were determined for interaction analysis. The southern Chinese replication dataset included 388 cases and 543 controls independent of the viral genome-wide scan. This analysis found that the interaction between HLA-A*11:01 and high-risk 85841G EBV was stronger than that observed with the high-risk subtype (BALF2-CCT or 7048C) alone (Supplementary Table 11). Individuals lacking HLA-A*11:01 and infected with high-risk 85841G strains showed a significantly elevated NPC risk attributable to interaction (RERI = 9.23, 95% confidence interval (CI): 4.11–20.61, P interaction = 4.01 × 10 −3 ), whereas no such effect was observed for 85841A-carrying strains (Fig. 2d). The Singapore dataset, comprising 223 cases and 204 controls, provided an additional replication in a geographically distinct population. Despite a smaller sample size, a significant HLA-A*11:01 × high-risk 85841G interaction was observed (RERI = 10.11, P interaction = 0.04; Fig. 2e). Collectively, across discovery and two independent replication datasets in southern China and Singapore, these results demonstrate that the host × EBV genetic interaction in NPC risk is primarily driven by EBV EBNA3B variant 85841G and HLA-A*11:01, with consistent evidence across populations (Fig. 2c–e and Extended Data Fig. 4a).

Despite a suggestive interaction between HLA-A*02:07 and high-risk EBV in modulating NPC risk, its lower allele frequency (13% for A*02:07 versus 28% for A*11:01) limited power for a genome-wide viral interaction scan (simulated power: 72% for A*02:07 versus 91% for A*11:01; Methods; Supplementary Note 1). Nevertheless, a targeted analysis revealed that the potential A*02:07 × EBV interaction effect, if present, is likewise restricted to high-risk 85841G strains (Extended Data Fig. 4b). The specific viral variants that drive any interaction with A*02:07 remain to be determined in future studies.

85841G encodes an A*11:01 T cell epitope

The repertoire of host HLA alleles determines which viral peptides are presented to T cells, shaping adaptive immune responses that control EBV replication. Previous studies have reported that in populations in which HLA-A*11:01 is prevalent, EBV strains frequently carry escape mutations in two immunodominant EBNA3B epitopes, IVTDFSVIK (IVT; N9 epitope loss: K>N at position 9) and AVFDRKSDAK (AVF; N4 epitope loss: D>N at position 4), which disrupt A*11:01 binding and impair T cell responses 30 , 31 , 32 , 33 , 34 . In our dataset, we found that high-risk EBV lineages had near-universal loss of the AVF and IVT epitopes (Extended Data Fig. 5 and Supplementary Table 12). Therefore, A*11:01 probably confers protection in carriers infected with high-risk 85841G EBV through alternative epitopes.

SNP 85841, the top variant interacting with HLA-A*11:01, encodes a non-synonymous change (85841A>G resulting in Q900R) in EBNA3B. To explore whether this variant might create a novel HLA-A*11:01-restricted epitope, we scanned all 9–11-mer peptides spanning the 85841G variant. A candidate peptide, VVILENVSR, with arginine encoded by 85841G, showed strong predicted binding to A*11:01 (percentile rank of 0.8% by NetMHCpan 4.1 (ref. 35 ); Supplementary Table 13). By contrast, the 85841A-encoded peptide, VVILENVGQ, had no predicted binding (percentile rank of 9.6%; Fig. 3a). We performed an enzyme-linked immunosorbent assay (ELISA) using a peptide-exchangeable HLA-A*11:01 monomer and observed specific binding of VVILENVSR, but not VVILENVGQ, to A*11:01, whereas neither peptide showed detectable binding to HLA-A*02:07 (Extended Data Fig. 6a,b). Using ELISA with positional substitution peptides, we further found that the C-terminal arginine (position 9) encoded by 85841G, rather than serine at position 8, was critical for HLA-A*11:01 binding, consistent with the known A*11:01-binding motif and structural predictions (Extended Data Fig. 6a–d). Consistently, using biolayer interferometry, we confirmed that the 85841G-encoded peptide bound HLA-A*11:01 with high affinity (dissociation constant (K d ) = 8.5 μM), whereas the 85841A-encoded peptide showed almost no binding (K d = 236.2 μM; Extended Data Fig. 6e).

Fig. 3: The 85841G-encoded EBNA3B peptide VVILENVSR binds to HLA-A*11:01 and elicits A*11:01-restricted cytotoxic CD8+ T cell responses. — **Fig. 3: The 85841G-encoded EBNA3B peptide VVILENVSR binds to HLA-A*11:01 and elicits A*11:01-restricted cytotoxic CD8 + T cell responses.**

We next investigated whether the 85841G-encoded peptide is functionally recognized by A*11:01-restricted CD8 + T cells. Using peripheral blood mononuclear cells from A*11:01 + individuals infected with high-risk 85841G EBV strains, we expanded polyclonal T cells with three EBNA3B peptide pairs (representing three epitopes and their epitope-loss variants): AVF versus AVF/N4, IVT versus IVT/N9 and VVILENVSR versus VVILENVGQ (Extended Data Fig. 7a). These T cells were then co-cultured with COS-7 cells engineered to express HLA-A*11:01 and β 2 -microglobulin (β 2 M). The COS-7 cell line naturally lacks endogenous human HLA expression and is commonly used as the antigen-presenting cell (APC) after HLA–β 2 M transfection 36 , 37 , 38 , 39 . IFNγ intracellular cytokine staining showed that the 85841G-encoded peptide VVILENVSR elicited markedly stronger CD8 + T cell activation than either AVF or IVT peptides and their corresponding epitope-loss peptides (AVF/N4 and IVT/N9), whereas the 85841A-encoded peptide VVILENVGQ did not activate peptide-expanded CD8 + T cell responses (Fig. 3b). Consistently, cytotoxicity assays showed specific killing of A*11:01-expressing target cells pulsed with VVILENVSR, but not with any of the other five peptides (Fig. 3c). Peptide titration further showed dose-dependent T cell activation and cytotoxicity for VVILENVSR, but not VVILENVGQ, when presented by A*11:01-expressing COS-7 APCs (Fig. 3d,e). This T cell response was restricted to HLA-A*11:01, as COS-7 APCs expressing HLA-A*02:07 or A*02:01 did not support robust VVILENVSR-specific activation or killing (Extended Data Fig. 7b,c). Furthermore, endogenous antigen-presentation experiments using A*11:01 + APCs expressing EBNA3B haplotypes indicate that only the high-risk 85841G EBNA3B-encoded peptide, but not the peptide from either high-risk 85841A or the low-risk EBNA3B haplotype, can be processed and presented to trigger VVILENVSR-specific CD8 + T cell cytotoxicity (Fig. 3f,g). Consistent with these findings, VVILENVSR-specific CD8 + T cells showed stronger IFNγ-producing responses to A*11:01 + lymphoblastoid cell lines (LCLs) transformed by the NPC-derived M81 EBV (high-risk 85841G + ) as well as Akata EBV (N1 lineage), which also carries the 85841G-encoded VVILENVSR epitope (Fig. 3h). By contrast, VVILENVGQ-expanded T cells did not show specific responses. These VVILENVSR-specific T cells also selectively lysed A*11:01 + LCLs transformed by M81 or Akata, with substantially reduced killing of A*11:01 + LCLs transformed by the infectious mononucleosis-derived B95-8 strain (85841A + ) or A*11:01 − LCLs (Fig. 3i). Of note, among HLA-A*11:01 + individuals infected with high-risk 85841G strains, VVILENVSR-HLA-A*11:01 tetramer + CD8 + T cells were reduced in NPC cases compared with healthy carriers (Fig. 3j), in line with previous reports of generally reduced EBV-specific T cell immunity in NPC 40 , 41 , 42 . Together, these results establish VVILENVSR as a functional HLA-A*11:01-restricted epitope generated by 85841G and show that it elicits specific CD8 + T cell responses against high-risk 85841G EBV-transformed B lymphocytes despite the loss of the EBNA3B AVF and IVT epitopes.

In healthy carriers, EBV intermittently reactivates from the latent EBV + B cell reservoir and sheds into saliva. We additionally examined salivary EBV load in healthy carriers as an in vivo readout linked to mucosal viral activity. Owing to the episodic nature of EBV mucosal shedding 43 , salivary EBV DNA levels showed large inter-individual variability and did not differ appreciably between high-risk and low-risk EBV infections overall. By contrast, among carriers infected with the HLA-interacting high-risk 85841G strains, the HLA-A background affected EBV load in opposite directions: HLA-A*11:01 was associated with lower levels of salivary EBV DNA, whereas A*02:07 was associated with higher levels (Fig. 3k and Supplementary Note 10). These associations were not evident in individuals infected with non-interacting strains (Fig. 3l and Extended Data Fig. 7d,e). Together, these functional assays and in vivo data link high-risk 85841G infection and the HLA-A background to differential EBV-specific T cell responses and viral activity.

EBV–HLA interaction stratifies NPC risk

As EBV migrated with humans during the out-of-Africa dispersal, serial founder events shaped the geographical stratification of both host and viral genomes 44 . We found marked spatial co-enrichment in southern China between high-risk 85841G EBV and its interacting HLA-A*11:01, and also enrichment of the established major NPC risk allele HLA-A*02:07 (refs. 5 , 8 , 11 , 12 , 45 ), although evidence for an independent A*02:07–EBV interaction remains suggestive. Both alleles show pronounced south–north frequency gradients (Fig. 4a): the protective allele A*11:01 reaches approximately 28% frequency in southern China, compared with approximately 17% in the north, whereas the risk allele A*02:07 is approximately 13% in the south but only about 7% in northern Chinese populations. These gradients parallel the increasing NPC incidence from northern to southern China. For EBV, the high-risk 85841A strains showing little evidence of interaction with HLA-A are common in regions of intermediate NPC incidence (for example, Fujian, Taiwan and Singapore). However, the high-risk 85841G EBV strains predominate in high-incidence regions, including Guangdong (46.63%), Guangxi (28.90%) and Hong Kong (26.95%), but sharply decline in intermediate-incidence regions in China (8.89%) and Singapore (14.71%), and are nearly absent in northern China and Northeast Asia, where NPC incidence is very low (Fig. 4a and Supplementary Table 14, all data from healthy donors). Thus, the high-risk 85841G lineage is geographically co-enriched with the HLA-A background defined by the A*11:01 and A*02:07 frequency gradients in southern China, mirroring the geographical pattern of NPC incidence.

**Fig. 4: Geographical co-enrichment of high-risk EBV and HLA-A alleles concentrates interaction-driven NPC risk in southern China.**

To dissect how host–virus interactions and their geographical co-enrichment translate into population-level disease burden, we stratified individuals by EBV subtype and HLA-A background, including A*11:01 (the interacting allele) and A*02:07 (the major NPC risk allele). We assessed NPC risk and population attributable fractions (PAFs) across strata (Fig. 4b,c and Supplementary Table 15). Here, PAF denotes the proportion of cases in the population attributable to each stratum 46 , 47 . We observed a stepwise increase in NPC risk across virus–host genotype combinations. The baseline group (approximately 23% of the population), carrying low-risk EBV strains and a protective HLA-A background (A*11:01 + and A*02:07 − ), accounted for only 16% (95% CI 8−23%) of cases. Infection with high-risk 85841A strains moderately increased risk (OR = 5−7; PAF = 11%), largely independent of host HLA-A background. By contrast, high-risk 85841G strains conferred higher risk, even among individuals with a protective HLA-A background (OR = 9, 95% CI 5−14; PAF = 20%, 95% CI 13−27%), but risk was further amplified among individuals with a susceptible HLA-A background (A*11:01 − or A*02:07 + ), reaching a joint effect of OR of 18.6 (95% CI 11.6−29.8). In this dual-risk subgroup, defined by both high-risk 85841G EBV and the susceptible HLA-A background, the total relative excess risk is approximately 17.6 (OR = −1), of which approximately 9.2 is attributable to interaction (RERI), approximately 0.7 to HLA-risk alone and approximately 7.6 to EBV-risk alone, with interaction accounting for 52.6% of the total excess risk in this stratum (Fig. 4b,c and Supplementary Table 15). Of note, this most susceptible dual-risk subgroup comprises approximately 20.5% of the southern Chinese population, approximately 35 million individuals, yet accounts for approximately 47% (95% CI 42−52%) of regional NPC burden (Fig. 4b,c and Supplementary Table 15).

We additionally estimated the corresponding absolute NPC incidence rates for each HLA–EBV stratum in endemic southern China, based on the regional population incidence and the observed HLA allele and EBV strain frequencies, using the corresponding OR estimates (Fig. 4c). The dual-risk subgroup reached an estimated incidence of approximately 57.0 per 100,000 person-years, compared with approximately 3.1 in the baseline stratum without either risk factor, approximately 5.3 in the HLA-risk-only stratum and approximately 16.5–26.4 in the EBV-risk-only strata. The EBV-risk-only strata include high-risk 85841A infection regardless of HLA-A background, as well as high-risk 85841G infection on the protective HLA-A background. These findings reveal a sharply stratified risk landscape driven by convergent enrichment of high-risk EBV and host susceptibility, with direct implications for precision prevention in endemic populations.

Recombination and natural selection shaped 85841G

The convergence between host and viral genetic risk factors raises a key evolutionary question about how the 85841G variant originated and became so frequent in southern China. To investigate this, we constructed an EBV phylogeny by integrating 954 newly sequenced genomes with 1,132 publicly available genomes from the NCBI database (including 1,852 type 1 genomes in total; Methods; Fig. 5, Extended Data Fig. 1, Supplementary Tables 16 and 17 and Supplementary Fig. 7). Mapping the 85841G variant onto this phylogeny revealed that although most 85841G-carrying strains were nested within the high-risk EBV clade, a distinct subset was located outside this clade, predominantly derived from northern China (denoted as N1; Fig. 5). This non-monophyletic distribution suggests a mosaic evolutionary history for the 85841G variant.

**Fig. 5: The evolutionary history of the EBNA3B 85841G variant.**

Looking into the high-risk EBV clade bearing the 85841G variant, we found two lineages: a clonally expanded lineage (C1) and a distinct subgroup outside of C1 (denoted as X1; Fig. 5). We previously estimated that the clonal lineage C1 originated approximately 4,000 years ago from a recombination event between S1-like EBV strains from southern China (S1) and N1-like strains from the north (N1) 44 (Fig. 5). Here, using multiple population genetic methods, including chromosome painting, we confirmed that the central genomic region spanning EBNA3B 85841G was introgressed from N1 into S1-like strains (Fig. 5 and Supplementary Note 11). Of note, C1 strains showed strong signatures of recent selective sweep, including markedly short phylogenetic tip lengths, high genome similarity and elevated extended haplotype homozygosity spanning EBNA3B (approximately 48–108 kb), with a strong estimated selection coefficient of 4.3% (95% CI 2.8–8.4%) 44 (Supplementary Figs. 8 and 9 and Supplementary Note 11). We note that, in the absence of ancient EBV genomes, the recombination date and the inferred selection coefficient may vary depending on model assumptions, whereas the recombination-derived mosaic structure of C1 and X1 and the evidence for clonal expansion of C1 in southern China are supported by multiple evolutionary analyses. This positive selection probably underlies its current dominance in southern China, where C1 strains account for approximately 31% of EBV infections and approximately 74% of high-risk strains (Supplementary Table 14).

Chromosome painting of the X1 lineage revealed a distinct recombination landscape, involving partially overlapping but distinct N1-derived segments that also contained the 85841G variant (Fig. 5 and Supplementary Note 11). These findings suggest that the 85841G variant was acquired from northern strains independently through at least two distinct recombination events. X1 showed an interaction trend with HLA-A*11:01 on NPC risk, in the same direction as the interaction observed for C1, although only C1 reached statistical significance (Extended Data Fig. 8). As previously reported, homoplasy arising from recurrent mutations or recombination events that break linkage can increase mapping power in pathogen GWAS 26 . This recurrence of 85841G across multiple recombined lineages and the consistent interaction with HLA-A*11:01 support 85841G as a causal factor in host−EBV genetic interactions shaping NPC susceptibility.

Discussion

Through a stepwise genome-to-genome interaction analysis, we pinpoint HLA-A*11:01 and EBNA3B 85841G as the primary host–EBV interaction shaping NPC risk (Supplementary Note 12). The long-established protective association of HLA-A*11:01 is shown to depend largely on infection with 85841G-carrying high-risk EBV strains 5 , 8 , 11 , 12 , 45 . The 85841G substitution encodes an A*11:01-restricted epitope VVILENVSR (VVI) that elicits CD8 + T cell responses against EBV + B cells transformed by the high-risk 85841G strain M81. Although EBNA3B is not thought to be expressed in latency II NPC tumour cells, it is expressed in infected B cells in latency III, an EBV transcriptional program that supports continuous LCL growth and is linked to viral replication and reseeding in tonsillar tissue 48 , 49 . These data support a model in which VVI-specific A*11:01-restricted CD8 + T cells restrain the EBV-infected B cell reservoir and limit reactivation, thereby reducing epithelial reseeding and the probability of NPC initiation (Extended Data Fig. 9). Together with previous evidence that the NPC-derived strain M81 exhibits enhanced epithelial tropism and reactivation 50 , these findings support a two-layer architecture of NPC risk: viral genotype sets a baseline gradient, and the host HLA background modulates risk in a strain-dependent manner (Supplementary Note 13).

In southern China, high-risk 85841G strains and HLA-A*11:01 are co-enriched, sharpening NPC risk stratification at the population level. In individuals lacking HLA-A*11:01, recognition of the VVI epitope is absent, and high-risk 85841G EBV infection is associated with substantially higher risk. Accordingly, high-risk 85841G infections occurring in individuals lacking A*11:01 capture approximately 40% of cases in southern China (Supplementary Table 18). This burden becomes further concentrated when we additionally incorporate the major risk allele HLA-A*02:07. The susceptible HLA-A background (A*11:01 − or A*02:07 + ) together with high-risk 85841G EBV infection defines the dual-risk subgroup that comprises 20.5% of the population but contributes nearly one-half (47%) of cases, representing a priority target for screening. Such concentration of disease burden, driven by a biologically interpretable host–virus interaction, contrasts with the more diffuse polygenic risk architecture commonly observed in complex, non-Mendelian diseases 51 . This interaction-driven risk stratification provides a rational basis for precision prevention, informing whom to screen and which EBV strains to target for vaccine design and T cell-based therapies. Reducing exposure to high-risk EBV strains could, in principle, prevent up to approximately 78% of cases, corresponding to approximately 32,000 fewer cases annually in southern China (Fig. 4c).

From an evolutionary perspective, EBV fitness is probably shaped primarily by persistence and transmission rather than oncogenicity. The expansion of the recombinant, 85841G-carrying C1 lineage is therefore unlikely to reflect direct positive selection on the A*11:01-restricted VVI epitope itself. Instead, recurrent loss of the immunodominant A*11:01-restricted EBNA3B epitopes AVF and IVT has been linked to immune-driven viral adaptation in multiple populations where A*11:01 is prevalent 30 , 31 , 33 , 34 . Selection on any advantageous features within the C1 recombination block, or through association with human admixture followed by demographic expansion, could have driven the rapid rise of the 85841G (and the VVI epitope), but the precise selection and co-evolutionary process remains unresolved.

The A*11:01 × EBV 85841G interaction is unlikely to be the sole determinant of NPC risk. Additional viral variation within the high-risk lineage has been linked to EBV replication, cell growth and immune modulation, and EBNA3B variation within the recombination block may further influence strain-specific oncogenicity or T cell recruitment via CXCL9 and CXCL10 regulation 50 , 52 , 53 , 54 , 55 , 56 . More broadly, these findings underscore that susceptibility to pathogen-driven diseases may arise from dynamic host–pathogen interplay rather than static host or pathogen polymorphism alone. As expanded paired host–EBV datasets become available, comprehensive genome-to-genome scans may uncover additional interactions beyond those captured here 57 , 58 . Emerging animal models of EBV infection may further elucidate the mechanisms of pathogenic host–EBV interactions and shared principles across EBV-associated diseases 56 .

Methods

Study participants

Participants from southern China were enrolled through two independent recruitments (sample sets 1 and 2). Participants in sample set 1 were participants in a population-based case–control study conducted in NPC-endemic regions of southern China (Guangdong and Guangxi provinces) between 2010 and 2014. The study design was previously described in detail 59 . In brief, 2,554 treatment-naive patients histologically confirmed with NPC were identified through a rapid case ascertainment system involving a network of local physicians. For the population-based control recruitment, 2,648 healthy control individuals were frequency-matched to cases by sex, 5-year age group and residential area, and were randomly selected from local population registries. Saliva DNA samples were available for 1,202 cases and 1,780 controls. In sample set 1, 747 patients with NPC and 1,251 healthy controls with age, sex, and available EBV genotype and host genome-wide genotyping were included in the discovery phase of the host genome-wide scan (Extended Data Fig. 1, Supplementary Table 1 and Supplementary Fig. 2). In sample set 2, an independent NPC case–control study was recruited from the same Chinese population in southern China. This study consisted of 883 cases and 1,537 controls of self-reported Chinese ancestry recruited between 2013 and 2022. For sample set 2, 644 patients with NPC and 880 healthy participants with age, sex, and available EBV genotype and host genome-wide genotyping were included in the validation phase of host genome-wide interaction study (Extended Data Fig. 1, Supplementary Table 1 and Supplementary Fig. 2). All study participants were recruited irrespective of EBV strain status, and no selection was performed on the basis of infection with any specific EBV subtypes.

The Singapore dataset was completely independently recruited in Singapore, comprising 226 cases with NPC and 209 controls. Among these, 223 cases with NPC and 204 healthy controls with available EBV genotype and HLA allele data were included in the replication of EBV subtype–host interaction analysis (Extended Data Fig. 1 and Supplementary Fig. 4).

Furthermore, to capture EBV genomic diversity across China for viral population genetic and phylogenetic analysis, we enrolled 163 healthy participants from different regions of China, including northern China (Shandong, Shanxi, Inner Mongolia, Xinjiang, Hebei and Liaoning), eastern China (Fujian) and southwestern China (Sichuan). Among them, 113 EBV whole-genome sequences passed quality control (see details below) and were included in subsequent EBV phylogenetic analyses (Extended Data Fig. 1, Supplementary Table 16 and Supplementary Fig. 7).

This study was approved by the institutional ethics committees of the Sun Yat-sen University Cancer Center, Guangzhou, China and the Genome Institute of Singapore, A*STAR, Singapore. Written informed consent was obtained from all participants.

Sample collection and processing

Saliva samples were collected into vials containing an equal volume of prepared lysis buffer (50 mM Tris, pH 8.0, 50 mM EDTA, 50 mM sucrose, 100 mM NaCl and 1% SDS) and stored at −80 °C. DNA was extracted from saliva using either the Chemagic STAR (Hamilton Robotics) or the QIAamp DNA Blood Midi Kit (Qiagen) according to the manufacturer’s instructions. The extracted DNA was subsequently used for human and EBV genotyping, as well as EBV whole-genome sequencing.

Peripheral blood samples were collected and processed within 24 h. Human peripheral blood mononuclear cells (PBMCs) were isolated by Ficoll density gradient centrifugation, cryopreserved in liquid nitrogen and subsequently used for T cell functional assays.

Data acquisition and processing

The majority of the paired data used for the host–EBV interaction analyses were newly generated in this study (Extended Data Fig. 1). For host genome genotyping, all data used in the genome-wide interaction scan (n = 3,522) were generated for this study, comprising case–control datasets for the host genome-wide interaction analysis with high-risk EBV status (step 1: 747 cases and 1,251 controls in discovery; and 644 cases and 880 controls in validation). For EBV whole-genome sequences, this study in total included 732 newly sequenced EBV genomes, together with 1,354 EBV genomes from NCBI databases. The EBV genome-wide interaction fine-mapping was conducted in the paired host–EBV genome dataset (n = 734), within which 154 out of 734 EBV genomes were previously deposited in the NCBI in our earlier work 6 , 44 .

Genotyping and quality control of human genetic variants

Extracted DNA was used for human genotyping using the Asian Screening Array (ASA) Chip (Illumina) or OmniZhongHua-8 Chip (Illumina). The genotype data from participants with complete demographic information and EBV genotyping information (BALF2-CCT SNPs 162215A>C, 162476T>C and 163364C>T) were subjected to sample and genotype quality control, as detailed below (see workflow and detailed quality control metrics in Supplementary Table 3 and Supplementary Fig. 2). Genotype quality control was performed using PLINK (v1.9) 60 . Samples were excluded based on the following criteria: (1) overall variant call rate ≤ 92%, (2) heterozygosity rate deviating beyond ± 3 standard deviations from the mean, or (3) relatedness > 20%. Variants were removed if they demonstrated: (1) call rate ≤ 95%, (2) minor allele frequency (MAF) ≤ 5%, or (3) significant deviation from Hardy–Weinberg equilibrium with thresholds of P < 1 × 10 −6 in controls and P < 1 × 10 −10 in cases. To detect population outliers and stratification, a MDS approach was employed using either the samples from this study alone or in combination with reference samples from the 1000 Genomes Project ( ) 61 , as illustrated in Supplementary Fig. 3. Following quality control, 1,998 participants from sample set 1, and 1,524 participants from sample set 2 were used for further analysis (Extended Data Fig. 1, Supplementary Table 3 and Supplementary Fig. 2).

HLA imputation and HLA typing

Following quality control of the human genotype data, variants within the HLA region (chromosome 6: 25000000–35000000, GRCh37/hg19) were extracted for imputation. HLA region imputation was conducted using SNP2HLA 23 with the Pan-Asian reference panel 24 , 25 . Post-imputation filtering retained variants with INFO > 0.5 (INFO is an imputation quality metric indicating the reliability of imputed genotypes) and MAF > 5% (Supplementary Table 4 and Supplementary Fig. 2). HLA fine-mapping results were visualized using the code proposed by Raychaudhuri et al. 62 .

To validate the accuracy of HLA imputation derived from our human genotype data, we performed high-resolution HLA typing on a randomly selected subset of 732 successfully imputed samples from sample set 2. In brief, genomic DNA was sheared to approximately 150–200 bp, followed by end-repair and adaptor ligation to construct sequencing libraries. Targeted enrichment of HLA loci was performed using the commercially available TargetSeq Human HLA Panel (iGeneTech). The enriched libraries were sequenced on the Illumina NovaSeq 6000 platform with 150-bp paired-end reads. Sequence data were processed using HLA-HD (v1.2.0) 63 for high-resolution HLA typing. Comparative analysis against HLA typing data in 732 samples demonstrated allele-level concordance rates of 92.9% at HLA-A, 88.5% at HLA-B, 95.9% at HLA-C, 88.9% at HLA-DQB1, 87.5% at HLA-DPB1 and 83.8% at HLA-DRB1 loci for the imputation approach.

EBV whole-genome sequencing, variant calling and principal component analysis

We first quantified EBV DNA in each sample by real-time PCR, targeting a fragment of the BALF5 gene (5′ primer: GGTCACAATCTCCACGCTGA; 3′ primer: CAACGAGGCTGACCTGATCC). Samples with an EBV cycle threshold (Ct) value below 31 were selected for subsequent EBV whole-genome sequencing (WGS), as EBV DNA with Ct values above 31 was insufficient for successful WGS (Supplementary Figs. 4 and 7).

Saliva DNA was prepared following the NadPrep EZ DNA Library Preparation Protocol (v2.2). In brief, DNA was fragmented, purified, end blunted, adaptor ligated and then amplified. The DNA libraries were subjected to hybrid capture using the EBV-targeting single-stranded DNA probes developed by Integrated DNA Technologies. After capture enrichment using the Nanodigmbio Hybridization Capture of DNA Libraries Protocol (v1.0), libraries were sequenced (paired-end 150 bp) using the Illumina NextSeq 500 platform.

Raw reads were trimmed and filtered using fastp 64 . The processed paired-end reads were aligned to the EBV B95-8 reference genome (RefSeq accession no. NC_007605.1) using the Burrows–Wheeler Aligner (v0.7.17) 65 . Duplicated reads were removed by Picard (v2.18.14) 66 . Samples with insufficient sequencing coverage (less than 90% coverage at 30× depth) or multiple EBV infections were excluded. Multiple infections were defined as samples with more than 10% biallelic or multi-allelic variants. In addition, VerifyBamID (v1.1.3) 67 was used to detect mixed infections, using a cut-off (Chipmix > 0.01681) determined from the average Chipmix values for simulated 1.5% mixed infections. A total of 734 samples (324 patients with NPC and 410 healthy donors) with available host genome information in sample sets 1 and 2 passed the EBV genome quality control criteria and were included in genome-wide EBV–HLA interaction analyses.

Following the GATK best practice workflows (v4.1.8.1), variants were identified after base quality score recalibration 68 . To avoid inaccurate calling, we removed low-quality variants using GATK recommended hard-filtering thresholds (threshold for SNPs: ‘QD < 2.0’, ‘QUAL < 30.0’, ‘SOR > 3.0’, ‘FS > 60.0’, ‘MQ < 40.0’, ‘MQRankSum < −12.5’ and ‘ReadPosRankSum < −8.0’; threshold for indels: ‘QD < 2.0’, ‘QUAL < 30.0’, ‘FS > 200.0’ and ‘ReadPosRankSum < −20.0’). After filtering for missingness < 10% and MAF > 10%, 1,942 variants were retained for further analysis (Supplementary Table 9 and Supplementary Fig. 4). Moreover, in the genome-wide EBV–HLA interaction analysis, stratification was performed based on disease status and the investigated HLA-A allele, and only EBV variants with MAF > 5% in each stratum were included in the analysis. The functional annotation of the EBV variants was performed using the SNPEff package (v5.0e) according to the reference genome (NC_007605.1, NCBI annotation, November 2013) 69 .

In the genome-wide EBV–HLA interaction analysis, EBV principal components were incorporated as covariates to adjust for viral population stratification. For principal component analysis (PCA), variants were filtered based on MAF > 10% and linkage disequilibrium pruning with a pairwise correlation r 2 > 0.5 within a 1,000-SNP sliding window with a 5-SNP step size. PCA was performed using PLINK (v1.9) 60 .

Genotyping of EBV variants by MassArray iPLEX

EBV genotypes, including the BALF2-CCT subtype SNPs 162215A>C, 162476T>C and 163364C>T, were determined using customized primers and the Agena Bioscience MassArray iPLEX platform, following the manufacturer’s protocol, for sample sets 1 and 2 (refs. 6 , 70 ). Using Agena Bioscience MassArray iPLEX platform, 37 EBV markers were genotyped for viral lineage inference, including BALF2-CCT SNPs 162215A>C, 162476T>C and 163364C>T, EBER2 SNP 7048A>C and EBNA3B SNP 84414G>T, in 406 cases and 597 controls from sample set 1 (see Supplementary Table 19). Imputation of SNP 85841 was performed using Beagle (v5.4) 71 , with a reference panel of 667 southern China EBV genomes (734 total, excluding 67 used for imputation accuracy evaluation). Among 67 samples with both EBV genotyping and WGS data, the imputation accuracy for SNP 85841 reached 94.03%.

Stepwise human–EBV genome-wide interaction analysis

The near-universal prevalence and large genome size of EBV pose major challenges for genome-to-genome interaction analyses, as adequately powered case–control studies would require prohibitively large samples with paired host–viral genomes. To address these limitations, we applied a stepwise analytical framework. In step 1, we scanned the human genome to identify variants showing interaction signals with the high-risk EBV subtype. In step 2, we fine-mapped the EBV genome to pinpoint viral variants driving the observed human–EBV interactions.

We adopted an approach for genome-wide interaction analysis recently developed by Zhu et al. 21 . This method uses two logistic regression models: one to estimate the overall genetic effects, and the other to estimate the interaction-independent genetic effects. By testing the difference between the overall genetic effect and the genetic effect excluding interaction, it allows robust detection of genome × EBV interaction effects. This stepwise genome-wide interaction analysis is detailed as follows.

Step 1: host genome-wide scan for interaction with the high-risk EBV

In this step, we tested whether the high-risk EBV subtype interacts with host genetic variants (additive coding) to influence NPC susceptibility. To evaluate this interaction, we calculated the difference between the overall and interaction-independent genetic effects, which is theoretically equivalent to testing the combined effect of genome × EBV interaction and mediation 21 . To disentangle interaction from mediation, we evaluated the association between genome and EBV and found no evidence of genetic mediation at a genome-wide significance threshold of 5 × 10 −7 (Supplementary Fig. 10). Thus, statistically significant loci identified in this analysis were interpreted as showing interaction effects.

In a human genome-wide scan for the interaction effects with high-risk EBV subtype, a GWAS on NPC was conducted to examine the genetic (G) contribution to NPC (Y) using a logistic regression model, adjusting for covariates (C) such as sex, age and the first four host MDS components to account for population structure:

(1)

Here, alpha 1 represents the overall effect of G. Then, a genome-wide G × EBV interaction analysis on NPC risk was modelled through the following logistic regression:

(2)

In this model, beta 1 and beta 3 correspond to the interaction-independent effect of G to differentiate from the overall genetic effect alpha 1 in equation (1) and the interaction effect of G × EBV, respectively. The terms G, EBV and G × EBV represent the host genotype, the infection of EBV subtype and their interaction, respectively.

Then, we tested the G × EBV interaction by the method proposed by Zhu et al. 21 :

(3)

where theta reflects the contribution of interaction-independent effect to overall effect. We estimated the causal effect theta by conducting a linear regression model that assessed the normalized interaction-independent effect’s contribution to the normalized overall effect, denoting as widehattheta:

(4)

where n is the sample size, and epsilon is the random noise. In the host genome study, theta in equation (4) was estimated by all host genetic variants. The genetic variants without G × EBV interaction will fall on the regression line but the variants with G × EBV interaction will depart from this line. Therefore, we searched the host genetic variants that deviate from this regression line to test the effect of G × EBV interaction.

Host genome-wide interaction analysis was conducted using the iterative Mendelian randomization and pleiotropy (IMRP) approach 72 . The statistic T interaction in equation (3) provides a statistical assessment of the significance of G × EBV interaction effects. To account for the effect of sample overlapping between the GWAS and interaction analysis, the IMRP method incorporated a correlation coefficient of standardized alpha 1 and beta 1, calculated from genome-wide variants lacking significant associations (P > 0.05). Following causal effect theta estimation, a genome-wide interaction test, conceptually analogous to pleiotropy testing in IMRP, was applied. The quantile–quantile plot indicates effective control of type 1 error in human genome-wide interaction analysis (Extended Data Fig. 2a–c). For conditional analyses, the GWAS and interaction analyses were performed by incorporating the top variant and the EBV–interaction term as covariates in two logistic regression models.

Step 2: EBV genome-wide scan for interaction with HLA-A*11:01

In this step, we tested whether EBV variants (binary coding) interact with HLA-A*11:01 (presence or absence) to influence NPC risk. To evaluate this interaction, an EBV GWAS was conducted to examine the EBV variant (G EBV ) contribution to NPC (Y) using a generalized linear mixed model. Specifically, to control for viral population structure, we constructed a viral genetic relatedness matrix (vGRM) from EBV variants after linkage disequilibrium-based pruning (sliding window of 1,000 SNPs, step of 1 SNP, r 2 < 0.4) using the GEMMA software (v0.98.5) 73 , following the recommended practice and previous work on highly linked genomes 6 , 74 , 75 . We also performed sensitivity analyses using alternative linkage disequilibrium-pruning threshold r 2 < 0.3 or 0.6, and the top-ranked interaction signal for EBV 85841G remained unchanged. We then fit a logistic mixed model in the GMMAT software (v1.4.2) 76 , using this vGRM as a random effect (Z) and including age, sex, the top four viral principal components and the top four host MDS components as fixed effects (C) to account for viral and host population structure (equation (5)):

(5)

Here alpha 1 EBV represents the overall effect of G EBV . The random effect Z was assumed to follow a multivariate normal distribution with mean zero and covariance matrix proportional to the viral genetic relatedness matrix, that is, Z sim N(0,sigma EBV 2 vGRM), to account for genome-wide EBV genetic similarity among samples. Then, a viral genome-wide G EBV × HLA interaction analysis on NPC risk was modelled through the following logistic mixed model:

(6)

In equation (6), beta 1 EBV and beta 3 EBV correspond to the interaction-independent effect of G EBV to differentiate from the overall EBV effect alpha 1 EBV in equation (5) and the interaction effect of G EBV × HLA, respectively. The terms G EBV, HLA and G EBV * HLA represent the EBV genotype, the HLA variant and their interaction, respectively.

Then, we tested the G EBV × HLA interaction by the statistic proposed by Zhu et al. 21 :

(7)

where theta EBV reflects the contribution of interaction-independent effect to overall effect. We estimated the causal effect theta EBV by conducting a linear regression model that assesses the normalized interaction-independent effect’s contribution to the normalized overall effect, denoting as hattheta EBV:

(8)

where n is the sample size and epsilon is the random noise. In the EBV genome study, theta EBV in equation (8) was estimated using independent EBV variants, which were selected based on a linkage disequilibrium (r 2 ) threshold of 0.3 across the whole EBV genome. We searched the EBV variants that deviate from this regression line to test the effect of G EBV × HLA interaction.

The EBV genome interaction test was performed by the IMRP 72 . The statistic T interaction in equation (7) is a test for evaluating the effect of G EBV × HLA interaction. After estimating the causal effect theta EBV, we performed the interaction test for all EBV variants. In the conditional analysis of the interaction test, the top EBV variant was included as covariates in two logistic mixed models.

To control for multiple testing and type 1 error, we used the method previously presented 27 . BRASS first fit a generalized linear mixed model that included host covariates (age, sex and host MDS components), viral principal components and vGRM to control the viral population structure and sample relatedness, excluding the viral SNP, host allele and their interaction (the null model). We computed the residuals and decorrelated them to obtain approximately exchangeable residuals. We then permuted these decorrelated residuals 10,000 times, and re-fit the interaction model on each replicate. This approach preserves the host–viral relatedness structure and accounts for clonality and relatedness to the viral genome while breaking only the association of interest. Using the minimum P value across the EBV genome in each permutation, we set the empirical genome-wide threshold at α = 0.05, yielding an empirical genome-wide threshold of 5.92 × 10 −4 for the EBV interaction scan.

Power simulation

To quantify power under this stepwise design, we performed simulations following a previous method 77 . In step 1, we resampled the discovery dataset 1,000 times and simulated NPC case–control status using the estimated main and interaction effects of the top interaction signal rs417162 (tagging HLA-A*11:01; r 2 = 0.73), high-risk EBV (BALF2-CCT) and their interaction. The power to rediscover the top interaction at the genome-wide significance threshold (P < 5 × 10 −8 ) was 83%. In step 2, we randomly resampled individuals and their EBV genomes from the paired host–EBV genome dataset and simulated NPC case–control status using the estimated effects of HLA-A allele, EBV SNP 85841 and their interaction. Across 1,000 simulations, the power to detect the HLA-A allele × 85841G interaction at the Bonferroni-corrected significant threshold for the EBV genome scan was estimated. The power to detect the A*11:01 × 85841G and A*02:07 × 85841G interactions was 91% and 72%, respectively. The results of power calculations are consistent with previous human G × G power analyses and support that the stepwise design is statistically efficient under current sample-size realities 77 , 78 , 79 , 80 .

Relative excess risk and attributable fraction due to interaction

RERI

RERI quantifies the magnitude to which the joint effect of two factors (here HLA and EBV) exceeds the sum of their individual effects, providing population-level metrics essential for evaluating absolute risk differences in preventive interventions 28 , 81 . The OR for NPC associated with the joint status of the presence of HLA risk variant (HLA = 1) and the high-risk EBV subtype (EBV = 1) was defined as OR 11. We estimated the interaction effect between the HLA variant and the EBV subtype infection on NPC risk as the relative excess RERI 28 . We first fitted the following logistic regression model for NPC:

(9)

Here, Y=1,0 represents the NPC case or control status, HLA=1,0 indicates the presence or absence of HLA variant, and EBV=1,0 denotes the presence or absence of EBV subtype, and C represents the set of covariates, as indicated in each table or figure legends. Because NPC is a rare disease, with an incidence of approximately 20 per 100,000 person-years in southern China 82 , 83 , we used OR as an approximation of the relative risk in the interaction analyses. Under the rare outcome assumption, then

(10)

The 95% CI for the relative excess RERI was analysed using bootstrap resampling (5,000 iterations) implemented in the R package boot (v1.3.30) or delta method. The P value for the relative excess RERI was calculated using the method proposed by VanderWeele et al. 81 , which uses the low-risk genotype and low-risk exposure as the reference group and applies a one-tailed test.

Attributable fraction due to interaction

The attributable fraction due to interaction quantifies the proportion of disease risk, associated with a given exposure, that can be explained by its interaction with another factor. We estimated the proportion of HLA-associated risk that is attributable to its interaction with high-risk EBV strains (Extended Data Fig. 2h). The OR of NPC (Y = 1) associated with the joint status of HLA risk genotype (HLA = 1) and the high-risk EBV subtype (EBV = 1) was defined as OR 11, which was calculated with a logistic regression model of NPC status (Y) on the HLA variants (HLA), the EBV subtypes (EBV) and their interactions (HLA × EBV), adjusting for age, sex and the top four host MDS components. The NPC risk due to an HLA risk genotype was calculated with OR 10 -1, whereas the NPC risk due to a high-risk EBV subtype was calculated with OR 01 -1. The increased NPC risk due to the coexistence of both HLA risk genotype and high-risk EBV subtype was OR 11 -1. The frequency of the high-risk EBV subtype in healthy controls is denoted as P(EBV=1), whereas the frequency of HLA risk genotype in healthy individuals is represented as P(HLA=1). Accordingly, the proportion of HLA-associated risk attributable to interaction was calculated as frac(OR 11 -OR 10 -OR 01 +1)P(EBV=1)(OR 10 -1)+(OR 11 -OR 10 -OR 01 +1)P(EBV=1) 84 . The 95% CI for the attributable fraction due to interaction was estimated via bootstrap resampling (5,000 iterations).

Population attributable fraction and NPC incidence estimation

PAF

PAF estimates the attributable fraction for a binary outcome (here NPC Y = 1) under the hypothetical scenario where the risk factor (X) is eliminated from the population and corresponds to the attributable fraction of disease outcome explained by the risk factor. The PAFs explained by the effects of HLA-A alleles (presence or absence), EBV subtypes (presence or absence) or their combinations (Fig. 4b,c, Extended data Fig. 2h and Supplementary Tables 15 and 18) were estimated using logistic regression models adjusting for covariates (C), including age, sex and the top four human MDS components, with the R package AF (v0.1.5) 85 .

(11)

In equation (11), gamma 1 corresponds to the effect of the risk factor (X; HLA-A alleles or EBV subtypes). Because NPC is a rare disease 82 , 83 , the relative risk can be approximated by the OR. The attributable fraction of NPC explained by the risk factor can be approximated by

(12)

NPC incidence estimation across EBV and HLA strata

Age-standardized incidence rates (ASRs; per 100,000) for NPC in Asia were sourced from: (1) the 2020 China Cancer Registry Annual Report 86 , (2) the 2021 Hong Kong Cancer Registry 87 , (3) the WHO International Agency for Research on Cancer databases 3 , and (4) peer-reviewed literature 88 . HLA-A allele frequencies were obtained from three resources: our dataset, the NyuWa Genome Resource and published studies 89 , 90 , 91 , 92 . To evaluate population burden across strata defined by EBV strains and HLA-A background in southern China, we estimated stratum-specific NPC incidence and preventable cases. NPC incidence rates for each stratum were derived based on stratum-specific odds ratio to the overall regional NPC incidence of approximate 20 per 100,000 person-years 82 , 83 , under the rare-disease approximation. We further translated the PAF into the estimated number of attributable cases under the counterfactual assumption 46 , 47 that exposure to high-risk EBV strains could be eliminated. For each stratum i, preventable cases were computed as Preventable cases i = PAF i × C South , where C South = 51,000 × 0.80 reflects an annual national NPC burden of approximately 51,000 cases 93 and approximately 80% occurring in southern China.

Peptide–HLA-binding assay

Prediction of HLA-A allele-binding affinity

Peptide–HLA-A*11:01-binding or A*02:07-binding affinities were predicted by NetMHCpan 4.1 (ref. 35 ). HLA–peptide-binding predictions were determined based on percentile ranks for 9–11-mer peptides spanning the EBV SNP 85841A>G-encoded Q900R variant in EBNA3B from high-risk 85841G and non-high-risk EBV strains. Peptides with normalized percentile ranks of less than 2% were defined as binders to the respective HLA-A allele (Supplementary Tables 20 and 21).

AlphaFold3-based structure prediction

The structures of peptide–HLA-A*11:01 complexes for the low-risk peptide VVILENVGQ (85841A encoded) and the high-risk peptide VVILENVSR (85841G encoded) were predicted using the AlphaFold3 public web server 94 . Inputs comprised the HLA-A*11:01 heavy chain, β 2 M and the EBNA3B peptide sequence, using the HLA-A*11:01–β 2 M scaffold derived from Protein Data Bank 6JOZ as a structural ref. 95 . Molecular graphics were generated in PyMOL (v3.1.6.1) 96 .

Biolayer interferometry assay

Peptides with a purity of more than 95% were synthesized by GenScript Biotech. Biotinylated UV-exchangeable HLA-A*11:01–RVFA-J-SFIK and HLA-A*02:07–LLDSD-J-ERL peptide–MHC monomers were purchased from BetterGen. Biolayer interferometry was performed using an Octet R8 (Sartorius), and data were analysed with Octet Analysis Studio (v12.2) software. After UV irradiation for 20 min on ice to degrade the original A*11:01-bound peptide RVFA-J-SFIK, biotinylated HLA-A*11:01 was loaded onto SSA Biosensors (Sartorius 18-5057) at a concentration of 10 μg ml −1 for 480 s. Coated sensor tips were dipped into dilution buffer for a 60-s baseline measurement, then into peptide solutions at serially diluted concentrations for a 90-s association phase, and finally into dilution buffer for a 120-s dissociation phase. Peptide concentrations are indicated in each figure. The equilibrium K d was determined by steady-state analysis of data at equilibrium for each peptide concentration, a method suitable for protein–small-molecule interactions exhibiting fast on and off rates and low signal intensity. The HLA-A*11:01-restricted epitope IVTDFSVIK served as the positive control (K d = 54.7 μM), whereas the peptide LLWTLVVLL, which cannot bind to A*11:01, was used as the negative control (K d = 242.7 μM).

HLA-A peptide exchange ELISA

First, 20 µl diluted peptide (400 µM in PBS) and 20 µl HLA-A monomer (0.2 mg ml −1 ) were added to a microcentrifuge tube. UV-mediated peptide exchange was performed on ice for 30 min under 365-nm UV irradiation. The Flex-T Human Class I Peptide Exchange ELISA Kit (BioLegend) was used to evaluate the efficiency of UV-activated peptide exchange according to the manufacturer’s instructions. DMSO was used as the blank control. For the HLA-A*11:01 monomer, IVTDFSVIK was used as the positive control and CLGGLLTMV as the negative control. For the HLA-A*02:07 monomer, CLGGLLTMV was used as the positive control and IVTDFSVIK as the negative control.

T cell functional assay

EBV production

EBV strains B95-8, Akata and M81 were produced by the B95-8 cells, CNE2-Akata cells and CNE2-M81 cells, respectively 50 , 97 , 98 , 99 , 100 . In brief, B95-8 cells were resuspended at a density of 2 × 10 5 cells per millilitre in RPMI 1640 medium with low-serum condition (2% fetal bovine serum (FBS)), 100 U ml −1 penicillin and 100 µg ml −1 streptomycin to trigger the EBV lytic cycle for 2 weeks. To induce CNE2-Akata and CNE2-M81 cells, 20 ng ml −1 12-O-tetradecanoylphorbol 13-acetate (TPA; Beyotime) and 2.5 mM sodium butyrate (NaB; Sigma-Aldrich) were added in RPMI 1640 medium with 10% FBS, 100 U ml −1 penicillin and 100 µg ml −1 streptomycin for 12 h. Then, the medium was replaced with fresh complete medium, and cells were cultured for 3 days. Cell supernatants were collected and filtered by 0.45-µm membrane filter (Millipore). EBV was concentrated by centrifugation at 50,000g for 2.5 h and resuspended in RPMI 1640 medium. Virus stocks were stored at −80 °C until use.

Cell line construction

COS-7 cell lines, each overexpressing a single HLA–β 2 M complex (A*11:01, A*02:07 or A*02:01) and luciferase, were gifted by the Li laboratory 36 . In addition, expression plasmids encoding EBV EBNA3B haplotypes from the high-risk 85841G, high-risk 85841A and low-risk strains were individually transiently transfected into the HLA-A*11:01 and luciferase co-expressing COS-7 cell line using polyethylenimine (PEI; Polysciences) at a DNA:PEI ratio of 1:3, followed by incubation for 24 h.

To establish EBV-infected LCLs, PBMCs were resuspended in B95-8–EBV-containing, Akata–EBV-containing or M81–EBV-containing solutions at a density of 5 × 10 5 cells per millilitre, respectively 101 . After 12 h, cells were centrifuged and cultured at 48-well plates in RPMI 1640 medium with 10% FBS, 1% (v/v) GlutaMax (Gibco), 1 mM sodium pyruvate (Gibco), 2 µg ml −1 cyclosporin A (CsA; Yeasen), 100 U ml −1 penicillin and 100 µg ml −1 streptomycin. Cell culture medium was replaced weekly (50% volume) with fresh complete medium containing 2 µg ml −1 CsA. After approximately 4–6 weeks, the LCLs were stable and proliferated continuously, and were used for subsequent experiments. Successful establishment of LCLs was confirmed by EBER fluorescence in situ hybridization (EBER-FISH) and flow cytometric analysis of B cell markers.

Peptide-specific T cell expansion

Peptide-specific T cells were generated as previously described 102 , 103 . In brief, PBMCs (1 × 10 6 cells) were cultured in 24-well plates containing complete RPMI 1640 medium, consisting of RPMI 1640 medium supplemented with 10% FBS, 100 U ml −1 penicillin, 100 µg ml −1 streptomycin, 10 mM HEPES (Gibco), 1 mM sodium pyruvate (Gibco), 1% (v/v) non-essential amino acids (Gibco), 1% (v/v) GlutaMax (Gibco) and 50 µM β-mercaptoethanol (Sigma-Aldrich). On days 1 and 5, PBMCs (5 × 10 5 cells) were incubated with single peptide (10 µM) or a peptide pool (2 µM per peptide) for 2 h at 37 °C in 5% CO 2 . Peptide-pulsed PBMCs were washed once and then mixed with the autologous PBMC cultures. After 3 days of cell culture, cultures with 10 U ml −1 recombinant human IL-2 (PeproTech), 10 ng ml −1 recombinant human IL-7 (Bio-Techne) and 10 ng ml −1 recombinant human IL-15 (PeproTech) were added, and half of the medium was replaced every 2 days thereafter. On day 13, expanded T cells were washed and rested overnight in cytokine-free medium.

T cell activation assay (intracellular cytokine staining)

HLA-A-expressing COS-7 cell lines (used as APCs) were pulsed with peptides in 96-well U-bottom plates and incubated for 2 h at 37 °C. After incubation, cells were washed twice and then resuspended in complete RPMI 1640 medium for subsequent activation and cytotoxicity assays. Peptide-expanded T cells were co-cultured with peptide-pulsed APCs, EBV EBNA3B-transduced APCs or LCLs in complete RPMI 1640 medium supplemented with anti-human CD28 (1 µg ml −1 ; CD28.2, BD Biosciences) and anti-human CD49d (1 µg ml −1 ; 9F10, BD Biosciences). Co-cultures were set up in 96-well U-bottom plates at an effector-to-target ratio of 1:5 (5 × 10 4 effector cells and 2.5 × 10 5 target cells per well) and incubated at 37 °C for 2 h. Brefeldin A (10 µg ml −1 ) and monensin (1 µM) were then added, followed by incubation at 37 °C for 6 h. Cells were subsequently kept overnight at 4 °C. Surface staining was performed for 1 h at 4 °C using the following reagents diluted at 1:100, including Fixable Viability Dye eFluor 506 (eBioscience), APC/cyanine7 anti-human CD45 (HI30, BioLegend), PE-Cy7 anti-human CD3 (SP34-2, BD Biosciences) and PerCP-Cy5.5 anti-human CD8 (RPA-T8, BD Biosciences). Cells were then washed, fixed and permeabilized for 1–5 h, followed by intracellular staining with Alexa Fluor 647 anti-human IFNγ (BD Biosciences; Supplementary Fig. 11a). Intracellular cytokine staining was performed using a CytoFLEX LX flow cytometer (Beckman Coulter). Data were analysed using CytExpert (Beckman Coulter) and FlowJo (v10; TreeStar).

T cell cytotoxicity assay

Peptide-expanded T cells were co-cultured with peptide-pulsed APCs, EBV EBNA3B-transduced APCs or LCLs in complete RPMI 1640 medium in 96-well U-bottom plates at an effector-to-target ratio of 5:1 (2.5 × 10 4 effector cells and 5 × 10 3 target cells per well). Co-cultures were incubated for 24 or 48 h at 37 °C in 5% CO 2 . For luciferase-expressing COS-7 APCs, T cell-mediated cytotoxicity was assessed by measuring the luciferase activity of the remaining viable target cells. After co-incubation, plates were centrifuged at 1,000 rpm for 5 min and washed twice to remove dead cells and cellular debris. Cells were then resuspended in culture medium, mixed with an equal volume of Bio-Lite Plus detection reagent (Bio-Lite Plus Luciferase Assay System, Vazyme) and incubated at room temperature for at least 3 min. Luminescence was measured using a multimode microplate reader (Spark 10 M, Tecan). Percent specific lysis was calculated using the following formula: Specific lysis (%) = 100 × (1 – test luminescence/maximum luminescence).

T cell-mediated cytotoxicity against LCL targets was assessed by flow cytometric analysis of target cell death. After co-incubation of LCLs and peptide-expanded T cells, cells were harvested and washed twice with PBS and then incubated with BV421 anti-human CD19 (HIB19, BioLegend) antibody for 30 min at 4 °C. Then, cells were washed with binding buffer (BioLegend) three times and stained with APC Annexin V Apoptosis Detection Kit (BioLegend) together with SYTOX AADvanced (Invitrogen) according to the manufacturers’ instructions (Supplementary Fig. 11b). LCLs co-cultured without peptide-expanded T cells served as the baseline control. Samples were analysed on a CytoFLEX LX flow cytometer (Beckman Coulter). Percent specific lysis was calculated as Specific lysis (%) = 100 × (1 – cell survival in test group/cell survival in control group).

IFNγ ELISPOT assay

The IFNγ enzyme-linked immunospot (ELISPOT) assay was conducted following the manufacturer’s instructions (Human IFNγ ELISPOT Set, BD Biosciences). In brief, ELISPOT plates were coated with the capture antibody (NA/LE purified anti-human IFNγ) overnight at 4 °C. After washing once, plates were blocked with complete RPMI 1640 medium for 2 h at room temperature. Peptide-expanded T cells (1 × 10 6 cells per millilitre) and single peptide (2 µM) or peptide pools (1 µM per peptide) in complete RPMI 1640 medium were plated and incubated at 37 °C with 5% CO 2 for 18–22 h. DMSO was used as a negative control. After incubation, the detection antibody (biotinylated anti-human IFNγ) was added and incubated at room temperature for 2 h, followed by incubation with streptavidin–horseradish peroxidase for 1 h at room temperature. The 3-amino-9-ethylcarbazole substrate (AEC Substrate Set, BD Biosciences) was then added to each well for 15–45 min, and deionized water was used to stop the substrate reaction. Spots were enumerated using an AID ELISPOT Reader.

Tetramer staining

PBMCs from patients with NPC or healthy controls were stained for 1 h at 4 °C with a PE-conjugated HLA-A*11:01 tetramer loaded with VVILENVSR (MBL), AVFDRKSDAK (BetterGen) or IVTDFSVIK (BetterGen), and with an APC-conjugated HLA-A*11:01 tetramer loaded with ATAAAAAAK (MBL) as a negative control. Cells were subsequently surface-stained with Fixable Viability Dye eFluor 506 (eBioscience), FITC anti-human CD3 (HIT3a, BioLegend) and PerCP-Cy5.5 anti-human CD8 (RPA-T8, BD Biosciences; Supplementary Fig. 11c). Flow cytometry was performed using a CytoFLEX LX flow cytometer (Beckman Coulter). Data were analysed using CytExpert (Beckman Coulter) and FlowJo (v10; TreeStar).

EBV population genetic, phylogenetic and selection analysis

EBV genome alignment

High-quality EBV WGS data were obtained from 388 cases and 453 controls from southern China, and 113 healthy controls from other areas in China, after EBV genome quality control. In addition, we retrieved raw EBV WGS fastq data for 176 individuals from published studies in the NCBI databases. Applying the same quality control criteria, EBV WGS data from 118 of these individuals were retained and included in subsequent analyses (for details, see Supplementary Table 16 and Supplementary Fig. 7). Variant calling was performed following GATK best practices workflow (see the section ‘EBV variant calling’ for details). Subsequently, consensus genome sequences were generated by integrating sample-specific variants into the reference genome, with ‘N’ indicating a missing value. The proportion of reads supporting the major allele was assessed for each site. If this proportion exceeded 60%, the major allele was assigned as the consensus base. If the supporting read proportion for the major allele was below 60%, this ambiguous site was masked with N to mitigate potential artefacts. A total of 1,072 consensus genome sequences were generated for the following multiple sequence alignment.

In addition, we retrieved 1,130 assembled EBV genomes (more than 100 kb) from the NCBI database (taxon ID: 10376; accessed July 2022), excluding any sequences that overlapped with the 1,072 consensus genomes generated from raw fastq data in this study. As raw sequencing data were unavailable for these public genomes, only the FASTA-format assemblies were used. Metadata, including geographical origin, EBV type and diseases of carriers, were extracted for all sequences (Supplementary Table 16). Together with the 1,072 consensus sequences generated in this study, all 2,202 genomes were included in the subsequent multiple sequence alignment.

Multiple sequence alignment was performed using MAFFT (v7.490) 104 with the ‘–keeplength’ parameter to preserve the original length of the reference genome (171,823 bp) and establish homologous alignment. To ensure high alignment quality, we first masked repetitive regions (annotated in NC_007605.1, 20.7% of the EBV genome) and low-coverage regions defined as more than 50% N characters across all alignments in a 100-bp window. Subsequently, sequences with missing values exceeding 25% were excluded from further analysis. Following these quality control steps, a total of 2,086 EBV genomes were retained for downstream analyses.

Tree construction

EBV has been classified into two major types, type 1 and type 2, which diverged early in evolution and exhibit distinct genomic features, particularly in the EBNA2, EBNA3A, EBNA3B and EBNA3C gene regions, with type 1 representing the predominant lineage globally. To account for this deep evolutionary split, we first distinguished type 1 and type 2 strains before downstream analyses. We performed PCA based on genome-wide variation from 2,086 EBV strains. Variants were identified based on the multiple sequence alignment, and variants with high levels of missing data (more than 10%) or a MAF below 5% were removed using VCFtools (v0.1.13) 105 . PCA was subsequently performed with PLINK (v1.9) 60 and identified two genetically distinct EBV clusters. These clusters were primarily differentiated by variation within the EBNA2 and EBNA3 genes, corresponding to the EBV type 1 and type 2 lineages (Supplementary Fig. 12). Given the clear separation and the global predominance of type 1, subsequent analyses were restricted to this group (n = 1,852).

The maximum likelihood of the phylogenetic relationship was inferred using IQ-TREE 2 (v2.2.0) 106 . The optimal model of nucleotide substitution was selected using the ModelFinder function within IQ-TREE 2, based on the Bayesian Information Criterion 107 . The model TVM + F + I + R7 was identified as the best fitting and was used for final tree construction, with statistical support assessed using 1,000 non-parametric bootstrap pseudoreplicates. To evaluate the robustness of phylogenetic inference to recombination, we performed a recombination-masked analysis by excluding the inferred N1-derived recombinant region (approximately 54–130 kb) from the multiple sequence alignment and reconstructing the maximum likelihood tree using the same IQ-TREE 2 workflow.

The phylogeny was rooted using an African-derived EBV strain (KP968262.1), as African lineages are known to occupy basal positions in global EBV trees, consistent with the co-evolutionary history of EBV and its human host originating in Africa. Rooting was performed using the root function in R package ape (v5.8.1) 108 . All phylogenetic trees were subsequently visualized and annotated using R packages ape (v5.8.1), treeio (v1.30.0), ggtree (v3.14.0) and ggtreeExtra (v1.16.0) 108 , 109 , 110 , 111 .

Chromosome painting analysis

To investigate the ancestry composition and detect potential admixture signals in EBV strains, we applied a chromosome painting approach using ChromoPainter, embedded within the fineSTRUCTURE (v4.1.0) 112 . ChromoPainter uses a hidden Markov model to reconstruct each recipient haplotype as a mosaic of genomic segments copied from the donor haplotypes. A total of 132 EBV genomes from northern China (N1) and 116 from southern China (S1) were designated as donors to paint 643 recipient genomes (C1 and X1).

As ChromoPainter requires phased haplotypes, we first converted filtered VCF files into ChromoPainter-compatible.phase format using PLINK (v1.9) 60 and a companion script supplied by ChromoPainter (plink2chromopainter). ChromoPainter was first run in estimation mode (-a 0 0) using all samples to infer the effective population size (Ne) and mutation rate (μ). In the subsequent painting step, these estimated parameters were applied, with the number of expectation-maximization iterations set to 10, and the -k parameter set to 20, defining the expected number of chunks copied from donors.

Genomic similarity and recombination analysis

Genetic similarity between consensus genome sequence of 85841G-carrying high-risk strains and other strains on the phylogeny was accessed using a sliding window approach (1,000-bp window size and 100-bp step size). To detect recombination, we generated consensus genome sequences for each focal clade (for example, S1, N1, X1 and C1) using EMBOSS (v6.6.0.0) 113 , and conducted recombination analysis on these consensus sequences using RDP5 (ref. 114 ). Genetic differentiation between clades was evaluated by calculating the fixation index (FST, Weir and Cockerham’s method), implemented in VCFtools (v0.1.13) 105 , with the same window parameters.

Extended haplotype homozygosity test

In healthy individuals from southern China, extended haplotype homozygosity spanning EBV variants (that is, SNP 85841, SNP 84414, SNP 84462, SNP 7048 and SNP 163364) was estimated using the R package rehh (v3.2.2) 115 .

登录APP查看全文