Emergence of human G2P[4] rotaviruses containing animal derived gene segments in the post-vaccine era

The introduction of Rotarix into the Belgian immunization program in June 2006 coincided with an increase of the relative prevalence of G2P[4] strains. However, the genetic composition of these persistent G2P[4] strains has not been investigated. Therefore, we have investigated the NSP4 gene of 89 Belgian G2P[4] strains detected between 1999 and 2013, covering both pre- and post-vaccination periods. The NSP4 genes were divided over seven separate clusters of which six were more closely related to animal than to human strains. The NSP4 genes that clustered more closely to animal DS-1-like strains were isolated after 2004–2005 and were found throughout multiple seasons. Complete genome sequencing of 28 strains identified several other gene segments that clustered more closely to animal than to human DS-1-like strains. These findings suggest that frequent interspecies reassortments may have played a role in the spread of G2P[4] rotaviruses in the post-vaccination period in Belgium.

Italy and the USA revealed the co-circulation of multiple lineages for several gene segments, some of which are likely of animal origin [10][11][12] . In addition, the circulation of animal-human reassortant viruses possessing an animal NSP4 gene segment has been increasingly reported [13][14][15][16] . Animal-human reassortant viruses have the potential to change viral fitness and promote viral transmission, especially in an immune naïve population. Previously, this was observed for G9 and G12 Wa-like rotaviruses, possessing VP7 genes that are thought to have originated from porcine rotaviruses and have efficiently spread globally in the human population in less than two decades 17 .
In Belgium, Rotarix was introduced in 2006 and RotaTeq in 2007. Approximately 85% of all administered rotavirus vaccines in Belgium is Rotarix. The introduction of Rotarix coincided with an increased relative prevalence of G2P [4] rotaviruses, which was sustained throughout multiple seasons 18 . Although it was previously shown that the vaccine effectiveness of Rotarix against DS-1-like G2P [4] strains was slightly lower than against Wa-like strains [19][20][21] , we hypothesize that also changes in the DS-1-like rotaviruses themselves, may have influenced the sudden increase in the relative prevalence of DS-1-like rotaviruses in Belgium. To determine to what extent reassortment had occurred in Belgian DS-1-like rotaviruses we analyzed 89 NSP4 genes and 28 complete genomes of Belgian DS-1-like rotaviruses that were isolated before and after vaccine introduction.

Results
Viral dataset. Based on a preliminary phylogenetic analysis of VP7, a selection was made of Belgian G2P [4] rotaviruses collected between 1999 and 2013. In total, 89 strains were selected and the complete NSP4 gene was sequenced. NSP4 sequences were obtained from every rotavirus season between 1999 and 2013 except for the 2002-2003 rotavirus season. Despite repeated attempts no reliable sequences could be obtained for the few available strains from 2002-2003, possibly due to low viral load.
Interspecies reassortment of NSP4 is common in Belgian G2P [4] rotaviruses. The 89 Belgian NSP4 sequences were compared with DS-1-like strains deposited in GenBank ( Fig. 1) and were found in seven distinct clusters, which were between 10.1% and 15.1% different on the nucleotide level. A little more than half of the G2P [4] strains (n = 47; indicated by red circles in Fig. 1) were found in a cluster together with other typical human DS-1-like strains isolated from various places around the world. The remaining 42 Belgian G2P [4] strains were either clustering more closely to bovine or bovine-like human strains or were clustering relatively distantly from any other known strain (indicated by circles in different shades of blue). The largest of these clusters (animal cluster 3) comprised 31 Belgian G2P [4] strains, which were all detected between 2008 and 2013. They showed on average only 1.9% nucleotide diversity, indicating sustained human-to-human transmission instead of independent interspecies transmission events. These strains were most closely related to environmental and human G2P [4] strains from Brazil, Germany, Australia and the USA. Clustering slightly more distantly to these, there were two animal derived human strains from Hungary, RVA/Human-wt/HUN/BP1879/2003/G6P [14] and RVA/ Human-wt/HUN/BP1362/2005/G3P [9], which are thought to be of bovine(-like) and feline origin, respectively. Finally, clustering even more distantly, at approximately 6% nucleotide difference, there were bovine and canine rotaviruses isolated in Turkey and Italy, suggesting that this large cluster of human E2 strains may have originated in animals.
The third largest cluster (animal cluster 2) contained five closely related Belgian G2P [4] strains detected between 2006 and 2010 that were probably the result of human-to-human transmission events and were intermingling with other human G2 and G12 strains isolated from the USA, Thailand, India, Pakistan and Australia. Close to this cluster at 1.0-1.6% nucleotide distance a group of bovine rotaviruses isolated in India was found as well as two bovine-like human P [14] strains, RVA/Human-wt/TWN/04-97s379/2008/G8P [14] and RVA/ Human-wt/IND/N-1/2009/G6P [14] (Fig. 1).
Based on the phylogenetic analysis, Belgian NSP4 genes were classified according to their most likely origin: human (n = 47; 52.8%) or animal (n = 42; 47.2%). The temporal distribution of the Belgian bovine and human NSP4 genes showed that before the introduction of Rotarix in Belgium bovine NSP4 genes were only sporadically detected (only 2 out of 17 strains; 11.7%) (Fig. 2). The two bovine-like NSP4 genes detected before vaccine introduction (animal cluster 6) formed a single cluster and probably represented a single interspecies reassortment event. After vaccine introduction bovine-like NSP4 genes constituted 56% of all NSP4 genes (40 out of 72 strains) in our sample set, including both clusters were human-to-human transmission had occurred (animal clusters 3 and 4).
Complete genome analysis reveals zoonotic reassortments in gene segments other than NSP4. Twenty-eight strains were selected for complete genome sequencing based upon the genetic diversity The phylogenetic tree of VP7 showed that Belgian strains were found in three clusters together with contemporary human G2P [4] and G2P [6] strains isolated worldwide (Fig. 4). One Belgian strain, RVA/Human-wt/BEL/ BE67/2009/G2P [4], was more distantly related to the other Belgian G2P [4] rotaviruses (95.2-96.5% similarity on the nucleotide level). None of the Belgian strains clustered closely to the G2 component of the RotaTeq vaccine (6.2-7.8% difference). The VP4, NSP1 and NSP3 genes of the 28 Belgian strains clustered closely together with other contemporary human DS-1-like rotaviruses although for every gene segment one or two exceptions could be identified (Figs 4 and 5). In particular in the NSP1 tree, strain RVA/Human-wt/BEL/BE45/2007/G2P [4] was found to be 7.1-8.5% divergent from other Belgian DS-1-like NSP1 genes. The Belgian DS-1-like strains  belonging to the H2 genotype showed only little diversity (<2.4%) (Fig. 6). The two strains belonging to the H3 genotype, BE38 and BE41, were identical and closely related with RVA/Rhesus-tc/USA/PTRV/1990/G8P [1], which has previously been characterized as an interspecies transmission from a ruminant to a monkey 4,23 , and the bovine-like human strain RVA/Human-wt/DNK/11H64158/2008/G8P [14]. Even though NSP5 is not a genetically diverse gene, the RotaTeq NSP5 was 2.2% different from BE38 and BE41.
The phylogenetic trees of VP2 and NSP2 displayed a similar clustering pattern (Figs 7 and 8) with 26 out of 28 strains clustering with typical human DS-1-like strains and BE38 and BE41 clustering with bovine, bovine-like human or feline-like human strains. The RotaTeq vaccine strain was 90.3-90.0% and 93.2-93.0% similar to BE38 and BE41 for VP2 and NSP2, respectively.
Taken together, the phylogenetic analysis for each gene segment revealed a wide variation in the genetic composition of Belgian DS-1-like rotaviruses. Twelve Belgian G2P [4] strains possessed a completely human DS-1-like genotype constellation, whereas 10 Belgian rotaviruses were most closely related to human DS-1-like strains for ten out of eleven gene segments. Four Belgian strains contained two animal derived gene segments. For BE38 and BE41 the majority of gene segments were more closely related to bovine than to human strains (Fig. 3).

Discussion
A relatively high prevalence of G2P [4] strains has been observed in several countries where Rotarix is the mainly used rotavirus vaccine 18,24,25 . A slightly lower vaccine effectiveness has been shown against G2P [4] rotaviruses when compared to Wa-like rotaviruses (G1, G3, G4 and G9) in Belgium, which may explain the relative increase in the prevalence of G2P [4] rotaviruses in Belgium 19,26 . We show here that the vaccine introduction in Belgium is also temporarily associated with substantial changes in the genetic makeup of G2P [4] rotaviruses in terms of the large scale transmission of bovine-like NSP4 gene segments in the human rotavirus population (Fig. 1).
Currently, it is unclear whether the emergence of rotaviruses possessing a bovine-like NSP4 gene is a cause, a consequence or completely unrelated to the increased prevalence of G2P [4] rotaviruses after vaccine introduction. Rotaviruses containing a bovine-like NSP4 gene have been detected in other countries with rotavirus vaccination programs, such as Brazil, Australia, the USA and several countries in Sub-Saharan Africa 10,12,15,27 , but also in countries were no rotavirus vaccines are universally used like the Democratic republic of the Congo, Italy, Russia and Thailand 11,14,16 . It is yet unknown for how long these bovine NSP4 gene segments have been circulating in the human rotavirus population and to elucidate this, detailed large scale analyses with time-stamped animal and human DS-1-like NSP4 genes are needed. However, human rotaviruses carrying a bovine-like NSP4 gene have been detected in humans only relatively recently, suggesting that at least the two largest bovine-like NSP4 clusters observed in this study (Fig. 1) are the result of worldwide viral migration rather than of a putative selection pressure caused by the usage of the Rotarix vaccine within Belgium. Therefore, we hypothesize that it is most likely that the increased relative prevalence of G2P [4] rotaviruses as a result of universal Rotarix vaccination in Belgium also promoted the spread of G2P [4] rotaviruses with a bovine-like NSP4. Such a selective sweep has previously been reported for the VP7 gene segment of human G2P [4] strains and suggests that this is probably an important mechanism for G2P [4] rotaviruses to antigenically adapt 10 . The results of this study suggest that for NSP4 these selective sweeps may also occur in a similar way as for VP7. The fact that animal NSP4 genes are globally reported suggests that animal NSP4 genes may not pose a selective disadvantage over typically human NSP4 genes. However, adaptive immunity is less well understood for NSP4 than for VP7 and VP4 and clearly defined NSP4 antigenic epitopes are lacking. It seems likely that selective sweeps are not only restricted to NSP4 and VP7, but could potentially involve other gene segments as well.
Interspecies reassortment between bovine-like and human rotaviruses was frequently detected in Belgian G2P [4] rotaviruses, as we have found evidence for at least 8 independent interspecies reassortments in 7 different gene segments (VP6, VP1-3, NSP2, NSP4 and NSP5). For at least three clusters (two clusters in NSP4 and one in VP3) evidence was found for sustained human-to-human transmission as similar bovine-like genes segments were detected in multiple rotavirus seasons (up to five consecutive rotavirus seasons). However, our data also indicated that human-to-human transmission is rather rare for most of these human-animal reassortants, suggesting that most interspecies reassortments are probably dead-end infections. Although the three human-to-human transmission clusters are unlikely to have originated in Belgium, the interspecies reassortant dead-end infections may have originated in Belgium as they often do not cluster closely to any other known strains from the global collection (Figs 1 and 6-8). Alternatively, this could be due to undersampling or the lack of sufficient complete genome data of DS-1-like and animal rotaviruses.   Twenty-eight Belgian rotaviruses were selected for complete genome sequencing and this revealed that interspecies reassortment was not only restricted to the NSP4 gene segment, although reassortment was less often observed in gene segments other than NSP4. As the NSP4 gene encodes an enterotoxin, it be interesting to see if there were differences in clinical severity in patients infected with RVA strains with a bovine derived NSP4 gene product. Unfortunately, this information was not present. It has been shown that bovine rotaviruses are infecting humans on a regular basis. Compared to other animal rotaviruses, particularly bovine rotaviruses with a P [14] VP4 have been detected relatively frequently in humans 9,28,29 . As thus far no G2P [4] rotaviruses have been found in ruminants, humans are then most likely a mixing vessel in which interspecies reassortment occurs. It is unknown whether reassortment between human and bovine rotaviruses is restricted to certain gene segments, but it seems likely that gene segments sharing the same genotype in human and ruminant rotaviruses (VP1-VP3, VP6, NSP2 and NSP4) may reassort without great loss in fitness. For the other gene segments this is less clear and additional experiments using in vitro generated reassortants are needed to answer this question. However, we found two strains, BE38 and BE41, which contained six bovine gene segments, including an H3 NSP5. These two strains were genetically almost identical, isolated from the same rotavirus season and most likely epidemiologically linked to each other, indicating that a rotavirus with an H3 NSP3 is able to successfully spread from one human to another.
One of the main limitation of this study is the way samples were collected. As we did not anticipate to find such a great diversity in the NSP4 gene and to maximize genetic diversity in our selection, samples were collected based on the genetic diversity of VP7. This makes it difficult to extrapolate our findings to the Belgian G2P [4] rotavirus population as a whole. However, we do believe that our study provides a good framework for future studies to investigate the prevalence of animal gene segments in DS-1-like rotaviruses in different geographical regions and in countries with and without universal rotavirus vaccination program. Our G2P [4] samples originated from multiple hospitals across Belgium, but the impact of this on genetic diversity is probably low as Belgium is a small country. However, little is known about the impact of geographical variation on rotavirus genetic diversity at such a small scale.
Previously it has been reported that the RotaTeq vaccine strain reassorted with wild-type Wa-like strains, creating wild-type Wa-like strains possessing RotaTeq NSP2 gene segments 30 . In our dataset we did not find any evidence of reassortment with RotaTeq, despite the fact that RotaTeq is used in Belgium (approximately 15% of all administered rotavirus vaccines is RotaTeq). However, the high number of reassortants in our dataset suggest that the chance of reassortment between RotaTeq and human DS-1-like strains is not negligible, especially since the process of vaccination itself provides for ample reassortment opportunities. Therefore rotavirus surveillance involving the sequencing of complete DS-1-like genomes is especially needed in countries with a universal rotavirus vaccination program.

Methods
Stool samples were collected at multiple hospitals in Belgium in the framework of the Belgian National Reference Center activities for rotavirus at the University Hospital Leuven, organized and approved by the WIV (Wetenschappelijk Instituut voor Volksgezondheid; Scientific Institute for Public Health) in Belgium (https:// nrchm.wiv-isp.be/nl/ref_centra_labo/default.aspx). All experiments were performed according to the guidelines and regulations of the WIV. Stool samples were collected from patients hospitalized with acute rotavirus gastroenteritis, and the G-genotype and P-genotype of these samples was determined. To minimize redundant sequencing of identical NSP4 genes a preliminary phylogenetic analysis of G2 genotypes was conducted, which was used to select 89 samples that reflected the genetic diversity of the G2-genotype in every rotavirus season. Based on the phylogenetic diversity of the obtained 89 NSP4 sequences, 28 strains were selected for complete genome sequencing. In this way, our sample set is likely to reflect the overall genetic diversity of Belgian G2P [4] samples between 1999-2013, but caution should be taken with quantitative extrapolation as the samples were not randomly selected.
Stool samples were diluted once in PBS and RNA was extracted using the Qiagen Viral RNA mini kit (Qiagen) according to manufacturer's recommendations. The NSP4 gene segment was amplified with a one-step RT-PCR kit (Qiagen) using forward primer: 5′ -GGC TTT TAA AAG TTC TGT TCC-3′ and reverse primer: 5′ -GGW YAC RYT AAG ACC RTT CC-3′ , and the following conditions: 50 °C for 30 min, 95 °C for 15 min, followed by 35 cycles of amplification (30 sec at 94 °C, 30 sec at 45 °C, and 3 min at 72 °C), with a final extension of 10 min at 72 °C. The PCR amplicons were purified with EXO-SAP-it ™ (Affymetrix) and subsequently sequenced using the same forward primer as used for the RT-PCR with the ABI PRISM ™ BigDye Terminator Cycle Sequencing Reaction kit (Applied Biosystems).
The complete genome sequences for eleven strains were determined using the primers as described in Supplementary Table S1. For 17 other complete G2P [4] genomes next-generation sequencing technology was used. For this, viral RNA was first amplified using the Whole transcriptome kit (WTA2, Sigma-Aldrich) and sequence libraries were prepared using the KAPA library preparation kit (Kapa Biosystems) and sequenced in paired-end mode on a HiSeq 2500 ™ (Illumina) for 151 cycles. Sequences were trimmed for quality and adapters with Trimmomatic v0.32 31 and assembled using SPAdes v3.1 32 . Rotavirus contigs were identified using Blastn 33 . Sequences have been deposited in GenBank and are available under accession numbers KR705141-KR705509. Maximum likelihood phylogenetic trees with 500 bootstrap replicates were constructed using a generalized time-reversible (GTR) model allowing for invariant sites in Mega 6 34 . For all the trees, reference strains were used, as well as a selection of contemporary strains which were genetically closely related to the Belgian G2P [4] strains under investigation.