Uncovering the First Atypical DS-1-like G1P[8] Rotavirus Strains That Circulated during Pre-Rotavirus Vaccine Introduction Era in South Africa

Emergence of DS-1-like G1P[8] group A rotavirus (RVA) strains during post-rotavirus vaccination period has recently been reported in several countries. This study demonstrates, for the first time, rare atypical DS-1-like G1P[8] RVA strains that circulated in 2008 during pre-vaccine era in South Africa. Rotavirus positive samples were subjected to whole-genome sequencing. Two G1P[8] strains (RVA/Human-wt/ZAF/UFS-NGS-MRC-DPRU1971/2008/G1P[8] and RVA/Human-wt/ZAF/UFS-NGS-MRC-DPRU1973/2008/G1P[8]) possessed a DS-1-like genome constellation background (I2-R2-C2-M2-A2-N2-T2-E2-H2). The outer VP4 and VP7 capsid genes of the two South African G1P[8] strains had the highest nucleotide (amino acid) nt (aa) identities of 99.6–99.9% (99.1–100%) with the VP4 and the VP7 genes of a locally circulating South African strain, RVA/Human-wt/ZAF/MRC-DPRU1039/2008/G1P[8]. All the internal backbone genes (VP1–VP3, VP6, and NSP1-NSP5) had the highest nt (aa) identities with cognate internal genes of another locally circulating South African strain, RVA/Human-wt/ZAF/MRC-DPRU2344/2008/G2P[6]. The two study strains emerged through reassortment mechanism involving locally circulating South African strains, as they were distinctly unrelated to other reported atypical G1P[8] strains. The identification of these G1P[8] double-gene reassortants during the pre-vaccination period strongly supports natural RVA evolutionary mechanisms of the RVA genome. There is a need to maintain long-term whole-genome surveillance to monitor such atypical strains.


Full-Genome Constellation Analysis
Whole-gene sequences of the 11 genes of strains, RVA/Human-wt/ZAF/UFS-NGS-MRC-DPRU1971/ G1P [8] and RVA/Human-wt/ZAF/UFS-NGS-MRC-DPRU1973/G1P [8], were determined, and their genotype constellations were revealed as G1-P[8]-I2-R2-C2-M2-A2-N2-T2-E2-H2 (Table 1). The sizes of full-length segments 1 to 11 and their respective open reading frames (ORFs) for the two study strains were determined ( Table 1). The ORF sequences for all the 11 genes of these two South African atypical G1P [8] strains were deposited in GenBank under accession numbers MT163245-MT163266. Color codes indicate genogroup attribution. Green color represents the genotype associated with the Wa-like genogroup, while red color represents the genotype belonging to the DS-1-like genogroup. The nomenclature of the RV strains indicates RV group, species where the strain was isolated, name of the country where the strain was originally isolated, common name, year of isolation, and genotypes for genome segments four and nine as proposed by the Rotavirus Classification Working Group (RCWG) [12]. ORF = open reading frame.

Neutralization Epitopes
Atypical G1P[8] strains Figure 2. Alignment of antigenic residues in VP7 between the strains contained in Rotarix ® and RotaTeq ® and wild type G1 strains. Antigenic residues are divided in three epitopes (7-1a, 7-1b, and 7-2). Amino acids that differ between Rotarix ® and RotaTeq ® are indicated in boldface. Sky blue colored residues are residues that are different from both Rotarix ® and RotaTeq ® , green colored residues are different from Rotarix ® , and brown colored residues are different from RotaTeq ® . Amino acid changes that have been shown to escape neutralization with monoclonal antibodies are indicated with a black dot. Atypical G1P [8] and countries of detection are indicated on the left side of the figure. South Africa atypical G1P [8] strains are in boldface characters.

Analysis of the VP4 Neutralization Epitopes
The VP4 spike protein is cleaved by trypsin into two distinct structural proteins, VP8* and VP5* [2]. Analysis of the two South African study strains' VP4 sequences showed a conserved trypsin cleavage site (arginine) at positions 230, 240, and 581 [39]. Furthermore, the neutralization epitopes in the VP8* and the VP5* regions were analyzed. The VP8* region has four (8-1 to 8-4) neutralization epitopes, while VP5* has five (5-1 to 5-5) (Figure 4) [40]. Comparison of the two South African P [8] strains relative to the Rotarix ® and the RotaTeq ® P [8] sequences displayed 32 and 35 identical amino acid residues, respectively, spanning the VP4 antigenic epitopes (Figure 4). Amino acid differences between the two P [8] study strains and the P[8] component of vaccine strains were only identified in 8-1, 8-2, and 8-3 VP8* epitopes. Five amino acid differences (E150D, N195G, S125N, S131R, and N135D) were identified in the study strains in relation to Rotarix ® P [8] strain, while two amino acid differences (E150D and D195G) were identified relative to P [8] strain of RotaTeq ® (Figure 4). Analysis with VP4 epitopes of other atypical G1P [8] strains identified similar amino acid residues with the exception of position 113 in the 8-3 epitope, whereby asparagine was observed in the study strains while other atypical strains had either an aspartate or serine at this position ( Figure 4). Further analysis of the study strain's VP4 neutralization epitopes with corresponding VP4 neutralization epitopes of globally selected P[8]-lineage-III strains identified two amino acid differences (S146G and N113D) ( Figure 4). Pathogens 2020, 9, x FOR PEER REVIEW 8 of 18

Discussion
This study described the first pre-vaccine era atypical reassortant G1P [8] strains, whose outer gene segments (VP7 and VP4) expressed a Wa-like genotype, whereas the backbone genes expressed a DS-1-like genotype constellation. Analysis of the whole-genome constellation showed that genetic reassortment mechanism generated the DS-1-like G1P [8] strains. Rotavirus reassortment events are mainly facilitated by the segmentation inherent in the RV genome [2], which can generate rare or novel RV strains and hence contribute to the vast RVA diversity [22]. Wa-like and DS-1-like intergenogroup reassortment events involving G1P [8] and DS-1-like genotype constellation have been described recently in six countries: Brazil, Japan, Philippines, Thailand, Vietnam, and Malawi [23][24][25][26][27][28][29][30]. According to literature, viable atypical reassortant strains can occur under natural conditions

Discussion
This study described the first pre-vaccine era atypical reassortant G1P [8] strains, whose outer gene segments (VP7 and VP4) expressed a Wa-like genotype, whereas the backbone genes expressed a DS-1-like genotype constellation. Analysis of the whole-genome constellation showed that genetic reassortment mechanism generated the DS-1-like G1P [8] strains. Rotavirus reassortment events are mainly facilitated by the segmentation inherent in the RV genome [2], which can generate rare or novel RV strains and hence contribute to the vast RVA diversity [22]. Wa-like and DS-1-like intergenogroup reassortment events involving G1P [8] and DS-1-like genotype constellation have been described recently in six countries: Brazil, Japan, Philippines, Thailand, Vietnam, and Malawi [23][24][25][26][27][28][29][30]. According to literature, viable atypical reassortant strains can occur under natural conditions involving Wa-like G1P [8] or G3P [8] outer capsid genes expressing DS-1-like genetic background [41][42][43][44]. This study identified two DS-1-like G1P [8] strains in the course of the ongoing AEVGI whole-genome characterization of South African RVA strains. While reported in low frequencies and limited settings in Brazil (1.6% during 2013-2017 seasons) [22], Thailand (0.4% during 2012-2014 seasons) [27], and Vietnam (14% during 2012/2013 season) [28,29], 31-62% of these DS-1-like G1P [8] strains accounted for RVA positive strains circulating across selected regions in Japan [26], and 40% of randomly sampled post-vaccine samples were reported in Malawi [24]. Such atypical reassortant strains have the potential to predominate in circulation. G1P [4] strains suggested to have emerged from intergenogroup reassortment events accounted for 41% of RVA strains circulating in the peak months of the 2001 RVA season in Detroit, USA [45], whereas a surge in G3P [4] strains also presumed to have emerged from intergenogroup reassortment events were detected in Brazil at 36% [46] and in Ghana at 64% [47]. The two South African study strains were identified during the pre-RVA vaccination period in South Africa in contrast to the previously reported atypical strains found during the post-RVA vaccination period. This implies that the reassortment events that led to the emergence of the South African atypical G1P [8] strains may not necessarily be driven by vaccine-induced selective pressure but by natural evolutionary processes of RVA genome.
Vaccine escape mutants can result due to mutations occurring in well-known VP7 neutralization epitope regions [48]. Host-antigen binding interactions involving human G1 strains are significantly impacted by mutations occurring at positions 94, 97, 147, and 291 [48]. The identified N94S substitution involving the substitution of asparagine (N) with a serine(S), which are both polar non-charged amino acid residues [49], may not significantly alter the overall morphology of the protein surface. However, since asparagine is usually N-glycosylated, there is a likely loss of glycosylation site, which could have a wide-ranging impact on the immunogenicity of the 7-1a epitope [50]. The D97E amino acid substitution involving polar negatively charged residues, aspartate (D) and glutamate (E), is likely to be a silent nucleotide change [49]. Similarly, a K291R substitution involving lysine (K), an amphipathic polar amino acid, and arginine (R), a positively charged amino acid, is unlikely to have a far-reaching structural effect on the VP7 protein surfaces. However, M217T substitution was identified resulting in substitution of methionine, a non-polar residue to threonine, a polar amino acid residue, which could likely result in significant changes in biochemical properties of VP7 [49]. This M217T substitution was also present in earlier strains as well as some post-vaccine strains that were included for analysis, and the role it plays in driving epidemiological fitness of G1 strains is not fully resolved. In the host cell, trypsin-like proteases cleave the VP4 spike protein into two structural domains (VP8* and VP5*) [2]. Four surface-exposed antigenic epitopes (8-1 to 8-4) have been described in the VP8* region, while five antigenic epitopes (5-1 to  in the VP5* region have been documented [40]. The amino acid changes E150D, N195G, S125N, and N135D that were observed relative to the vaccine strains were conservative. However, a S131R substitution that resulted in a change in polarity might play a role in escape of host immunity [49]. Another amino acid substitution R131S resulting in a change in charge from positively charged amino acid to non-charged amino acid that was identified when comparison was made against globally selected lineage-III VP4 strains might impact vaccine escape effect [49].

Ethics Approval
The study was approved under ethics number UFS-HSD2018/0510/3107 by the Health Sciences Research Ethics Committee (HSREC) of the University of Free State, Bloemfontein, South Africa. The patient identities and demographics were de-linked from their unique laboratory identifiers to ensure confidentiality.

Sample Collection
Rotavirus positive stool samples from children under five years of age treated for gastroenteritis at Dr. George Mukhari Hospital, Pretoria North, South Africa and conventionally genotyped as G1P [8] were sourced from archival storage (2002 to 2017) of South Africa Medical Research Council-Diarrheal Pathogens Research Unit (MRC-DPRU), a WHO Rotavirus Regional Reference Laboratory (WHO-RRL) in Pretoria, South Africa. The two stool samples that were later genotyped as DS-1-like G1P [8] strains were collected from 6-month female and 12-month male children on 15 and 16 May 2008, respectively, from Soshanguve, Pretoria.

Extraction and Purification of Double-Stranded RNA
The extraction of RV ds-RNA was conducted by utilizing a previously described method [51], albeit with modifications (UFS-NGS unit extraction SOP). Briefly, a pea size (~100 mg) sample of stool was added to 200 µL of phosphate-buffered saline (PBS) solution, pH 7.2 (Sigma-Aldrich ® , St Louis, MO, USA). The solution was mixed by pulse-vortexing for five seconds. A 1 mL volume of TRI-Reagent ® -LS (Molecular Research Center, Inc, Cincinnati, OH, USA) was added and let to stand for five minutes. Phase separation was achieved by addition of 270 µL of chloroform (Sigma-Aldrich ® , St Louis, MO, USA). Afterward, centrifugation for 13,000 revolutions per minute (RPM) was performed for 20 min at 4 • C in a temperature-controlled microcentrifuge (Eppendorf microcentrifuge 5427R, Hamburg, Germany). A volume of 1 mL isopropanol (Sigma-Aldrich ® , St Louis, MO, USA) was added to the supernatant, and centrifugation was performed at 13,000 RPM for 30 min at room temperature. The supernatant was poured off, and the tubes were let to dry for 10 min, after which 95 µL of elution buffer (EB) from the MinElute Gel extraction kit (Qiagen, Hilden, Germany) was added. A 30 µL volume of of 8M LiCl2 (Sigma, St. Louis, MO, USA) was added, and the solution was precipitated for 16 h at 4 • C in a water bath in a Tupperware box. The MinElute gel extraction kit (Qiagen, Hilden, Germany) was used to purify the extracted RNA according to manufacturer's instructions, and 1% 0.5 X TBE agarose gel stained with Pronasafe (Condalab, UK) electrophoresis was used to verify the integrity and the enrichment of dsRNA, which was visualized on a G:Box Syngene UV transilluminator (Syngene, Cambridge, UK).

Synthesis and Purification of Complementary DNA (cDNA)
The Maxima H Minus Double-Stranded cDNA Synthesis Kit (Thermo Fischer Scientific, Waltham, MA) was utilized to synthesize cDNA from the extracted viral RNA. Briefly, denaturation at 95 • C for 5 min of the extracted RNA was performed followed by addition of 1 µL of 100 µM Random Hexamer primer. Incubation was performed in a thermocycler at 65 • C for five minutes. The First-Strand Reaction mix (5 µL) and the First Strand Enzyme Mix (1 µL) were added, and the solution was incubated at 25 • C for 10 min followed by 2 h at 50 • C, and then the reaction was terminated by heating at 85 • C for 5 min. A volume of 55 µL of nuclease-free water, 20 µL of 5X Second Strand Reaction Mix, and 5 µL of Second Strand Reaction Mix was then added. The solution was then incubated at 16 • C for 60 min, after which the reaction was stopped by adding 6 µL 0.5M EDTA. A volume of 10 µL RNAse I was then added, and the synthesized cDNA was incubated for five minutes at room temperature. Subsequently, the MSB ® Spin PCRapace (Stratec) Purification Kit was used to purify the synthesized cDNA.

DNA Library Preparation and Whole-Genome Sequencing
The Nextera ® XT DNA Library Preparation Kit (Illumina, San Diego, California, US) was utilized to prepare DNA libraries by following manufacturer's instructions. Briefly, the genomic DNA was tagmented by using the Nextera ® transposome enzyme, and the tagmented DNA was subsequently amplified using a limited-cycle PCR program. The DNA libraries were cleaned-up using AMPure XP magnetic beads (Beckman Coulter, Pasadena, CA, USA) and 80% freshly prepared ethanol. The quantity of the DNA was determined using Qubit 2.0 fluorometer (Invitrogen, Carlsbad, CA, USA), and the quality of the libraries and the fragment sizes was assessed using Agilent 2100 BioAnalyzer ® (Agilent Technologies, Waldbronn, Germany) by following the manufacturer's specified protocol. The Illumina MiSeq ® sequencer (Illumina, San Diego, CA, USA) was utilized to perform paired-end nucleotide sequencing (301 × 2) for 600 cycles by using a MiSeq Reagent Kit v3 at the University of the Free State-Next Generation Sequencing (UFS-NGS) Unit, Bloemfontein, South Africa.

Genome Assembly
Geneious Prime ® software, version 2019.1.1 (Biomatters, https://www.geneious.com/; [52]) was used for genome assembly. Briefly, for use with the reference mapping tools integrated in Geneious Prime version 2019.1.1, the default medium sensitivity parameter was selected to generate contigs from the FASTQ files data generated by the Illumina MiSeq ® instrument. Complementary RV genome assembly was also performed using an in-house genome assembly pipeline and CLC Genomics Workbench 12 (https://www.qiagenbioinformatics.com/).

Determination of Rotavirus Whole-Genotype Constellations
The genotype of each gene segment was determined using Rota C, v 2.0 [13], an online server for genotyping RVA strains. This was used to generate the full genotype constellations for each RV strain.

Phylogenetic Analyses
Complete sequences for each gene segment were aligned and sequence comparisons performed as described previously [53][54][55]. Multiple sequence alignments were implemented utilizing the MUSCLE package in Molecular Evolutionary Genetics Analysis (MEGA) 6 software ( [56]; http: //www.megasoftware.net/). Upon alignment, the DNA Model Test program in MEGA 6 was used to determine the evolutionary model that best fits each gene sequence datasets. The models identified as best fitting with the sequence data for the indicated genes using the Corrected Akaike Information Criterion (AICc) were as follows: GTR+G+I (VP7, VP4, VP6, VP1, VP2, VP3, NSP1, NSP2, NSP3) and HKY+G+I (NSP4 and NSP5). These models were utilized in maximum-likelihood trees' construction using MEGA 6 with 1000 bootstrap replicates to estimate branch support. Genetic distance matrices were prepared using the p-distance algorithm of MEGA 6 software [56]. In addition to the two whole-genome sequences of the strains in this study, other cognate sequences were acquired from GenBank ( [57]; http://www.ncbi.nlm.gov/genbank). Further phylogenetic analysis by geographical regions Africa (Eastern Africa, Southern Africa, and West Africa), Asia, Americas, Europe, and Oceania was also performed. mVISTA software was used to visualize the comparative sequence similarities of concatenated whole-genome of genetically related strains [58].

Conclusions
Whole-gene analyses showed that the South African DS-1-like G1P [8] strains were generated involving locally circulating G2P [6] strains by acquiring the VP7 and the VP4 outer capsid proteins of locally co-circulating G1P [8] strains. Similar to their pre-vaccine era detection in Vietnam and Philippines, the identification of these atypical DS-1-like G1P [8] strains during the pre-vaccine period in South Africa, as opposed to their detection during post-vaccination era in selected settings in Brazil (Sao Paulo and Goias in 2013), Japan (Okayama, Aichi, Akita, Kyoto, and Osaka Prefectures in 2012), Thailand (Phetchabun and Sukhothai in 2013), and Malawi (Blantyre in 2013/2014), suggests that they originated from natural evolutionary processes of RVA genome. Whole-genome surveillance of RVA genotypes is imperative to understand the occurrence rate, the mechanisms that drive emergence of such atypical strains, and their epidemiological fitness as well as to assess the effect of vaccine selective pressure in shaping the antigenic landscape of RVA strains.