Static and Evolving Norovirus Genotypes: Implications for Epidemiology and Immunity

Intra-genotypic diversity

To investigate the diversity and evolutionary differences of the distinct norovirus genotypes, more than 2000 sequences of the gene (ORF2) encoding the VP1 were retrieved from GenBank, with 101 and 1909 genes from GI strains and GII strains, respectively (Table 1). Individual sequences were vetted with an online norovirus typing tool that follows a widely-used universal classification and nomenclature system for norovirus genotypes [25], and an effort was made to verify the date of occurrence with the supporting documentation. Genotypes with 10 or more complete (or nearly complete) ORF2 sequences were selected for further phylogenetic analysis at both the nt and amino acid (aa) level, and included 16 out of the 31 current GI and GII genotypes (Table 1). We defined an intra-genotypic variant as a group of strains (≥ 2) that clustered together in the phylogenetic tree and that showed <5% difference in their nt or aa sequences, but ≥5% difference compared to other strains. Most genotypes segregated into 1 to 5 phylogenetic variants when nt sequences were analyzed (Table 1), with the exception of genotype GII.4 that displayed at least 10 different variants (Table 1 and Fig 1A). The number of variants was lower in seven genotypes (GI.1, GI.4, GI.6, GII.2, GII.12, GII.13 and GII.14) when aa sequences were used for tree reconstruction and distance analyses (Table 1). To link the different intra-genotypic variants to phenotypic characteristics in the VP1, we focused on the analysis of aa sequences. Interestingly, non-GII.4 genotypes presented variants with strains that have been detected many years apart (mean: 24.9 years, standard deviation: 12.97 years) while having only a few differences in their aa sequence (Table 1). An example is the GII.6 genotype, with each of the three variants (A-C) containing strains that differed by only a few aa residues (≤1.2%) but that were detected up to 41 years apart (Fig 2A). In contrast, the GII.4 genotype was comprised of variants that were present in the human population from 3 to 8 years (mean: 5.3 years, Fig 1A).

We developed an algorithm to illustrate visually the relationship between amino acid diversity and time among strains from a given genotype. The algorithm generated a heat map in which each square represents the number of strains with such aa difference in their VP1 and plotted against the timespan of detection. Analysis of 16 genotypes with sufficient data in GenBank revealed two distinct patterns of variant evolution; one in which the number of aa differences accumulated continually over time (Fig 1B), and one in which the number of aa differences remained relatively constant over time, regardless of the timespan between strains (Fig 2B). The first pattern (Fig 1B) was related to a constantly changing “evolving genotype” represented by the GII.4, while the second pattern (Fig 2B) was a highly conserved or “static genotype” represented by the 15 other genotypes with sufficient data for analysis (Table 1). The static genotypes resolved readily into distinct intra-genotypic variants, with one exception. The diversity plot of the GII.12 genotype strains displayed a subtle accumulation of differences over a time period of approximately 20 years; however, those differences were not associated with different intra-genotypic variants (S1 Fig). A larger number of sequences over a longer period will be helpful for defining variant diversity in GII.12 and other static genotypes.

During 2014–2015, a sharp increase in the number of gastroenteritis outbreaks was reported in Asia [22, 34, 43] that were associated with the emergence of a new variant of the genotype GII.17, i.e. variant C or Kawasaki_2014 [35, 43]. Our phylogenetic analysis of this genotype revealed four distinct variants, with one of the variants having strains that spanned over 37 years with a low level of sequence diversity consistent with a static pattern (S2 Fig). In contrast, the emerging strains showed multiple substitutions at the nt and aa sequence level [22, 35][44], which led to diversification into two separate phylogenetic variants (C and D). These differences appeared to accumulate over time, with predominant strains that circulated during the 2014–2015 season (variant D) differing by 5.4±1.1% from those in the 2013–2014 season (variant C) in their aa sequence.

To gain insight into the evolutionary differences noted at the aa level, the nt rate of evolution and nonsynonymous substitutions (dN)/synonymous substitutions (dS) ratio were calculated for each of the genotypes included in this analysis. The nt rates of evolution were similar among all the norovirus genotypes (range: 5.40x10^-3–2.23x10^-4 nt substitutions/site/year). However, differences were seen in the dN/dS ratios, with genotypes GII.4 and GII.17 presenting slightly higher values than any other norovirus genotype (Table 1). Note that in five different genotypes (GI.3, GI.4, GII.4, GII.13, and GII.17), the dN/dS ratio in the P2 encoding region was at least two times higher than the complete ORF2 (encoding VP1) dN/dS ratio. The dN/dS ratio has been used to inform the evolutionary pressures on a gene: a dN/dS > 1 (higher number of non-synonymous mutations) indicates positive evolution as the phenotype is changing due to pressures of the environment (e.g. immune responses), while a dN/dS < 1 indicates a purifying selection (also known as negative selection), where the new phenotype is mostly deleterious and eliminated from the population [45–47]. Despite the small differences in the dN/dS ratios, purifying selection (dN/dS < 1) is strongly acting at the VP1 protein level in all norovirus genotypes, including GII.4.

Intra- and inter-host evolution of norovirus genotypes

We next examined whether the different patterns of diversification observed in ORF2 extended to other regions of the genome. We analyzed by NGS full-length norovirus genomes from immunocompetent individuals infected in different settings. The first set of samples was from a child who was consecutively infected with three different genotypes (GII.4, GII.6 and GII.17) over a 3-year period [18, 35]. Although in each episode the child resolved the symptoms within ~ 72 hours, viral RNA was detected in stools for weeks after onset of symptoms. The full-length genome sequences were compared within the first 3 weeks for GII.4 and GII.6 viruses and the first 2 weeks for GII.17 viruses. A total of up to 78,000 reads/site (mean: 14688, standard deviation: 5675) were obtained for each sample by NGS (S1 Table). All consensus sequences were identical to the reference (day 1 [d1]); however, mutations ranging from 5 to 50% of the total reads (S1 Table and Fig 3) were found at later time points for each virus. In GII.4 viruses, 12 nt mutations arose in the subpopulations of the sample collected at d14, with nine of them being non-synonymous mutations. Three aa mutations were located in the P domain, the capsid surface-exposed region of the VP1 protein, with two present in antigenic sites A (E368A) and E (S412N) [48] (S3 Fig). By d21, the GII.4 virus had acquired 20 nt mutations in its subpopulations, with most of them (19/20) new mutations as compared to the d14 sequence. Nine out of the 20 nt mutations changed the aa sequence, with 4 of them mapping to the P domain, near or at antigenic sites (S3 Fig). In contrast, GII.6 and GII.17 viruses presented only four and two substitutions in their subpopulations at d21 and d14, respectively, with no evidence of the accumulation of mutations over time (S2 Table).

Although a large amount of norovirus is shed in stool, the infectious dose in natural transmission is likely low [41]. Thus, during inter-host transmission events noroviruses may undergo an initial reduction in the number of replicating viruses, creating a bottleneck effect. To compare inter-host evolution in individuals involved in outbreaks, we analyzed samples from outbreaks that occurred in the state of Maryland where the causative agents were identified as GII.6 (a hospital outbreak in 1971) [49] or GII.4 noroviruses (nursing home outbreaks in the 1987–1988 winter season) [32]. The samples and their dates of collection are indicated in Fig 4. Comparison of the NGS sequences with the outbreak consensus sequence revealed only a few substitutions, ≤ 5 nt and ≤ 2 aa, among samples from the same outbreak (Fig 4A and Fig 4B). However, when the consensus sequences from the GII.4 outbreaks were compared, a progressive accumulation of mutations (up to 86 nt and 16 aa) were detected in a period of three months, with no aa substitutions detected in the ORF2 (Fig 4B). To confirm these observations, we compared 151 genomes from GII.4 viruses (variant Den Haag) detected during three epidemic seasons in Japan [50, 51], and showed the accumulation of nt and aa substitutions over time (S4 Fig). In both sets of samples, the ORF2 (encoding VP1) acquired fewer amino acid substitutions as compared with ORF1 and ORF3, and thus maintained the VP1 phenotype for the GII.4 variant circulating in that given season. The data from full-length genome analyses were consistent with those from analyses of ORF2 in the GenBank database: different patterns of evolution exist among the norovirus genotypes in an acute outbreak setting.

Norovirus re-infection and genotype clustering

To reconcile our observations on the different mechanisms of diversification and data on re-infection and epidemiology of noroviruses, we investigated whether additional relationships might exist among the genotypes from the two major genogroups. Our phylogenetic tree constructed with representative strains from each genotype (strains from each lineage described here were included) showed clustering among certain genotypes (e.g. GI.3, GI.7, GI.8 and GI.9), while others appeared as single genotypes (e.g. GI.1, GII.3, GII.6; Fig 5A, S3 Table). The genotypic clustering was reproducible with a second phylogenetic methodology (S5 Fig). We designated each of the separate branches as groups A-L (Fig 5A), and the deduced aa sequences showed an approximate cut-off value of ≥20% aa differences between groups (Fig 5B). In a review of data from research groups that have documented norovirus re-infections and determined the genotype for each infection [17, 18, 38, 52–54], we observed that the pattern of re-infection might be consistent with the new grouping system as a predictor of antigenically-distinct strains. To test the hypothesis that these groups, provisionally designated here as “immunotypes,” might play a role in norovirus immunity, we developed a matrix that recorded the data from each of the consecutive re-infection cases with documented norovirus genotyping [17, 18, 22, 38, 52–54]. For example, a child consecutively infected with a GII.4 (“immunotype” G), GII.6 (“immunotype” H), and GII.17 (“immunotype” J) norovirus would count as one individual for the cell of the matrix that compares immunotype G and H, and as one individual for the cell that compares immunotype H and J. A matrix was constructed using re-infection data available from 116 children and 2 adults (Fig 5C, S6 Fig). Overall, the majority of re-infections occurred with viruses from different immunotypes, with re-infection rare from strains within an immunotype. A notable exception was immunotype G, which is comprised of genotypes GII.4 and GII.20. Re-infection of eight children (as shown by the 8 individuals in the black cell) was documented to have occurred with different variants of GII.4 viruses [17, 38].

Ethics statement

Stool specimens from the child were obtained with the written informed consent of the parent, and enrollment in National Institutes of Health (NIH) clinical study NCT01306084. Archival stool samples stored in the Laboratory of Infectious Diseases Calicivirus Repository were waived as exempt from IRB review by the NIH Office of Human Subjects Research and Protection (OHSRP 11833). Epidemiological information relating to the sample collection has been published elsewhere [18, 32, 35, 49].

Data mining and genomics

The full-length (or nearly full-length) ORF2 sequences (encoding for VP1) from each of the 31 genotypes described for GI and GII were retrieved from GenBank (accessed on March 2015) for analyses. Alignments were performed with Clustal W as implemented in MEGA v6 [75]. Sequences from each genotype were aligned separately to minimize the presence of insertions or deletions (indels), which can arise when different genotypes are compared. Phylogenetic trees were constructed using Kimura 2-parameter as method of nt substitution and Neighbor-Joining as algorithm of reconstruction as implemented in MEGA v6 with default settings. Phylogenetic trees that used aa sequences were reconstructed using a Poisson method of aa substitution. Bootstrap analyses were used to support the clustering of the variants. Information on the strains used for the phylogenetic analyses is provided in S4 Table. Evolutionary rates (nt substitutions/site/year) for each genotype were estimated using the ORF2 sequences and the Bayesian Markov Chain Monte Carlo (MCMC) approach as implemented in the BEAST package [76]. For each set of data the General Time Reversible (GTR) model with gamma rate distribution and invariable sites parameter was used and the MCMC was run for a sufficient number of generations to reach convergence of all parameters. All evolutionary rates were calculated using strict clock model and coalescent constant size tree prior, except for genotypes GI.4, GI.6, GII.14 and GII.16, which reached convergence using Bayesian Skyline and random local clocks. Selection pressures acting in the VP1 sequences were investigated by estimating the mean rate of nonsynonymous substitutions (dN) and synonymous substitutions (dN) and the dN/dS ratio as implemented in MEGA v6. The nearly full-length genome sequences from 151 GII.4 viruses detected in Japan during 2006–2009 [50, 51] were downloaded from GenBank and analyzed using MEGA v6 and Prism software (GraphPad Prism version 7).

Heat map plots of genotypic diversity

To visualize the aa substitutions within each genotype, a Python script (available upon request) was developed to calculate the number of aa differences and the isolation year differences between two individual strains. Isolation years were extracted from strain descriptions. The difference values were added into a matrix where the y-axis represents the isolation year differences and the x-axis the amino acid differences. Note that some cells will present more than one comparison, since strain pairs presenting the same number of aa differences and the same year difference, despite the years detected, will be included in the same cell. Heat map plots were calculated for each genotype using GraphPad Prism version 7 (GraphPad Software, La Jolla California USA), with the values representing the number of strains compared.

Full-length genome amplification

A platform was developed to analyze the plasticity of norovirus genotypes at the full genome level. Briefly, viral RNA was extracted from 10% (w/v) stool suspensions using the MagMax Viral RNA Isolation Kit (Ambion, California, USA) following manufacturer’s recommendations. Complementary DNA was synthesized from the viral RNA using the Tx30SXN primer (GACTAGTTCTAGATCGCGAGCGGCCGCCCTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT [77]) at 5μM final concentration, and the Maxima H Minus First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, California, USA) following manufacturer’s recommendations except that only 0.1 μL of Enzyme Mix was used per reaction. Amplification of the full-length genome was performed using 5 μl of the RT reaction, a set of primers that target the conserved regions of the 5’- and 3’-end of GII noroviruses (GII1-35: GTGAATGAAGATGGCGTCTAACGACGCTTCCGCTG, and Tx30SXN), and the SequalPrep Long PCR Kit (Invitrogen, California, USA) following manufacturer’s recommendations. Amplicons were excised from an agarose gel and purified with the QIAquick Gel Extraction Kit (Qiagen, California, USA).

Next-generation sequencing

Ion Torrent libraries were prepared by using 300–500 ng of full-length genome PCR amplicons following standard Ion Torrent library prep protocol. DNA was fragmented followed by the introduction of ligation barcode adapters. Adapted-ligated libraries were amplified using 13 PCR cycles, and size selected from agarose gels. Final libraries were quantified by Qubit (Invitrogen, California, USA), Bioanalyzer (Agilent), and qPCR. Libraries were normalized to 1nM, pooled at an equal molar ratio, and loaded onto a 318 v2 Chip in an Ion OneTouch2 machine. The sample from Ion OneTouch2 was transferred to an Ion OneTouch ES and then to an Ion PGM for sequencing with a 400bp kit (Life Technologies, California, USA). Ion Torrent sequence reads were de-multiplexed, and each individual set of reads was aligned to reference sequences using Bowtie2 and SAMtools [78, 79]. Aligned reads were visualized in the Integrative Genomics Viewer (IGV) [80] for single nt polymorphisms (SNPs) identification. Consensus sequence for each full-length genome was calculated using IGV. Read coverage (reads/nt position) was calculated using the genomecov command from BEDTools [81]. Sequence analyses were performed using MEGA v6 and Sequencher 5.4 (Gene Codes Corporation, Michigan, USA). The consensus sequence was calculated using default settings in Sequencher v5.4, and genomic sequences determined in this study were deposited into GenBank under Accession numbers KY424328 through KY424350. All other relevant data are within the paper and its Supporting Information files.