Human norovirus (HuNoV) belongs to the genus Norovirus in the family Caliciviridae and is a major causative agent of acute gastroenteritis in humans worldwide (Green, 2013; Robilotti et al., 2015). The large genetic divergence of HuNoV results in two genogroups and many genotypes (Green, 2013; Robilotti et al., 2015; Vinjé, 2015). Previous epidemiological studies revealed that GI.3, GII.2, GII.3, GII.4, GII.6, GII.12, and GII.17 were prevalent among HuNoV genotypes (Mathijs et al., 2011; Hoa-Tran et al., 2013; Vega et al., 2014a; Kumazaki and Usuku, 2015; Thongprachum et al., 2016). Importantly, GII.4 has caused pandemics of acute gastroenteritis in people of all ages in various countries during the past 10 years (Bull et al., 2006; Siebenga et al., 2009; Robilotti et al., 2015). Notably, multiple GII.4 variants strains having different antigenicity suddenly emerged and caused pandemics of gastroenteritis in several countries (Debbink et al., 2012a,b, 2013; Lindesmith et al., 2012a; Kumazaki and Usuku, 2015; Allen et al., 2016; Qiao et al., 2016; Parra et al., 2017). For example, in the 2006/07 season, the GII.4 variant emerged and spread rapidly in a large number of areas, including Japan, leading to pandemics in people of all ages (Tu et al., 2007; Motomura et al., 2008; Siebenga et al., 2009; Lam et al., 2012). The GII.4 variant 2006b, a representative strain (Den Haag 2006b strain), caused pandemics in five seasons (Hoa-Tran et al., 2013; Kumazaki and Usuku, 2015; Sato et al., 2017). In the 2006/07 season in Japan, it was estimated that over 3 million children and adults presented with acute gastroenteritis due to infections with the GII.4 2006b (Umeda, 2007). Furthermore, in the 2012/13 season, the GII.4 variant 2012, a representative strain (Sydney 2012 strain), also suddenly emerged and caused pandemics of gastroenteritis similar to that caused by the GII.4 2006b strain in the USA, Europe, and Japan (Siebenga et al., 2009; Allen et al., 2014; Vega et al., 2014a; Kumazaki and Usuku, 2015; Vinjé, 2015; Parra et al., 2017).

The HuNoV major capsid (VP1) protein is considered to be closely associated with the infectivity and antigenicity of these strains (Debbink et al., 2012a). Many previous reports suggested that the HuNoV VP1 gene rapidly evolved, resulting in a large divergence of antigenicity (Allen et al., 2008; Bok et al., 2009; Yang et al., 2010; Debbink et al., 2012a; Zakikhany et al., 2012). Earlier findings also indicated that the rapid evolution of the VP1 gene in HuNoV is strongly associated with gastroenteritis pandemics caused by HuNoV (Robilotti et al., 2015).

Various bioinformatic techniques for investigating viral genes may contribute to a better understanding of the evolution of the viruses (Andernach et al., 2013; Añez et al., 2013; Sun et al., 2017; Wei and Li, 2017). Of these, phylogenetic analyses using the Bayesian Markov chain Monte Carlo (MCMC) method can characterize evolutionary processes by estimating the date of divergence and evolutionary rates using sequences with a known sample-collection date for various organisms (Drummond and Bouckaert, 2015). Additionally, the Bayesian skyline plot (BSP) allows prediction of the transition of the infectious population size for pathogens (Hall et al., 2016). Conformational epitope analyses using in silico approaches also enable us to identify and map binding sites of host antibodies on three-dimensional structures of proteins (Potocnakova et al., 2016). For HuNoV, a large number of studies on viral evolution have been performed using bioinformatics (Giammanco et al., 2014; Kimura et al., 2015; Kobayashi et al., 2015, 2016; Lu et al., 2016; Parra et al., 2017; Siqueira et al., 2017; Tohma et al., 2017). Previous reports described the evolutionary rates of some HuNoV genotypes using the Bayesian MCMC method and linear regression analyses, which suggested differences in such rates among the genotypes (Parra et al., 2017). In addition, the blockade epitopes on the VP1 protein in some HuNoV genotypes were identified by amino acid variation analyses and molecular biological techniques using the virus-like particle (VLP) (Lindesmith et al., 2014). Previous phylogenetic approaches also showed that GII.4 continuously created genetically different clusters, which resulted in the emergence of many variants (Fioretti et al., 2014; Qiao et al., 2016). However, the differences of evolutionary features among all GII.4 variants, including recently discovered strains, have not been comprehensively elucidated (Bok et al., 2009; Siebenga et al., 2010; Qiao et al., 2016). It may be important to understand why GII.4 has frequently produced escape from previously acquired immunity and repeatedly caused large outbreaks. Moreover, a recent study constructed a system for predicting the HuNoV genotype composition in the peak season using a unique algorithm, namely, NOROCAST (Suzuki et al., 2016). Thus, these bioinformatics technologies may contribute to our understanding of the evolution of the HuNoV GII.4 capsid gene, including the changes of antigenicity, for advancing the prevalent prediction systems for HuNoV. In the present study, using various bioinformatic techniques, we compared the molecular evolution of all GII.4 variants strains using globally collected full-length HuNoV GII.4 VP1 gene sequences.

In the present study, we show the molecular evolution of the capsid VP1 gene in globally collected HuNoV GII.4, including recent strains, during a period of approximately 40 years. Our new findings can be summarized as follows: (1) The evolutionary pattern, including evolutionary rates, phylogenetic distance, and infectious population size, may differ between GII.4 variants. (2) The evolutionary rates and constraints on amino acid mutations are likely distinct between domains within GII.4 VP1. (3) There may be few epitopes on the exterior surface of the VP1 structure other than those positions on the P2 domain, as previously identified. (4) Prevalent GII.4 variants specifically undergo selective pressure that leads to new variant that can escape the immune system of the host.

The Bayesian MCMC method is used to characterize evolutionary processes in terms of date of divergence and the rates of evolution for organism genes with high evolutionary rates, such as RNA viruses (Drummond and Bouckaert, 2015). This method has already been used to study molecular evolution not only for noroviruses but also for various other RNA viruses (Añez et al., 2013; Vijaykrishna et al., 2015). In particular, it is important to estimate nucleotide evolutionary rates for HuNoV since previous reports showed that the viruses could evolve with the distinct rates among the genotypes depending on difference of polymerase fidelity and immune pressure for each genotype (Bull et al., 2010). Estimation of evolutionary rates using the method may contribute to predict the frequency of emergence of variants that can escape from previously acquired host immunity in the HuNoV genotypes. However, the results of the method are affected by the sequences included in the dataset. In this study, the time-scaled evolutionary phylogenetic tree showed that the common ancestor of GII.4 and GII.20 dated back to about 175 years ago. Furthermore, during the past 40 years, the GII.4 strains evolved and formed seven clusters including 14 known variants (Figure 1). Next, we estimated the HuNoV GII.4 VP1 gene evolutionary rate as 7.68 × 10⁻³ substitutions/site/year. Our observations are generally consistent with previous observations overall, with GII.4 having a higher evolutionary rate than other genotypes such as GII.7 (2.3 × 10⁻³ substitutions/site/year), GII.3 (4.16 × 10⁻³ substitutions/site/year), and GII.2 (1.31 × 10⁻³ substitutions/site/year) (Bok et al., 2009; Bull et al., 2010; Siebenga et al., 2010; Boon et al., 2011; Fioretti et al., 2014; Qiao et al., 2016; Mizukoshi et al., 2017; Parra et al., 2017). In addition, we estimated the different evolutionary rates among GII.4 variants. The rates among domains (shell, P1, and P2 domains) in the VP1 gene also tended to be distinct, with that of the P2 domain being the highest (Figure 2). A previous report showed that GII.4 had replication rates of polymerase that differed among the variants (Bull et al., 2010). In addition, the prevalence of GII.4 differed among the variants (Vinjé, 2015), which may produce differences of herd immunity pressure for each GII.4 variant. We also performed similarity plot analyses in the GII.4 variant strains. Our data showed that the similarities of the shell domains were relatively high (85–100%), followed by a lower level in the P1 domains, and the lowest level in the P2 domains (75–95%), which is consistent with previous reports that showed similarities in the shell, P1, and P2 domains of around >85, >80, and >75%, respectively (Eden et al., 2013; Doerflinger et al., 2017). Next, we calculated the phylogenetic distance among the GII.4 strains. Notably, the values of the intra-variant distance were different. Of them, Osaka 2007 variant strains tended to have high nucleotide variability, while Asia 2003 probably contained fewer nucleotide substitutions than the other variants (Figure 4). These results suggest that GII.4 has uniquely evolved among these variants having different evolutionary rates.

Earlier molecular biological and bioinformatic studies identified or estimated the blocking epitopes in the P2 domain within the GII.4 VP1 structure (Debbink et al., 2012b, 2013; Lindesmith et al., 2012a,b; Parra et al., 2012; Giammanco et al., 2014). These studies were mainly focused on the P2 domain and based on the amino acid variations on the exterior surface in the domain that arose over time. However, other conformational epitopes may also be contained in the shell and P1 domains. Since amino acid sequences of the shell and P1 domains are conserved than those of the P2 domain, it may be possible to use the epitopes of the shell and P1 domains for the development of diagnostic tools for HuNoV, such as immunochromatography kits. Additionally, antibodies associated with antigenicity in positions different from previously identified epitopes may be produced by the B-cell response in the host. To the best of our knowledge, conformational epitope analyses of the whole VP1 structure, including the shell and P1 domains, of various GII.4 variants using mathematical methods have not yet been conducted. The epitopes predicted in our data were close to or overlapped with those positions identified by previous in vitro experiments using VLP and monoclonal antibodies (Lindesmith et al., 2012a). Furthermore, our prediction identified few epitopes, including in shell and P1 domains, other than the previously reported positions. This result may suggest difficulty in acquiring antibodies broadly responding to various GII.4 variants that can escape from previously acquired herd immunity. To the best of our knowledge, this is the first observation of this kind in this field. This study suggests the need for additional clarification of the immunology of the GII.4 VP1 protein. However, the present study was performed only in silico. To obtain a better understanding of norovirus biology, further studies including in vitro epitope-mapping experiments for all antibodies produced by the host should be performed. Next, our analyses revealed that many common epitopes among the GII.4 variants strains were present in the P2 domain. Moreover, a number of substitutions of amino acids were detected in the conformational epitopes (Table 2). Other in vitro studies using the HuNoV GII.4 VLP established that the antigenicity changes of the GII.4 were associated with the amino acid substitutions in the P2 domains of GII.4 VP1 proteins (White, 2014). Additionally, previous in vitro and bioinformatic studies predicted contact residues specific to antibodies in epitopes A (aa294-298, aa362, 368), B (aa333, aa382), C (aa340, aa376), D (aa393–395), and E (aa407, aa412–413). It has also been confirmed that epitope A is an antigenic region, whereas epitope D is associated with HBGA binding sites. Moreover, epitope C has been predicted to be an evolving epitope and was confirmed using a monoclonal antibody. Epitope E is associated with the potential blockade of Farmington Hills 2002 variant. However, there are no reports on epitope B based on in vitro experiments. Previous bioinformatic studies suggested that the five epitopes are associated with the production of antibodies against various GII.4 strains on the exterior surface in the P2 domain (Debbink et al., 2012b, 2013; Lindesmith et al., 2012a,b; Parra et al., 2012; Giammanco et al., 2014; Garaicoechea et al., 2015). Thus, based on our results and those of other researchers, we propose that previously identified epitopes are robustly associated with the antigenicity of HuNoV GII.4.

We also estimated the selective pressures in the present GII.4 VP1 gene. The results obtained by the SLAC, FEL, and IFEL methods revealed positive selection for the whole GII.4 genotype at aa6, aa9, and aa534, which are located in the shell and P1 domains, whereas no positive selection sites were identified in the P2 domain (Table 3). However, separating for each GII.4 variant, positive selection sites in three prevalent GII.4 variants strains were predicted in the P2 domain (aa393 and aa412 in Den Haag 2006b, aa294, and aa376 in New Orleans 2009, and aa393 in Sydney 2012) (Table 4). These sites mostly overlapped with the epitopes that were identified by our prediction and previous in vitro studies (Table 2). Positive selection implies that non-synonymous substitutions facilitate the survival of organisms, such as escape from the immune response of the host (Kryazhimskiy and Plotkin, 2008). Moreover, a previous report predicted an association with positive selection of efficient RNA translation via secondary RNA structures (Siebenga et al., 2010). These results indicate that the prevalent GII.4 variants specifically undergoes strong selective pressure that may produce immune escape variants due to the numerous infectious populations with herd immunity. However, the importance of the positive selection in the shell and P1 domains requires further investigation by in vitro studies. Previous studies also estimated the positive selection sites in the GII.4 VP1 gene (Lindesmith et al., 2008; Bok et al., 2009; Siebenga et al., 2010; Vega et al., 2014b). For example, using the SLAC method, Bok et al. (2009) predicted six sites of positive selection (aa6, aa9, aa15, aa47, aa395, and aa534). However, the threshold for statistical significance in that study was higher than that in our study (p < 0.25). Siebenga et al. (2010) identified eight positive selection sites using FEL and IFEL (aa6, aa9, aa47, aa352, aa372, aa395, aa407, and aa534), which is partially consistent with our results. The partial discrepancy between our data and those in previous reports might have resulted from the difference in the datasets between studies. Moreover, our results revealed many sites under negative selection. In each domain, negative selection sites accounted for approximately 60–80%. The proportion of negative selection sites tended to be distinct between the domains in the VP1 gene (Table 5). However, low numbers of strains may result in no or few negative selection sites in some GII.4 variants. Negative selection exhibits non-synonymous substitutions reduce the survival of organisms. Namely, the sites should preserve synonymous substitutions to maintain the survival (Suzuki, 2004). These results suggest that a number of mutations in the VP1 gene may be under negative and neutral selection.

Next, we mapped the conformational epitopes, positive selection, and amino acid substitutions on the VP1 structures in GII.4 variants. Most putative epitopes and substitutions were located on the exterior surface of the P2 domain in all variants. The epitopes overlapped with or were proximal to the substitution sites. Moreover, positive selection sites for shell and P1 domains existed at the exterior surface of the structure that was distal to the predicted epitopes. However, the positive selection sites for Den Haag 2006b, New Orleans 2009, and Sydney 2012 variants were located on the exterior surface of the epitopes (Figure 5). These results indicate that HuNoV GII.4 has evolved with frequent alteration of antigenicity. In addition, the circulation of GII.4 may also be associated with mechanisms other than immune escape.

In addition, we estimated the fluctuation in the size of infectious populations in HuNoV GII.4. The population size increased around 1994–1996 and 2005–2006. Moreover, Den Haag 2006b, New Orleans 2009, and Sydney 2012 variants accumulated infectious populations in 2005–2007, 2004–2010, and 2006–2011, respectively (Figure 6). These results suggest that GII.4 has a capacity to repeatedly produce highly epidemic variants. A previous report showed evidence of higher affinity to a broad range of HBGA types in the Den Haag 2006b strains compared with that in other variants (de Rougemont et al., 2011). Thus, the prevalence of HuNoV GII.4 may be associated with multiple factors, such as the capacity to bind to host receptors and immune evasion.

Overall, this study utilized limited numbers of GII.4 VP1 sequences. For example, there were some biases among GII.4 variants in the amount of original data deposited in GenBank. Additionally, sequences with >98.9% identity were omitted. These biases may have affected the results of the bioinformatic analyses.

In conclusion, our findings show that the evolutionary pattern of GII.4 may be distinct among various variants. Moreover, our data may contribute to developing candidates of an efficient vaccine and prevalent prediction systems for HuNoV. In the near future, the accumulation of data from molecular epidemiological studies with continuous surveillance will be needed.

TM, KK, and HK designed this study. TM, KN, YMa, NN, AR, TS, AY, MK, YMo, YS, and NS analyzed the data. TM, KN, YMa, and HK wrote the manuscript. KK and HK supervised this study. All authors read and approved the manuscript.