Whole-Genome Analyses Identifies Multiple Reassortant Rotavirus Strains in Rwanda Post-Vaccine Introduction

Children in low-and middle-income countries, including Rwanda, experience a greater burden of rotavirus disease relative to developed countries. Evolutionary mechanisms leading to multiple reassortant rotavirus strains have been documented over time which influence the diversity and evolutionary dynamics of novel rotaviruses. Comprehensive rotavirus whole-genome analysis was conducted on 158 rotavirus group A (RVA) samples collected pre- and post-vaccine introduction in children less than five years in Rwanda. Of these RVA positive samples, five strains with the genotype constellations G4P[4]-I1-R2-C2-M2-A2-N2-T1-E1-H2 (n = 1), G9P[4]-I1-R2-C2-M2-A1-N1-T1-E1-H1 (n = 1), G12P[8]-I1-R2-C2-M1-A1-N2-T1-E2-H3 (n = 2) and G12P[8]-I1-R1-C1-M1-A2-N2-T2-E1-H1 (n = 1), with double and triple gene reassortant rotavirus strains were identified. Phylogenetic analysis revealed a close relationship between the Rwandan strains and cognate human RVA strains as well as the RotaTeq® vaccine strains in the VP1, VP2, NSP2, NSP4 and NSP5 gene segments. Pairwise analyses revealed multiple differences in amino acid residues of the VP7 and VP4 antigenic regions of the RotaTeq® vaccine strain and representative Rwandan study strains. Although the impact of such amino acid changes on the effectiveness of rotavirus vaccines has not been fully explored, this analysis underlines the potential of rotavirus whole-genome analysis by enhancing knowledge and understanding of intergenogroup reassortant strains circulating in Rwanda post vaccine introduction.


Introduction
Group A rotaviruses (RVA) are the most common cause of severe virus-induced diarrhea in infants and children less than five years of age in both developed and developing countries, accounting for over 128,500 deaths globally in 2016, of which 104,733 occurred in sub-Saharan Africa [1]. Children living in low-and middle-income countries have a higher disease burden due to poor living conditions, inadequate sanitation and limited supply of clean drinking water, among other factors [2,3]. Prior to the introduction of the vaccine RotaTeq ® (RV5, Merck & Co. Inc., Kenilworth, NJ, USA) in Rwanda, approximately 3500 RVA related mortality was reported in children annually, accounting for 8.8% of all mortality rate in children less than five years of age [4]. RotaTeq ® was introduced into the Rwandan National Immunization Program in May 2012 with a vaccine coverage of 99% within a year after introduction [5,6]. During the first three years post-RotaTeq ® introduction, the proportion of total diarrheal hospitalization due to rotavirus declined by 25-44% among children less than five years in the Eastern Province of Rwanda suggesting a positive vaccine impact on the population [7].
Genome reassortment between human RVA strains have been occasionally reported among circulating strains globally [24][25][26][27][28][29][30]. Such strains are suspected to influence the low effectiveness of rotavirus vaccines in regions with a significant strain diversity [31]. The RVA strains derived from intergenogroup reassortment are normally influenced by co-infection with multiple RVA strains, possibly generating strains with novel genome constellations [8,25,[32][33][34]. However, such strains are reported to exhibit decreased evolutionary fitness compared to strains having a pure Wa-like or pure DS-1-like genotype constellations and thus transmit less rapidly across human populations [26,29]. In an effort to contribute to the understanding of dynamics of intergenogroup reassortment of RVA strains, the present study reports on five intergenogroup reassortant strains observed postvaccine introduction in Rwanda between 2013 and 2015 using a whole-genome sequencing approach.

Ethics Statement
This study was reviewed and approved by the Health Sciences Research Ethics Committee of the UFS, Bloemfontein on the 15 October 2019 and assigned an ethics number (UFS-HSD2019/1601/2810). The diarrheal stool samples were collected as a routine diagnostic clinical sample when the parents brought their child to the hospital for clinical management, requiring no written informed consent. The archived RVA-positive samples were anonymized and utilized for strain characterization.

Sample Collection
Stool samples (n = 158) were collected from hospitalized children presenting with acute gastroenteritis between 2011 and 2016 in Rwanda as part of the WHO/AFRO supported rotavirus surveillance program. These samples were retrieved from "African stool repository" established to archive stool samples as part of the WHO African surveillance network, maintained at the Diarrheal Pathogens Research Unit, WHO Rotavirus Regional Reference Laboratory (Pretoria, South Africa. Five samples had the reassortant genome constellations, which meet the criteria for this study. Two samples were collected from unvaccinated children aged 14 months (female) and 36 months (male) from the South and East Province of Rwanda, respectively. The remaining three samples were collected from children (males) who received all three doses of the RotaTeq ® vaccine, one aged 12 months (from the North Province) and two aged 24 months (both from the East Province). The samples were whole-genome sequenced and analyzed at the UFS-NGS Unit. The Rwandan strains described in this study were deposited in GenBank database under accession numbers MT163179-MT163266.

Preparation of Purified dsRNA and cDNA for Rotavirus Whole-Genome Sequencing
The dsRNA was extracted using a protocol previously described by Nyaga et al. [35], and purified using the Qiagen MinElute gel extraction kit (Qiagen, Hilden, Germany). The quality of the purified dsRNA was thereafter verified by 1% agarose gel electrophoresis prior to quantification using a BioDrop-µLITE spectrophotometer (Biodrop, Cambridge, UK). For the cDNA synthesis, the Maxima H Minus Double-Stranded cDNA Synthesis Kit (ThermoFisher Scientific, Waltham, MA, USA) was used in accordance to manufacturer's instructions with minor modifications. The modification included the denaturation of the dsRNA at 95 • C as a replacement for 65 • C for five minutes prior to synthesis of the first strand for 2 h at 50 • C in a thermocycler (Merck, Darmstadt, Germany).

DNA Library Preparations and Whole-Genome Sequencing
The Nextera ® XT DNA library preparation kit (Illumina Inc., San Diego, CA, USA) was used to construct DNA libraries for whole-genome sequencing, following the manufacturer's instructions. The libraries were dual-barcoded using Nextera index kit (Illumina Inc.) and purified using AMPure XP magnetic beads (Beckman Coulter, Indianapolis, IN, USA) while simultaneously selecting for~300 bp DNA fragments and removing short library fragments. Fragment sizes were confirmed using a High Sensitivity DNA Kit (Agilent Technologies, Waldbronn, Germany) and run on a Bioanalyzer 2100 (Agilent Technologies). The libraries were quantified on a Qubit 3.0 Fluorimeter (Invitrogen, Carlsbad, CA, USA). Whole-genome sequencing was performed on an Illumina MiSeq (Illumina Inc., San Diego, CA, USA) platform for 600 cycles (301 × 2 paired ends) at a final concentration of 8pM and a PhiX (20pM) spike of 20% was used as a positive control.

Genome Constellations
In order to gain insight into the genetic variability among the study strains and their genetic relatedness with selected global reference RVA strains, the full-genome sequences of 158 Rwandan samples sequenced from the pre-(2011) and post-(2012-2016) vaccination period were determined. From that dataset pool, only five samples from the pool possessed reassortant constellations ( Table 1). The full or near full-length nucleotide sequences reads were assembled for all the 11 gene segments of the five surveillance strains. The length size of contigs and number of reads post-assembly are given ( Table 1)  Wa-like genogroups are represented by the green color, DS-1like is represented by the red color and AU-1-like is represented by the yellow color.

The VP4 and VP7 Antigenic Region Analyses
Amino acid changes at the VP7 and VP4 antigenic epitopes of circulating RVA strains and RVA vaccine strains can affect the ability of antibodies to neutralize virus infectivity and could also undermine vaccine effectiveness [40]. In this regard, we compared the antigenic differences between RotaTeq ® vaccine and the Rwandan study strains by analyzing the amino acid sequence of the VP7 and VP4. The VP7 gene contains two structurally defined antigenic epitope regions: 7-1 and 7-2 made up of 29 amino acid residues [41,42]. The 7-1 epitope is further subdivided into 7-1a and 7-1b. The VP7 epitopes of the Rwandan G4 strain was compared to the G4 VP7 protein of strain RVA/Vaccine/USA/RotaTeq-BrB-9/1996/G4P75, which showed 27 amino acid differences distributed across the VP7 epitope regions (Figure 1A,B). Only two residues in position 190 (within the 7-2 region) and 291 (within the 7-1a region) were conserved between the RotaTeq ® vaccine strains and the Rwandan study strain ( Figure 1A). Amino acid substitutions from uncharged polar molecules to charged polar molecules were observed in three positions, T96D, T217E and S221D. Furthermore, amino acid substitution from charged polar molecule to uncharged polar molecule was also observed in two positions E97T and D211T. The polarity of the other amino acids was not affected by either of the substitutions. In the presence of trypsin, the VP4 is cleaved into two domains, the VP8* (8-1 to 8-4) and the VP5* (5-1 to 5-5) made up of 37 amino acid residues in the antigenic epitope regions [40,43]. The VP4 epitopes of the Rwandan strains was compared to the P [8] VP4 protein of strain RVA/Vaccine/USA/RotaTeq-WI79-4/1992/G6P1A8 (Figure 2). The Rwandan study strains differed from the RotaTeq ® vaccine strain in only three positions, E150D and D195G at 8-1 epitope region and L388I at 5-1 epitope region, while the rest of the residues in the epitope region were conserved (Figure 2). At position 150, the amino acid changed from a glutamic acid to aspartic acid. While the change at position 195 was from an aspartic acid (charged polar molecule) to a glycine molecule (nonpolar molecule) and the change at position 388 was from leucine to isoleucine.

Phylogenetic Analysis of the VP7 Gene of G4, G9 and G12
Phylogenetic trees were constructed for each of the 11 genome segments of the five Rwandan RVA reassortant strains in comparison with global reference strains from Gen-Bank (Figure 3, Figure 4 and Figure S1). The three Rwandan G12P [8] strains clustered in G12 lineage III, and shared 87.6-99.3% nucleotide (nt) similarity with other lineage III G12 strains ( Figure 3 and Figure S2). The G12 strains clustered together and shared 99.8-100% nt similarity amongst themselves. On the other hand, the Rwandan G9 strain clustered in G9 lineage III distantly from globally circulating G9 strains (81.7-92.1% nt similarity). The G4 strain clustered in G4 lineage I with circulating G4 global strains and exhibited 86.0% nt similarity with a RotaTeq ® vaccine strain.

Phylogenetic Analysis of the VP4 Gene of P[4] and P[8]
The VP4 gene of the two Rwandan P [4] strains and three Rwandan P [8]

Discussion
In the present study, we analyzed the whole-genome of five Rwandan rotavirus strains that have undergone intergenogroup reassortment identified as part of WHO supported on-going rotavirus sentinel surveillance. All the G12P [8] strains were from vaccinated children while the G9P [4] and G4P [4] strains were from non-vaccinated children. Bányai et al. [31] stated that most atypical RVA strains are result of natural intergenogroup reassortment between the Wa-like and DS-1-like strains due to the segmented nature of RVA. The emergence of these intergenogroup reassortant strains may be attributed to either, the lack of RNA polymerase proofreading ability or co-infection of multiple strains from various human RVA strains [8,15,[44][45][46]. Co-infections have been reported in high frequencies in several RVA strains across Africa [21,22,35,47].
The detection of the unusual G9P [4] strain (RVA/Human-wt/RWA/UFS-NGS:MRC-DPRU566/2013/G9P [4]) in Rwanda is noteworthy as it is more prevalent in South-East Asia, Japan and Central America and detected at very low frequency (2%) in Africa [23,[48][49][50][51][52][53][54]. Phylogenetic analysis revealed that the five Rwandan strains clustered mostly with strains circulating globally, with few exceptions that clustered separately. This suggests a direct importation of these variants from abroad rather than local emergence through multiple reassortment events between locally circulating strains. It is evident that the reassortment events of the five Rwandan strains mostly occurred with contemporary human rotavirus strains as they did not show sufficient evidence of animal/human reassortment. Furthermore, RVA/Human-wt/RWA/UFS-NGS:MRC-DPRU8020/2015/G12P [8] and RVA/Human-wt/RWA/UFS-NGS:MRC-DPRU9995/2015/G12P [8] exhibit a high similarity amongst each other across the genome and both human and RotaTeq ® vaccine strains in the VP2, NSP2, NSP4 and NSP5 gene segment. This finding suggests a reassortment between both human and RotaTeq ® vaccine strains might have also transpired. This observation is consistent with the findings of Rose et al. [55], who have reported that such reassortment events are expected considering the attenuation of RotaTeq ® vaccines and the segmented nature of the RVA genome [55].
In the construction of the RotaTeq ® vaccine in the 1980s, the G4 and P [8] components were included in the composition of this pentavalent vaccine [18]. When comparing the G4 and P [8] components of the Rwandan strains and the RotaTeq ® vaccine strains, we observed that the Rwandan G4 strain clustered distinctly to the RotaTeq ® vaccine strain in G4-lineage I while the Rwandan P [8] strains clustered in P [8] lineage III distantly from the RotaTeq®vaccine strain in P[8]-lineage II. This phenomenon can be attributed to the constant evolutionary changes that rotaviruses undergo. Hence the currently circulating RVA strains may cluster in VP7 and VP4 lineages different from the RVA vaccine strains [40]. Such changes may influence a varying selective pressure against these VP7 and VP4 lineages ultimately reducing vaccine effectiveness over time. The VP1, VP3, NSP1 and NSP3 genes of strain RVA/Human-wt/RWA/UFS-NGS:MRC-DPRU6235/2014/G4P [4] were phylogenetically distinct from strains circulating in other places thus suggesting that these gene segments may be unique.
Comparison of the deduced amino acid sequences of the G4P [4] and G12P [8] Rwandan strains and the RotaTeq®vaccine strain showed that none of the Rwandan strains were identical to the RotaTeq®vaccine strain. Amino acid substitutions were observed throughout the VP7 epitope regions excluding position 291 and 190, while the VP4 exhibited only three amino acid substitutions at position 150, 195 and 388. The amino acid substitution at position 96, 97, 211, 217 and 221 on the VP7 epitope region suggests a radical change in polarity. McDonald et al. [29] suggested that such substitutions do not change the genotype specificity of the rotavirus strain. However, they may influence the binding of neutralizing antibodies thus affecting the viral fitness through selection pressure. Despite the distant clustering between the Rwandan P [8] strains and the vaccine strain, only three changes were observed in the VP4 neutralizing epitope regions. Antigenic variations between rotavirus strains and the vaccine strains is often implicated in the decreased effectiveness of rotavirus vaccines in low-income countries such as Rwanda [56][57][58]. There is no evidence that the amino acid changes on the VP4 epitope regions are due to vaccination. The changes could have occurred as a natural evolutionary process as they were observed in other strains circulating globally over the years.

Conclusions
The detection of five intergenogroup reassortant Rwandan rotavirus strains from the whole-genome analysis further emphasizes the ubiquitous nature and diversity of RVA strains in circulation. Whether vaccine introduction is responsible for the observed reassortment events in vaccinated children or not, remains unknown as several natural factors can be attributed to the evolution of these RVA strains. The amino acid substitutions observed in the antigenic regions in the neutralizing epitope of the VP7 and VP4 proteins of the Rwandan strains when compared with the RotaTeq ® vaccine strain related to the reduced effectiveness of rotavirus vaccines in Rwanda as well as other low-income countries in the region. Continuous surveillance of rotavirus, using the whole-genome sequencing analysis, is very important to monitoring the impact of vaccine pressure, on the circulating rotavirus strains in African countries.

Study Limitations
The children in this study represents a small portion of the population of children that tested positive for rotavirus in Rwanda as part of a WHO whole-genome surveillance pilot study.