Molecular evolutionary analysis of novel NSP4 mono-reassortant G1P[8]-E2 rotavirus strains that caused a discontinuous epidemic in Japan in 2015 and 2018

In the 2010s, several unusual rotavirus strains emerged, causing epidemics worldwide. This study reports a comprehensive molecular epidemiological study of rotaviruses in Japan based on full-genome analysis. From 2014 to 2019, a total of 489 rotavirus-positive stool specimens were identified, and the associated viral genomes were analyzed by next-generation sequencing. The genotype constellations of those strains were classified into nine patterns (G1P[8] (Wa), G1P[8]-E2, G1P[8] (DS-1), G2P[4] (DS-1), G3P[8] (Wa), G3P[8] (DS-1), G8P[8] (DS-1), G9P[8] (Wa), and G9P[8]-E2). The major prevalent genotype differed by year, comprising G8P[8] (DS-1) (37% of that year’s isolates) in 2014, G1P[8] (DS-1) (65%) in 2015, G9P[8] (Wa) (72%) in 2016, G3P[8] (DS-1) (66%) in 2017, G1P[8]-E2 (53%) in 2018, and G9P[8] (Wa) (26%) in 2019. The G1P[8]-E2 strains (G1-P[8]-I1-R1-C1-M1-A1-N1-T1-E2-H1) isolated from a total of 42 specimens in discontinuous years (2015 and 2018), which were the newly-emerged NSP4 mono-reassortant strains. Based on the results of the Bayesian evolutionary analyses, G1P[8]-E2 and G9P[8]-E2 were hypothesized to have been generated from distinct independent inter-genogroup reassortment events. The G1 strains detected in this study were classified into multiple clusters, depending on the year of detection. A comparison of the predicted amino acid sequences of the VP7 epitopes revealed that the G1 strains detected in different years encoded VP7 epitopes harboring distinct mutations. These mutations may be responsible for immune escape and annual changes in the prevalent strains.


Introduction
Rotavirus A (RVA) is the major cause of gastroenteritis in infants and young children worldwide.In 2017-2018, RVA caused an estimated approximately 200,000 deaths globally in children under five years of age (Cohen et al., 2022).Although most RVA-related deaths occur in developing countries, RVA also imposes a substantial burden in developed countries, including Japan (Tate et al., 2016;Burnett et al., 2020;Tsugawa et al., 2021;Hallowell et al., 2022).Two live attenuated RVA vaccines were introduced in Japan in 2011 (Rotarix ® , GlaxoSmithKline Biologicals, Rixensart, Belgium) and in 2012 (RotaTeq ® , Merck & Co., Inc., Kenilworth, NJ) and included in the routine vaccination program since October 2020.RVA vaccines are effective; however, selective pressure on a vaccine may induce an epidemic strain shift.To monitor Japanese RVA epidemic strains, since 2012, our laboratory has conducted RVA surveillance based on full-genome analysis (Fujii et al., 2014(Fujii et al., , 2019(Fujii et al., , 2020)).
In this study, we performed comprehensive molecular epidemiological research based on nearly complete genome analyses of RVAs isolated in Japan between 2014 and 2019.Notably, in 2015 and 2018, we detected 42 strains of G1P[8]-E2 (G1-P[8]-I1-R1-C1-M1-A1-N1-T1-E2-H1), a novel NSP4 mono-reassortant strain.Consequently, we focused on G1P[8]-E2 strains and conducted a detailed analysis of their genetic evolutionary trends.This is the last comprehensive surveillance data collected before the COVID-19 pandemic and before the initiation of routine RVA vaccination in Japan.

Sample collection
As shown in the map in Supplementary Figure S1, surveillance was conducted at 12 sentinel hospitals in Hokkaido Prefecture (Hakodate Municipal Hospital, Hokkaido Medical Center, Iwamizawa Municipal General Hospital, JCHO Sapporo Hokushin Hospital, Japan Red Cross Urakawa Hospital, NTT Medical Center Sapporo, Rumoi City Hospital, Sunagawa City Medical Center, Steel Memorial Muroran Hospital, Takikawa Municipal Hospital, Tomakomai City Hospital, and Yakumo General Hospital) between 2014-2019.Patients under 15 years of age who were admitted to these hospitals with acute gastroenteritis were screened for the presence of the RVA antigen using commercially available immunochromatographic tests.Analyses were performed on RVA-positive stool specimens collected from patients.The samples were suspended in phosphate-buffered saline (at approximately 10%) and stored at −20°C until required for further analyses.All patients or patient guardians from which samples were collected provided their informed consent.The study protocol was approved by the Ethics Committee of Biomedical Science in the National Institute of Infectious Diseases, Japan (approval number: 956).TRIzol ® LS Reagent was added to 80 μL of stool suspension and homogenized by vortexing.After a 5-min incubation at room temperature, 320 μL of ethanol was added, and the mixture was loaded directly onto the provided spin column.The column was centrifuged at 12,000 × g for 1 min and washed with prewash and wash buffer.The purified RNA was then eluted with 40 μL of DNase/RNase-free water.

cDNA library and nucleotide sequencing
Next,-generation sequencing was performed as previously described (Dennis et al., 2014;Doan et al., 2016;Fujii et al., 2019).Briefly, a 200-bp fragment library was constructed for each sample using the NEBNext Ultra RNA Library Prep Kit for Illumina v1.2 (New England Biolabs, Ipswich, MA, United States), according to the manufacturer's instructions.Library purification was performed using Agencourt AMPure XP magnetic beads (Beckman Coulter, Pasadena, CA, United States) according to the manufacturer's instructions.DNA concentrations were determined on a Qubit 2.0 fluorometer using the Qubit HS DNA Assay (Invitrogen, Carlsbad, CA, United States).A 151-cycle paired-end-read sequencing run was conducted on a MiSeq desktop sequencer (Illumina, San Diego, CA, United States) using the MiSeq Reagent Kit v2 (300 cycles).Sequence data were analyzed using CLC Genomics Workbench Software v7.0.3 (CLC Bio, Aarhus, Denmark), and sequences of representative strains were deposited in the GenBank/EMBL/DDBJ databases under Accession Nos.LC750855-LC752086, LC763123-LC763221.

Genotyping and phylogenetic analysis
Genotypes of the 11 genome segments were determined using the Rotavirus A Genotyping Tool v 0.1 on Rijksinstituut voor Volksgezondheid en Milieu (RIVM: https://www.rivm.nl/mpf/typingtool/rotavirusa/) and Nucleotide BLAST on the National Center for Biotechnology Information (NCBI: https://blast.ncbi.nlm.nih.gov/Blast.cgi).Near-full-length genome sequences were aligned with reference sequences using MEGA 7 software (Tamura et al., 2013) and MAFFT multiple-sequence alignment software, v 7.0 (Katoh et al., 2009).The best substitution models were selected based on the corrected Akaike information criterion value, as implemented in MEGA 7. Phylogenetic trees were constructed using the maximum likelihood estimation method with 1,000 bootstrap replicates.Representative reference sequences detected before 2019 were selected from GenBank for the construction of phylogenetic trees and subsequent evolutionary analysis.Only sequences with lengths >90% for each segment were used.Sequences of almost identical strains (such as strains detected in the same area and same year) were excluded.

Bayesian evolutionary analysis with BEAST
Evolutionary rates and the Times of Most-Recent Common Ancestor (TMRCAs) were calculated by the Bayesian Markov Chain Monte Carlo (MCMC) method implemented in BEAST v1.8.1 (Drummond et al., 2012).The best-substitution models used for BEAST analyses were calculated using the MEGA 7 software.A strict clock and coalescent exponential growth model (Drummond et al., 2002) was employed.MCMC runs were carried out for 200 million generations to achieve convergence, with sampling conducted every 1,000 steps.Convergence was assessed from the effective sample size after a 10% burn-in using the Tracer software v1.6. 1 Only parameters with an effective sample size of >200 were accepted.Maximum clade credibility (MCC) time-scaled trees were annotated with Treeannotator and viewed using FigTree v1.4.4. 2

Genotype distribution in Hokkaido, Japan
Between 2014 and 2019, a total of 489 RVA-positive stool specimens were identified and analyzed using next-generation sequencing.The genotype constellations of these strains were classified into nine patterns (G1P S1).The genotype distributions for each year are shown in Table 1.In Japan, RVA gastroenteritis cases occur mainly from January to June, with a peak in March and April.The same trend was observed in our present study.The major prevalent genotype differed by year as follows: G8P in 2018, and G9P[8] (Wa) (26%) in 2019.The G2P[4] (DS-1) and G9P[8] (Wa) genotypes were detected continuously in all six years, and the G8P[8] (DS-1) genotype was detected across five years except in 2018.The G1P[8] (DS-1) genotype was prevalent until 2015 but was subsequently (since 2016) replaced by G3P[8] (DS-1), which possesses a similar genetic background.The G1P[8]-E2 and G9P[8]-E2 genotypes are novel and unusual NSP4 mono-reassortant strains.The G1P[8]-E2 genotype initially spread (to a certain extent) in 2015 but was not detected in 2016 or 2017.This genotype subsequently re-emerged in 2018 but was not detected in 2019.In contrast, the G9P[8]-E2 genotype was continuously detected in 2018 and 2019.The genotype distribution for each hospital is shown in Supplementary Figure S2.

Phylogenetic tree of VP7 gene in G1 strains
As the first step in clarifying the relationship between the new NSP4 mono-reassortant G1P[8]-E2 strain and other major G1 strains, a phylogenetic tree of the G1-type VP7 gene segment was constructed (Figure 1).This tree incorporated sequences from 168 strains, including the 99 detected in this study, 25 other recent Japanese strains, and 44

Evolutionary rate of RVA genome
Next, to investigate the evolutionary history of the Japanese RVA strains, the evolutionary rate was calculated for each genotype of each gene (Supplementary Table S2).The number of sample sequences used for the calculation and best substitution models are shown in Supplementary Table S2.For this analysis, only human wild-type RVA strains were used to calculate the natural rate of evolution in the human population.Vaccine, tissue-cultured, and animal-derived strains were excluded.The estimated mean evolutionary rates ranged between 5.48 × 10 −4 (R1-type VP1) -1.18 × 10 −3 (E2-type NSP4) nucleotide substitutions per site per year.The rates in genotype-2 (i.e., DS-1-like genotypes) were higher than those of genotype-1 (i.e., Wa-like genotypes) in many genes (VP6, VP1, VP2, VP3, NSP2, NSP3, NSP4, and NSP5).

Time-scaled MCC trees
Focusing on the G1P[8]-E2 strains, the time-scaled MCC trees were constructed for each genotype for each gene.The numbers of sequences and models used for each tree are listed in Supplementary Table S2.In the VP7 G1 MCC tree (Figure 2), the Japanese G1 strains detected in this study were concentrated into four clusters, as shown in the phylogenetic tree (Figure 1).The G1P [8] strains detected in Japan between 2011-2019 were sorted into different clusters each year, indicating that the same strains were unlikely to be continuously prevalent in a given region.However, all the detected G1P[8]-E2 strains shared a recent common ancestor according to the TMRCAs of 2013.2 (95% highest posterior density (HPD) interval: 2012.1-2014.2) (Figure 2; Supplementary Table S3).

Amino acid sequences of VP7 epitope
To investigate the cause of the discontinuity in the G1 epidemic strains, the amino acid sequences of the VP7 epitopes of representative Japanese strains from each year were compared (Figure 7) (Aoki et al.,

Discussion
Our study analyzed the nearly complete genome sequences of RVA strains detected in Hokkaido, Japan, over the course of six years  1; Supplementary Table S1).Each of these novel strains achieved nationwide prevalence for 2-3 years, after which each strain was detected in several smaller outbreaks, as reported in the Infectious Agents Surveillance Report 3 and previous studies (Fujii et al., 2019(Fujii et al., , 2020;;Tsugawa et al., 2021).However, G1P[8] (DS-1) strains have become almost undetectable in recent years and have been replaced by equine-like G3P[8] (DS-1) strains with a similar genetic background.In contrast, G2P[4] and G9P[8] strains remained prevalent in Japan throughout the study 3 https://www.niid.go.jp/niid/en/iasr.htmlperiod.G8P[8] (DS-1) strains were continuously detected since 2014 (except for 2018).However, the recent emergence of two classes of NSP4 mono-reassortant strains, G1P[8]-E2 and G9P[8]-E2, is remarkable.The G9P[8]-E2 strains are widely prevalent in Japan (Fujii et al., 2020;Fukuda et al., 2022).Phan et al. reported the detection of a single G1P[8]-E2 strain in Hokkaido, Japan, in 2018 (Phan et al., 10.3389/fmicb.2024.1430557Frontiers in Microbiology 09 frontiersin.org2020), but only limited data were provided.This is the first study that provides evidence that the G1P[8]-E2 strain caused large-scale epidemics in Hokkaido, Japan in 2015 and 2018.The G1P[8]-E2 strains were detected in Tomakomai, Muroran, and Iwamizawa in 2015, as well as in Tomakomai, Muroran, Iwamizawa, Sapporo, Hakodate, Sunagawa, and Urakawa in 2018.However, no G1P[8]-E2 strains were detected in any part of Hokkaido in 2016 or 2017, and to our knowledge, none have been detected outside Hokkaido to date.In Japan, RVA strains similar to those found in China and Southeast Asia are often detected.On the other hand, the genotype distribution seems to vary from country to country, with G9P[8] continuing to be the major prevalent strain in China (Tian et al., 2021), and G4P[6] and  G8P [6] being prevalent in Korea (Lee et al., 2018).Further studied are required to understand the causes of these differences.
The observation of prevalent G1P[8]-E2 strains in discontinuous years is interesting.However, it remains unknown where this virus was present during the intervening two-year interval.This virus could have been maintained among older children, adults, and vaccinators with asymptomatic infections, as suggested by several studies (Gunawan et al., 2019;Bennett et al., 2021).The discontinuity in the prevalence of these epidemic strains may be attributed to developed herd immunity.Most patients with RVA gastroenteritis are infants; once a particular strain is prevalent in an area, many children acquire immunity to that strain.Consequently, other strains with different antigenicities gain prevalence in the following season(s).As generations shift, immunity to the previous strain decreases, and other strains become more prevalent.The same mechanism may be responsible for the annual changes in the major epidemic strains observed in this study (G8P[8] 1).In Japan, RVA vaccines were introduced in November 2011 and were incorporated into routine vaccination programs in October 2020.The vaccination rate during the study period was estimated to be approximately 50-70% (Tsugawa et al., 2021).Therefore, the prevalence of RVA strains in Japan reflects a complex situation involving two major factors: herd immunity acquired by vaccination and natural infection with wild strains.Further studies are needed to clarify the differences in the quality of cross-reactivity among immune responses acquired by vaccination and natural infection with wild strains.While previous studies have elucidated upon the development of herd immunity induced by RNA vaccines (Paulke-Korinek et al., 2011;Mast et al., 2015;Ito and Higashigawa, 2021), the G12 strain continued to be predominant in the United States even after vaccine introduction.It has been suggested that one reason may be due to a lower level of vaccine efficacy against G12 strains (i.e., different VP7 genotypes) (Esona et al., 2021).Although difficult to evaluate, Ianiro et al. ( 2016) also noted that the genetic shift of G1 strains in Italy could be attributed herd immunity.
In this study, G1 RVA was detected each year, except in 2016, and phylogenetic analyses showed that recent G1 strains were segregated into many clusters (Figures 1, 2).The prevalent strains were sorted into a single cluster in 2018 and 2019, however, in 2014 and 2015, they were sorted into multiple clusters.As the years progressed, epidemic strains were often replaced by strains from other clusters.This observation suggests that epidemic strains may shift in response to the effects of herd immunity, even among strains of a specific genotype.The data obtained by comparing VP7 epitope sequences (Figure 7) support this hypothesis.Previous studies have shown that amino acid substitutions at positions 94, 97, 147, and 291 of the VP7 gene significantly affect the antigenic recognition of human G1 strains (Coulson and Kirkwood, 1991).Compared to the vaccine strains, the Japanese G1 strains have several mutations.Mutation patterns differed depending on the year of detection.Mutations N94S, S123N, M217T, and K291R have been reported in strains isolated from other countries (Novikova et al., 2020;Zhou et al., 2020), suggesting that these substitutions have already been reported in recent G1 epidemic strains.Additionally, in this study, the VP7 gene harbored the L148F mutation in 2012 and 2014, the D100E mutation in 2015 and 2019, the T242M mutation in the G1P[8] (DS-1) strain in 2012-2015, and the N147D mutation in the G1P[8]-E2 strain in 2015 and 2018.These substitutions may be responsible for immune escape and annual changes in the prevalent strains.Mwangi et al. previously reported the presence of the N147D mutation in South African strains following vaccine introduction and showed that this substitution destabilized the structure of the VP7 protein (Mwangi et al., 2020).However, strains with the same mutations were observed in Thailand (LC086752) and Vietnam (KX362880) and were closely related to the G1P[8]-E2 strains detected in this study (Figure 1).Thus, it is likely that viruses encoding VP7 with the N147D substitution retain sufficient stability to spread.Therefore, these South African, Thai, and Vietnamese strains are thought to be derived from the same parental strain as the  S3).
The epitope sites shown here were identified using a limited number of strains.Epitopes that have not yet been identified may also be involved in antigenicity; however further studies are required to clarify this.In this study, the evolutionary rate of each gene was calculated for each genotype (Supplementary Table S2), and the resulting values ranged from 5.48 × 10 −4 (for R1-type VP1) to 1.18 × 10 −3 (for E2-type NSP4) /site/year.For eight genes (VP1, VP2, VP3, VP6, NSP2, NSP3, NSP4, and NSP5), the evolutionary rates tended to be higher in genotype-2 viruses than in genotype-1.Similar trends were reported in other studies (Liu et al., 2022).This difference may reflect high sequence diversity due to the emergence and prevalence of new DS-1like strains in recent years.Previous studies have suggested that reassortment frequently occurs in the NSP4 segment (Agbemabiese et al., 2016;Fujii et al., 2019).The emergence of G1P[8]-E2 and G9P[8]-E2 may be attributed to the characteristics of the NSP4 segment, which is easily substitutable, and the increasing presence of strains bearing DS-1-like genotypes.NSP4 is a known multifunctional protein associated with enterotoxin activity, RNA replication, interaction with viroplasm, intracellular calcium regulation, and pathogenicity (Estes and Greenberg, 2013).NSP4 is a potential target for protective immunity, but its proportion in the total immune response to RVA is unknown.Since the E2 type of the NSP4 gene is also carried by conventionally prevalent G2, G8, and equine-like G3 strains, it is unlikely that the emergence of G1P[8]-E2 alone will dramatically change the quality of herd immunity.
This study demonstrates the complexity of RVA epidemic strains in Japan.Because many new reassortant strains have emerged in recent years, it has become increasingly difficult to predict which strains will be prevalent in the future.Our findings will be useful in predicting future epidemic strains in areas where RVA remains prevalent.Since the SARS-CoV-2 pandemic in 2020, the outbreak of RVA has dramatically decreased in Japan, as reported in the Infectious Diseases Weekly Report (IDWR; https://www.niid.go.jp/niid/en) and previous studies (Fukuda et al., 2021;Hibiya et al., 2022).This shift is considered to have been due to travel, social restrictions, and non-pharmaceutical interventions (Dapper et al., 2022).After resuming social activities, RVA strains are speculated to once again readily spread domestically and internationally, as noted in other articles (Lappe et al., 2023;Wang et al., 2023).Continuous molecular epidemiological studies based on full genome analysis will continue to be critical for detecting the emergence of new epidemic RVA strains.Moreover, it is important to study the cross-reactivity of the protective immune responses against recent epidemics.

FIGURE 2
FIGURE 2 Maximum clade credibility (MCC) tree of the VP7 gene.Bayesian time-scaled tree for the G1-type VP7 gene was constructed.The representative Japanese strain groups are indicated by square brackets, and their type and year of detection are shown on the right.Estimated ages are indicated with arrows for representative nodes.The strains detected in this study are colored according to the group: blue, G1P[8] (Wa); purple, G1P[8] (DS-1); light blue, G1P[8]-E2.