A comprehensive analysis of genetic diversity of EBV reveals potential high-risk subtypes associated with nasopharyngeal carcinoma in China

Abstract Epstein–Barr virus (EBV), a widespread oncovirus, is associated with multiple cancers including nasopharyngeal carcinoma (NPC), gastric cancer and diverse lymphoid malignancies. Recent studies reveal that specific EBV strains or subtypes are associated with NPC development in endemic regions. However, these NPC specific subtypes were only identified in a portion of infected individuals due possibly to the limited samples size studied or the complicated population structures of the virus. To identify additional high-risk EBV subtypes, we conducted a comprehensive genetic analysis of 22 critical viral proteins by using the largest dataset of 628 EBV genomes and 792 sequences of single target genes/proteins from GenBank. The phylogenetic, principal component and genetic structure analyses of these viral proteins were performed through worldwide populations. In addition to the general Asia-Western/Africa geographic segregation, population structure analysis showed a ‘Chinese-unique’ cluster (96.57% isolates from China) was highly enriched in the NPC patients, compared to the healthy individuals (89.6% vs. 44.5%, P < 0.001). The newly identified EBV subtypes, which contains four Chinese-specific NPC-associated amino acid substitutions (BALF2 V317M, BNRF1 G696R, V1222I and RPMS1 D51E), showed a robust positive association with the risk of NPC in China (Odds Ratio = 4.80, 20.00, 18.24 and 32.00 for 1, 2, 3 and 4 substitutions, respectively, Ptrend <0.001). Interestingly, the coincidence of positively selected sites with NPC-associated substitutions suggests that adaptive nonsynonymous mutation on critical proteins, such as BNRF1, may interact with host immune system and contribute to the carcinogenesis of NPC. Our findings provide a comprehensive overview of EBV genetic structure for worldwide populations and offer novel clues to EBV carcinogenesis from the aspect of evolution.


Introduction
Epstein-Barr virus (EBV) belongs to the gamma-herpesviruses family and is ubiquitous in adults across the world. To some extent, it can be considered as the most successful oncovirus. After primary infection usually occurred before teenagers, EBV stays largely dormant in its life-long persistent infection. It can be reactivated under certain condition (Xu et al. 2012;He et al. 2019), which may be related to a wide spectrum of malignant diseases, including epithelial tumors like nasopharyngeal carcinoma (NPC), gastric cancer and diverse lymphoid malignancies such as Burkitt's lymphoma (BL), Hodgkin's lymphoma (HL) and NK/T-cell tumors (Young and Rickinson 2004). These EBVassociated malignancies show a great disparity of geographical distribution. For example, NPC is extremely high prevalent in South China and Southeast Asia, while BL is highly endemic in Africa and EBV-associated NK/T-cell lymphoma is relatively prevalent in East Asia, especially in Japan and Korea (Shannon-Lowe and Rickinson 2019).
Several factors such as environmental agents, behavior factors and human genetic susceptibility have been associated with the substantial geographic disparity of NPC. The role of EBV genetic diversity on NPC development, which has long been speculated, however, was only depicted recently with the identification of several geographically specific high-risk EBV subtypes for NPC. For example, individuals infected with EBV subtype BALF2-CCT showed 11 times higher risk in NPC development than the carriers of low-risk BALF2-ATC in South China (Xu et al. 2019). In Hong Kong population, the high-risk EBV subtypes are identified by the genetic variation in EBV-encoded small RNA locus (Hui et al. 2019) as well as the genome-wide single nucleotide variant profiles (Lam et al. 2020). Other risk EBV subtypes appear to be associated with replication-competent since the identified mutations are located in genes that promote EBV replication, including BZLF1, EBNA1 and BRLF1 (Dheekollu et al. 2017;Bristol et al. 2018;Zhang et al. 2018). However, these NPC-specific EBV subtypes were only identified in a portion of infected individuals due possibly to the limited study samples or the concurrence of geographical variants that are not disease associated (Lin et al. 2019). For instance, a previous identified risk variant for NPC (e.g.G155391A) in China was found to be scarce in NPC from Indonesia (Correia et al. 2018), indicating the complexity of virus genetic populations in the different areas (Wegner et al. 2019). Thus, it is important to clarify the genetic population structure with sufficient EBV strains from both the controls and cases before exploring NPC associated variants/subtypes.
To this end, we retrieved all available unique EBV genomes (N ¼ 628) from GenBank including 231 cancer cases and 397 controls. Since the high-risk EBV subtypes identified to date seem all related with the virus replication, we, therefore, selected 22 genes that are involved in the EBV replication from the EBV genome (Kenney and Mertz 2014;Al Moustafa et al. 2018) for the genetic diversity analysis. These genes include (1) EBNAs, LMP1, LMP2 and BARF1, which inhibit viral replication and function in carcinogenesis of EBV-associated cancers (Chakravorty et al. 2019), such as activating NF-jB signaling (LMP1), promoting Lyn/ Syk kinases (LMP2A) and suppression of apoptosis (BARF1) (Tsang et al. 2019); (2) BZLF1 and BRLF1, two viral immediateearly proteins which are associated with EBV reactivation (Morales-Sanchez and Fuentes-Panana 2018; Wu et al. 2018);and (3) genes that are vital to persistent EBV infection and carcinogenesis through recently pan-cancer analysis (Chakravorty et al. 2019), including A73, BALF2, BALF3, BALF5, BARF0, BZLF1, BRLF1, EBNA1, EBNA3B, EBNA3C and RPMS1. These viral genes were expressed in nearly all cancer type, indicating their functional roles in cancer pathophysiology (Hu et al. 2016;Borozan et al. 2018;Chakravorty et al. 2019;Peng et al. 2019). In this study, we conducted phylogenetic and principal component analyses for all selected EBV genes and identified high-risk EBV subtypes that are associated with NPC in China. Positive selection for candidate proteins provide further support for these NPC associated high-risk subtypes.

Data sources and samples
A total of 781 sequences of entire or partial (>100 kb) EBV genome (released before Feb 2019) were downloaded from GenBank (https://www.ncbi.nlm.nih.gov/genbank). The information including diseases/health condition, geographic origin, sample type was also retrieved. Sequences lacking the information of diseases/health condition and geographic origin were excluded. If multiple isolates were from identical cell lines or individuals, the more complete sequences were included. A total of 22 EBV proteins involving in viral replication and carcinogenesis were identified through literature mining (Supplementary Table S1). The complete coding sequences of these targeted genes were extracted. An additional set of EBV genomic sequences from China released between Feb and Oct in 2019 were acquired to confirm the cancer-related amino acid (AA) substitutions (Supplementary Table S2). The complete coding sequences of single target gene were obtained from GenBank (Supplementary Table S3). Duplicates in the combined dataset were excluded. In addition, a set of 244 sequences for BRLF1 from our previous study were included, which contains 80 NPC cases, 80 controls from Guangdong province, and 84 healthy persons from Shanxi province, China (Zhang et al. 2018).

Phylogenetic and principal component analyses of EBV protein sequences
The coding sequences for each gene were aligned using MUSCLE following the instruction of MEGA-X and translated into AA sequences. The repeated regions in EBNA1 and EBNA3C were removed from the sequence alignments. The sequence from a chimaera of B95-8 and Raji EBV was defined as prototypic reference (Accession number NC007605). Maximum-likelihood phylogenetic tree was generated for each protein by IQ-TREE version1.6 with 1,000 bootstrap value using the best-fit model (Hoang et al. 2018;Zhou et al. 2018). Visualization and annotation of the trees were performed using treeio and ggtree packages in R (Wang et al. 2020;Yu et al. 2017). Principal component analysis (PCA) was conducted for each protein on AA substitutions (excluding repeated regions) using adegenet package (Jombart 2008).

Population structure analysis
Population structure analyses were conducted for different set of proteins. The analyses were performed using variable locus excluding incidental mutations (mutated AA ¼ 1) by Structure software version 2. 3 (Hubisz et al. 2009) with a burn-in period of 10,000 and 10,000 MCMC simulations at five iterations. Estimation for the most likely number of putative clusters follows the direction raised by Evanno et al. (2005) using STRUCTURE HARVESTER web v0.6.94. The results were visualized using DISTRUCT (Rosenberg 2004).

Association analysis and identification of NPCassociated variants
The association analyses between single AA substitution and cancers or geographic origin were conducted using treeWAS (https://github.com/caitiecollins/treeWAS) in R (version 3.6.3) (Collins and Didelot 2018;San et al. 2019). Rare substitutions (minor frequency <0.05) were excluded before the analysis.
To further identify the variants associated with certain cancer such as NPC while not with geographical area, AA substitutions which were major in cancer-related cluster (>50%), while minor in other geographically distinct clusters (<20%) were extracted based on the proportions assigned to each cluster in population structure analyses.
In the present study, a novel mutation (BNRF1 V1222I, loc5399) was validated via nested PCR and Sanger sequencing using mouthwashes from 108 NPC cases and 179 controls in our EPI-NPC-2005 project (Xue et al. 2018). All samples from the cases were obtained before any therapy. The cases and controls were matched by age and gender.

Positive selection analysis
Positive selection analyses of viral proteins in East Asian strains (including the isolates from China, South Korea and Japan) were performed using the HyPhy software (Weaver et al. 2018) on the Datamonkey webserver (www.datamonkey.com). Four different methods were employed to detect the positively selected sites: Single-Likelihood Ancestor Counting, Fixed Effects Likelihood, Mixed Effects Model of Evolution, and Fast, Unconstrained Bayesian Approximation for inferring selection. Positively selected site was inferred and reported if it was predicted by any one of the methods.

Characteristics of the target protein sequences
The final dataset of EBV genome contains 628 sequences that were released from January 2006 to 2019 (Supplementary Table  S2). Nearly 40 per cent sequences were from China (N ¼ 247), including 75 from NPC patients or cell lines (C666, M81), 16 from gastric cancer, 4 from lung cancer and 152 from healthy controls. A total of 164 sequences were from other East Asia countries, mainly Japan and South Korea. The isolates from lymphoma, gastric cancer and non-malignant disorders such as IM, PTLD and CAEBV were 45, 5 and 114, respectively. The sequences from Europe and America account to 54 and 43, respectively, most of which were from non-cancerous diseases (N ¼ 41 for Europe; N ¼ 30 for America). There were 57 sequences from Africa, including 34 from BL, 22 from non-cancerous individuals or cell lines and 1 from NPC. A total of 11 sequences were from Papua New Guinea in Oceania, with one from BL and ten from non-cancerous isolates ( Table 1).

Strong geographic segregation was corroborated in most viral proteins
The phylogenetic tree without branch lengths was built for each protein to detect potential virus population. For most of the proteins (18/22), we observe a strong geographic segregation between Asians and Western populations, which is consistent with a recent study using entire genome sequences of EBV (Wegner et al. 2019). The isolates from China tend to aggregate with those from Japan and South Korea, while sequences from western countries were clustered to those from Africa. The isolates from Southeast Asia were clustered with those from Papua New Guinea (Fig. 1, Supplementary S1-S3). Similar conclusion can be drawn from the PCA results (Fig. 2, Supplementary S1-S3). The isolates from Europe, America and Australia were assigned to the group of "Western world" in the ensuring analysis. This remarkable distinction between Asians and Western populations was also found in the genetic structure analysis (Fig. 3).
Although strong geographic segregation was detected among EBV strains, no cancer-specific clade was observed based on these trees. The EBV sequences from cancer patients or cell lines clustered together with those from non-cancerous isolates in the same geographic area (Fig. 1, Supplementary S1-S4). When the AA substitutions in analyzed proteins were compared between isolates from China and other populations, after correcting for population structure, fewer substitution were associated with geographical areas compared to the NPC-associated ones, especially in EBNA3B, EBNA3C and LMP1 (Supplementary  Table S4), indicating the genetic diversity between these populations were mainly due to the clonal population structure.

Phylogenetic analyses and PCA revealed diverse patterns among different proteins
Although most of the proteins showed similar geographic segregation, inconsistent patterns were observed among different viral proteins. The most distinct difference lies in the relationship between isolates from China and other districts. According to the phylogenetic analysis, twenty-two proteins were assigned into four sets as follows: i. The clustered Chinese set: this set exhibited a cluster composed of the isolates from China, while largely excluding those from other Asians. Take protein BNRF1 as an example, 52.70 per cent of the isolates from China can be delineated as a separate cluster apart from most other Asians (Figs. 1a,2a). Other proteins assigned to this set included: BALF2, BALF4, BALF5, BHRF1, LMP2 and RPMS1 ( Supplementary Fig. S1). ii. The clustered Asian set: six proteins included into this set are BALF3, BNLF2A, BNLF2B, EBNA1, LF1 and LMP1. The typical cluster (BALF3) showed that the isolates from China were mingled with the strains from other Asian isolates, but the segregation of Asia-Western/Africa remains conspicuous (Figs. 1b, 2b; Supplementary Fig. S2). iii. The partitioned Western set: this was represented by the protein BARF0 (Figs. 1c, 2c). The isolates from Western world and Africa were partitioned into different clusters and some were close to those from Asians. This set also included four other proteins: BRLF1, BZLF1, EBNA3B, and EBNA3C ( Supplementary Fig. S3).  iv. The ambiguous set: it contains four proteins, namely A73, BARF1, BILF1 and LF2. The Asia-Western/Africa segregation were ambiguous in these trees ( Supplementary Fig. S4), possible due to the relative conservativeness of the sequences among these genes.

Genetic structure analyses identified a Chinesespecific cluster enriched in the NPC patients
To identify the possible cancer associated with EBV, genetic structure analyses were further performed using the gene sets assigned above. Agreed with phylogenetic results, a 'Chinesespecific' cluster emerged from the first gene set analysis, in which the majority of the Chinese isolates (56.02%) were assigned. Particularly, in cluster A1, 96.57 per cent of isolates were from China (Fig. 3a). For the second gene set (clustered Asian set), 55.50 per cent isolates in cluster B1 were from China (Fig. 3b). In the plot of the third gene set (partitioned Western set), no Chinese-specific cluster was observed (Fig. 3c). The Western populations and Africans were assigned similarly by all three gene sets analysis, with corresponding percentages of 92.38 per cent, 78.22 per cent and 92.30 per cent for cluster A5, B3 and C5, respectively (Fig. 3).
Interestingly, the Chinese-specific cluster A1 was significantly more prevalent in NPC patients than in controls in China (89.6% vs. 44.5%, P < 0.001, Fig. 3a, Supplementary Table S5). However, the other relatively 'district-specific' clusters (the cluster A5, B3 and C5 in the Western populations and Africans) did not exhibit significant divergence between lymphomas and corresponding controls.

NPC risk EBV subtype defined by four substitutions conferred a much higher risk for NPC
To check whether specific EBV subtypes were related with NPC, we performed association study in Chinese-specific cluster A1 using treeWAS, which corrects the confounding effects, such as clonal population structure and recombination, and therefore enhances the statistical power of the association analysis (Collins and Didelot 2018;San et al. 2019). Four AA substitutions (BALF2 V317M, BNRF1 G696R, V1222I and RPMS1 D51E) were significantly associated with NPC ( Supplementary Fig. S5), with Odd ratios (ORs) range from 5.38 to 21.10 (all P < 0.001, Table 2). These results were confirmed by newly released EBV genomes from 162 NPC patients and 38 controls in China, with the ORs range from 2.94 to 8.36 (all P < 0.01, Table 2). Particularly, the   BNRF1 V1222I mutation, which yields the highest OR of 21.10, was further validated using our independent samples from EPI-NPC-2005 project. The adjusted effect size was 5.55 [95% confidence interval (CI): 3.22-9.55, P < 0.001, Supplementary Table  S6). When these substitutions were used as a combined indicator, a significant positive effect was observed in the original dataset (P trend < 0.001). The ORs were 4.80, 20.00, 18.24 and 32.00 for the EBV isolates containing 1, 2, 3 and 4 AA substitutions, compared to the wild-type. These effects remained significant in additional Chinese sequence set, in which the ORs were 21.67 (95% CI: 2.51-187.16) and 7.99 (95% CI: 3.11-20.50) for EBV isolates with three and four AA substitutions, respectively (Table 2).
3.6 The NPC-associated AA substitutions in BNRF1, BALF2 and RPMS1 are under positive selection The entire coding sequences of proteins BNRF1, BALF2 and RPMS1 were tested for positive selection and the results are shown in Table 3. It is intriguing to find that many NPCassociated mutations, including BALF2 V317M and BNRF1 V1222I, were under positive selection. In BNRF1, there were fifteen positively selected sites and seven were found to be associated with the NPC. In BALF2, three positively selected sites (317, 700 and 1093) also associated with NPC (3/8). While one site (50) in RPMS1 was associated with NPC (1/2) ( Table 3, Supplementary S4). It is worth noting that not all detected positive sites are the determinant variants for NPC in China. For example, nonsynonymous mutations at BALF2 codon 1093 and BNRF1 codon 797 showed positive association with NPC in China, but their proportions in non-endemic areas of NPC (East Asia, the Western world and Africa) achieved to 36-96 per cent, suggesting these could not be the determinant mutations associated with Chinese NPC. The substitutions at other seven sites showed inverse association with NPC in China. And most of these substitutions had relatively higher frequencies in nonendemic areas, which was consistent with the former result.

Discussion
In this study, we identified potential NPC-risk EBV subtypes in China through the extensive genetic variation analysis. With the rapid development of sequencing technologies and accumulating viral sequences, inferring the EBV inter-and intrapopulations genetic structure across the world are only possible in recent years. Here, we performed, by far, the largest EBV population structure analysis, and detected its geographic segregation, as well as the disparity of the mutations among different proteins. Importantly, a set of proteins generated a 'Chinesespecific' genetic cluster that was enriched in NPC. From this cluster, four NPC-associated AA substitutions located on three proteins (BALF2, BNRF1 and RPMS1) were identified. A combined indicator, which contains the NPC risk subtype, yields the highest risk estimation (OR) of 32.00. In addition, nearly half of the positively selected sites on these proteins were overlapped with the NPC-associated substitutions, implying that adaptive mutations of EBV proteins in Chinese population may drive the development of NPC.
A geographic segregation was observed for most of the viral protein in this study. The isolates from East Asians clustered distantly to those from the Western populations and Africans. Similar patterns have been reported by the studies using 127-270 entire genome sequences of EBV (Palser et al. 2015;Chiara et al. 2016;Borozan et al. 2018;Bridges et al. 2019;Wegner et al. 2019;Xu et al. 2019). Such highly structured virus population may relate with the host immune selection. As suggested in recent study, the genetically homogenous host populations (like Chinese and other Asians, and among Western populations) may force the virus to evolve along similar trajectories (Wegner et al. 2019). Long-time hybridization with immigrants, recombination between individuals and genetic drift also contribute to its distinct yet connected population structure (Chiara et al. 2016). Interestingly, our results highlight the subtle divergence of genetic structure among different viral proteins. About one third of the proteins like BNRF1 and BALF2 generate a 'Chineseunique' cluster, in which a large proportion of isolates were clustered separately apart from other Asians. One recent study also showed the different tree topologies among nine EBV genes using 188 sequences (Zanella et al. 2019), but the geographic discrepancy has not been mentioned. The theory of co-evolution of proteins under host immune selection may give clues to this divergence. Different proteins of EBV may induce varied levels and different types of host adaptive immune response (Taylor et al. 2015). Researches on non-cancerous carriers of EBV found that lytic proteins, mainly BZLF1 and BRLF1, and latent proteins would normally induce frequent and intensified CD4 þ T-cell responses, while the viral proteins such as BZLF1, BRLF1, BALF2, BALF4, BALF5 and EBNA3s induce a relatively large CD8 þ T-cell response. Similarly, a small-size but stable CD8 þ T-cell response could be induced by BNRF1 (Taylor et al. 2015). The nonsynonymous mutations in these immunogenic proteins showed strong linkage disequilibrium with each other, indicating co-evolution of these proteins (Wegner et al. 2019). Thus, we speculate that the distinction of immune selection pressure represented by MHC molecules from various populations may result in the population structure disparity among viral proteins, as indicated in HIV studies (Moore et al. 2002;Bhattacharya et al. 2007).
A Chinese-unique cluster A1 is highly enriched in the NPC patients but also exists in the controls (89.6% vs. 44.5%), which is consistent with previous studies in identification of NPCassociated variants (Hui et al. 2019;Xu et al. 2019). These potential risk subtype of EBV may partly contribute to the carcinogenesis. Here, we identified four NPC-associated AA substitutions located on BALF2, BNRF1 and RPMS1 from the cluster A1. These substitutions remain significantly associated with NPC in Chinese after correction of confounding factors, such as the population structure and recombination, indicating the highly potential of the causal association. Although two variates (BALF2 V317M and RPMS1 D51E) had been reported (Feng et al. 2015;Xu et al. 2019), the newly identified BNRF1 variant (V1222I) yields the highest OR. Importantly, the combination of these four variants generated the highest risk estimation of 32.00 (95% CI: 9.18-111.49), indicating that EBV subtype containing these risk mutations may contribute to carcinogenic virulence in Chinese population. It is noteworthy that the genetic variation of BALF2 V317M and BNRF1 V1222I showed the strong evidence of positive selection, which possibly caused by the viral evolution to avoid host immune surveillance. In other words, the interaction between host immune system and virus shapes potential NPC risk subtypes and both contributed to the carcinogenesis. Indeed, the human leucocyte antigen genes have shown strong association with NPC (Bei et al. 2012). Thus, future combined research of host genetic susceptibility and viral subtype with careful consideration of possible founder effects may provide sufficient evidence of carcinogenicity as well as more precise prediction for NPC.       All the positively selected codon for the three proteins were present, apart from the four NPC-associated amino acid substitutions reported in Tables 2 and 3. The positively selected codon was identified by at least one method with P-value <0.1 or FUBAR posterior probability >0.9. SLAC, single-likelihood ancestor counting; FEL, fixed effects likelihood; MEME, mixed effects model of evolution; FUBAR, fast, unconstrained Bayesian approximation for inferring selection; NA, not applicable; NPC, nasopharyngeal carcinoma.
* P-values are estimated using chi-square tests. Bold text indicates statistically significant difference between NPC cases and controls with P-value less than 0.05.
# Bold values indicate positively selected codon inferred by different methods.
The protein BNRF1 exhibit largely coincidence between the NPC-associated AA substitutions (7/12) and its positively selected sites. There seems to be two tendencies, BNRF1 V1222I showing positive association with NPC only presented in China, while five others enrolling lower risks was relatively more frequent in non-endemic areas. BNRF1 is a major tegument protein which plays the roles in the inhibition of host intrinsic defenses against viral reactivation, virus entry and host chromosomal instability (Shumilov et al. 2017;Tsai et al. 2011). We speculate that some adaptive mutations on BNRF1 may have effect on the development of NPC in Southern China. Further experimental investigations are needed in order to validate the roles of these mutations in NPC.
There are some limitations in this study. First, the detailed ethnic information of the host was not available. Some sequences from immigrants and multiethnic areas such as Southeast Asians, may be misallocated, although the EBV sequences from these areas were relatively less (N ¼ 30). The results of PCA showed distinct cluster for these isolates from other Asians, indicating different origins. Second, in this study, no lymphomaassociated cluster in any populations was found, which may be due to the mixture of pathological subtypes and small sample size. Lastly, the phylogenic tree construction was limited by the available sequences, and the trees without detailed branch length cannot reflect evolutionary order and distance. Therefore, phylogenetic analyses were used only as a genetic classification method and the conclusions of genetic diversity between populations were drawn.
In summary, we have proposed three distinct patterns of population structure for different EBV proteins on the basis of the general Asia-Western/Africa segregation. These potential NPC risk EBV subtypes were identified from a Chinese-unique cluster, indicating geographic-and disease-specific feature concurrently. The evidence of positive selection suggests that viral adaptive mutation on some critical proteins like BNRF1 may involve in NPC development. Our findings provide a comprehensive overview of EBV population structure worldwide, implicating that viral disease-specific mutations should be investigated jointly with the population-specific variations. More importantly, our results provide the novel insights into the EBV carcinogenesis in NPC from the aspect of viral adaptive evolution and acknowledge the significance of combined research on host immunity-related genes and EBV.