**Rotavirus A in Brazil: Molecular Epidemiology and Surveillance during 2018–2019**

**Meylin Bautista Gutierrez, Alexandre Madi Fialho, Adriana Gonçalves Maranhão, Fábio Correia Malta, Juliana da Silva Ribeiro de Andrade, Rosane Maria Santos de Assis, Sérgio da Silva e Mouta, Marize Pereira Miagostovich, José Paulo Gagliardi Leite and Tulio Machado Fumian \***

Laboratory of Comparative and Environmental Virology, Oswaldo Cruz Institute, Oswaldo Cruz Foundation, Avenida Brasil 4365, Rio de Janeiro 21040-900, Brazil; meylin.gutierrez@ioc.fiocruz.br (M.B.G.); amfialho@ioc.fiocruz.br (A.M.F.); adriana.maranhao@ioc.fiocruz.br (A.G.M.); fabio.malta@ioc.fiocruz.br (F.C.M.); juliana@ioc.fiocruz.br (J.d.S.R.d.A.); rmsassis@ioc.fiocruz.br (R.M.S.d.A.); mouta@ioc.fiocruz.br (S.d.S.e.M.); marizepm@ioc.fiocruz.br (M.P.M.); jpgleite@ioc.fiocruz.br (J.P.G.L.)

**\*** Correspondence: tuliomf@ioc.fiocruz.br or fumiantm@gmail.com; Tel.: +55-(21)-25621817

Received: 12 May 2020; Accepted: 7 June 2020; Published: 27 June 2020

**Abstract:** Rotavirus A (RVA) vaccines succeeded in lowering the burden of acute gastroenteritis (AGE) worldwide, especially preventing severe disease and mortality. In 2019, Brazil completed 13 years of RVA vaccine implementation (Rotarix™) within the National Immunization Program (NIP), and as reported elsewhere, the use of Rotarix™ in the country has reduced childhood mortality and morbidity due to AGE. Even though both marketed vaccines are widely distributed, the surveillance of RVA causing AGE and the monitoring of circulating genotypes are important tools to keep tracking the epidemiological scenario and vaccines impact. Thus, our study investigated RVA epidemiological features, viral load and G and P genotypes circulation in children and adults presenting AGE symptoms in eleven states from three out of five regions in Brazil. By using TaqMan®-based one-step RT-qPCR, we investigated a total of 1536 stool samples collected from symptomatic inpatients, emergency department visits and outpatients from January 2018 to December 2019. G and P genotypes of RVA-positive samples were genetically characterized by multiplex RT-PCR or by nearly complete fragment sequencing. We detected RVA in 12% of samples, 10.5% in 2018 and 13.7% in 2019. A marked winter/spring seasonality was observed, especially in Southern Brazil. The most affected age group was children aged >24–60 months, with a positivity rate of 18.8% (*p* < 0.05). Evaluating shedding, we found a statistically lower RVA viral load in stool samples collected from children aged up to six months compared to the other age groups (*p* < 0.05). The genotype G3P[8] was the most prevalent during the two years (83.7% in 2018 and 65.5% in 2019), and nucleotide sequencing of some strains demonstrated that they belonged to the emergent equine-like G3P[8] genotype. The dominance of an emergent genotype causing AGE reinforces the need for continuous epidemiological surveillance to assess the impact of mass RVA immunization as well as to monitor the emergence of novel genotypes.

**Keywords:** acute gastroenteritis; rotavirus A; incidence; genotyping; Brazil

#### **1. Introduction**

Acute gastroenteritis (AGE) remains as a major cause of mortality in children under five years old worldwide [1,2]. Among the AGE-causing pathogens, rotavirus A (RVA) is one of the leading agents, responsible for approximately 200,000 deaths per year among children <5 years old in developing countries [3–5]. Regarding severe disease, RVA accounts for around 20% and 40% of all AGE-hospitalization in countries with and without RVA vaccines implemented, respectively [6,7]. Currently, four World Health Organization (WHO)-prequalified live-attenuated oral RVA vaccines are available internationally—Rotarix™, RotaTeq™, Rotavac™, and RotaSiil™—and over 100 countries have introduced one of these vaccines into their national immunization program [8] (https://www.who. int/immunization/diseases/rotavirus/en/).

Rotaviruses belong to the *Reoviridae* family, genus *Rotavirus*. While nine rotaviruses species have been described (A–I), RVA is by far the most important species infecting humans worldwide [9,10]. The non-enveloped triple-layered viral particle has 70–75 nm in diameter with 11 segmented double-stranded RNA (dsRNA) genes, encoding for six structural (VP1-VP4, VP6, VP7) and depending on the strain, five or six non-structural proteins (NSP1-NSP5 or NSP6) [11]. Genetically, RVA is classified into G- and P-types, based on nucleotide sequence of genomic segments coding VP7 and VP4 proteins (binary classification), and currently there have been described 36 G- and 51 P-types [12]. Although many G and P combination would be possible to emerge, a few genotypes (G1P[8], G2P[4], G3P[8], G4P[8], G9P[8], and G12P[8]) have prevailed worldwide causing the majority of RVA infections in children [13–15].

Brazil has implemented the Rotarix™ vaccine in the National Immunization Program (NIP) in March 2006, which led to a significant reduction of diarrhea-associated mortality and hospitalization [16–18]. Linhares et al. [19] demonstrated the higher effectiveness of Rotarix™ among Brazilian infants aged up to 12 months and decreasing in older children. Concerning the genotype distribution in Brazil after the introduction of Rotarix™, G2P[4] was by far the most prevalent genotype detected until 2010. From 2011 onwards, a gradual decrease in the prevalence of G2P[4] was observed, being replaced by G3, G9, and G12 harboring a P[8]-type [20–23]. Nevertheless, unusual RVA genotypes have been frequently detected, such as: G3[P6], G12[P6], G8P[4], and G8P[6] and more recently the equine-like G3P[8] [17,23,24]. Similarly, recent studies from other countries have reported the detection of rare RVA genotype combination [25–28].

It has been demonstrated that the distribution of RVA genotypes over the years is characterized by natural and cyclical genotype fluctuations [20,29,30]. However, the selective pressure due to mass RVA vaccination could favor specific G and P combinations [9,31]. Therefore, the new and dynamic epidemiological scenario reinforces the need to continuously document RVA prevalence in AGE cases, molecular epidemiology and the potential emergence of unusual genotypes.

Our study investigated RVA prevalence, features and the molecular characterization of G and P genotypes among patients with AGE from three regions (Southern, Southeastern and Northeastern) in Brazil, 2018–2019. RVA was detected and quantified by quantitative RT-PCR (RT-qPCR) from diarrheic stool samples received from eleven Brazilian states, and G and P genotypes were determined by multiplex one-step RT-PCR or sequencing.

#### **2. Materials and Methods**

#### *2.1. Stool Collection and Ethics Statements*

This study included stool samples that were collected between January 2018 and December 2019 from children and adults with symptoms of AGE, characterized as ≥three liquid/semi liquid evacuations in a 24 h period. Inpatients and outpatients diarrheic stool samples were collected from eleven states from three regions of Brazil: Southern, Southeastern, and Northeastern. Samples were systematically sent together with clinical-epidemiological records to the Regional Rotavirus Reference Laboratory–Laboratory of Comparative and Environmental Virology (RRRL–LVCA). The laboratory is part of the ongoing national network for AGE surveillance and coordinated by General Coordination of Public Health Laboratories, Brazilian Ministry of Health.

This study is approved by the Ethics Committee of the Oswaldo Cruz Foundation (FIOCRUZ), number CAAE: 94144918.3.0000.5248. The surveillance is performed through a hierarchical network in which samples are provided by medical request in hospitals and health centers, monitored by the Brazilian Unified Health System (SUS). Patients' data were maintained anonymously and securely.

#### *2.2. Viral RNA Extraction*

Viral RNA was purified from 140 µL of clarified stool suspension (10% *w*/*v*) prepared with Tris-calcium buffer (pH = 7.2). Samples were subjected to an automatic nucleic acid extraction procedure using a QIAamp® Viral RNA Mini kit (QIAGEN, CA, USA) and a QIAcube® automated system (QIAGEN), according to the manufacturer's instructions. RVA RNA was eluted in 60 µL of the elution buffer AVE. The isolated RNA was immediately stored at −80 ◦C until the molecular analysis. In each extraction procedure, RNAse/DNAse-free water was used as negative control.

#### *2.3. RVA Detection and Quantification*

RVA was detected and quantified by using a TaqMan®-based quantitative one step PCR (RT-qPCR) with primers and probe targeting the conserved NSP3 segment, according to Zeng et al. (2008). Briefly, RT-qPCR reactions were performed with 5 µL of the extracted RNA in a final volume of 25 µL using the SuperScript™ III Platinum™ One-Step qRT-PCR Kit (ThermoFisher Scientific, Invitrogen Division, Carlsbad, CA, USA) in the Applied Biosystems® 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). NSP3 primers and probe final concentrations used were 0.8 and 0.5 µM, respectively. The thermal cycling conditions were carried out as follows: RT step at 55 ◦C for 30 min, an initial denaturation step at 95 ◦C for 10 min and 40 cycles of PCR amplification at 95 ◦C for 15 s and 60 ◦C for 1 min. Samples that crossed the threshold line showing a characteristic sigmoid curve were regarded as positive. All runs included negative and non-template controls, and a standard curve with serial dilutions (106–10<sup>1</sup> ) of double-stranded DNA fragments (gBlock® Gene Fragment, Integrated DNA Technologies, Iowa, USA) containing the RVA NSP3 target region to ensure the correct interpretation of the results throughout the study. RVA viral loads were expressed as genome copies per gram (GC/g) of stool.

#### *2.4. Genotyping and Sequencing*

RVA-positive samples obtained by RT-qPCR were G- and P-genotyped using a one-step multiplex RT-PCR. The reactions were performed using the Qiagen One Step RT-PCR kit (Qiagen), using forward conserved primers VP7uF or VP4uF and specific reverse primers for G types G1, G2, G3, G4, G9, and G12, or P types P[4], P[6], P[8], P[9], and P[10] as recommended by the Centers for Disease Control and Prevention, USA. The G- and P-genotypes were assigned based on different amplicon sizes [base pairs (bp)] using agarose gel analysis. Sanger sequencing was also used to characterize the nucleotide (nt) sequence of specific strains, such as non-typeable samples or the equine-like G3, using consensus primers directed to the conserved regions within the VP4 and VP7 genes. The amplicons fragments of 876 bp and 881 bp for VP4 and VP7, respectively, were purified using the ExoSAP clean-up kit (ThermoFisher Scientific) and sent to the FIOCRUZ Institutional Platform for DNA sequencing (PDTIS). All primers used for RVA genotyping were based on previously studies [32–34].

#### *2.5. Phylogenetic Analysis*

Chromatogram analysis and consensus sequences were obtained using Geneious Prime (Biomatters Ltd., Auckland, New Zealand). RVA genotypes were confirmed in terms of closest homology sequence using Basic Local Alignment Search Tool (BLAST). Phylogenetic trees were constructed using the maximum likelihood method and the Kimura two-parameter model (2000 bootstrap replications for branch support) in MEGA X v. 10.1.7 [35], with RVA reference sequences obtained from the National Center for Biotechnology Information (NCBI) database. Nucleotide sequences obtained from clinical samples were submitted to NCBI GenBank (accession numbers: MT386419 to MT386453).

#### *2.6. Statistical Analysis*

Statistical analyses were performed using GraphPad Prism software v. 8.4.1 (GraphPad Software, San Diego, CA, USA). As appropriate, Mann–Whitney U test, Chi-squared or Fisher test was used to

assess significant difference between RVA detection rates, years of collecting samples and age groups, as well as to compare RVA viral load according to different age groups. A *p* value < 0.05 was considered to be statistically significant.

#### **3. Results**

#### *3.1. Rotavirus A Epidemiology*

During the two-year period of this study (2018–2019), a total of 1536 stool samples were collected from symptomatic inpatients with AGE (1161 and 375 from children and adults, respectively). Overall, we detected RVA in 12% of samples (n = 185), 10.5% in 2018 and 13.7% in 2019. We observed a slight increase in RVA incidence in 2019, but without statistical significance (*p* = 0.053). Except for three months in 2018 (April, June, and December), RVA circulated year-round, with monthly detection rates varying from 1.6% to 36.7% in May 2018 and September 2019, respectively (Figure 1A). In relation to seasonal patterns, we observed higher RVA circulation during winter/spring months, especially marked in Southern region states (Figure 1B,C), whilst RVA detections were lowest in autumn months.

**Figure 1.** Monthly distribution of tested acute gastroenteritis samples, rotavirus A (RVA)-positive samples and RVA detection rates in Brazil (**A**), Northeastern and South-eastern states (**B**), and Southern states (**C**), during 2018–2019.

In regard to regional analysis, higher RVA prevalence was observed in the Northeast region (18.7%) compared to Southeastern and Southern regions (3.4% and 12.5%, respectively). Comparing the two year of the study, RVA detection rates were higher in 2019 for the three regions, but only with statistical significance in Southeastern region (*p* = 0.022). Table 1 shows detailed analysis by regions and states. It is interesting to note that the two states of Southern region (Santa Catarina and Rio Grande do Sul) accounted for almost half of the AGE cases and RVA-positive samples (Figure 2).


**Table 1.** Number of tested and rotavirus-positive fecal samples through laboratory-based surveillance by region and state in Brazil during 2018 and 2019.

**Figure 2.** Map of Brazil highlighting the eleven states with sentinel surveillance service attended by the Rotavirus Regional Reference Laboratory, IOC, FIOCRUZ. Number of tested samples (**A**) and number of RVA-positive samples (**B**).

≤ Most of stool samples received were from children less than five years old, representing 72.1% (1108/1536) of the AGE cases. RVA detection rate was significantly higher among children aged between 24 and 60 months (18.8%) compared to the other age groups, where detection rates varied from 9.3% to 12.1% (Table 2). We also analyzed RVA viral load (GC/g of stool) among different age groups. The median values of RVA viral loads varied from 4.2 to 6.8 log<sup>10</sup> GC/g among the different age groups. RVA-positive samples showed viral load values statistically lower in AGE cases among children ≤6 months compared to older patients (*p* < 0.05) (Figure 3).

**Table 2.** Number of tested and rotavirus-positive fecal samples through laboratory-based surveillance by age group in Brazil during 2018–2019.


\* *p*-values were calculated between the age group of >24–60 and each other. All other combinations were not statistically different.

≤ ≤ ≤ **Figure 3.** Rotavirus A (RVA) viral load expressed as log<sup>10</sup> genome copies per gram of stool (log<sup>10</sup> GC/g) among different age groups in Brazil, 2018–2019. Box-and-whisker plots show the first and third quartiles (equivalent to the 5th and 95th percentiles), the median (the horizontal line in the box), and range of log<sup>10</sup> GC/g values. \* *p* ≤ 0.05; \*\* *p* ≤ 0.01; \*\*\*\* *p* ≤ 0.0001.

#### *3.2. RVA Genotyping*

A total of 186 RVA-positive samples were subjected to G and P genotyping by one-step multiplex RT-PCR. From these, 167 samples (89%) were successfully genotyped; 80 from 2018 and 87 from 2019. We characterized seven different RVA genotypes circulating during this study: G3P[8], G3P[6], G9P[8], G1P[8], G2P[6], G12P[6], and G6P[8]. G3P[8] was detected year-round and was by far, the most prevalent genotype, accounting for 83.8% (n = 67) of genotyped samples in 2018 and 65.5% (n = 57) in 2019 (Figure 4). Two other usual RVA genotypes were detected, but in lower prevalence—G1P[8] detected in one sample in 2018 and 2019, and G9P[8] detected in two and eight samples from 2018 and 2019, respectively. We also detected unusual G/P combinations, especially in 2019, as follows: G3P[6] in 6.3% of samples from 2018; G6P[8], G12P[6], and G2P[6] in 13.8%, 4.6% and 1.2% of samples from 2019 (Figure 4). G or P non-typed (NT) samples (GNTP8, GNTP6, and G3P[NT]) accounted for 5.4% of samples, and were represented mostly by samples with low RVA viral load (high Ct values).

In addition to RT-PCR genotyping, we sequenced some of the RVA-positive samples in order to get detailed information of the circulating strains and their respective lineages. We successfully obtained 22 and 21 consensus sequences of VP7 and VP4 genes, respectively. Phylogenetic analysis of the VP7 gene confirmed the characterization of Brazilian strains belonging to G3 and G6. Eighteen G3 strains from both years and from the three Brazilian regions were sequenced. From these, 94.4% (n = 17) of sequences clustered within the lineage 1, represented by equine-like G3P[8] strains. Our sequences were genetically related to previously detected equine-like G3P[8] strains from Brazil (KX469400) and other countries, such as Germany (KY000546), Slovakia (MN203563), Dominican Republic (MG652313), and Japan (LC47366). One G3 sequence clustered within lineage 3 that comprises the Wa-like G3P[8] group. The Brazilian Wa-like G3 sequence was closely related to strains from Brazil (KJ454454), Argentina (KJ583190 and KJ583201) and Hungary (JQ693568), with nt similarity varying from 98.4

to 99.8% (Figure 5A). The four G6 strains sequenced in our study, harboring a P[8]-type, clustered within lineage 1 showing moderate nt identity (97.8–98.1%) with G6P[8] strains detected in Bulgaria (KM590371 and KM590373) and with G6P[9] strains from Germany (KX880436) and Italy (KC152917). None of our G6 sequences clustered within the G6 lineage 3, that comprises human-bovine reassortant strains (Figure 5A).

**Figure 4.** Rotavirus A (RVA) genotypes distribution in Brazil, 2018 (**A**) and 2019 (**B**). Bi-monthly genotypes circulation during the two-year of study (**C**).

Phylogenetic analysis of 21 sequences of VP4 gene, demonstrated that, except for one, all P[8] Brazilian strains harboring two different G-types (G3 and G6) grouped into lineage 3. The 20 strains were closely related (99.2–99.6% of nt similarity) to P[8]-3 Brazilian strains isolated in 2016 (KX469415 and MH569765) and strains from other countries, such as USA (MF997038), Japan (LC477395), Spain (KU550282), Australia (KU059769), and Italy (MK158257). One strain was characterized into P[6] lineage 1, and was closely related (99.5–99.7% of nt similarity) to strains detected in Argentina (KJ583199), Iraq (JX891397), and China (MG78835) (Figure 5B).

**Figure 5.** Phylogenetic analyses based on VP7 and VP4 nucleotide (nt) sequences of circulating Brazilian rotavirus strains. Strains obtained in this study are marked with a black filled circle and names contain the register number, state, and collection date (M/Y). Reference strains were downloaded from GenBank and labeled with their accession number followed by country, register number, year, and genotype. Neighbor-joining phylogenetic trees of VP4 (**A**) and VP7 (**B**) were constructed with MEGA X software and bootstrap tests (2000 replicates) based on the Kimura two-parameter model. Bootstrap values above 70% are given at branch nodes.

#### **4. Discussion**

In this study, we provide laboratory-based RVA national surveillance in eleven states from three regions in Brazil, during 2018–2019. We tested 1536 AGE stool samples and found an overall RVA-positivity of 12%. RVA detection rates were higher during winter/spring months and among children aged 24–60 months. By far, G3P[8] was the most frequently detected genotype, and showed a year-round circulation.

Despite the development of vaccines, RVA are still a major cause of severe AGE in infants worldwide [5]. Here, we detected RVA in 10.5% and 13.7% of samples from 2018 and 2019, respectively. In Brazil, after Rotarix™ implementation, different studies have investigated RVA circulation among AGE cases. A study from the Enteric Diseases Laboratory at Adolfo Lutz Institute, one of the three Brazilian Reference Laboratory for RVA surveillance, reported annual RVA prevalence varying from 9.9% to 25.3% during 2013–2017, with AGE samples from five states in the Midwestern and part of the Southeastern and Southern regions [24]. A previous study from our group demonstrated an overall RVA positivity of 20.8% among children up to 12 years old between 2006 and 2017, with annual detection rates varying between 5% to 35% [20]. Studies conducted at Evandro Chagas Institute, the national and regional reference center for RVA surveillance in Northern Brazil, demonstrated RVA positivity rates of 33% in samples from six states from North Brazil, 2011–2012 [23], and 24.2% in samples collected between June 2012 and June 2015 [36]. However, it is worth mentioning that both studies involved children hospitalized for severe AGE. In Argentina, RVA positivity decreased from 26.8% to 13.6% comparing the pre- and post-vaccination periods [37]. Other studies performed elsewhere have described RVA detection rates varying from 8.4% to 23.2% [38–42].

RVA seasonality has been well defined, especially for temperate climate countries, where RVA peaks during dry and cold months. In tropical areas, RVA circulates year-round without marked peaks of infections [43,44]. In Brazil, we observed a year-round RVA circulation without marked seasonality, but high detections rates of RVA was observed during winter/spring months, in agreement with other studies [42,45]. RVA highest detection rate was observed in September 2019 (36.7%) in line with findings observed over a 21-year period in Brazil [20], and also with Luchs et al. [46] that demonstrated the peak of RV incidence in September during a five-year RVA surveillance study (2007–2012) in Brazil. As a continental-size country, we analyzed separately, RVA circulation in Southern states in comparison with Southeastern and Northeastern states (Figure 1B,C). We observed more clear peaks of RVA infections in winter/spring months (June 21st to December 20th) in Southern states (Rio Grande do Sul—RS, and Santa Catarina—SC) compared to Southeastern and Northeastern Brazil. This could be explained as both RS and SC states are in a subtropical area, characterized by different climate pattern compared to the other states. A three-year study conducted in Vietnam to access RVA epidemiology in AGE cases also demonstrated varied seasonally positivity, with different RVA-detection peaks among the three regions analyzed—North, Central, and South [47]. The fact that RVA usually peaks in September in Brazil, observed here and by others [20,46], is an important information to authorities to prepare strategies to reduce AGE impacts in the health system.

Regarding RVA infections among different age groups, we observed a significantly high positivity rate among children aged >24 and 60 months compared to other age groups. This shifting in the age of children more affected by RVA illness (older children) has been observed, especially in countries that have introduced RVA mass vaccination. Our data are consistent with previous findings reported from Brazil [20] and the USA [48,49]. In contrast, countries where RVA vaccines are yet to be introduced into national immunization programs, have reported the majority of RVA positive children (~90%) within the first 2 years of life [44,47]. By analyzing RVA shedding among the age groups, we found a statistically lower viral load among children less than six months (Figure 3). We believe that this lower viral load could be mostly explained by the passive protection mediated by breast milk maternal antibodies [50], but also by the higher effectiveness and prompt immune response generated by Rotarix™ after the oral doses administered at the age of 2 and 4 months [17]. However, this second hypothesis alone could not explain the high viral load among children aged >6 and 12 months. In addition, high Ct

values could indicate less severe disease [51]. In that study, authors demonstrated that the severity of diarrhea, determined by the Vesikari score, was significantly and negatively associated with Ct values of children stool samples.

Regarding RVA genotype characterization, we successfully identified G- and P-types in 89% of positive samples, by one-step multiplex RT-PCR and sequencing. By far, G3P[8] was the most prevalent genotype in both years. The phylogenetic analysis of the VP7 gene revealed that the majority of the Brazilian strains sequenced (94%) belong to the equine-like G3 genotype (G3-1). Moreover, all the P[8] strains sequenced clustered within the P[8]-3 lineage. This P[8]-3, harboring a G12-type, was the dominant strain in Brazil in 2014, detected in 75% of genotyped samples [52].

Emergent equine-like DS-1-like G3P[8] RVA strains were firstly identified in children with AGE in Australia in 2013 [53]. From 2013 onwards, the equine-like G3P[8] DS-1-like genotype has spread and become endemic worldwide [54–60]. In Brazil, the first evidence of the circulation of equine-like G3P[8] date from 2015, when Luchs et al. [24] detected the reassortant RVA strain in a touristic city of Southern Brazil, Foz do Iguaçu, that borders Argentina and Paraguay. Subsequently, these novel viruses quickly spread to other states in Brazil, being the most prevalent genotype in 2017 (66.2%). The occurrence of DS-1-like G3P[8] RVA strains was also reported in Amazon region, Northern Brazil in 2016 [60]. In the previous study from our group, we demonstrated the increase of G3P[8] from 2015, peaking in the last year of the study—2017. However, it was not investigated whether they belonged to the DS-1-like RVA group [20]. More recently, countries such as Australia, Italy and Pakistan, have demonstrated the high prevalence of the emergent equine-like G3P[8] genotype [42,45,61].

Atypical genotypes G3P[6], G6P[8], G2P[6], and G12P[6] were also detected as minor genotypes in our study. The phylogenetic analysis of the VP7 gene demonstrated that Brazilian G6-1 strains were closely related to strains circulating in Bulgaria and Italy [62,63]. The genotype G12P[6] characterized in our study has been frequently detected in Nepal, with detection rates of 46.4% in 2013 and 36% in 2014, among AGE cases in children less than five years of age [64,65]. Unexpectedly, we did not detected the former dominant G2P[4] genotype. In Brazil, after Rotarix™ implementation in March 2006, this genotype has been the most frequently detected until 2015 [66], however, the recently low prevalence of G2P[4] viruses could be explained a cyclical pattern of circulation along with the herd induced homotypic immunity and depletion of the susceptible population [20].

A major strength of our study is that we included data from eleven states, representing around 100 million inhabitants (almost half of Brazilian population). Albeit, this could be considered as a major limitation as well, since the variability in reporting and collecting AGE cases by states generates surveillance biases. Another limitation is that important RVA genes, such as VP6 and NSP4, were not characterized. Nevertheless, future studies approaching a more complete genetic characterization of G3P[8] strains, as well as unusual genotypes detected here (G3P[6] and G12P[6]) will be performed, in order to monitor RVA genotypes spread and evolution over time.

In conclusion, we found a 12% of RVA-positivity in AGE cases from Brazil, and according to global trends, the equine-like G3P[8] was the dominant genotype in 2018 and 2019. The constant shifting of RVA genotypes circulation and the potential emergence of unusual/reassortant strains reinforces the importance and the need for continuous country-based epidemiological and molecular surveillance programs.

**Author Contributions:** Conceptualization, M.B.G., M.P.M., J.P.G.L., and T.M.F.; methodology, M.B.G., A.M.F., A.G.M., F.C.M., J.d.S.R.d.A., R.M.S.d.A., and S.d.S.e.M.; formal analysis, M.B.G. and T.M.F.; investigation, M.B.G., F.C.M., J.d.S.R.d.A., R.M.S.d.A., and T.M.F.; writing—original draft preparation, M.B.G. and T.M.F.; writing—review and editing, T.M.F., M.P.M., J.P.G.L., A.M.F., A.G.M., F.C.M., J.d.S.R.d.A., R.M.S.d.A., and S.d.S.e.M.; supervision, T.M.F.; project administration, T.M.F.; M.P.M., and J.P.G.L.; funding acquisition, T.M.F. and J.P.G.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by The Brazilian National Council for Scientific and Technological Development (CNPq) (Programa Universal grant no. 424376/2016-4) and (Programa Estratégico de Apoio à Pesquisa em Saúde—PAPES VII grant no. 401795/2015-2—FIOCRUZ/CNPq), Carlos Chagas Filho Foundation for Research Support of the State of Rio de Janeiro (FAPERJ) (grant no. 202.796/2019—Jovem Cientista do Nosso Estado, TMF; grant no. 202.779/2018—Cientista do Nosso Estado, JPGL) and PAEF-2—Oswaldo Cruz Institute. Further support was given by CGLab. We are grateful for support from the Coordination for the Improvement of Higher Education Personnel (CAPES—PrInt-Fiocruz-CAPES program).

**Acknowledgments:** We would like to thank the Rotavirus Surveillance System, coordinated by the Coordenação Geral de Laboratórios de Saúde Pública (CGLab), Brazilian Ministry of Health, and the State Central Laboratories involved in the study.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

### **Multiple Introductions and Predominance of Rotavirus Group A Genotype G3P[8] in Kilifi, Coastal Kenya, 4 Years after Nationwide Vaccine Introduction**

**Mike J. Mwanga <sup>1</sup> , Jennifer R. Verani 2,3, Richard Omore <sup>4</sup> , Jacqueline E. Tate <sup>3</sup> , Umesh D. Parashar <sup>3</sup> , Nickson Murunga <sup>1</sup> , Elijah Gicheru <sup>1</sup> , Robert F. Breiman <sup>5</sup> , D. James Nokes 1,6 and Charles N. Agoti 1,7,\***


Received: 20 October 2020; Accepted: 20 November 2020; Published: 24 November 2020 -

**Abstract:** Globally, rotavirus group A (RVA) remains a major cause of severe childhood diarrhea, despite the use of vaccines in more than 100 countries. RVA sequencing for local outbreaks facilitates investigation into strain composition, origins, spread, and vaccine failure. In 2018, we collected 248 stool samples from children aged less than 13 years admitted with diarrheal illness to Kilifi County Hospital, coastal Kenya. Antigen screening detected RVA in 55 samples (22.2%). Of these, VP7 (G) and VP4 (P) segments were successfully sequenced in 48 (87.3%) and phylogenetic analysis based on the VP7 sequences identified seven genetic clusters with six different GP combinations: G3P[8], G1P[8], G2P[4], G2P[8], G9P[8] and G12P[8]. The G3P[8] strains predominated the season (*n* = 37, 67.2%) and comprised three distinct G3 genetic clusters that fell within Lineage I and IX (the latter also known as equine-like G3 Lineage). Both the two G3 lineages have been recently detected in several countries. Our study is the first to document African children infected with G3 Lineage IX. These data highlight the global nature of RVA transmission and the importance of increasing global rotavirus vaccine coverage.

**Keywords:** gastroenteritis; rotavirus; G3[P8]; phylogenetics; equine-like

#### **1. Introduction**

Following progressive introduction of rotavirus vaccines into national immunization programs (NIP) of more than 100 countries since 2006, a significant decline of rotavirus group A (RVA) disease burden has occurred [1,2]. However, despite these successes, RVA remains a leading cause of diarrhea morbidity and mortality [3,4], resulting in an estimated 128,500 deaths annually among under-5-year-olds, a majority occurring in low-income settings [5]. Consistently, licensed oral RVA vaccines have underperformed in low-income settings compared with high-income settings [6,7]. After monovalent Rotarix® vaccine was introduced into Kenya's NIP in July 2014, with doses given at 6 and 10 weeks of life, a multi-site case-control study found an overall 2-dose vaccine effectiveness of only 64% (95% confidence interval (CI): 35–80%) in under-5-year-olds [8]. In England, the same vaccine showed effectiveness of 77% (95% CI: 66–85%) [9].

In humans, RVA immunity is partly conferred by neutralizing antibodies directed against the VP4 (protease-sensitive) and VP7 (glycoprotein) viral capsid surface proteins that define P and G types, respectively [10]. These two viral proteins are highly diverse, with up to 36 different G and 51 different P types recorded to-date [11], some of which predominantly infect non-human animal species [12]. Among other factors, the higher number of co-circulating GP genotypes in low-income settings has been proposed to be a potential contributor to rotavirus vaccine underperformance [6].

Currently, there are four licensed and WHO pre-qualified RVA vaccines; all live attenuated and administered orally, but with different strain compositions. These are monovalent Rotarix® (G1P[8]), pentavalent RotaTeq® (5 reassortant viruses; G1, G2, G3, G4 and G6 genotypes in combination with P[8]), monovalent ROTAVAC® (G9P[11]) and pentavalent ROTASIIL® (5 reassortant viruses; G1, G2, G3, G4 and G9). All four vaccines were shown to be largely cross-protective against heterotypic strains in both clinical trials and following vaccine implementation in several settings [6,13]. Paradoxically, post-vaccine rollout, outbreaks caused by strains heterotypic to the vaccine in use have been sometimes reported in countries, occurring in patterns seeming to be influenced by the vaccine regimen in use [14–16].

Recent genotyping studies of RVA have found increased proportions of G2P[4], G3P[8] and G12P[8] genotypes in rotavirus vaccinating countries [14,16–18]. These genotypes appeared to play only a minor role in the pre-vaccine era; thus, their increasing prevalence is consistent with increased capacity in escaping vaccine immunity [12,19]. Furthermore, there have been several reports of human infection with equine-like G3 viruses suggestive of greater human vulnerability to antigenically novel RVA strains [20–29]. At the Kenya Medical Research Institute (KEMRI)—Wellcome Trust Research Programme (KWTRP), we have maintained a RVA surveillance at Kilifi County Hospital (KCH), located in rural coastal Kenya since 2009 [30]. The aim of the current analysis was to determine the genetic relatedness of the strains that were in circulation in the 2018 RVA season in Kilifi, their origins, global phylogenetic context, and role in the local sub-optimal vaccine performance.

#### **2. Results**

#### *2.1. Study Population Characteristics*

Between January and December 2018, 384 children aged less than 13 years were admitted to KCH with diarrhea as one of their illness symptoms. Of these, 208 (54.2%) were Kilifi Health and Demographic surveillance system (KHDSS) area residents (Figure S1). A stool sample was obtained from 248 (64.6%). The main reasons for non-sampling were death (*n* = 13), discharge or transfer before sample collection (*n* = 22), consent refusal (*n* = 52), or other (*n* = 16). Among study eligible children (*n* = 384), the distribution of the sampled and not sampled children differed significantly across age strata (*p* = 0.002) and discharge outcome (*p* < 0.001), Table 1. The distribution of the sampled and not sampled children were similar across sexes and by rotavirus vaccine eligibility status. The majority of the eligible participants were aged less than 2 years (68.2%) and were age eligible to have received one or two doses of rotavirus vaccine (83.6%). By EIA testing, RVA was detected in 55 children (22.2%), Figure 1a, 32 (58.1%) of which were KHDSS area residents. Fifty-one (92.7%) of the RVA positive children were age eligible to have received two doses of the RVA vaccine. Of these, the vaccination status was known for 36 (70.6%), of which 29 (80.6%) were confirmed to have received two doses of Rotarix® vaccine while the remainder (19.4%) received one dose, Table 1.


**1.**AcomparisonofdemographiccharacteristicsofchildrenwithdiarrheaadmittedtoKilifiCountyHospital(KCH)thatsampledthose who were

**Value \***

0.008

0.073

 1.000

0.254

0.063

0.209

*Pathogens* **2020**, *9*, 981

¶ SD stands for standard deviation; δ IQR stands for interquartile range; **\$** *p* value for comparison of sampled and not sampled groups; **\*** *p* value for comparison of RVA positive and negative groups.

One or 2 dose eligible but status unknown

Outcome (*<sup>n</sup>* = 379)

Died

Alive

 38 (10.0) 341 (90.0)

 13 (5.3) 234 (94.7)

 25 (18.9) 133 (81.1)

 none

 15 (29.4)

> <0.001

 63 (40.1)

12 (6.3) 180 (93.8)

0.194

 126 (39.3)

 78 (37.5)

 48 (42.5)

**Figure 1.** Summary of rotavirus group A (RVA) surveillance in Kilifi County Hospital (KCH) in 2018 and identified genotypes. Panel (**a**) sample flowgram from patient recruitment to VP4 and VP7 genotyping results for the RVA positives. Panel (**b**) monthly cases of diarrhea in children aged less than 13 years recorded at KCH in 2018 (grey bars) compared with monthly proportions of RVA positive samples (black dashed line on the secondary axis). Panel (**c**) the number of RVA positive samples by month in 2018 and by the GP genotype. The black circle size is proportional to the number of samples (the smallest indicates one sample and the largest is 13 samples). Panel (**d**) genotypes identified in children according to rotavirus vaccination status.

#### *2.2. Characteristics of the RVA Infections and the Infected Children*

RVA prevalence was higher in female compared to male children admitted with diarrhea (29.8% vs. 15.7%, *p* = 0.008), Table 1. RVA was detected in all months of 2018 except January and February Figure 1b. Diarrhea cases peaked in June while RVA prevalence peaked in August (50% of all collected samples were RVA positive). Sequencing and GP typing was successful for 48 (87.3%) of the 55 RVA-positive samples. Five G types (G1, G2, G3, G9 and G12) and two P types (P[4] and P[8]) were identified in the successfully sequenced samples. From these, six GP combinations were identified, namely: G3P[8] (*n* = 37, 77.1%), G1P[8] (*n* = 6, 12.5%), G2P[4] (*n* = 2, 4.2%), G2P[8] (*n* = 1, 2.1%), G9P[8] (*n* = 1, 2.1%) and G12P[8] (*n* = 1, 2.1%). The G3P[8] and G1P[8] strains were the only genotypes detected for > 2 months while the other four genotypes were detected sporadically (1–2 months), Figure 1c. The distribution of the infecting genotype (summarized as G3P[8] versus non-G3P[8]) did not differ significantly by sex, patient age, vaccination status or discharge outcome, Table 2 and Figure 1d.


**2.**CharacteristicsofchildrenwhomwereinfectedwithrotavirusG3P[8]versusthosewhomwereinfectedwith

*Pathogens* **2020**, *9*, 981

#SD stands for standard deviation, δ IQR stands for interquartile range.

#### *2.3. Genetic Diversity in the Sequenced Viruses*

For the VP4 segment, a 579 nt long region (~25%) was recovered for 47 viruses (88.5%) while for the VP7 segment, a 644 nt long region (~65%) was recovered for 48 viruses (87.3%). One virus (KEN/KLF0879/2018), genotyped G9P[8], yielded a significantly shorter VP4 fragment relative to the other viruses (<500 nt) due to low quality sequencing data and was excluded from subsequent analyses. Consistent with the greater number of assigned G types (*n* = 5) compared to P types (*n* = 2) types, the range of pairwise nt differences was much greater in the VP7 (up to 203 nt differences) compared to VP4 segment (up to 87 nt differences), Figure 2a,b, respectively. A multi-modal distribution of nt differences was observed for both VP4 and VP7 segments. A total of 328 (~51%) and 141 (~24%) SNP positions were identified in the sequenced VP7 and VP4 fragments, respectively. Of the 48 sequenced samples, 22 (45.8%) yielded unique VP7 sequences while 17 (36.2%) gave unique VP4 sequences.

δ

**Figure 2.** Genetic diversity in the sequenced RVA positives from Kilifi County Hospital (KCH). Panel (**a**) shows the distribution of pairwise nt differences in the sequenced portion of VP7 (644 nt long) of 48 RVA positives. Panel (**b**) shows the distribution of pairwise nt differences in the sequenced portion of VP4 (579 nt long) of 47 RVA positives.

#### *2.4. Molecular Genetic Clusters*

Using the range of pairwise nt differences observed in first modal distribution for the VP7 (0 to 20 nt differences, i.e., >97% nt similarity) to define a molecular genetic cluster, seven G clusters were assigned (named Clu\_1-7). Members of a cluster were universally of same G type. All G type sequences identified to be of the same type formed a single cluster except G3P[8] that occurred in three clusters, named Clu\_3/G3P[8], Clu\_4/G3P[8] and Clu\_5/G3P[8]. The temporal pattern of the assigned clusters is shown in Figure 3a. Most of the high incidence months (April to August) had multiple genetic clusters co-circulating, except for July, which had a single G3P[8] cluster. The reconstructed phylogenetic relationship between strains of the different G and P types sequenced is shown in Figure 3b,c. The VP7 phylogeny showed segregation of the seven clusters we identified from the pairwise nt difference analysis. The VP4 phylogeny showed less clear-cut phylogenetic clustering with respect to the assigned genetic clusters. The two phylogenies were not entirely congruent, a feature suggestive of reassortment in the local strains. The minimum spanning networks reconstructed for both the VP7 and VP4 sequences are shown in Figure 3d,e. Viruses in the same genetic cluster consistently had four or less nt differences to the closest next virus within the same genetic cluster.

**Figure 3.** Temporal and genetic relatedness of the sequenced Kilifi rotaviruses. Panel (**a**) number of RVA positive samples by molecular genetic cluster and month. The circle sizes are proportional to the number of samples (the smallest indicates one sample and the largest is 13 samples). Panel (**b**) shows a Maximum Likelihood (ML) tree of the Kilifi 48 VP7 sequences. Panel (**c**) shows an ML tree of the Kilifi 47 VP4 sequences. Panel (**d**) shows the reconstructed POPART minimum spanning network from the 48 VP7 sequences. The vertexes represent the sequenced VP7 haplotypes. The size of the vertex is proportional to the number of haplotypes (identical sequences) and is colored by the assigned molecular genetic cluster. The numbers shown on the edges represent the number of nucleotide changes from one vertex (haplotype) to the next. Panel (**e**) same as panel (**d**) above but for the Kilifi 47 VP4 sequences.

#### *2.5. Spatial Distribution of the Kilifi G3 Genetic Clusters*

A few viruses in different VP7-based genetic clusters had identical VP4 sequences and we explored if these were spatially clustered. Twenty-eight of the 48 genotyped samples were from KHDSS area residents. The geographical distribution of all diarrhea admissions and the RVA positives by genetic cluster is shown in Figure S1. Cases of the predominant Clu\_3/G3P[8] strains came from only a few locations although it appeared that road access (especially the Malindi-Mombasa highway) may have played a role in influencing which patients were turning up at KCH due to easier access.

#### *2.6. Global Genetic Context of the Kilifi 2018 G3 Strains*

A total of 338 G3 sequences from 26 countries fully met the criteria for inclusion as comparison data, including 39 previously collected in Kenya. The phylogeny derived from the combined Kilifi and global G3 viruses is shown in Figure 4a while Figure 4b shows the phylogenetic relatedness of all previous G3 sequences of RVA sampled in Kenya (5 locations including Kilifi).

**Figure 4.** Global phylogeny derived from nucleotide sequences of G3 strains sampled between 2012–2018. (**a**) The phylogenetic tree reconstructed from 375 VP7 sequences of G3 type (338 collated from GenBank sampled across 26 countries including 39 from Kenya, and 37 G3 viruses sequenced in the current study) to determine the lineage and global context of the Kilifi sequences. The countries included were Australia, Belarus, Brazil, China, Dominican Republic, Ethiopia, Hungary, India, Indonesia, Italy, Japan, Kenya, South Korea, Kuwait, Nigeria, Pakistan, Peru, Russia, Spain, Taiwan, Thailand, USA, Uganda and Vietnam. The taxa for Kenya G3 sequences are provided by filled circles colored green and with the assigned Kilifi clusters names indicated next to the branches containing these sequences. Panel (**b**) a phylogeny of all Kenya G3 sequences (*n* = 76). The different colors of the filled circle symbols indicate the Kenya taxa distinguished by their location of sampling. The names assigned to the Kilifi clusters are indicated next to the nodes leading to their branches as similarly shown in panel (**a**).

A majority of the global viruses fell within two of nine previously identified G3 lineages [25]; Lineage I and equine-like G3 lineage (named Lineage IX). The Kilifi G3 sequences had representation in both these two lineages: Lineage I (*n* = 35, 94.6%) and equine-like G3 Lineage (*n* = 2, 5.6%). Viruses of the genetic cluster Clu\_4/G3P[8] clustered with the equine-like G3 Lineage while the Kilifi G3 Lineage I viruses separated into two groups that corresponded to the Clu\_3/G3P[8] cluster (*n* = 30) and the Clu\_5/G3P[8] cluster (*n* = 5).The distribution of the pairwise nt differences in the compiled global G3 sequences dataset, like for the Kilifi G3 viruses, showed a multi-modal distribution (figure not shown). The first major trough was observed at 27 nt differences.

On applying the threshold used to identify the local molecular genetic clusters (>97% genetic similarity) on the global G3 dataset, 18 clusters were identified (Table S1). Of these, eight were singletons, six comprised of between 2 and 3 members and the remaining four clusters had 10, 47, 116 and 181 members. All the Kilifi G3 viruses fell in the three clusters that had the highest membership overall, Table S1. For each of the three Kilifi G3 genetic clusters we explored their closest genetic relative in the global dataset by network reconstructions (Figure 5). For the Kilifi Clu\_3/G3P[8] the closest similar sequences were from India (G3P[8] collected in 2016) and Singapore (G3P[8] collected in 2016) that had 2 nucleotide differences Figure 5a. For the Kilifi Clu\_4/G3P[8] (the equine-like G3 Lineage) the closest relative was from Taiwan (G3P[8] collected in 2016) with zero nucleotide difference in the sequenced region Figure 5b. For the Kilifi Clu\_5/G3P[8] the closest relatives were from Kenya (G3P[6] collected in 2014) and Uganda (G3P[6] collected in 2013) that had zero and 2 nucleotide difference, respectively, Figure 5c. Overall, within these three major global G3 genetic clusters, clustering by country was common.

**Figure 5.** Haplotype network showing relationships of the identified global G3 lineages that included Kilifi viruses. Panel (**a**) shows the network for Lineage I cluster viruses that included the Kilifi Clu\_3/G3P[8] strains. The vertices represent the VP7 haplotypes. The size of the vertex is proportional to the number of haplotypes (identical sequences) and is colored by the country of sampling. The numbers shown on the edges represent the number of nucleotide changes from one vertex (haplotype) to the next. Panel (**b**) and (**c**) have the same description as panel (**a**) above but represent Lineage IX (equine-like G3) cluster that included Kilifi Clu\_4 G3P[8] and the Lineage I cluster that included Kilifi Clu\_5 G3P[8] sequences, respectively.

#### **3. Discussion**

Four years after Kenya introduced Rotarix® vaccine into its NIP, multiple RVA GP genotypes circulated during the 2018 season in Kilifi, Kenya, with the G3P[8] genotype predominating at 67.2%. At this study site, the preceding two years (2016 and 2017) were dominated by the G2P[4] and G1P[8] genotypes, respectively, with only six cases of G3P[8] detected from September 2009 to December

2017 [30] and an additional three partially genotyped G3P[x] detected in 2013 [31]. The G3P[8] strains are partially heterotypic to the monovalent Rotarix® vaccine, which is comprised of an attenuated G1P[8] strain. During 2018, this local G3P[8] predominance is consistent with the previously documented season-to-season spatial-temporal fluctuations in the prevalence of RVA genotypes [12], hypothesized to be driven by the prevailing population-level immunity derived from natural infections and the use of vaccines [14].

Vaccination records were available for 70.6% of the children with an RVA positive test. Of these, 92.7% were age eligible to have received the two doses of Rotarix® vaccine and, in that subgroup, the vast majority (80.6%) had indeed received the full 2-dose series. However, overall, the vaccination status of these children did not appear to predict either their RVA diagnosis result or the infecting GP genotype. These findings, albeit from a single season and site, suggest that for these children who acquired an RVA infection despite one or two-dose vaccination, host factors rather than viral characteristics or vaccine composition may explain the vaccine failures. A follow-up study is planned.

At least seven distinct genetic clusters constituted the 2018 coastal Kenya RVA season. The VP7 sequences showed greater genetic diversity and provided a better phylogenetic resolution compared to the VP4 sequences. Each of the identified G types corresponded to a single genetic cluster except G3 viruses that segregated into three genetically distinct clusters. Strikingly, some samples with different G types yielded identical VP4 sequences, indicating that some of the children may have been infected by reassortant viruses or harbored mixed infections [25]. Our analyses improve understanding on the recent composition and transmission patterns of local RVA seasons, providing insight into the design of final stretch RVA control strategies following vaccine introduction.

Several recent studies have reported the increased proportion of G3P[8] strains, e.g., in Australia [14], Japan [32], Thailand [28], Indonesia [29], Pakistan [33], Dominican Republic [25], Brazil [34], Spain [20], Mozambique [24], Malawi [35] and Botswana [36]. The global G3 sequences available from GenBank showed extensive genetic diversity. The significance of this diversity in relation to human immune recognition should be investigated. Notably, recent years have also observed the emergence and global spread of a new G3 lineage named equine-like G3, of putative equine origin, assigned G3 Lineage IX [25]. Strains of G3 Lineage IX were first detected in 2013 in Japan and have since been widely detected in several other countries (Australia [21], Taiwan (unpublished data in GenBank), Indonesia [29], Thailand [28], USA [26], Dominican Republic [25], Brazil [34], Italy [23], Germany [27], Hungary [22] and Spain [20]). Our study is the first to document African children infection with the G3 Lineage IX. Continued surveillance to monitor whether this particular strain becomes endemic in Kenya and the wider Africa continent in the face of increased RVA vaccine coverage is important to optimize RVA vaccine-mediated control. Notably, recent studies in Botswana [36], Mozambique [24], Malawi [35] and Ethiopia [37] reported increased prevalence of G3 type viruses but sequencing data from these studies are not yet available.

Based on sequence data deposited in GenBank, the predominant Kilifi G3 cluster (Clu\_3/G3P[8]) was the second most common genetic cluster globally. The closest sequences were from Singapore and India, both countries that did not yet have RVA vaccine in their NIP in 2018. The second most prevalent Kilifi G3 genetic cluster was Clu\_5/G3P[8]. Notably, this cluster has not been detected frequently around the globe and the closest genetic links were Kenyan strains collected in Kiambu County (Central province) in July and August 2014 [38], Kilifi in 2017, and strains from Ethiopia (collection date: April 2016 [39]) and Uganda (collection date: January 2013 [40]), neighboring countries which included RVA vaccines in their NIP in 2013 and 2018, respectively. Although the Kilifi Clu\_4/G3P[8] (equine-like G3 Lineage) was the least prevalent locally, it was the most prevalent globally. The closest relatives to the Kenyan strains were from Taiwan, a country yet to introduce RVA vaccination.

This study had some limitations. First, the sequence data from the cohort represents a single site and one season. Second, we only sequenced portions of the VP4 and VP7 segments. Whereas these data were adequate to assign genotypes, lineages and estimate the number of genetic clusters, whole genome sequences provide a better resolution in examining reassortment events, evolution in

internal genes and studying genetic clusters [18,25,41]. Third, to determine the origin and pathways of spread of the imported genetic clusters, background sequence data from more countries and including populations neighboring coastal Kenya would have been ideal. Unfortunately, sequence data in public sequence databases to facilitate such phylogeographic analysis are currently limited. Fourth, the absence of significant epidemiological data for some variables e.g., vaccine status for ~30% of the RVA positive children and geographic origin for children from outside the KHDSS area limited our analyses.

In conclusion, the finding that >20% of diarrheal stools from children admitted to KCH with diarrhea in 2018 were RVA positive highlights that RVA is still a significant contributor to severe childhood diarrhea in coastal Kenya, despite the introduction of Rotarix® into Kenya's NIP in 2014. The cross-continent detection of the emerging equine-like G3 viruses and other typical human G3 strains demonstrates the global nature of RVA transmission. Strikingly, strains found circulating in the Kilifi population were most closely related to strains circulating in countries that were yet to introduce RVA vaccines into their NIP. This observation reminds of the global connectedness regarding pathogen movement and emphasizes the importance of vaccinating all eligible populations across the world, as failure to do so builds a reservoir for strains that continue to seed transmission in vaccinated populations. Identifying factors responsible for RVA vaccine underperformance in low-income settings is a priority research area that may support efforts to further reduce RVA burden. Our study did not ascertain that viral genetic diversity is a contributor to the vaccine underperformance in this setting. Studies investigating the relationship between RVA vaccine immunogenicity and infant characteristics, such as malnutrition, age at first RVA dose, concomitant receipt of oral polio vaccine (OPV), enteric co-infections and enteric dysbiosis may provide better insight into RVA vaccine performance characteristics.

#### **4. Materials and Methods**

#### *4.1. Study Population and Location*

KCH is the main referral hospital in Kilifi County (population size ~1.5 million people). The major economic activities in the county are subsistence farming, fishing and tourism [42]. An area around KCH (~900 km<sup>2</sup> with a population of ~300,000 people) is monitored by the KWTRP and is known as the KHDSS area [42], Figure S1. A high proportion of the patients seeking care at the KCH are KHDSS area residents [42]. Vaccination data of admitted children were collected using an electronic registry [8,43,44].

In the current analysis, stool samples were collected from eligible and consented pediatric patients admitted to KCH between January and December 2018 (the surveillance period), as part of the ongoing rotavirus surveillance program [8,31,43]. All children aged <13 years old admitted with diarrhea (defined as passing three or more watery stools in the last 24-h) were eligible for inclusion [8,31,43]. Following a review of demographic and clinical data collected by a clinical staff, parents or caregivers of eligible children were approached for consent, and a single stool sample was collected. The samples were immediately transferred into a cool box with ice blocks before transportation to the KWTRP for RVA testing and long-term storage at −80 ◦C.

#### *4.2. Specimen Laboratory Processing*

RVA in the stool samples was detected using ProSpecT™ enzyme immunoassay (EIA) kit (Oxoid, Basingstoke, UK) following the manufacturer's instructions. RVA positive samples were amplified in the VP4 and VP7 segments using One-step Reverse Transcriptase PCR Kit (Qiagen, Valencia, CA, USA) using previously published primers [45,46]. Successful amplification of the target regions was confirmed by the presence of expected bands (VP4: 660 bp and VP7: 881 bp) following gel electrophoresis of the PCR products. Products from successful PCRs were purified using GFX DNA purification kit (GFX-Amersham, Amersham, UK) and sequenced bi-directionally (both in forward and

reverse directions) using Big Dye Terminator 3.1 (Applied Biosystems, Foster City, CA, USA) chemistry. The primers used during PCR amplification were used for sequencing on an ABI Prism 3130xl Genetic Analyzer (Applied Biosystems, Foster City, CA, USA).

#### *4.3. Genotyping and Phylogenetic Analysis*

The sequence reads were assembled using Sequencher v5.4.6 (Gene Codes Corp Inc., Ann Arbor, MI, USA). Nucleotide (nt) sequence alignments were prepared using MAFFT v7.222 and visualized using Aliview v1.8. G and P genotypes were determined using Virus Pathogen Resource (ViPR) online classification tool [47]. The best nt substitution model for the alignments were determined IQ-Tree v1.6.6 [48]. Phylogenetic trees were reconstructed using the maximum likelihood (ML) method in RaxML v8.2.12 [49] and MEGA v7 [50]. Support for the tree branching patterns was evaluated by 1000 bootstrap iterations.

#### *4.4. Genetic Clusters*

Molecular genetic clusters were defined from the distribution of pairwise nt differences of VP7 segment sequences. Pairwise nt differences were determined using pairsnp (https://github.com/ gtonkinhill/pairsnp/). Viruses within the same molecular genetic clusters were those which pairwise nt differences occurred within the first modal distribution. Using this threshold, clusters were identified using the USEARCH algorithm [51]. Single nucleotide polymorphic (SNP) positions in alignments were assessed using parseSNP [52]. The minimum spanning networks between the RVA positive patients were reconstructed using POPART v1.70 program [53].

#### *4.5. Comparison Dataset*

The phylogenetic context of the locally predominant genotype in global RVA populations was investigated by co-analysis with similar G type strains sequence data deposited in GenBank. The search in GenBank was conducted in October 2020. The criteria for comparison data inclusion were (i) detection in a human stool/rectal swab specimen, (ii) sequence fully overlapping with the VP7 region sequenced for the Kilifi viruses, (iii) information on country and date of sampling available and (iv) sample collected in 2012–2018. G3 sequences collected previously from around Kenya including Kilifi were included in the analysis.

#### *4.6. Statistical Analysis*

Numerical data were analyzed in STATA v15.1. Continuous variables were summarized using various measures of dispersion. Differences between groups were assessed using a t-test or Wilcoxon rank-sum test. Binary data were summarized using proportions and comparison between groups made using either χ <sup>2</sup> or Fisher's exact test (depending on group sample size). The 95% CI were presented for proportions and standard deviation for means. A *p*-value of <0.05 was considered significant.

#### *4.7. Data Availability*

Partial sequences for the VP7 and VP4 segments reported in this work have been deposited to GenBank database under the sequence accession numbers MN194408-MN194485 for VP7 and MN194325-MN194364 for VP4.

#### *4.8. Ethical Statement*

Before sample collection informed written consent was obtained from the child's parent or guardian. The Scientific Ethics Review Unit (SERU) board that sits at KEMRI, Nairobi, approved the study protocols (SERU#3049).

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-0817/9/12/981/s1, Figure S1: Geographic origin distribution of sampled children who presented with diarrhea symptoms at KCH

and were Kilifi Health Demographic Surveillance System (KHDSS) area residents; Table S1: The global distribution of the identified G3 global genetic clusters.

**Author Contributions:** Conceptualization, C.N.A., D.J.N., M.J.M. and R.F.B.; methodology, D.J.N., C.N.A. and M.J.M.; software, M.J.M. and C.N.A.; validation, C.N.A., E.G. and M.J.M.; formal analysis, C.N.A. and M.J.M.; investigation, C.N.A., M.J.M. and D.J.N. and E.G., resources, R.F.B., J.E.T., U.D.P. and D.J.N.; data curation, N.M., M.J.M. and C.N.A.; writing—original draft preparation, M.J.M. and C.N.A.; writing—review and editing, R.O., J.R.V., R.F.B., J.E.T., U.D.P. and D.J.N.; visualization, C.N.A. and M.J.M.; supervision, D.J.N. and C.N.A.; project administration, E.G., R.O., M.J.M., N.M. and D.J.N.; funding acquisition, D.J.N., R.O., J.R.V., J.E.T., U.D.P., R.F.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by Gavi, The Vaccine Alliance through Emory University, to the Rotavirus Immunization Program and Evaluation in Kenya (RIPEK) and the Wellcome Trust (grant numbers 203077 and 102975). CNA is supported through the DELTAS Africa Initiative [DEL-15-003]. The DELTAS Africa Initiative is an independent funding scheme of the African Academy of Sciences (AAS)'s Alliance for Accelerating Excellence in Science in Africa (AESA) and supported by the New Partnership for Africa's Development Planning and Coordinating Agency (NEPAD Agency) with funding from the Wellcome Trust [107769/Z/10/Z] and the UK government. Views expressed in this publication are those of the authors and not necessarily those of Gavi, AAS, NEPAD Agency, Wellcome Trust or the UK government.

**Acknowledgments:** We thank the participants who provided samples for analysis, and the Clinical and Laboratory staff at Virus Epidemiology and Control group, KEMRI Wellcome Trust Research Programme for collection and processing the samples. This work is published with permission from Director KEMRI.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

**Disclaimer:** The findings and conclusions in this report are those of the authors and do not necessarily represent the official position of the U.S. Centers for Disease Control and Prevention (CDC) or the authors' affiliated institutions.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*

## **Epidemiological Trends of Five Common Diarrhea-Associated Enteric Viruses Pre- and Post-Rotavirus Vaccine Introduction in Coastal Kenya**

**Arnold W. Lambisia 1,2,\* ,**† **, Sylvia Onchaga 1,**† **, Nickson Murunga <sup>1</sup> , Clement S. Lewa <sup>1</sup> , Steven Ger Nyanjom <sup>2</sup> and Charles N. Agoti 1,3**


Received: 8 July 2020; Accepted: 10 August 2020; Published: 15 August 2020

**Abstract:** Using real-time RT-PCR, we screened stool samples from children aged <5 years presenting with diarrhea and admitted to Kilifi County Hospital, coastal Kenya, pre- (2003 and 2013) and post-rotavirus vaccine introduction (2016 and 2019) for five viruses, namely rotavirus group A (RVA), norovirus GII, adenovirus, astrovirus and sapovirus. Of the 984 samples analyzed, at least one virus was detected in 401 (40.8%) patients. Post rotavirus vaccine introduction, the prevalence of RVA decreased (23.3% vs. 13.8%, *p* < 0.001) while that of norovirus GII increased (6.6% vs. 10.9%, *p* = 0.023). The prevalence of adenovirus, astrovirus and sapovirus remained statistically unchanged between the two periods: 9.9% vs. 14.2%, 2.4% vs. 3.2 %, 4.6% vs. 2.6%, (*p* = 0.053, 0.585 and 0.133), respectively. The median age of diarrhea cases was higher post vaccine introduction (12.5 months, interquartile range (IQR): 7.9–21 vs. 11.2 months pre-introduction, IQR: 6.8–16.5, *p* < 0.001). In this setting, RVA and adenovirus cases peaked in the dry months while norovirus GII and sapovirus peaked in the rainy season. Astrovirus did not display clear seasonality. In conclusion, following rotavirus vaccine introduction, we found a significant reduction in the prevalence of RVA in coastal Kenya but an increase in norovirus GII prevalence in hospitalized children.

**Keywords:** viral diarrhea; real-time PCR; rotavirus vaccination; Kenya

#### **1. Introduction**

In the year 2016 alone, approximately 300,000 children aged <5 years succumbed to diarrhea in sub-Saharan Africa [1]. Viral pathogens including rotavirus group A (RVA), adenovirus (type 40/41), astrovirus, norovirus (genogroup GI and GII) and sapovirus are among the top causative agents of severe diarrhea globally [2,3]. Understanding their epidemiological patterns such as prevalence, incidence, seasonality, clinical severity and infection age distribution in local settings is essential for designing and prioritizing interventions. Historically, RVA has been the single most important cause of severe childhood diarrhea, responsible for ~38% (95% CI: 4.8–73.4%) of hospital cases (<5 years) pre-vaccine introduction [4]. However, RVA prevalence has been rapidly declining since 2009 and was approximately 23% (95% CI: 0.7–57.7%) in 2016, in settings where the rotavirus vaccine was in use [4]. Due to the shared ecological niche and the apparent decline of all-cause gastroenteritis-associated

hospital admissions, it has been hypothesized that rotavirus vaccination has likely impacted the epidemiology of the other enteric viruses [5]. However, there are contradicting reports on the specific impact of rotavirus vaccination on the prevalence of the individual enteric viruses—for example, norovirus [6,7]. This has not been adequately examined in African populations where diarrhea burden is highest. Kenya began rotavirus vaccination in July 2014 using the monovalent Rotarix® (RV1), derived from G1P[8] strain, administered at 6 and 10 weeks of life. RV1 vaccine coverage in Kenya has increased over time since 2014 but is varied by age group, number of doses and geographic region in Kenya [8]. Within Kilifi County, coastal Kenya, coverage in 2017 in <1-year-olds was 73% (at least one dose) vs. 65% (complete two doses), while in <12–24 month-olds, it was 86% (at least one dose) vs. 84% (complete two doses) [9].

The KEMRI/Wellcome Trust Research Programme (KWTRP) has been running surveillance of RVA since 2002 in children admitted to the Kilifi County Hospital (KCH). The current study screened archived diarrheal samples from KCH, spanning both the pre- and post-rotavirus vaccine introduction periods in Kenya for RVA, astrovirus, adenovirus (all serotypes), sapovirus and norovirus (only GII) using real-time reverse-transcription polymerase chain reaction (RT-PCR) approach. We update on the prevalence of these viral diarrheal agents and their seasonal patterns pre and post introduction of the rotavirus vaccination program in Kenya.

#### **2. Results**

#### *2.1. Study Population Characteristics*

Out of 2156 children aged <5 years who presented with diarrhea at KCH during the four selected years (2003, 2013, 2016 and 2019), 1397 (64.8%) provided a stool sample; see Table 1. Overall, the demographic characteristics of the eligible children sampled, and eligible children not sampled, differed in age strata distribution (*p* = 0.001) and discharge outcome (*p* < 0.001); see Table 1. The main reasons for failure to sample eligible children were as follows: death (*n* = 21, 2.8%), discharge or transfer before sample collection (*n* = 296, 40.0%), consent refusal (*n* = 315, 41.5%) or other (*n* = 127, 16.7%). Among the sampled cases, 984 (70.4%) had a specimen available and tested by real-time RT-PCR for the five enteric viruses, and these were included in subsequent analysis. The median age of the sampled participants was significantly higher for the post-vaccine introduction period compared to pre-vaccine introduction period (*p* < 0.001); see Table 2.




**Table 1.** *Cont.*

SD means standard deviation; IQR means interquartile range. Not sampled: sample was not collected due to lack of consent, time-up, death and others. # Discharge outcome data for three subjects were missing.



SD means standard deviation; IQR means interquartile range; RVA means rotavirus group A. Values given are the counts and percentages are provided in brackets. # Discharge outcome for two subjects was missing. Disease Severity Was Calculated Using the Vesikari Clinical Severity Scoring System Manual [10].

#### *2.2. Overall Virus Detection*

Of the 984 samples analyzed, at least one of the viruses was detected in 401 samples (40.8%) at the real-time RT-PCR cycle threshold (Ct) value of <35.0. The lower the Ct value, the higher the virus titer in the sample. The detection frequency differed significantly for adenovirus (*p* = 0.001) and sapovirus (*p* < 0.001) pre- and post-rotavirus vaccine introduction when the Ct cut-off value was gradually lowered (<30, <35, <40), unlike for RVA, astrovirus and norovirus GII; see Figure 1. All our subsequent analyses were undertaken at Ct value <35.0 Single infections were detected in 354 specimens (36.0%) and included RVA (*n* = 149, 42.1%), adenovirus (*n* = 91, 25.7%), norovirus GII (*n* = 75, 21.2%), sapovirus (*n* = 20, 5.7%) and astrovirus (*n* = 18, 5.1%).

**Figure 1.** Detection frequency of RVA, adenovirus, norovirus GII, astrovirus and sapovirus at different cycle threshold (Ct) cutoffs for children under 5 years of age admitted to KCH Kenya with diarrhea symptoms. The error bars represent 95% confidence interval for the proportions. Proportions were compared using chi-square test. RVA stands for rotavirus group A, ADV stands for adenovirus, NOR stands for norovirus GII, ASV stands for astrovirus and SAP stands for sapovirus.

#### *2.3. Patterns Pre-Post Vaccine Introduction*

RVA showed a significant decrease (23.3% vs. 13.8%, *p* < 0.001) in prevalence while norovirus GII showed a significant increase (6.6% vs. 10.9%, *p* = 0.02) post-vaccine introduction compared to pre-vaccine introduction; see Table 3. There were no significant changes in the prevalence of astrovirus (*p* = 0.585), adenovirus (*p* = 0.053) and sapovirus (*p* = 0.133) pre- and post-RVA vaccine introduction (chi-squared (χ 2 ) test); see Table 3. Notably, norovirus GII had a gradual increase in prevalence across the four years, from 6.7% (95% CI: 3.8–10.9%) to 12.4% (95% CI: 8.8–16.7%); see Figure 2. RVA was the most commonly detected virus across all years, except in year 2019, in which adenovirus had the highest prevalence; see Figure 2. χ – –

**Table 3.** Comparison of the prevalence of viral detection in children under 5 years of age admitted to KCH Kenya with diarrhea symptoms pre- and post-rotavirus vaccine introduction.


**Figure 2.** Prevalence of RVA, adenovirus, norovirus GII, astrovirus and sapovirus in 2003, 2013, 2016 and 2019 in children under 5 years of age admitted to KCH Kenya with diarrhea symptoms. The error bars represent 95% confidence interval for the proportions. Proportions were compared using chi-square test. Abbreviations used for viruses as in Figure 1.

Notably, RVA and sapovirus cases in the post-vaccine introduction period had statistically significant lower and higher median Ct values, respectively, compared to the pre-vaccine period (Wilcoxon, *p* value < 0.001); see Figure 3. This was not observed for the other three screened viruses preand post-rotavirus vaccine introduction. The median age of the RVA positive cases was significantly higher for the post-vaccine introduction period (14.0 months) compared to the pre-vaccine introduction period (10.4 months) (Wilcoxon, *p* < 0.001). A similar shift was not observed for the other viruses; see Figure 4.

**Figure 3.** Distribution of Ct values among cases under 5 years of age admitted to KCH Kenya with diarrhea symptoms pre- and post-vaccine introduction. RVA stands for rotavirus group A, ADV stands for adenovirus, NOR GII stands for norovirus GII, ASV stands for astrovirus and SAP stands for sapovirus.

**Figure 4.** Distribution of age in months among cases under 5 years of age admitted to KCH Kenya with diarrhea symptoms pre- and post-vaccine introduction.

#### *2.4. Virus Coinfections (i.e., Two or More Viruses in a Single Specimen)*

These were detected in 47 specimens (4.8%). In 583 specimens (59.2%), none of the targeted viruses was detected. The prevalence of coinfections pre-vaccine was 4.4% (95% CI: 2.7–6.7%), while in the post-vaccine introduction period, this value was 5.7% (95% CI: 3.9–8.0%), *p* = 0.454. RVA and astrovirus were the most common coinfections in the pre-vaccine introduction period (*n* = 6), while in the post-vaccine introduction period, it was RVA and adenovirus (*n* = 15); see Table 4.

**Table 4.** Coinfections pre- and post-rotavirus vaccine introduction. RVA stands for rotavirus group A, ADV stands for adenovirus, NOR GII stands for norovirus GII, ASV stands for astrovirus and SAP stands for sapovirus.


Abbreviations used for viruses as in Figure 3.

#### *2.5. Circulating RVA Genotypes Pre- and Post-Vaccine Introduction*

G1P[8] was the predominant RVA genotype pre vaccine introduction. However, in the post-vaccine introduction period, the predominant genotypes were G2P[4] (2016) and G3P[8] (2019); see Table 5.


**Table 5.** Frequency of RVA genotypes detected in coastal Kenya pre- (2003 and 2013) and post- (2016 and 2019) vaccine introduction.

#### *2.6. Seasonality of the Detected Viruses*

We constrained this analysis to the years 2013, 2016 and 2019, where >70% of the eligible patients had been analyzed. Pre-vaccine introduction (in 2013), for RVA, there were two peak months, in June and September. However, post-vaccine introduction (in 2016 and 2019), there was only a single peak month for RVA in September and August, respectively. For norovirus GII, cases were observed throughout the year, with peak months varying from year-to-year, in July, April and June in 2013, 2016 and 2019, respectively. Similarly, adenovirus cases appeared to occur throughout the year, with two peak months in 2013 (June and September) and one peak month in 2016 and 2019 (August for both). For sapovirus and astrovirus, we observed less than five cases monthly between January and August and no cases in the last quarter of each the three years; see Figure 5.

**Figure 5.** The frequency of detection of RVA, adenovirus, norovirus GII, astrovirus and sapovirus by month in children under 5 years of age admitted to KCH Kenya with diarrhea in 2013, 2016 and 2019.

#### *2.7. Primer*/*Probe Mismatches with Contemporary Sequences*

3′ 3′ — Nucleotide mismatches were observed in either or both the primers and probes and the viral target sequences for all the viruses except for norovirus GII; see Figure 6. The RVA forward primer had a G-A and A-G mismatches at positions 12 and 15, respectively. Adenovirus had two mismatches in the forward primer (C-G and G-A), three mismatches in the probe (C-T, C-T and T-C) and two mismatches in the reverse primer (T-C and C-T), and none of them were within five bases of the 3 ′ end. Mismatches within the sapovirus primer/probe binding sites were pronounced in sapovirus genogroup V and included six mismatches in the forward primer, three mismatches in the probe and two mismatches in the reverse primer. Some of the mismatches were within five bases of the 3 ′ end (forward primer: C-G, probe: T-C, reverse primer: A-C and T-C). Astrovirus primers and probe did not have pronounced mismatches present in all the sequences—rather, they had mismatches in individual sequences; see Figure 6.

**Figure 6.** *Cont*.


**Figure 6.** The primers and probes target sites for RVA, adenovirus and norovirus GII, sapovirus and astrovirus were aligned using MAFFT v.7.31313 and the alignments were trimmed to the region of the primer and probe target sites. Nucleotide differences between the expected primer and probe target sites and the viral sequences were identified and highlighted. Dots indicate identity with primer or probe sequences.

#### **3. Discussion**

We observed a significant decrease in the prevalence of RVA in the post-vaccine introduction period in KCH, concurring with findings of a recent multi-site study in Kenya that reported RVA vaccine effectiveness of ~64% (95% CI: 35–80%) and a reduction in rotavirus-associated hospital admissions two years post-vaccine introduction of ~80% (95% CI: 46–93%) [9,11]. Note that Kenya rotavirus vaccine coverage was considered medium in 2018 (70–79%) [12]. Our pre- and post-vaccine introduction analysis observed a significant increase in the prevalence of norovirus GII in KCH post-rotavirus vaccine introduction, as similarly observed in the United States, Nicaragua and Bolivia following RVA vaccine introduction [13–15]. It is unclear if this has been driven by an established biological interaction between these two viruses or that this reflects natural norovirus GII fluctuation in prevalence across multiple years.

The shift in the predominant genotypes pre- and post-vaccine introduction from G1P[8] to G2P[4] in 2016 and G3P[8] in 2019 in our setting has also been described elsewhere, e.g., in Belgium, Madagascar and Ethiopia [16–18]. G3P[8] was the predominant genotype in this setting in 2018 [19] and it continued being the dominant genotype in 2019. Although these dominant post-vaccine genotypes are either partially or fully heterotypic to the Rotarix G1P[8] strain, in their surface exposed immunodominant proteins, there is not enough evidence yet to directly attribute their increased incidence to vaccine introduction [20]. Additional analysis will help to bring better understanding on the reason behind their dominance.

Despite RV vaccine introduction in Kilifi, Kenya, no significant difference was observed in the discharge outcome for all causes of diarrhea pre- and post-rotavirus vaccine introduction. We suggest two explanations for this. Firstly, the majority of the children who were eligible to be in this study and died did not have a sample collected to determine their RVA and other enteric pathogens' status. Secondly, inpatient mortality of children treated for diarrhea in Kilifi County Hospital has been previously found to be predicted by a positive HIV test, bacteremia and poor nutritional status [21]. This may have not changed pre- or post-introduction of rotavirus vaccination.

RVA Ct values were decreased in post-vaccine samples compared to pre-vaccination years. This was despite RVA disease severity remaining unchanged between the two periods. Different extraction methods were used to process the samples between 2003, 2013 and 2016, 2019. However, according to Liu et al., the difference in the extraction methods for enteric pathogen studies is not significant, except for norovirus GII, which showed a higher Ct value with kits targeting RNA purification alone compared to those targeting total nucleic acid (TNA) (difference within 1 Ct value). Different extraction kits were used in this study because raw stool samples from 2003 to 2016 were already destroyed following a directive by the WHO in 2016 that was part of the larger global polio eradication effort.

It has been previously noted the introduction to rotavirus vaccines may result in the shift of diarrhea disease burden to slightly older age groups [20]. Our study found a significant increase in the median age of diarrhea cases post-vaccine introduction (12.5 months) compared 11.2 months pre-introduction. This in part may be explained by the higher immunity at both individual and population levels against rotavirus that wanes as children grow older.

On local seasonality patterns, in each year, a peak month(s) of occurrence was observed for RVA, norovirus GII, sapovirus and adenovirus but not astrovirus. The Kilifi area has a tropical climate with two rainy seasons; the main rains usually peak in May (up to July) while the short rains usually peak in November (can run from October to December). RVA and adenovirus appeared to peak in the dry months while norovirus GII and sapovirus peaked in the rainy season. Similar patterns in the seasonality of RVA, adenovirus, norovirus GII and sapovirus have been observed elsewhere [22–25]. The seasonality of astrovirus is not well described.

The performance of qPCR assays can be impacted by mismatches within the last five bases at the 3 ′ end of primers and probe or/and the number of mismatches being more than five in the primers and probe [26,27]. The mismatches observed in the primer and probe binding sites of adenovirus, astrovirus and sapovirus may have impaired the real-time PCR function by blocking the amplification or increasing the quantification cycles. Consequently, this may have impacted the estimated frequency of detection of these viruses. Unlike for RVA, the magnitude of the mismatches in qPCR function could have been shown better using recent local sequences of the other viruses.

This study had limitations: firstly, we did not analyze healthy children in the community to inform on the background prevalence of the five viruses in our study population. Secondly, the adenovirus assay was not specific to type 40/41 alone; thus, some of the adenoviruses detected may not be associated with diarrhea. Thirdly, a significant number of eligible cases were not sampled, including those who died before sampling. This potentially biased prevalence of the screened pathogens in the study population. Fourthly, extracting TNA from samples after many years of storage could lead to lower Ct values due to deterioration. Finally, the seasonality of examined pathogens will be best described if we examine more years.

In conclusion, we found a significant decline in the prevalence of rotavirus in hospitalized children in coastal Kenya after rotavirus vaccine introduction. This finding reinforces evidence of the continued benefit of rotavirus vaccination in this setting. Concomitantly, there has been a surge in norovirus GII prevalence, but the factors driving this increase are unclear and will require future investigation. The observation that the screened viruses peak at different times of the year also would benefit further investigation in order to understand drivers of their transmission and inform the design of effective intervention measures.

#### **4. Materials and Methods**

#### *4.1. Study Site and Population*

This study was undertaken at KCH, a referral hospital serving the Kilifi County population, which is majorly a rural population. We utilized stool specimens collected during routine surveillance of rotavirus in children with diarrhea as one of their illness symptoms, aged below five years and admitted to KCH [9,11]. Diarrhea was defined as observation of three or more loose stools in the preceding 24-h period. In this study, we selected two pre-vaccine years (2003 and 2013) and two post-vaccine years (2016 and 2019) for analysis. A stool specimen was collected from children who met the diarrhea case-definition following parental or guardian consent. The study protocol was approved by the Scientific and Ethics Review Unit (SSC#2861 and SERU#CGMRC/113/3624) based at KEMRI, Nairobi, Kenya.

#### *4.2. Laboratory Methods*

Irrespective of their previously determined rotavirus status, TNA were extracted from 0.2 g of 2003 and 2013 specimens (or 200 µL if liquid) using the cador Pathogen 96 QIAcube HT Kit (Qiagen, Manchester, UK). For 2016 and 2019 specimens, TNA were extracted using QIAamp Fast DNA Stool Mini kit (Qiagen, Manchester, UK) as per the manufacturer's instructions. Fecal specimens from the post-vaccine period (0.2 mg or 200 µL) were subjected to bead beating prior to TNA extraction and collected in a 200 µL of elution buffer [28].

The TNA extracts were screened for the five viruses by a two-step real-time RT-PCR assay [29]. First, cDNA was synthesized in a total volume of 20 µL using random hexamers and 5µL of TNA using the Omniscript Reverse Transcriptase kit (Qiagen, Manchester, UK), as per the manufacturer's instructions. Two µL of the cDNA was henceforth used for real-time RT-PCR in a total volume of 20 µL using the QuantiFast RT-PCR Kit (Qiagen, Manchester, UK) and run on the ABI 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). Primers and probes were adopted from previously published work [30]. The presence of nucleotide mismatches in the primer and probe binding sites was investigated by aligning the primers/probes to genomic sequences deposited in GenBank from 2010 to 2019, using MAFFT v.7.313 [31]. The adenovirus probe/primer pair used in this study detected adenovirus serotypes beyond type 40/41. We used three Ct cut-off values (<40.0, <35.0

and <30.0) to define positive samples. Samples that were positive for RVA in 2003, 2013, 2016 and 2019 were processed for RVA genotyping using VP4 and VP7 RT-PCR, followed by either dideoxy sanger sequencing, as described elsewhere [19], or next-generation sequencing on the Illumina Miseq platform [32].

#### *4.3. Statistical Analysis*

All statistical analyses were performed using R version 3.6.1 [33]. Prevalence was defined as the proportion of these viruses in a hospital-admitted diarrhea patient population during the study period in Kilifi, Kenya. Means and medians of continuous variables were compared using a Kruskal Wallis and Wilcoxon rank-sum test, respectively. Binary data were summarized using proportions and comparisons between groups made using χ 2 statistics. A *p* value of <0.05 was considered statistically significant. Diarrhea severity in RVA positive cases pre- (year 2013) and post- (years 2016 and 2019) was assessed using the Vesikari Clinical Severity Scoring System Manual [10], with a modification in the treatment parameter. If the participant was given oral rehydration therapy or intravenous fluid therapy, they received a score of one or two, respectively.

**Author Contributions:** Conceptualization, A.W.L. and C.N.A.; methodology, A.W.L. and C.N.A.; formal analysis, A.W.L., N.M. and C.N.A.; investigation, A.W.L., S.O., C.S.L. and C.N.A.; resources, D.J.N. and C.N.A.; data curation, A.W.L. and N.M.; writing—original draft preparation, A.W.L.; writing—review and editing, A.W.L., S.O., N.M., C.S.L., S.G.N. and C.N.A.; visualization, A.W.L. supervision, S.G.N. and C.N.A., project administration, C.N.A.; funding acquisition, C.N.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was funded by the Wellcome Trust (102975,203077) The authors Arnold Lambisa, Sylvia Onchaga and Charles Agoti were supported by the Initiative to Develop African Research Leaders (IDeAL) through the DELTAS Africa Initiative (DEL-15-003). The DELTAS Africa Initiative is an independent funding scheme of the African Academy of Sciences (AAS)'s Alliance for Accelerating Excellence in Science in Africa (AESA) and supported by the New Partnership for Africa's Development Planning and Coordinating Agency (NEPAD Agency) with funding from the Wellcome Trust (107769/Z/10/Z) and the UK government. The views expressed in this publication are those of the authors and not necessarily those of AAS, NEPAD Agency, Wellcome Trust or the UK government. This paper is published with the permission of the Director of KEMRI.

**Acknowledgments:** We thank all the study participants for their contribution of study samples, their parents/guardians, members of the viral epidemiology and control research group (http://virec-group.org/) and colleagues at the KEMRI Wellcome Trust Research Programme for their useful discussions during the preparation of the manuscript. We are grateful to James Nokes of KEMRI-Wellcome Trust for his comments and suggestions on the presentation of this work This paper is published with the permission of the Director of KEMRI.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Genotype Diversity before and after the Introduction of a Rotavirus Vaccine into the National Immunisation Program in Fiji**

**Sarah Thomas 1,\*, Celeste M. Donato 1,2 , Sokoveti Covea <sup>3</sup> , Felisita T. Ratu <sup>3</sup> , Adam W. J. Jenney 4,5,6 , Rita Reyburn 4,6, Aalisha Sahu Khan <sup>3</sup> , Eric Rafai <sup>3</sup> , Varja Grabovac <sup>7</sup> , Fatima Serhan <sup>8</sup> , Julie E. Bines 1,2,9,† and Fiona M. Russell 4,6,†**


**Abstract:** The introduction of the rotavirus vaccine, Rotarix, into the Fiji National Immunisation Program in 2012 has reduced the burden of rotavirus disease and hospitalisations in children less than 5 years of age. The aim of this study was to describe the pattern of rotavirus genotype diversity from 2005 to 2018; to investigate changes following the introduction of the rotavirus vaccine in Fiji. Faecal samples from children less than 5 years with acute diarrhoea between 2005 to 2018 were analysed at the WHO Rotavirus Regional Reference Laboratory at the Murdoch Children's Research Institute, Melbourne, Australia, and positive samples were serotyped by EIA (2005–2006) or genotyped by heminested RT-PCR (2007 onwards). We observed a transient increase in the zoonotic strain equine-like G3P[8] in the initial period following vaccine introduction. G1P[8] and G2P[4], dominant genotypes prior to vaccine introduction, have not been detected since 2015 and 2014, respectively. A decrease in rotavirus genotypes G2P[8], G3P[6], G8P[8] and G9P[8] was also observed following vaccine introduction. Monitoring the rotavirus genotypes that cause diarrhoeal disease in children in Fiji is important to ensure that the rotavirus vaccine will continue to be protective and to enable early detection of new vaccine escape strains if this occurs.

**Keywords:** rotavirus; Fiji; Rotarix; genotype; equine-like G3P[8]

#### **1. Introduction**

Rotavirus is the most common cause of severe diarrhoea in children under 5 years of age worldwide. In 2016, rotavirus was responsible for 258 million episodes of diarrhoea and was attributed to ~128,500 deaths in children under 5 years, with the majority occurring in countries in Asia and Africa [1]. Genotyping of rotavirus strains underpins global rotavirus surveillance. The binomial classification of rotavirus genotypes is based on the outer capsid proteins VP7 and VP4 that define G and P genotypes, respectively [2]. There are 36 G

**Citation:** Thomas, S.; Donato, C.M.; Covea, S.; Ratu, F.T.; Jenney, A.W.J.; Reyburn, R.; Sahu Khan, A.; Rafai, E.; Grabovac, V.; Serhan, F.; et al. Genotype Diversity before and after the Introduction of a Rotavirus Vaccine into the National Immunisation Program in Fiji. *Pathogens* **2021**, *10*, 358. https:// doi.org/10.3390/pathogens10030358

Academic Editor: Thirumalaisamy P. Velavan

Received: 8 February 2021 Accepted: 12 March 2021 Published: 17 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

types and 51 P types described in humans and various animal species to date; however, the most common rotavirus genotypes observed in humans are the VP7 genotypes: G1, G2, G3, G4 and G9 and the VP4 genotypes: P[4] and P[8], representing three quarters of all genotypes causing human disease [3,4]. Previously uncommon genotypes including G12 and equine-like G3P[8] genotypes are increasingly being identified as a cause of rotavirus disease globally [5,6].

Fiji is a Pacific Island Nation with a population of approximately 837,271 [7]. Although designated as an upper middle-income country, it was estimated that prior to the COVID-19 pandemic, 24% of the population were living in poverty [7]. The child under-5-year mortality rate in Fiji was reported as 25.7 deaths per 1000 live births in 2019 [8]. Rotavirus was a major cause of diarrhoea-related hospitalisations in Fiji prior to rotavirus vaccine introduction, detected in 52% (2006) and 60% (2007) of children less than 5 years hospitalised with acute diarrhoea, with an annual incidence estimated at 486 per 100,000 children less than 5 years [9]. Due to this burden of rotavirus gastroenteritis, Fiji introduced a rotavirus vaccine (Rotarix, GlaxoSmithKline, Belgium) into the National Immunisation Program in October 2012. Rotarix is a monovalent vaccine containing a single, human, G1P[8] strain that is administered in a two-dose schedule at 6 and 14 weeks of age. The uptake of Rotarix in Fiji was prompt, reaching 85% coverage by 2013 and 99% coverage in eligible infants from 2014 onward [10]. The introduction of rotavirus vaccines in Fiji has been highly successfully resulting in an 82% reduction in rotavirus diarrhoea related hospitalisations in children less than 5 years of age [11].

The aim of this study was to describe the pattern of rotavirus genotype diversity from 2005 to 2018, specifically to describe any changes in genotype patterns that may have occurred following the introduction of the rotavirus vaccine in Fiji in 2012.

#### **2. Results**

#### *2.1. Study Samples*

During the study period 2005–2018, a total of 1504 stool samples was collected and sent to the WHO Rotavirus Regional Reference Laboratory (RRL) at the Murdoch Children's Research Institute (MCRI). Of these, 1208 samples had sufficient data available on the date of collection and stool volume to enable analysis. Of the 1208 samples, a total of 576 were confirmed as rotavirus positive and proceeded to genotype characterisation (Figure 1). Thirty-four samples were not genotyped due to laboratory error, comprising 1 sample from 2010 and 33 samples from 2011, and were subsequently excluded from further analysis. The remaining 542 samples were proceeded with for further analysis (Figure 1).

#### *2.2. Genotype Distribution and the Impact of Vaccine Introduction*

In the pre-vaccine period (2005–2012), 58% (479/827) of samples received were confirmed as rotavirus positive, compared to only 18% (63/347) of samples in the post-vaccine era (2013–2018) (Table 1). These values may be affected by sampling changes over the study period. Between 2005 and 2009, only positive samples were received; between 2010 and 2016, all positive and negative samples were received; and from 2017 onward, all positive and 10% of all negatives were sent to MCRI. Overall, between 2005–2018, G1P[8] was the most commonly detected genotype (*n* = 157, 29%), with both G2P[4] (*n* = 155, 29%) and G3P[8] (*n* = 144, 27%) detected at similar frequencies, followed by G12P[8] (*n* = 33, 6%) (Table 1). Other genotypes including G2P[8], G3P[6], G8P[8], G9P[8] and G12P[4] as well as mixed or partially typed samples were infrequently detected (*n* = 1–5, 0.2–3%). However, marked differences were observed following vaccine introduction. Prior to vaccine introduction, genotype dominance varied across years, with G3P[8] dominant in 2006 (*n* = 74, 94%) and 2009 (*n* = 26, 59%), G2P[4] dominant in 2008 (*n* = 30, 40%) and 2010 (*n* = 88, 74%), and G1P[8] dominant in 2011 (*n* = 127, 85%) and 2012 (*n* = 4, 67%). However, the number of samples available for genotyping was low in 2005 and 2007 and no clear dominant genotype could be determined.

**Figure 1.** Consort diagram of samples included in this study.

Following vaccine introduction, there was a marked decrease in the number of rotavirus positive samples available for genotyping, reflecting the reduction in rotavirus disease observed (Table 1). The diversity of genotypes decreased following vaccine introduction (Figure 2) with some genotypes (G2P[8], G3P[6], G8P[8], G9P[8], G12P[4]) no longer detected. There was only one mixed genotype sample identified in the post-vaccine era, compared with 13 mixed genotype samples detected in the pre-vaccine era. The dominant genotype continued to vary annually following vaccine introduction, with G3P[8] dominant in 2013 (*n* = 6, 50%), G1P[8] in 2014 (*n* = 17, 71%), G12P[8] in 2017 (*n* = 11, 100%) and G3P[8] in 2018 (*n* = 6, 100%). The previously dominant G1P[8] disappeared 3 years after vaccine introduction, and G2P[4] strains were not detected after 2014 (Table 1). Emergence of the novel, equine-like G3P[8] reassortant strain, previously not detected in Fiji, was reported in the years following vaccine introduction. This equine-like G3P[8] was dominant for two consecutive years (2015–2016), accounting for 83% (*n* = 5/6) and 100% (*n* = 4/4) of samples genotyped. However, it was not detected in 2017 or 2018. G3P[8] re-emerged in 2018 after not being detected for 4 years.

*Pathogens* **2021**, *10*, 358


**Table 1.** Genotype distribution in samples received by the WHO Regional Reference Laboratory.

–

– **Figure 2.** Distribution of main genotypes in samples collected in the pre-vaccine (2005–2012) and post-vaccine (2013–2018) period. (**a**) Number of samples in each of the main genotype groups. (**b**) Proportion of samples identified in each of the main genotype groups, of the total number of samples genotyped.

–

#### **3. Discussion**

This is the first study in a low- or middle-income country in the Western Pacific Region to describe rotavirus genotypes following national rotavirus vaccine introduction. Prior to rotavirus vaccine introduction, G1P[8], G2P[4] and G3P[8] were the predominant genotypes causing rotavirus diarrhoea in children less than 5 years of age in Fiji. Genotype diversity decreased following rotavirus vaccine introduction in Fiji; with G2P[8], G3P[6], G8P[8], G9P[8] and G12P[4], which all represented minor genotypes in the pre-vaccine period, subsequently undetected in the vaccine era. Following rotavirus vaccine introduction, G2P[4] has not been detected since 2014 and G1P[8] has not been detected since 2015. This is in contrast to changes in genotype distribution observed in Australia following introduction of the Rotarix (GlaxoSmithKline, Rixensart, Belgium) and RotaTeq (Merck, Kenalworth, NJ, USA) vaccines. In Australia, although there was an overall decrease in the common genotypes (G1, G2, G3, G4 and G9) from 83% to 63% observed following rotavirus vaccine introduction, an increase in G2P[4] (pre-vaccine era 5%; post-vaccine era 21%) was observed in states and territories implementing the Rotarix vaccine and G1P[8] continued to be detected [5].

We found that equine-like G3P[8] was the dominant genotype in 2015 to 2016 but was not detected in the following years (2017 or 2018). An increase in novel zoonotic strains such as equine-like G3P[8] following the introduction of Rotarix has also been observed in other countries (Australia, Japan, Hungary and Brazil) [12–16]. The segmented rotavirus genome allows reassortment to occur both within and between human and animal strains if the human host is infected with two different rotavirus strains, thus giving rise to novel and unusual genotype combinations [17]. In Australia, an increase in G12P[8], equine-like G3P[8], G8, G10 and other zoonotic reassortant strains has also been observed following rotavirus vaccine introduction [5]. In Fiji, human G3P[8] was not detected during 2015 and 2016 when the equine-like G3P[8] was circulating, but this strain re-emerged two years later when equine-like G3P[8] was no longer detected. This is consistent with reports from Asia, Australia, Europe and the U.S. [5,17].

The G12 genotype was first identified in Fiji in 2008, with both G12P[4] and G12P[8] detected. These strains accounted for 32% (*n* = 24/75) of all rotavirus positive samples in 2008 but were not detected again until 2017 when all available samples (*n* = 11) were identified as G12P[8] (Table 1). The emergence of G12 following vaccine introduction has been observed in other countries but does not appear to be dependent on vaccine coverage. In Finland, a 9% increase in G12P[8] was observed five years after vaccine introduction, with a higher frequency of G12P[8] detected in vaccinated children (14%) than observed in unvaccinated children (7%) [18]. In Australia, a small G12P[8] outbreak was reported in 2005 prior to vaccine introduction; however, since vaccine introduction G12P[8] has become common, detected in 18% of samples from children less than 5 years with acute diarrhoea [5]. Similarly, G12, originally detected in Brazil in 2008 following vaccine introduction (2006), has emerged to be the most prevalent genotype (G12P[8]) in 87% of samples in 2014 [19].

A key strength of this study is the ability to observe genotypic changes over time, following the introduction of Rotarix into a national program in a Pacific nation associated with very high vaccine coverage. Monitoring rotavirus genotypes that continue to cause diarrhoea in children provides critical information regarding the ongoing effectiveness of the vaccine program and can assist in outbreak investigation. It also enables early identification of the emergence or importation of new strains that may have a public health impact. This is particularly relevant for Fiji as an island nation with an economy highly dependent on tourism where there is potential for importation of novel strains resulting in disease outbreaks.

This study aligns with data on the impact of rotavirus vaccines on rotavirus disease hospitalisations in children less than 5 years of age in Fiji. Fiji has been notable within the Pacific as a country that has introduced new vaccines based on local data and is committed to monitoring vaccine impact. No other Pacific nation participates in WHO rotavirus surveillance. Data from Fiji may assist in informing vaccine decisions of neighbouring countries in the region. A limitation of this study is that it can only report on samples received for analysis by the WHO Rotavirus Regional Reference Laboratory. Despite attempts, not all children admitted to hospital with diarrhoea have a stool sample collected and sent for analysis. Following introduction of a rotavirus vaccine, the number of children hospitalised with rotavirus disease has dramatically decreased; as a result, the number of stool samples available to provide comparisons of genotypic distribution between the pre-vaccine and post-vaccine era has been impacted. As stool collection is still requested for hospitalised patients with acute diarrhoea in Fiji, it is unlikely that there is a bias impacting on stool collection between the period before and after introduction of the rotavirus vaccine.

The variation in the proportion of rotavirus negative samples reported reflects differences in the rotavirus detection status of stool samples submitted to MCRI for genotypic analysis over the 14-year surveillance period (Table 1). The lower proportion of rotavirus negative samples early in the surveillance period (2005–2009) has limited impact on the

outcome of this paper given the focus is the period following the introduction of the rotavirus vaccine.

In 2015, there was a marked increase in the number of negative samples tested. From 2010 to 2013, there was an increase in typhoid detection as the result of six typhoid outbreaks in Fiji, along with outbreaks of both Zika virus and Chikungunya virus, which were both initially detected in 2015, all which may have led to an increased number of negative samples being sent to MCRI for analysis during this time period [20–22]. Being negative samples only, this also would have had minimal impact on the rotavirus distribution observed in this study.

The effect of age on rotavirus detection in the stool following introduction of rotavirus vaccines in Fiji has recently been reported [11]. Due to the success of the rotavirus vaccination program, the ability to compare age related differences in genotype distribution in samples from the pre- and post-vaccine eras has been impacted by the limited number of samples available for analysis in the post-vaccine era (pre-vaccine era *n* = 462; post-vaccine era *n* = 50). This decline in number of available rotavirus positive samples was not likely to be due to a lack of sampling due to the ongoing surveillance program operating in Fiji.

In this study we report changes in the pattern of rotavirus genotypes causing diarrhoea in children in Fiji since rotavirus vaccine introduction. We observed a transient increase in the zoonotic strain equine-like G3P[8] and a reduction in the previously dominant G1P[8] and G2P[4]. A decrease in detection of rotavirus strains G2P[8], G3P[6], G8P[8] and G9P[8] was also detected in the years following rotavirus vaccine introduction. Monitoring rotavirus genotypes provides key information regarding the ongoing effectiveness of the vaccine program and can assist in outbreak investigation. It also enables early identification of the emergence or importation of new strains that may have a public health impact.

#### **4. Materials and Methods**

#### *4.1. Population and Study Sites*

Fiji has participated in the WHO Global Rotavirus Surveillance Program since 2006, monitoring rotavirus disease burden and rotavirus genotype diversity associated with hospitalisations in children less than 5 years of age. The samples from children hospitalised with acute diarrhoea are sent to the WHO Rotavirus Regional Reference Laboratory at Murdoch Children's Research Institute (MCRI) according to the case definitions and methods defined for the WHO Global Rotavirus Surveillance Program [11]. Two major hospitals admitting children with gastroenteritis in Fiji participated in this study. The Colonial War Memorial Hospital (CWMH) is Fiji's largest general hospital and is the main referral centre for the greater Suva area with approximately 34,920 children under 5 years of age. The Savusavu District Hospital is a secondary health inpatient and outpatient facility serving mainly a semiurban and rural population with an estimated 6563 children under 5 years of age. Samples were received from inpatients at Savusavu district hospital only in 2013 and 2014. Details on rotavirus surveillance in Fiji has previously been described [9,11].

#### *4.2. Participants*

A prospective rotavirus surveillance program was established in 2005 to effectively capture rotavirus detected in faecal samples from children less than 5 years with acute nonbloody diarrhoea in Fiji. Acute diarrhoea was defined as 3 or more loose, nonbloody stools within a 24-h period for <14 days. Eligible participants were identified by checking admission data and children's wards daily, parental/guardian consent was obtained, and stool was collected within 48 h of admission. Once rotavirus positivity was determined, demographics and clinical information were obtained via medical records.

Ethics approval for these studies was obtained from the Fiji National Research Ethics Review Committee (number 2013-40) and from the University of Melbourne Human Research Ethics Committee for the initial study surveillance in Colonial War Memorial Hospital from 2005–2012 (Ethics ID:050546X) and Savusavu from 2010–2012 (Ethics ID:0931282); during this period written informed consent was obtained from participants' parents. From

June 2012 onward, the Ministry of Health and Medical Services considered this public health surveillance and no longer required written consent.

#### *4.3. Genotyping*

Faecal specimens were collected and stored at 4–8 ◦C prior to being transported to the Fiji Centre for Communicable Disease Control in Suva for rotavirus antigen testing via the ProSpecT Rotavirus test, a commercial rotavirus enzyme immunoassay (EIA) (Thermofisher Scientific, Waltham, MA, USA) as per manufacturer's instructions. Stool samples were then stored at −70 ◦C. De-identified rotavirus specimens were transported on dry ice to the WHO Rotavirus Regional Reference Laboratory at the Murdoch Children's Research Institute, Parkville, Australia. Sample selection for shipment to MCRI varied during the surveillance program. Between 2005 and 2009 only stool samples that tested positive to rotavirus in the Fiji laboratory were sent to MCRI, with the negative samples reflected in Table 1 having been identified by EIA at MCRI. Between 2010 and 2016, stool samples were sent to MCRI for EIA and RT-PCR genotyping irrespective of whether they were rotavirus positive or negative via EIA conducted in Fiji. From 2017 onward, all rotavirus positive samples and 10% of rotavirus negative samples by EIA in Fiji were sent to MCRI for further analysis. All samples were retested at MCRI to confirm rotavirus positivity prior to proceeding with genotypic analysis. Rotavirus positivity (or negativity) was confirmed using the ProSpecT Rotavirus test, (EIA) (Thermofisher Scientific, Waltham, MA, USA) as per the manufacturer's instructions. Stool samples that tested positive or equivocal for rotavirus antigen were further characterised to determine the G and P genotype. Samples from 2005 and 2006 were routinely serotyped using an in-house monoclonal antibody based serotyping EIA. This EIA consisted of a panel of monoclonal antibodies specific to the VP7 outer capsid protein of group A rotavirus serotypes G1, G2, G3, G4 and G9 [23]. Prior to 2007, P-typing and RT-PCR were not routinely performed. From 2007 onward, rotavirus G and P genotypes were determined by heminested multiplex RT-PCR assay. All samples collected prior to 2007 have retrospectively been characterised by heminested RT-PCR G and P genotyping. In brief, viral RNA was extracted from 20% (*w*/*v*) faecal extracts in a virus dilution buffer (0.01 M Tris-HCL [pH7.5], 10.5 mM CaCl, 145 mM NaCl) using the QIAamp Viral RNA mini extraction kit (QIAGEN, Hilden, Germany) according to the manufacturer's instructions.

The One-step RT-PCR kit (QIAGEN) was used to perform first round PCR, using VP7 primers VP7F and VP7R and VP4 primers VP4F and VP4R [24,25]. Second round genotyping PCR was performed using AmpliTaq DNA Polymerase with Buffer II (Applied Biosystems, Foster City, CA, USA), with specific G and P oligonucleotide primers for G typing (G1, G2, G3, G4, G8 and G9) or P typing (P[4], P[6], P[8], P[9], P[10] and P[11]) as described previously [26]. Amplified products were run on a 1.5% or 2% agarose gel for G and P types respectively, and genotypes were determined based on amplicon band size. PCR non-typeable samples were determined by Sanger sequencing. Strains including equine-like G3, G12 and unusual or uncommon strains were unable to be genotyped using standard primers. VP7 or VP4 amplicons from first round PCR products were purified for sequencing using the Wizard SV Gel and PCR Clean up System (Promega, Madison, WI, USA) as per manufacturer's protocol. Purified DNA with oligonucleotide primers (VP7F/R or VP4F/R) were sent to the Australian Genome Research Facility (AGRF) Melbourne and sequenced using an ABI PRISM BigDye Terminator Cycle Sequencing Reaction Kit (Applied Biosystems, Foster City, CA, USA) in an Applied Biosystems 3730xl DNA analyser. Sequencher version 4.10.1 (Gene Codes Corporation, Ann Arbor, MI, USA) was used to edit the sequences. BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi (accessed on 10 October 2020)) and RotaC version 2.0 (http://rotac.regatools.be (accessed on 10 October 2020)) [27] were used determine the genotype of each sample.

#### *4.4. Data Analysis*

Samples were excluded if there was no date of stool collection available, if there was insufficient sample to process, if the sample was not confirmed as rotavirus positive by EIA at MCRI, or if samples were rotavirus positive by EIA at MCRI but genotype could not be determined. To describe the impact of rotavirus vaccine introduction on genotype distribution, samples were grouped into pre-vaccine (2005–2012) and post-vaccine (2013–2018) eras according to the date of collection. Analysis is by descriptive observations and comparisons between the pre-vaccine and post-vaccine introduction eras.

**Author Contributions:** Conceptualization, S.T., J.E.B. and F.M.R.; data curation, S.T. and C.M.D.; formal analysis, S.T., C.M.D. and J.E.B.; funding acquisition, J.E.B. and F.M.R.; investigation, S.T., C.M.D., S.C., F.T.R., A.W.J.J., R.R., A.S.K., E.R., V.G., F.S., J.E.B. and F.M.R.; methodology, S.T. and C.M.D.; project administration, J.E.B. and F.M.R.; resources, J.E.B. and F.M.R.; software, S.T. and C.M.D.; supervision, J.E.B. and F.M.R.; visualization, S.T. and C.M.D.; writing—original draft, S.T.; writing—review and editing, C.M.D., S.C., F.T.R., A.W.J.J., R.R., A.S.K., E.R., V.G., F.S., J.E.B. and F.M.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** The identification, collection, transportation and laboratory analysis of surveillance samples was funded as part of the World Health Organization Rotavirus Regional Surveillance activity. The World Health Organization provided funds to set up the initial rotavirus surveillance (2005–2007) and an additional World Health Organization grant [Registry File No. V27-181-188] for data collected between 2006 and 2013. A Merck investigator grant (IISP ID#:35248) funded the early work in Savusavu (2009–2010). Surveillance data collected between 2014 and 2018 by the Department of Foreign Affairs and Trade of the Australian Government and Fiji Health Sector Support Program (FHSSP). FHSSP is implemented by Abt JTA on behalf of the Australian Government. The Murdoch Children's Research Institute is supported by the Victorian Government's Operational Infrastructure Support program. Funders had no role in study design, data collection, data analysis, interpretation, writing of the report. C.M.D. is supported through the Australian National Health and Medical Research Council with an Early Career Fellowship (1113269). F.M.R. is supported through Australian National Health and Medical Research Council with an Early Career Fellowship and Translating into Practice Fellowship.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Fiji National Research Ethics Review Committee (number 2013-40) and from the University of Melbourne Human Research Ethics Committee for the initial study surveillance in Colonial War Memorial Hospital 2005–2012 (Ethics ID:050546X) and Savusavu from 2010–2012 (Ethics ID:0931282).

**Informed Consent Statement:** Informed consent was obtained from participant's parents between 2005 and 2012; from June 2012, the Fiji Ministry of Health and Medical Services considered this public health surveillance and no longer required written consent.

**Data Availability Statement:** Data is contained within the article.

**Acknowledgments:** We gratefully acknowledge Rachel Devi; Kathryn Bright; Beth Temple; Lisi Tikoduadua; Joe Kado; E. Kim Mulholland; Kimberley K. Fox; and Mike Kama for assistance in collection of stool samples in Fiji. We acknowledge Josephine Logronio, WHO Western Pacific Regional Office. We gratefully acknowledge Carl D. Kirkwood, Nada Bogdanovic-Sakran, Huy Tran and the Enteric Disease Group MCRI for their assistance within the laboratory.

**Conflicts of Interest:** C.M.D. has served on a rotavirus advisory board for GSK (2019); all payments were paid directly to an administrative fund held by Murdoch Children's Research Institute. J.E.B. is lead for the Rotavirus Vaccine Program at Murdoch Children's Research Institute that aims to develop an affordable rotavirus vaccine, RV3-BB. J.E.B. is Director of the Australia Rotavirus Surveillance Program that receives funding from the Australian Commonwealth Department of Health and Aging and GlaxoSmithKline. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. All other authors declare no conflict of interest.

#### **References**


### *Article* **Rotavirus Strain Distribution Before and After Introducing Rotavirus Vaccine in India**

**Tintu Varghese <sup>1</sup> , Shainey Alokit Khakha <sup>1</sup> , Sidhartha Giri <sup>1</sup> , Nayana P. Nair <sup>1</sup> , Manohar Badur <sup>2</sup> , Geeta Gathwala <sup>3</sup> , Sanjeev Chaudhury <sup>4</sup> , Shayam Kaushik <sup>5</sup> , Mrutunjay Dash <sup>6</sup> , Nirmal K. Mohakud <sup>7</sup> , Rajib K. Ray <sup>8</sup> , Prasantajyoti Mohanty <sup>8</sup> , Chethrapilly Purushothaman Girish Kumar <sup>9</sup> , Seshadri Venkatasubramanian <sup>9</sup> , Rashmi Arora <sup>10</sup>, Venkata Raghava Mohan <sup>11</sup> , Jacqueline E. Tate <sup>12</sup> , Umesh D. Parashar <sup>12</sup> and Gagandeep Kang 1,\***

	- drrajib2007@gmail.com (R.K.R.); prasantij53@gmail.com (P.M.)

**Abstract:** In April 2016, an indigenous monovalent rotavirus vaccine (Rotavac) was introduced to the National Immunization Program in India. Hospital-based surveillance for acute gastroenteritis was conducted in five sentinel sites from 2012 to 2020 to monitor the vaccine impact on various genotypes and the reduction in rotavirus positivity at each site. Stool samples collected from children under 5 years of age hospitalized with diarrhea were tested for group A rotavirus using a commercial enzyme immunoassay, and rotavirus strains were characterized by RT-PCR. The proportion of diarrhea hospitalizations attributable to rotavirus at the five sites declined from a range of 56–29.4% in pre-vaccine years to 34–12% in post-vaccine years. G1P[8] was the predominant strain in the pre-vaccination period, and G3P[8] was the most common in the post-vaccination period. Circulating patterns varied throughout the study period, and increased proportions of mixed genotypes were detected in the post-vaccination phase. Continuous long-term surveillance is essential to understand the diversity and immuno-epidemiological effects of rotavirus vaccination.

**Keywords:** rotavirus diarrhea; rotavirus genotyping; Rotavac vaccine

**Citation:** Varghese, T.; Alokit Khakha, S.; Giri, S.; Nair, N.P.; Badur, M.; Gathwala, G.; Chaudhury, S.; Kaushik, S.; Dash, M.; Mohakud, N.K.; et al. Rotavirus Strain Distribution Before and After Introducing Rotavirus Vaccine in India. *Pathogens* **2021**, *10*, 416. https://doi.org/10.3390/pathogens 10040416

Academic Editors: Julie Bines and Celeste Donato

Received: 22 February 2021 Accepted: 18 March 2021 Published: 1 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Rotavirus is the leading etiology of acute gastroenteritis in children under 5 years old worldwide, causing high mortality, especially in middle- and low-income countries. India accounts for 22% of the total global rotavirus mortality [1]. In India, 40% of all diarrhea-related hospitalizations among children under 5 years of age is caused by group A rotavirus [2].

The genome of group A rotavirus is composed of 11 double-stranded RNA segments, of which the VP7 and VP4 genes coding for the outer capsid proteins are used for the classification of the virus into G and P types, respectively. Studies have been conducted across the globe to understand the natural evolution of rotavirus and its relevance in the context of vaccine introduction. Globally, G1P[8], G2P[4], G3P[8], and G9P[8] are the most common genotypes associated with rotavirus diarrhea [3]. However, it is hypothesized that large-scale vaccination may exert pressure on circulating strains, leading to possible changes in strain circulation.

In 2009, the World Health Organization (WHO) recommended the inclusion of rotavirus vaccines in the national immunization program of all countries. Currently, four live-attenuated oral vaccines are prequalified by WHO, which includes Rotarix (Glaxo-SmithKline Biologicals, Rixensart, Belgium), RotaTeq (Merck & Co., Inc., West Point, PA, USA), Rotavac (Bharath Biotech, India), and Rotasiil (Serum Institute of India PVT. LTD., Pune, India). These rotavirus vaccines differ in their genotypic composition, with Rotarix and Rotavac being the monovalent vaccines and RotaTeq and Rotasiil being the pentavalent vaccines [4]. Rotarix and RotaTeq vaccines have been available on the market since 2006 and are currently used by nearly 90 countries in their immunization programs [4]. Early studies reported a decline in G1P[8] and the emergence of G2P[4] after Rotarix vaccination [5,6], while others showed no change [7]. The emergence of G9P[8] and G12P[8] was reported with the use of the RotaTeq vaccine [8]. However, such changes were also observed in other countries without rotavirus vaccination [9,10]. Hence, the vaccine impact on the circulating pattern of rotavirus strains is not clearly understood.

In India, the indigenously developed Rotavac vaccine, based on the human-bovine reassortant neonatal attenuated 116E strain, is a monovalent vaccine with the genotypic composition G9P [11]. It was introduced to the Universal Immunization Program (UIP) in April 2016 in a phased manner [11]. Other rotavirus vaccines like Rotateq and Rotarix were available in the private sectors for immunization before nationwide rotavirus vaccine implementation. India is the first Asian country to introduce rotavirus vaccines to the national immunization schedule, and currently, the Rotavac vaccine is used only in India and a few smaller countries [4]. The National Rotavirus Surveillance Network was established in India in 2005 to generate data on disease burden and monitor the trends of circulating genotypes [12,13]. This study describes the reduction in rotavirus prevalence and temporal trends in rotavirus strain distribution before and after Rotavac vaccine introduction in five sites in India.

#### **2. Results**

#### *2.1. Prevalence of Rotavirus Diarrhea*

Between September 2012 and June 2020, 8499 children under 5 years of age were enrolled in the surveillance study at the five sites. The details of enrollment and rotavirus testing are summarized in Table 1.


**Table 1.** Enrollment and rotavirus testing details from 5 surveillance sites (September 2012–June 2020).

The proportion of diarrhea hospitalizations attributable to rotavirus at the five sites declined from a range of 56–29.4% in pre-vaccine years to 34–12% in post-vaccine years. The maximum annual positivity rate was in 2014 (46.2%), and the minimum was in 2019 (13.3%). The positivity rates declined steadily after vaccine implementation and were more marked towards the later years with higher vaccine coverage (Figure 1). The maximum reduction in rotavirus diarrhea was seen in Tirupati (72.1%), the site with maximum vaccine coverage, compared to a 32.5% reduction in Vellore, which was the last to introduce the vaccine and hence had the lowest overall vaccine coverage among the five sites.

**Figure 1.** Impact of rotavirus vaccine after its introduction into the universal immunization programme in India, prevaccination and post-vaccination introduction surveillance comparison data from study sites at Rohtak (**A**), Tandak (**B**), Tirupati (**C**) Bhubaneswar (**D**), Vellore (**E**), and all the sites combined (**F**).

#### *2.2. Rotavirus Genotype Distribution in India*

During the study period, genotyping was performed for 76.04% of the samples. The proportion of positive samples tested by genotyping PCR was greater in the postvaccination period (97.02%) compared to the pre-vaccination period (64.97%), when the protocol changed for genotyping of a subset of samples.

G1P[8] was the most common strain (49.5%) in the pre-vaccine period. The other common genotypes were G2P[4] (8%), G9P[4] (7.5%), G9P[8] (4.5%), and G12P[6] (3.8%). Conversely, G3P[8] (44.3%) was the most common genotype in the post-vaccine period, with G1P[8] (15.4%), G2P[4] (7.4%), G9P[4] (4.9%), and G1P[6] (3.7%) being the next most common genotypes (Figure 2). Marked yearly changes were seen among the circulating strains. Circulation of G9P[8] peaked during the year 2013, while G12P[6] increased in 2014/2015. Some reassortant strains like G1P[4], G2P[6], G2P[8], G3P[4], G3P[6], and G4P[6] were occasionally reported during the study period (Figure 2).

**Figure 2.** Data represented as the proportion of a specific genotype compared to the total genotype results. Uncommon genotype: <1% of total results; Mixed genotypes: those with >1 G or P-type; Untypables: those with either G or P untyped.

> The genotype distribution also varied across the sentinel sites in North India (Tanda and Rohtak) and South India (Vellore, Tirupati, and Bhubaneswar). G1P[6] was seen predominantly in northern sites, while G9P[8] and G12P[6] were seen in southern sites

during the pre-vaccination period. In the post-vaccination period, the major circulating strains remained the same in northern sites, with G3P[8] topping the list. G3P[8] emerged in the southern sites as well, with a decline in G9P[8] and G12P[6] (Figure 3).

**Figure 3.** Comparison of genotype distribution between Northern sites (Tanda and Rohtak) and Southern sites (Vellore, Tirupati and Bhubaneswar). The major genotypes are compared during pre-vaccination (September 2012–April 2016) and post-vaccination period (May 2016–June 2020). The mixed genotype infections were excluded from the analysis.

> An increased prevalence of G3P[8] and decreased prevalence of G1P[8] were noted in the post-vaccination period compared to the pre-vaccination period. G1P[8] peaked during the year 2014 (62.6%) and has declined steadily since then. G3P[8] started appearing in 2015 and was the predominant genotype in the following years. No novel strains were detected during the post-vaccination period. Mixed genotype infections occurred in a higher proportion in the post-vaccine period (17.4%) compared to the pre-vaccine period (6.4%). G1 (33%) was the most common G-type found in mixed infections, mainly in combination with G12 (10.8%) and G3 (9%). Similarly, P[8] (93.4%) was the most common P-type in mixed infections, along with P[4] (55.6%) and P[6] (37.8%).

#### **3. Discussion**

Pre- and post-introduction surveillance at five sites in India indicate that vaccination is impacting severe rotavirus gastroenteritis. The overall prevalence of rotavirus in children with hospitalized gastroenteritis decreased after vaccine introduction, reaching 13.3% by the third year post-vaccine introduction, indicating the effectiveness predicted by clinical trials and modeling [14,15]. The maximum reduction rate was seen in Tirupati (72.1%), and the minimum was observed in Vellore (32.45%), which are the sites with maximum and minimum vaccine coverage, respectively.

During the study period, from 2012 to 2020, the major genotypes were G1P[8], G3P[8], G2P[4], G9P[4], G9P[8], G12P[6], and G1P[6], which include some reassortants that are not common in other parts of the world. There was marked temporal fluctuation, with G9P[8] detected at a high frequency in 2013/2014, only to disappear by 2019/2020, while G12P[6] was high in 2014/2015. Our findings are consistent with surveillance data from India and neighboring countries that also saw the emergence of G12 strains [16]. We also noted a geographic variation, with G12P[6] and G9P[8] seen more in the southern sites in 2012-2016 and G1P[6] observed more in the northern sites. These findings are in agreement with other studies conducted in the northern and southern parts of India [17–20]. The genotypic pattern in the northern sites had a rise in G3P[8] in the post-vaccine period compared to the pre-vaccine period, along with the disappearance of G9P[8] and the emergence of G12P[6]. However, the southern sites had a greater proportion of G3P[8] in the post-vaccination period, with a decline in both G12P[6] and G9P[8]. Variation in the geographic and temporal trends of rotavirus strains emphasizes the importance of multicentric studies.

Changes in genotype distribution and increased diversity are seen with other rotavirus vaccines. In Brazil, G2P[4] emerged as the major strain, while no change in genotype distribution was seen in Kenya after Rotarix introduction [7,21]. An increase in G3P[8] strain prevalence was seen in the United States after RotaTeq introduction [21,22]. In our study, G1P[8] was the predominant strain in the pre-vaccine period, coinciding with other studies conducted during this period [17,23,24], which declined thereafter with the emergence of G3P[8]. However, the rise in G3P[8] in 2017/2018 is likely to be a natural fluctuation rather than the effect of the vaccine, as there was a similar trend seen in other countries without rotavirus immunization [25,26]. There are currently ongoing efforts to examine rotavirus vaccine effectiveness against diseases caused by specific strains, which will help further address this issue.

Globally, the rate of mixed rotavirus infections is similar to our findings [23,27]. Mixed rotavirus infections can facilitate the evolution of novel strains by genetic reassortment between the segmented genes of rotavirus, eventually increasing its diversity. Other studies have reported an increased frequency of unusual and novel strains in the post-vaccine surveillance period [24,28]. In our study, reassortant strains including G1P[4], G2P[6], G3P[6], G3P[4], and G4P[6] were occasionally seen, with no specific increase in the postvaccination phase. Whole-genome sequencing and phylogenetic analysis will help to identify possible reassortment of rotavirus genes and detect mutational events. This will help us in tracking the virus evolution over time, which might give us more insight into the drivers of viral strain circulation and the impact of vaccines.

To conclude, our study showed a reduction in rotavirus diarrhea across five sites in India after Rotavac vaccine introduction. Changes in circulating strains with an increased rate of mixed infections were also seen in the post-vaccine period. In our study from 2016, additional methods were used for the genotyping of samples that remained untyped with standard laboratory protocols. Some of the differences in genotypes before and after vaccination introduction may have been caused by the change in genotyping methods. Due to the short period of surveillance, it is difficult to determine whether the changes were due to natural strain variations or vaccine pressure. Continued surveillance is warranted to determine the long-term effects of rotavirus vaccination.

#### **4. Materials and Methods**

#### *4.1. Study Sites*

Active hospital-based surveillance for diarrhea was established in five sentinel sites consisting of major referral hospitals from September 2012 to June 2020. The hospitals included were Christian Medical College (Vellore, Tamil Nadu), Sri Venkateshwara Medical College (Tirupati, Andhra Pradesh), Hi-Tech hospital (Bhubaneswar, Odisha), Pt. Bhagwat Dayal Sharma Postgraduate Institute of Medical Sciences (Rohtak, Haryana), and Rajendra Prasad Government Medical College (Tanda, Himachal Pradesh).

#### *4.2. Sample Collection and Laboratory Testing*

Sample collection and laboratory methods are detailed in the study protocol [2]. In brief, children under 5 years of age hospitalized with diarrhea were enrolled in the study. A stool sample, vaccination card copy, and case report form with clinical and demographic details were collected from each child. Samples were stored at the appropriate temperature until transported to CMC, Vellore, which served as the main testing laboratory. All testing was done as per the modified WHO generic protocol for rotavirus surveillance [29]. Stool samples were screened for rotavirus VP6 antigen using a commercial enzyme immunoassay (EIA). All EIA positive samples were further characterized by reverse transcription-polymerase chain reaction (RT-PCR) for VP7 (G Type) and VP4 (P Type) genes. In brief, RNA was extracted from 20% fecal suspension using the QIAamp Viral RNA Mini Kit (Qiagen). Complementary DNA (cDNA) synthesized by reverse transcription using Moloney murine reverse transcriptase enzyme (Superscript II MMLV-RT, Invitrogen) and random primers (Invitrogen) were used as templates for VP7 and VP4 typing by a hemi-nested multiplex PCR using published primers [30,31]. For the samples collected in the post-vaccination period, additional typing methods were used if they remained untyped with standard laboratory testing protocols [32]. The negative samples by genotyping PCR were confirmed for rotavirus positivity by VP6 PCR [32]. The untyped samples and unusual rotavirus strains were sequenced by the Sanger sequencing method.

**Author Contributions:** Conceptualization, G.K., U.D.P., and J.E.T.; methodology, G.K., S.G., S.A.K., T.V.; software, S.V.; validation, G.K., S.G., C.P.G.K.; formal analysis, T.V.; investigation, all authors; resources, U.D.P., J.E.T., G.K.; data curation, S.V., V.R.M.; writing—original draft preparation, T.V.; writing—review and editing, G.K.; visualization, T.V., N.P.N., S.G.; supervision, S.G., R.A., G.K.; project administration, G.K.; funding acquisition, G.K., U.D.P., J.E.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Indian Council of Medical Research [5/8-1(189)/TF/2011- 12-ECD] and the Bill and Melinda Gates Foundation [OPP1053989]. SG was supported on TW0007392.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Institutional Review Board of Christian Medical College, Vellore, and all participating institutions.

**Informed Consent Statement:** Written informed consent was obtained from parents/guardian of all children enrolled in this study.

**Data Availability Statement:** Since the study is continuing in some settings, the data are still being generated and have not yet been placed in a public repository. The data analyzed during the period reported will be made available on request after de-identification.

**Acknowledgments:** We would like to thank the families of infants and children who participated in the surveillance study, the sentinel site study teams, and Christian Medical College laboratory team.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Disclaimer:** The findings and conclusions in this report are those of the authors and do not necessarily represent the official position of the Centers for Disease Control and Prevention.

#### **References**


#### *Article*

### **Molecular Epidemiology of Rotavirus A Strains Pre- and Post-Vaccine (Rotarix**®**) Introduction in Mozambique, 2012–2019: Emergence of Genotypes G3P[4] and G3P[8]**

**Eva D. João 1,2,\* , Benilde Munlela 1,3 , Assucênio Chissaque 1,2 , Jorfélia Chilaúle <sup>1</sup> , Jerónimo Langa <sup>1</sup> , Orvalho Augusto 4,5, Simone S. Boene 1,3 , Elda Anapakala <sup>1</sup> , Júlia Sambo 1,2 , Esperança Guimarães 1,2, Diocreciano Bero <sup>1</sup> , Marta Cassocera 1,2, Idalécia Cossa-Moiane <sup>1</sup> , Jason M. Mwenda <sup>6</sup> , Isabel Maurício 2,7, Hester G. O'Neill <sup>8</sup> and Nilsa de Deus 1,9**


Received: 1 July 2020; Accepted: 14 August 2020; Published: 19 August 2020

**Abstract:** Group A rotavirus (RVA) remains the most important etiological agent associated with severe acute diarrhea in children. *Rotarix*® monovalent vaccine was introduced into Mozambique's Expanded Program on Immunization in September 2015. In the present study, we report the diversity and prevalence of rotavirus genotypes, pre- (2012–2015) and post-vaccine (2016–2019) introduction in Mozambique, among diarrheic children less than five years of age. Genotyping data were analyzed for five sentinel sites for the periods indicated. The primary sentinel site, Mavalane General Hospital (HGM), was analyzed for the period 2012–2019, and for all five sites (country-wide analyses), 2015–2019. During the pre-vaccine period, G9P[8] was the most predominant genotype for both HGM (28.5%) and the country-wide analysis (46.0%). However, in the post-vaccine period, G9P[8] was significantly reduced. Instead, G3P[8] was the most common genotype at HGM, while G1P[8] predominated country-wide. Genotypes G9P[4] and G9P[6] were detected for the first time, and the emergence of G3P[8] and G3P[4] genotypes were observed during the post-vaccine period. The distribution and prevalence of rotavirus genotypes were distinct in pre- and post-vaccination periods, while uncommon genotypes were also detected in the post-vaccine period. These observations support the need for continued country-wide surveillance to monitor changes in strain diversity, due to possible vaccine pressure, and consequently, the effect on vaccine effectiveness.

**Keywords:** rotavirus type A; Mozambique vaccine surveillance; G3 genotype; Rotarix

#### **1. Introduction**

Group A rotavirus (RVA) remains the most important etiological agent associated with severe acute diarrhea in children worldwide [1–3]. In 2016, RVA was estimated to cause more than 128,000 deaths among children younger than five years throughout the world, with more than 104,000 deaths occurring in sub-Saharan Africa [3].

RVA is a non-enveloped, double-stranded RNA virus. The segmented genome has 11 gene segments which encode six structural viral proteins (VP1, VP2, VP3, VP4, VP6, and VP7) and six non-structural viral proteins (NSP1, NSP2, NSP3, NSP4, and NSP5/6) [4–6]. The viral capsid is composed of three concentric layers which encapsulate the 11-segmented genome. The outer layer is composed of the viral spike protein, protease-sensitive VP4, and glycoprotein VP7. A dual typing system for RVA is based on the gene segments encoding VP4 (P genotypes) and VP7 (G types). The rotavirus classification-working group has identified 36 G and 51 P genotypes globally in humans and in the young of many mammalian and avian species [7–10]. Six G types (G1, G2, G3, G4, G9, G12) and 3 P types (P[8], P[4], P[6]) predominate globally [11–14], although in Africa and Asia genotypes, such as G5, G6, and G8, are also described as important [15]. The six most frequently reported G/P combinations associated with infections in humans worldwide are G1P[8], G2P[4], G3P[8], G4P[8], G9P[8], and G12P[8] [10–14,16].

In 2009, the World Health Organization (WHO) recommended the introduction of rotavirus vaccines in national immunization programs worldwide and particularly in countries with a high under-five mortality rate associated with diarrhea [17]. The WHO has coordinated the Global Network of Rotavirus surveillance (GNRS) since 2006 to support countries with evidence-based decision-making [10]. Mozambique has actively participated in WHO rotavirus surveillance since 2016. Continuous surveillance of circulating genotypes, as well as the monitoring of disease burden, is important to evaluate the effectiveness of rotavirus vaccines.

Before the introduction of rotavirus vaccines, a high rotavirus disease burden was reported in particular the southern Mozambican region. However, due to a lack of surveillance, no information was available from the center and northern regions of the country [18–20]. In the Global Enteric Multicenter Study (GEMS), which determined the burden and etiology of diarrhea in children under five years of age in four sub-Saharan African and three Asian countries, Mozambique had the highest attributable fraction (27.0%) of rotavirus-associated diarrhea among infants [20]. In Mozambique, the prevalence of rotavirus in under-five year old children from urban (Maputo City) and rural (Manhiça District) areas in 2012 and 2013 was higher than 40.0% [19]. A lower infection rate (24.0%) was, however, reported in 2011 in Gaza province, a rural area [18]. Data from the National Surveillance of Diarrhea also showed a high rotavirus infection rate of 40.2% and 38.3% in 2014 and 2015, respectively, before vaccine introduction in Mozambique [21]. The monovalent vaccine, *Rotarix*® (GlaxoSmithKline, Rixensart, Belgium), was introduced into the Expanded Program on Immunization of Mozambique in September 2015. Since then, the prevalence of rotavirus infections of 12.2% and 13.5% in 2016 and 2017, respectively, has been reported [21].

The evolution of RVA through the accumulation of point mutations, gene reassortment, recombination and interspecies transmission [5,22,23], call for rotavirus strain surveillance to elucidate the effect, if any, of rotavirus vaccine usage on the circulation of rotavirus genotypes in Mozambique. The main objective of the present study was to evaluate the distribution of rotavirus genotypes prior to (2012–2015) and following (2016–2019) rotavirus vaccine introduction in Mozambique, among diarrheic children less than five years of age.

#### **2. Results**

#### *2.1. Comparison of Rotavirus G- and P-Types in Mozambique Pre- and Post-Vaccine Introduction*

From May 2014 to December 2019, a total of 1736 diarrheal stool samples were collected in five sentinel sites as part of the National Surveillance of Diarrhea program in Mozambique. Of these stool samples, 468 tested positive for RVA by ELISA (27.0%) (Supplementary Table S1). A total of 94.0% (440/468) of these samples were genotyped, *n* = 245 from Maputo (HGM and HJM), *n* = 149 from Nampula (HCN), *n* = 34 from Quelimane (HGQ) and *n* = 12 from Beira (HCB) (Supplementary Table S2). During the pre-vaccine period (2014–2015) a total of 246 samples were genotyped and in the post-vaccine period (2016–2019) 194 samples (Supplementary Table S1). In total, 6.0% (28/468) were excluded from genotyping as an insufficient amount of sample was available.

For HGM, a total of 200 genotyped samples corresponded to the pre-vaccine period (2012–2015) and 43 to the post-vaccine period (2016–2019) (Supplementary Table S3). The samples from the pre-vaccine period also included 91 genotyped samples collected at HGM between 2012 and 2013 from a cross-sectional study [24] to extend the analyses for this particular site (Supplementary Table S3).

The analyses for HGM showed that G9 was the most prevalent G type (30.5%) in the pre-vaccine period (*n* = 200), but was significantly reduced to 9.3% during the post-vaccination period (*n* = 43). Similarly, G12 was also significantly reduced (from 18.5% to 2.3%) (Table 1). In contrast, during the pre-vaccination period, no G3 strains were detected; but during the post-vaccine period, the genotype was the most prevalent genotype (48.8%). Interestingly, a small increase in prevalence was observed for the G1 genotype, although this increase was not statistically significant (Table 1).


**Table 1.** Prevalence of G and P types at Mavalane General Hospital pre- and post-vaccine introduction in Mozambique (2012–2019).

1 It is not possible to calculate the Odds-ratio (OR) for cells with a value of 0; <sup>2</sup> Mix G: 2012–2015: G12G8 (2.0%), G12G9 (1.5%), G9G2 (1.5%); 2016–2019: G12G3 (2.3%); <sup>3</sup> x—refers to strains that were non-typeable for G; <sup>4</sup> x—refers to strains that were non-typeable for P; <sup>5</sup> Reference category: Pre-vaccine; Bold: The most prevalent genotypes per period.

P[8] was the most predominant P type in the pre-vaccine period (54.0%) (Table 1), as well as the post-vaccine period (51.2%). Only P[4] (37.2%) (Table 1) had a statistically significant increase during the post-vaccine period (*p* < 0.001). No mixed P types were detected during the post-vaccine period.

When all five sentinel sites (including HGM) were analyzed for the period of 2015–2019, a similar trend was observed for the G9 genotype. During the pre-vaccine period (*n* = 213), G9 was the

most prevalent G type at 49.3%, but a significant reduction for G9 (25.3%) was reported during the post-vaccine period (*n* = 194). The emergence of G3 was also observed, becoming the most prevalent genotype, although only at 26.3% (Table 2). In contrast, a reduction in the prevalence of the G1 genotype was observed (31.5% reduced to 21.6%) for all five sentinel sites.


**Table 2.** Prevalence of G and P types at five sentinel sites in Mozambique during surveillance pre- and post-vaccine introduction (2015–2019).

1 It is not possible to calculate the Odds-ratio (OR) for cells with a value of 0; <sup>2</sup> Mix G—2016–2019: G12G3 (0.5%), G2G1 (0.5%), G3G1 (2.1%), G9G3 (3.1%); <sup>3</sup> x—refers to strains that were non-typeable for G; <sup>4</sup> x—Refers to strains that were non-typeable for P; <sup>5</sup> Reference category: Pre-vaccine; Bold: The most prevalent genotypes per period.

During the pre-vaccine period, P[8] was the most frequently detected P genotype accounting for 85.4% of all genotypes detected (Table 2). However, this high frequency was significantly reduced in the post-vaccination period to less than half (39.2%). An increase in the detection of P[6] (19.1%) and P[4] (36.6%), from almost undetectable, were recorded during this period (Table 2).

Analyses of the data recorded for samples collected at the HGM, showed a slight increase in the odds ratio for G1 type from pre-vaccine to the post-vaccine period of 1.47 times (OR = 1.47 95CI = 0.59–3.44, *p* > 0.330), but a decrease in the odds ratio for genotypes G12 of 90.0% (OR = 0.10, 95CI = 0.003–0.66, *p* < 0.008) and G9 of 77.0% (OR = 0.23, 95CI = 0.01–0.69, *p* < 0.004), respectively (Table 1). Considering all the sentinel sites, a significant decrease was observed in the odds ratio for G1 genotype from pre-vaccine to the post-vaccine period of 40.0% (OR = 0.60, 95CI = 0.37–0.96, *p* < 0.030), as well as a reduction for G9 of 65.0% (OR = 0.35, 95CI = 0.22–0.54, *p* < 0.001) (Table 2).

A reduction for genotype P[8] from pre-vaccine to the post-vaccine period was also observed at HGM (11.0%, OR= 0.89, 95CI= 0.43–1.83, *p* > 0.740, Table 1), as well as for the country-wide sentinel sites (90.0% (OR = 0.10, 95CI = 0.06 to 0.16, *p* < 0.001, Table 2). In contrast, a significant increase in the odds ratio of genotype P[4] of 3.23 times (OR = 3.23, 95CI = 1.44–7.04, *p* < 0.001) was observed at the HGM (Table 1). Analyses for all the sentinel sites showed a high prevalence for P[4] during the post-vaccine period (36.6%) compared to the pre-vaccine period (0.5%).

#### *2.2. Comparison of G*/*P Genotype Combinations in Mozambique Pre- and Post-Vaccine Introduction*

At HGM, the most predominant combinations during the pre-vaccine period were G9P[8] (28.5%), G1P[8] (17.0%), G12P[6] (13.0%) and G2P[4] (10.0%), comprising a total of 68.5% of all genotypes analyzed (Table 3). During the post-vaccine period, G1P[8] (20.9%) was still one of the predominant combinations, although G3P[8] and G3P[4] strains were detected at 25.6% and 18.6%, respectively

(Table 3). A significant reduction in G9P[8] detection was observed following vaccine introduction (*p* < 0.001). Instead, the G9 genotype was now detected in combination with P[4] and P[6] both at a frequency of 4.7% (Table 3).


**Table 3.** G/P type combinations prevalent at Mavalane General Hospital pre- and post-vaccine introduction in Mozambique (2012–2019).

1 It is not possible to calculate the Odds-ratio (OR) for cells with a value of 0; <sup>2</sup> Other genotypes: 2012–2015: G12P[4] (0.5%), G2P[6] (1.0%), G2P[8] (0.5%), G8P[8] (0.5%); 2016–2019: G1P[4] (2.3%), G3P[6] (2.3%), G12P[4] (2.3%); <sup>3</sup> Mixed types: 2012–2015: G12G8P[4] (1.0%), G12G8P[6] (0.5%), G12G8P[6]P[4] (0.5%), G12G9P[6] (0.5%), G12G9P[8]P[6] (1.0%), G12P[8]P[6] (1.0%), G9G2P[4] (0.5%), G9G2P[6] (0.5%), G9G2P[8] (0.5%), G9P[8]P[4] (0.5%); 2016–2019: G12G3P[4] (2.3%); <sup>4</sup> Partial G/P types: 2012–2015: G12P[x] (1.0%), G2P[x] (1.0%),G9P[x] (1.5%), GxP[4] (1.0%), GxP[6] (0.5%), GxP[6]P[4] (0.5%), GxP[8] (4.0%), GxP[8]P[6] (0.5%); 2016–2019: GxP[4] (2.3%),GxP[8] (2.3%); <sup>5</sup> Reference category: Pre-vaccine; Bold: The most prevalent genotypes per period.

The most frequent G/P combinations observed for all the sites participating in the National Surveillance of Diarrhea program during the pre-vaccine period were G9P[8] and G1P[8] at 46.0% and 31.0%, respectively. These combinations comprised a total of 77.0% of all genotypes analyzed (Table 4).

In the post-vaccine period, G1P[8] remained the most frequent G/P combination, but at a reduced frequency of 20.6%. G2P[4] (at a slightly higher frequency) and G2P[6] (similar frequency as in 2015) were, again, detected in the post-vaccine period. Similar to the analysis for HGM, G3 in combination with P[4] (14.4%) and P[8] (9.8%) were detected during the post-vaccine period, together with G9P[4] (12.4%) and G9P[6] (8.8%). Mixed infections, as determined with RT-PCR, was detected for 6.2% of the samples (Table 4).

Analyses for HGM showed an increase in the odds for G1P[8] at 1.29 times (95CI = 0.50–3.07, *p* > 0.54), but a significant decrease in the odds ratio for G9P[8] at 94.0% (OR = 0.06, 95CI = 0.002–0.40, *p* < 0.001) (Table 3).

In contrast, a significant decrease in the odds ratio for all the sentinel sites was observed for G1P[8] at 42.0% (OR = 0.58, 95CI= 0.36–0.93, *p* < 0.020) and G9P[8] at 96.0% (OR = 0.04, 95CI= 0.02–0.10, *p* < 0.001) (Table 4).


**Table 4.** G/P type combinations prevalent at five sentinel sites in Mozambique during surveillance preand post-vaccine introduction (2015–2019).

1 It is not possible to calculate the Odds-ratio (OR) for cells with a value of 0; <sup>2</sup> Other genotypes: 2015: G12P[8] (0.9%), G1P[6] (0.5%); 2016–2019: G12P[4] (0.5%), G12P[8] (0.5%), G1P[4] (1.0%); <sup>3</sup> Mixed types: 2016–2019: G12G3P[4] (0.5%),G2G1P[8] (0.5%), G3G1P[8] (2.6%),G9G3P[6] (2.6%); <sup>4</sup> Partial G/P types: 2015: G9P[x] (3.3%),GxP[6] (0.5%), GxP[8] (7.0%); 2016–2019: G9P[x] (1.0%),GxP[4] (4.6%), GxP[6] (1.6%), GxP[8] (2.1%); <sup>5</sup> Reference category: Pre-vaccine; Bold: The most prevalent genotypes per period.

#### *2.3. Yearly Distribution of Rotavirus Genotypes at the Mavalane General Hospital (HGM) and National Surveillance Sites*

As reported before, G12P[6] (28.6%) and G2P[4] (23.1%) were the most predominant genotype combinations at HGM during 2012–2013 [24]. In 2014 and 2015, G1P[8] and G9P[8] with 84.8% and 73.7%, respectively, were detected at the highest frequencies. In 2016, during the post-vaccine period, the most frequent genotype was G1P[8] with 66.7%. The emergence of new genotypes was observed in 2016 (G3P[4]), which increased in 2017, to the most prevalent genotype (25.0%) followed by G1P[8] (18.8%) (Table 5). In 2018, G3P[8] and G3P[4] became the most prevalent genotype combinations with 36.4% and 27.3%, respectively. Finally, in 2019, only G3P[8] were detected at the HGM. No G1P[8] strains were, therefore, detected in 2018 and 2019 (Table 5).

Since data is available for only one year for all five participating sentinel sites during the pre-vaccine period, yearly analysis for the national surveillance sites are presented from 2015–2019. The results showed that in 2015 the most frequent G/P combination was G9P[8] (46.0%), followed by G1P[8] (31.0%). In 2016, G1P[8] was detected at the highest frequency (43.6%) (Table 6). Other genotype combinations, such as G2P[6] (17.9%), G9P[6] (12.8%), G9P[4] (7.7%), and G3P[4] (2.6%), were also observed in 2016 (Table 6). These results were comparable to those from HGM.

In 2017, G1P[8], as well as G9P[4], were detected at similar frequencies (19.2%), while G3P[4] was detected at 13.5% (Table 6). In 2018 and 2019, G3P[4] and G3P[8] became the most frequently detected genotype combination with 38.7% and 60.0%, respectively (Table 6). G3 was also observed in combination with P[4] (13.5%) and P[6] (1.9%) in 2017, whereas G1P[8] genotype was not detected in 2018, although this genotype was detected at 15.0% in 2019 (Table 6). The results reported for all the sentinel sites participating in the National Surveillance of Diarrhea program is comparable to that observed for HGM for the reporting period (2015–2019), except that G1P[8] was not detected in 2019 for HGM.


**Table 5.** Prevalence of G/P type combinations at Mavalane General Hospital in Mozambique by year.

<sup>1</sup> Other genotypes: 2012: G12P[4] (1.5%), G2P[8] (1.5%), G8P[8] (1.5%); 2015: G2P[6] (2.6%); 2016: G12P[4] (11.1%); 2017: G1P[4] (6.3%), G3P[6] (6.3%); <sup>2</sup> Mixed types: 2012: G12G8P[4] (3.0%), G12G8P[6] (1.5%), G12G8P[6]P[4] (1.5%), G12G9P[6] (1.5%), G12G9P[8]P[6] (3.0%), G12P[8]P[6] (3.0%), G9G2P[4] (1.5%), G9G2P[6] (1.5%),G9G2P[8] (1.5%),G9P[8]P[4] (1.5%); 2016: G12G3P[4] (11.1%); <sup>3</sup> Partial G/P types: 2012: G12P[x] (3.0%), GxP[6]P[4] (1.5%), GxP[6]P[4] (1.5%), GxP[8]P[6] (1.5%); 2013: G2P[x] (8.3%), GxP[4] (8.3%); 2014: GxP[6] (3.0%), GxP[8] (9.1%); 2015: G9P[x] (4.0%), GxP[8] (5.3%); 2017: GxP[4] (6.3%), GxP[8] (6.3%); Grey: The most prevalent genotypes per year.


**Table 6.** Prevalence of G/P type combinations at five sentinel sites in Mozambique during surveillance by year.

<sup>1</sup> Other genotypes: 2015: G12P[8] (0.9%), G1P[6] (0.5%), G2P[4] (0.5%); 2016: G12P[4] (2.6%), G2P[4] (5.1%); 2017: G12P[8] (1.0%), G1P[4] (1.9%); 2018: G2P[4] (3.2%); <sup>2</sup> Mixed types: 2016: G12G3P[4 (2.6%), G2G1P[8] (2.6%); 2017: G3G1P[8] (3.9%), G9G3P[6] (5.8%); <sup>3</sup> Partial G/P types: 2015: G9P[x] (3.3%), GxP[6] (0.5%), GxP[8] (7.0%); 2017: G9P[x] (1.9%), GxP[4] (6.7%), GxP[6] (1.9%), GxP[8] (2.9%); 2018: GxP[4](3.2%), GxP[6] (3.2%); 2019: GxP[4] (5.0%), GxP[8] (5.0%); Grey: The most prevalent genotypes per year.

#### *2.4. Geographical Distribution of Rotavirus Genotypes*

A variation in rotavirus genotypes between the five sentinel sites in Mozambique was observed (Supplementary Table S4).

In the pre-vaccine period (2015), it was observed that G1P[8] occurred in all regions included in this study, with the highest frequency (78.0%) detected in the northern region, at Nampula (HCN) (Supplementary Table S4). In contrast, the G9P[8] genotype combination was mostly detected in the southern region, Maputo (HGM and HJM) at 68.8%. Other uncommon genotypes, such as G2P[6], were mostly detected at Nampula at 10.2% but were not detected in Quelimane (HGQ) or Beira (HCB) (Supplementary Table S4). Similarly, in the post-vaccine period (2016–2019), the combination G1P[8]

was observed across the country. In 2016 at Maputo and Nampula, G1P[8] was the most prevalent genotype with 66.7% and 43.5%, respectively. The G1P[8] genotype was, however, also detected in Quelimane and Beira, which had small sample sizes.

In 2017 the genotype combination G3P[4] was the most prevalent (29.7%) in Maputo, while in Nampula and Quelimane G9P[4] and G9P[6] were the most prevalent at 32.7% and 35.7%, respectively. In 2018 and 2019, the G3 genotypes were predominantly detected in Maputo and Quelimane in combination with P[4] and P[8]. In Nampula and Beira, G3 was detected in combination with P[4] (Supplementary Table S4).

#### **3. Discussion**

Before rotavirus vaccine introduction in Mozambique, RVA surveillance studies focused in the southern region of the country [18–20]. Instituto Nacional de Saúde (INS) initiated national RVA surveillance in the southern region of Mozambique in 2014, which was expanded to other regions (center and north) in 2015. Following the country-wide introduction of *Rotarix*® in September 2015, its impact has been monitored and a substantial reduction in the prevalence of RVA infection rate to 12.2% and 13.5% in 2016 and 2017, respectively, was reported [21]. Since the country is vast, it is important to expand strain surveillance to include the entire country.

In the present analysis, rotavirus surveillance that form part of the National Surveillance of Diarrhea during 2014–2019, as well as data from a cross-section study at the HGM from 2012 and 2013, are reported [24].

During the surveillance at HGM (2012–2019), as well as country-wide sentinel sites (2015–2019), variations in the prevalence of genotypes in the pre- and post-vaccine periods were observed. Genotypes G9 and P[8] were consistently the most prevalent in the pre-vaccine period and in the post-vaccine period, genotypes G3 and P[8] were the most prevalent. However, the proportion of P[8] was reduced, and the prevalence of genotype P[4] increased. These results suggest that genotype prevalence can vary from year to year pre- or post-vaccination in Mozambique.

When comparing the most predominant G/P combinations before and after vaccine introduction at the HGM, G9P[8] was the most predominant genotype combination in the pre-vaccine period, while G1P[8] was the most prevalent genotype combination in the post-vaccine period. The country-wide surveillance also revealed a decreased odds ratio for G9P[8] after the introduction of the vaccine. However, this reduction was accompanied by the emergence of G9P[4] and G9P[6], especially in the northern part of Mozambique, after vaccine introduction. Finally, the emergence of G3P[4] and G3P[8] was also observed. These results showed that in this early phase of rotavirus strain surveillance, it is not clear whether these variations in genotype combinations between both periods were due to the rotavirus vaccine or simply natural variation in genotype frequency. Our results are consistent with previously published studies, as a number of countries from Africa, Europe and America reported a variation in the strain diversity between the two periods [16,25–30].

Countries that introduced the monovalent *Rotarix*® vaccine similar to Mozambique, reported a decline of genotype G1P[8] with a concurrent rise in other combinations in the post-vaccine period. For example, South Africa reported an increase in non-G1P[8] strains [25]. In contrast, in Malawi, the reduction of G1P[8] was not significant [27]. In Ghana, G1P[8] returned as one of the dominant strains in the fourth year post-vaccine introduction [26]. Other studies reported from England, Brazil, Belgium, Scotland, a decline in the proportion of G1P[8] with a rise in the proportion of heterotypic strains, such as G2P[4], was observed [28–31].

Additionally, Belgium reported a slightly lower vaccine effectiveness against G2P[4], and in Malawi, a lower vaccine effectiveness against G2 strains than G1 strains was reported [27].

In our analyses, HGM, with at least four years pre-vaccine data showed a slight increase of G1P[8] after vaccine introduction, although in the country-wide analyses the G1P[8] prevalence was reduced. This needs careful interpretation, due to the difference in the number of years in the pre-vaccine period, one of the limitations of this analysis.

Regarding the variation in the prevalence of some uncommon genotypes (e.g., G9P[4] G9P[6], G3P[4], G3P[6]) detected after vaccine introduction in Mozambique, it is important to mention that a number of studies in Africa [16,25,32] and Asia (India and Japan) also reported these uncommon genotypes before vaccine introduction in low frequency [33,34]. These uncommon genotypes, apart from G9P[6], were also observed in Ireland before vaccine introduction [35–37]. However, a study conducted in Ghana reported the emergence of G9P[4] at a low frequency only during the fourth rotavirus season after vaccine introduction [26].

The emergence of the genotype combinations G3P[4], detected in 2016, 2017, 2018, and G3P[8] in 2018 and 2019 was observed in Mozambique. These strains were also reported in the same period in Botswana after vaccine introduction in 2012 [38]. Botswana also reported an outbreak of G3P[8] in 2018 [39]. In addition, several countries reported G3 in combination with P[4] and P[8] during the 12th African Rotavirus Symposium 2019 [40–42]: Malawi (introduced vaccine in 2012, reported G3P[8] in 2018), South Africa (introduced vaccine in 2009, reported G3P[4] in 2015–2016), Kingdom of Eswatini (introduced vaccine in 2015, reported G3P[8] in 2018). These observations suggest that G3 strains were circulating in Southern Africa during 2015–2018, with a sharp increase in 2018. Around the world, the emergence of genotype G3P[8] and equine-like G3P[8] in 2013 in Australia and re-emergence of G3P[8] were observed in Brazil in the post-vaccine introduction [43–45]. The European Rotavirus Network (EuroRotaNet) reported 2017–2018 for the first time since inception, G3P[8] as the most prevalent strain [28].

Temporal variation of rotavirus strains was observed in Mozambique, in particular in the model site, Mavalane General Hospital (HGM), as data from a cross-sectional study that characterized rotavirus strains at the HGM from 2012 and 2013 [24], was combined with data generated at the same site as part of the National Surveillance program with its inception in 2014. As already mentioned, G12P[6] was the most predominant genotype in 2012, and in 2013, G2P[4] was the most prevalent [24]. In a similar time period, G12P[6]was also reported in the Manhiça District, while in 2011 in the Chókwè district, G12P[8] was the most prevalent genotype [18,24]. These results suggest circulation of G12 during 2011–2012 in southern Mozambique. The G12 genotype was detected at a prevalence of almost 20% in Sub-Saharan Africa during 2012–2013 [10,16]. In 2013 the G2P[4] was the predominant genotype in the Manhiça district [24] and also in South Africa in 2013 [46]. A shift in genotypes was observed in 2014 and 2015 when mostly G1P[8] and G9P[8] strains were detected.

In the post-vaccine period (2016–2019), G9P[8] was replaced by G1P[8] in 2016, while in 2017, G3P[4] was the most predominant followed by G1P[8]. In 2018 and 2019, no G1P[8] strains were detected; instead, the G3P[8] genotype was the most prevalent. The G3P[8] genotype combination is one of the most prevalent strains associated with human rotavirus infection globally [11–14]. However, G3P[4], which is considered an uncommon combination, was also detected. Studies published previously in Mozambique during the pre-vaccine period did not detect these strains. These temporal analyses clearly showed a yearly variation of rotavirus strains, complicating the assessment of vaccine introduction impact on changes in strain diversity [11–14]. These observations are further supported by data generated by the National Surveillance of Diarrhea that also showed a temporal variation of rotavirus strains and may rather represent the natural variation in rotavirus strains.

Evaluation of strains detected at the various sentinel sites between 2015–2019, showed that G1P[8] was detected at all sentinel sites, albeit at a variation in frequency. It is interesting to note that G9P[8] occurred mostly in Maputo (HGM and HJM) in the southern region of the country, while G9 in combination with P[4] and P[6] were observed mostly in the north, Nampula (HCN), and central region, Quelimane (HGQ). The occurrence of G2P[6] was mostly observed in Nampula. The emergence of G3 strains was, however, detected at all sites under surveillance suggesting that the occurrence of these strains was not location bound. Differences in the geographical distribution of genotypes within a country was previously reported [11].

Various challenges and limitations were experienced during the study. These include logistical issues, which led to a delay in the start of surveillance at some sentinel sites. The study was limited by its small sample size; therefore, it was not possible to perform in-depth temporal analyses by the site to access the genetic variability of strains. Furthermore, bias in strain diversity is possible since a low number of strains were characterized at some sentinel sites. Extended pre-vaccine genotyping data (four years) was available for only one sentinel site, whereas only one year genotyping data were available for the remainder of the sentinel sites.

Despite the circulation of diverse rotavirus strains and the emergence of some genotypes, the National Surveillance of Diarrhea reported a reduction in rotavirus prevalence during the early impact study of the rotavirus vaccine after vaccine introduction.

The whole genome characterization of rotavirus strains circulating pre- and post-vaccine introduction will be useful to evaluate any potential vaccine-induced selection of specific antigenic profiles. Moreover, with recent reports related to the emergence of double-reassortant G1P[8] on a DS-1–like genetic backbone [47–49], whole-genome characterization will be important for strains surveillance.

#### **4. Materials and Methods**

#### *4.1. Study Population and Stool Samples Collection*

RVA positive samples, as tested by Enzyme-Linked Immunosorbent Assay (ELISA), were included. Samples were obtained from children under five years of age suffering from moderate-to-severe acute and non-acute diarrhea. These samples were collected as part of an ongoing hospital-based diarrhea surveillance program, called the National Surveillance of Diarrhea (ViNaDia) that commenced in May 2014. Samples were included for this study up to December 2019. In addition, data from a cross-sectional study conducted at the Mavalane General Hospital (HGM) from January 2012 to September 2013 were also included in the analyses [24].

The National Surveillance of Diarrhea in children was led by the "Instituto Nacional de Saúde" (INS), started in May 2014 at the Mavalane General Hospital (HGM, first sentinel site) in the Maputo province (Figure 1). In March 2015, José Macamo General Hospital (HJM), also Maputo Province, and Nampula Central Hospital (HCN), in Nampula province in the northern region of the country were added. Surveillance was extended to two additional sentinel sites in June 2015: Beira Central Hospital (HCB) in Sofala Province and Quelimane General Hospital (HGQ) in the Zambézia province (Figure 1). Since 2016, Mozambique participates and actively report data to the WHO African Rotavirus Surveillance Network (ARSN). ARSN monitors rotavirus infection in children with severe acute watery diarrhea as part of a hospital-based sentinel-site surveillance program.

In the surveillance at HGM and HJM samples were collected and immediately transferred to the INS laboratory, while at HCB, HCN and HGQ, samples were collected and stored at −20 ◦C. Samples were transported on a weekly basis on dry ice to the INS laboratories located in Maputo City for testing and stored in −70 ◦C as previously described [21]. The cross-sectional study was conducted at the Centro de Investigação em Saúde de Manhiça (CISM). The sampling, testing procedures, clinical, socio-demographic information and characterization of rotavirus strains, as previously described [19,24].

#### *4.2. Ethical Approval*

The National Surveillance of Diarrhea in children protocol was reviewed and approved by the Mozambican National Committee on Bioethics for Health (CNBS) (reference N◦ : 348/CNBS/13; IRB00002657), as well as the rotavirus cross-sectional study (reference N◦286/CNBS/10; IRB00002657).

#### *4.3. Laboratory Testing*

#### 4.3.1. Rotavirus Detection and RNA Extraction

All samples analyzed, were tested for rotavirus using the commercial Enzyme-immuno-sorbent assay (ELISA) kit (Prospect, Oxoid Ltd., Hampshire, UK) following the manufacturer's instructions. Total RNA was extracted from ELISA-positive samples using the QIAamp Viral RNA protocol (QIAGEN, Hilden, Germany), and stored at −70 ◦C.

**Figure 1.** Map of Mozambique indicating the geographical location of study sites. Abbreviations for hospitals are indicated in red. HGM (Mavalane General Hospital), HJM (Jose Macamo General Hospital), HCB (Beira Central Hospital), HGQ (Quelimane General Hospital) and HCN (Nampula Central Hospital).

#### 4.3.2. Reverse Transcriptase (RT) and G/P Typing PCR

− − Extracted RNA (8 µL) was reverse transcribed using Con2/Con3 for the partial VP4-encoding gene (VP8\*, 876 bp) and sBeg9/End9 for the VP7-encoding gene. G genotypes were subsequently determined using a multiplex semi-nested PCR as described before [24]. Specific primers that identified the VP7-encoding gene with the following G genotypes: G1, aBT1; G2, aCT2; G3, aET3 or mG3; G4, aDT4; G8, aAT8; G9, aFT9, or mG9; G12, G12b; G10, mG10 in combination with the common primer RVG9 were used as described previously [50–52].

Similarly, Con3 was used in combination with specific primers that identify P genotypes: P[8], 1T-1D or 1T-1v; P[4], 2T-1; P[6], 3T-1; P[9], 4T-1, and P[10], 5T-1, P[11], mp11, P[14], P4943, as described previously [53–55]. The PCR product was analyzed using 2% agarose gel electrophoresis, stained with ethidium bromide and visualized under ultraviolet illumination.

#### *4.4. Data Management and Statistical Analyses*

The rotavirus vaccine, *Rotarix*®, was introduced in September 2015 in Mozambique. Therefore, the pre-vaccine period was considered to be before December 2015, due to logistical problems associated with vaccine introduction across the country.

The genotyping data from the primary sentinel site, Mavalane General Hospital (HGM), was analyzed separately from other sites because data at this site was available from 2012 and other sites from 2015.

Frequencies of identified genotypes are reported. To assess the magnitude of change in genotypes from the pre- to post-vaccine periods, unadjusted odds ratios (OR) and their 95% confidence intervals (95CI) were computed. In this analysis, the genotype was the dependent variable and time the predictor. All statistical analysis was conducted using Stata software version 15.0 (Stata Corp., College Station, TX, USA). A *p*-value of <0.05 was considered statistically significant.

#### **5. Conclusions**

This is the first report describing the circulation of rotavirus genotypes in three regions of Mozambique. A comparison between the pre- and post-vaccine introduction periods showed a shift in circulating genotypes following vaccine introduction. However, due to the short surveillance period, it is not clear if the observed changes were due to the introduction of the vaccine or a consequence of natural strain variation. In addition, the emergence of unusual strains, such as G3P[4] and G3P[8], was also observed, which support the need for continued country-wide surveillance to monitor changes, due to possible vaccine pressure, and consequently, the effect on vaccine effectiveness.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-0817/9/9/671/s1: Table S1: Total number of stool samples collected at sentinel sites in Mozambique during surveillance between May 2014 and December 2019; Table S2: Total number of stool samples collected per sentinel sites in Mozambique during surveillance between May 2014 and December 2019; Table S3: Total number of stool samples collected at Mavalane General Hospital during a cross-sectional study (2012–2013) and the National Surveillance of Diarrhea program (2014–2019); Table S4: Distribution of rotavirus genotypes between geographical regions.

**Author Contributions:** N.d.D. and E.D.J conceptualized the main project. E.D.J., B.M., J.C., J.L., A.C., E.A., J.S., E.G., D.B., M.C., I.C.-M. and S.S.B. performed investigations. Formal analysis and methodology was done by E.D.J. and O.A. Data curation was performed by A.C., M.C. and O.A. Writing of the original draft preparation was performed by E.D.J. Review and editing was done by E.D.J., N.d.D., I.M., J.M.M. and H.G.O. Visualizations was performed by A.C. Project administration was done by J.C. Validation was done by N.d.D., H.G.O. and I.M. Supervision and funding acquisition was done by N.d.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** The National Surveillance of Diarrhea was supported by European Foundation Initiative into African Research in Neglected Tropical Diseases (EFINTD) through a senior fellowship awarded to N.D. (grant number 89539), World Health Organization (WHO), Gavi, the Vaccine Alliance and Fundo Nacional de Investigação (FNI). E.D.J. Ph.D. was supported by Calouste Gulbenkian Foundation. B.M., A.C. and S.B. received a scholarship from the Deutsche Forschungsgemeinschaft (DFG; JO369/5-1).

**Acknowledgments:** We would like to thank the surveillance teams in Maputo, Nampula, Beira and Quelimane, as well as the parents that provided consent for the collection of stool samples and data from children. We are grateful for Centro de Investigação da Polana Caniço, INS, for providing the map illustration and to Adilson Bauhofer for assistance with data analyses.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 World Health Organization; Licensee MDPI, Basel, Switzerland. This is an open access article distributed under the terms of the Creative Commons Attribution IGO License (CC BY) license (http://creativecommons.org/licenses/by/3.0/igo/legalcode), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. In any reproduction of this article there should not be any suggestion that WHO or this article endorse any specific organisation or products. The use of the WHO logo is not permitted.

#### *Article*

### **Whole Genome Characterization and Evolutionary Analysis of G1P[8] Rotavirus A Strains during the Pre- and Post-Vaccine Periods in Mozambique (2012–2017)**

**Benilde Munlela 1,2,\* ,**† **, Eva D. João 1,3,\* ,**† **, Celeste M. Donato 4,5,6, Amy Strydom <sup>7</sup> , Simone S. Boene 1,2, Assucênio Chissaque 1,3 , Adilson F. L. Bauhofer 1,3 , Jerónimo Langa <sup>1</sup> , Marta Cassocera 1,3, Idalécia Cossa-Moiane 1,8 , Jorfélia J. Chilaúle <sup>1</sup> , Hester G. O'Neill <sup>7</sup> and Nilsa de Deus 1,9**


Received: 26 October 2020; Accepted: 3 December 2020; Published: 6 December 2020

**Abstract:** Mozambique introduced the Rotarix® vaccine (GSK Biologicals, Rixensart, Belgium) into the National Immunization Program in September 2015. Although G1P[8] was one of the most prevalent genotypes between 2012 and 2017 in Mozambique, no complete genomes had been sequenced to date. Here we report whole genome sequence analysis for 36 G1P[8] strains using an Illumina MiSeq platform. All strains exhibited a Wa-like genetic backbone (G1-P[8]-I1-R1-C1-M1-A1-N1-T1-E1-H1). Phylogenetic analysis showed that most of the Mozambican strains clustered closely together in a conserved clade for the entire genome. No distinct clustering for pre- and post-vaccine strains were observed. These findings may suggest no selective pressure by the introduction of the Rotarix® vaccine in 2015. Two strains (HJM1646 and HGM0544) showed varied clustering for the entire genome, suggesting reassortment, whereas a further strain obtained from a rural area (MAN0033) clustered separately for all gene segments. Bayesian analysis for the VP7 and VP4 encoding gene segments supported the phylogenetic analysis and indicated a possible introduction from India around 2011.7 and 2013.0 for the main Mozambican clade. Continued monitoring of rotavirus strains in the post-vaccine period is required to fully understand the impact of vaccine introduction on the diversity and evolution of rotavirus strains.

**Keywords:** rotavirus group A; G1P[8]; whole genome sequencing; Rotarix®; Bayesian analysis; Mozambique

#### **1. Introduction**

Rotavirus is one of the leading causes of diarrheal disease in children under five years of age [1,2]. Worldwide, the number of deaths due to rotavirus infection in children under five years of age in 2016 was estimated to be 128,500, of which 104,733 occurred in sub-Saharan Africa [2]. Rotavirus is a member of the *Reoviridae* family. The genome is comprised of 11 double-stranded ribonucleic acid (dsRNA) segments. The mature virus has an icosahedral capsid formed by three concentric protein layers. The 11 segments of the rotavirus genome encode 12 viral proteins: 6 structural proteins VP1-VP4, VP6 and VP7and 6 non-structural proteins (NSP1-NSP6) [3–6].

The gene segments encoding the external capsid proteins, VP7 and VP4, are used in a binary classification system defining G and P genotypes, respectively [5,7]. Currently, 36 G and 51 P genotypes have been described in humans and various animal species [7–10]. At least 73 combinations of human rotavirus group A (RVA) G/P genotypes have been described, of which the most common combinations are G1P[8], G2P[4], G3P[8], G4P[8], G9P[8] and G12P[8] [10,11]. However, the implementation of whole genome sequencing has led to comprehensive sequence-based classification of all RVA genes into genotypes, which are identified and differentiated according to particular cut-off values of nucleotide sequence identities [9,11]. Currently, 26 I (VP6), 22 R (VP1), 20 C (VP2), 20 M (VP3), 31 A (NSP1), 22 N (NSP2), 22 T (NSP3), 27 E (NSP4) and 22 H (NSP5) genotypes have been described [8]. The whole genome constellation of a strain can be described following the nomenclature Gx-P[x]-Ix-Rx-Cx-Mx-Ax-Nx-Tx-Ex-Hx. Two major genotype constellations have been designated for strains that commonly infect humans: Wa-like (I1-R1-C1-M1-A1-N1-T1-E1-H1) and DS-1-like (I2-R2-C2-M2-A2-N2-T2-E2-H2). A third constellation also observed in human strains, called AU-1-like (I3-R3-C3-M3-A3-N3-T3-E3-H3), has been shown to have a feline/canine origin [9,11].

Four live oral vaccines, namely Rotarix® (GlaxoSmithKline Biologics, Rixensart, Belgium), RotaTeq® (Merck & Co., Kenilworth, NJ, USA), Rotavac® (Bharat Biotech, Hyderabad, India) and Rotasiil® (Serum Institute of India Pvt. Ltd., Pune, India) have been prequalified by the World Health Organization [12,13]. Rotarix® and RotaTeq® have been introduced into the immunization programs of more than 100 countries [13]. Rotarix® is a monovalent vaccine containing a single human G1P[8] strain and is administered from the age of six weeks [13]. Prior to vaccine introduction in Mozambique, a high burden of rotavirus disease was reported in children under five years old. The rate of rotavirus infection in urban (Maputo City) and rural (Manhiça District) areas between 2012 and 2013 was 42.4% [6]. In 2011, a 24.0% infection rate was reported in the Gaza province, another rural area in southern Mozambique [14]. In both studies G1P[8] was detected at a low frequency [14,15]. Data from the National Surveillance of Diarrhea (ViNaDia) revealed a high rotavirus infection rate of 40.2% and 38.3% in 2014 and 2015, respectively [16]. The Rotarix® vaccine was introduced in Mozambique in 2015 with increasing vaccine coverage of 70% and 80% in 2016 and 2017, respectively [16,17]. Post vaccine introduction, the rotavirus infection rate was reduced to 12.2% and 13.5% in 2016 and 2017, respectively [16]. During ViNaDia surveillance, G1P[8] strains were consistently observed in the pre- (2012–2015) and post-vaccination period (2016–2019). However, in the post-vaccine period a decrease in G1P[8] strains was observed which coincided with the emergence of other non-G1P[8] genotypes such as G3P[4] and G3P[8] [18].

The whole genomes of Mozambican G2P[4], G8P[4], G12P[6] and G12P[8] RVA strains from the pre-vaccination period have been described [19,20]. However, there are no reports of the whole genome analyses of G1P[8] strains from Mozambique. To address this, the consensus sequences of 36 G1P[8] strains collected between 2012–2017 from vaccinated and non-vaccinated children were analyzed to investigate the diversity and evolution of G1P[8] strains.

#### **2. Results**

#### *2.1. Genome Constellation*

A total of 36 G1P[8] (12 from the pre-vaccine period and 24 from the post-vaccine period) strains were successfully sequenced with an average coverage ranging from 450.0 to 46060.5 per sequence (Supplementary Table S1). Complete open reading frames (ORFs) were obtained for 393 of the 396 genome segments analyzed. A partial ORF (99.0%) for segment four of RVA/Human-wt/MOZ/HCN0690/2015/G1P[8] was obtained, while two genome segments (encoding VP2 and VP3, respectively) of RVA/Human-wt/MOZ/HGM0059/2014/G1P[8] could not be determined as insufficient data were generated for these two segments (Supplementary Table S1). The genotype constellations were determined and all strains exhibited a Wa-like genetic backbone (G1-P[8]-I1-R1-C1-M1-A1-N1-T1-E1-H1). The nucleotide (nt) identities among Mozambican strains varied from 92.5–100.0% and the comparison between Rotarix® and the 11 genes of the Mozambican strains revealed 84.0–97.9% nt identity (Supplementary Table S2).

#### *2.2. Phylogenetic Analyses*

#### 2.2.1. Sequence Analyses of VP7 and VP4

The VP7 encoding sequences of the 36 Mozambican G1P[8] strains, collected between 2012 to 2017 from non-vaccinated and vaccinated children (Supplementary Table S3), were compared with human rotavirus sequences representing VP7 G1 lineages (I-VII) [21–28]. The Mozambican strains clustered into two distinct lineages, I and II (Figure 1a). The majority of Mozambican strains formed a highly conserved clade, and were closely related to various Indian strains circulating between 2012 and 2013. HGM0544 was moderately divergent to the rest of the strains in the clade sharing 99.2–99.7% nucleotide (nt) identity and 98.8–99.4% amino acid (aa) identity. Strains from the pre- and post-vaccine era were intermingled in lineage II. Only two Mozambican strains from this study clustered in the VP7 lineage I and were more diverse than the 34 strains clustering in lineage II. MAN0033, collected in a rural area in southern Mozambique before vaccine introduction, was closely related to Malawian strains from 2012 and to previously characterized Mozambican strains detected in 2011 [14]. HJM1646, collected in southern Mozambique after vaccine introduction, clustered distinctly and only shared 92.5–92.9% nt and 92.7–93.3% aa identity to the other Mozambican strains. HJM1646 clustered with contemporary Indian strains in a sub-lineage of African and global strains (Figure 1a).

The P[8] encoding sequences of the 36 Mozambican strains were compared with human rotavirus sequences representing the four lineages (I–IV) [21–27] (Figure 1b). The Mozambican strains clustered in the major P[8] lineage III. Similar to the VP7 tree, the majority of Mozambican P[8] sequences formed a highly conserved clade, and were closely related to various Indian strains circulating between 2012 and 2013. The P[8] encoding sequence of HJM1646 clustered with HGM0544, despite clustering in different VP7 lineages. These two strains were moderately divergent to the rest of the study strains in the Mozambican clade and clustered close to another Mozambican strain, RVA/Human-wt/MOZ/0060a/2012/G12P[8]P[14], which was previously detected in the Manhiça district in southern Mozambique [19]. MAN0033 clustered distinctly to the rest of the Mozambican strains sharing 95.7–96.2% nt and 98.3–98.7% aa identity and was closely related to contemporary Malawian strains isolated in 2012 (Figure 1b, Supplementary Table S1).

**Figure 1.** Phylogenetic trees based on the ORF (open reading frame) nucleotide sequence of the (**a**) VP7 and (**b**) VP4 genes of G1P[8] strains circulating in Mozambique and global strains obtained from GenBank. The trees were constructed based on the maximum likelihood method implemented in MEGA X [29], applying the best-fit nucleotide substitution model Tamura-3-parameter (T92+G+I) for VP7 and General Time Reversible (GTR-G) for VP4, determined by JModelTest [30]. Bootstrap values (1000 replicates) ≥70% are shown with DS-1 serving as an out-group (not shown in the final tree). Scale bar indicates genetic distance expressed as the number of nucleotide substitutions per site. Pre-vaccine Mozambican strains are indicated by blue squares, post-vaccine by red circles, the Rotarix® vaccine strain by a green triangle and Mozambican strains from previous studies [14,19] are indicated by black triangles. Lineages are defined from I-VIII for VP7 and I-IV for VP4 [21–23,25–28].

#### 2.2.2. Sequence Analyses of VP1-VP3 and VP6

Thirty-three Mozambican strains formed conserved, monophyletic clades that were observed in the VP1, VP2 and VP3 trees, closely related to RVA/Human-wt/IND/CMC00034/2013/G1P[8], and within lineages comprising of contemporary African and global strains (Figure 2a–c). In the VP6 tree, the 35 Mozambican strains clustered together but did not form a discrete monophyletic clade, and were closely related to contemporary Indian G1P[8] strains (Figure 2d). HJM1646 and HGM0544 clustered together in the VP1 tree were moderately divergent from the main Mozambican clade (Figure 2a). In the VP2 tree, HJM1646 clustered close to the monophyletic clade while HGM0544 fell within the Mozambican clade (Figure 2b). In the VP3 tree, these strains clustered together, distinct from the Mozambican clade, adjacent to previously characterized G12P[6] Mozambican strains (Figure 2c) [19]. MAN0033 fell within the same lineage as the main Mozambican clade, showing minor divergence in the VP1 and VP2 tree and more pronounced divergence in the VP3 tree, closely related to contemporary Malawian G1P[8] strains (Figure 2a–c). This strain clustered within a different lineage in the VP6 tree, closely related to the same group of Malawian G1P[8] strains and adjacent to the Mozambican G12P[6] strains (Figure 2d).

#### 2.2.3. Sequence Analyses of NSP1-NSP5/6

The conserved monophyletic clade, comprised of 33 Mozambican strains, was observed in the NSP1–NSP4 trees, with RVA/Human-wt/IND/CMC00034/2013/G1P[8] interspersed within the clade in the NSP2 tree (Figure 2e–h). HJM1646 and HGM0544 continued to show varied clustering patterns across the trees. In the NSP3 tree these strains clustered together and were divergent from the main Mozambican clade, clustering with Indian strains including RVA/Human-wt/IND/CMC00034/2013/G1P[8], and close to Mozambican G12P[6] strains (Figure 2g). HJM1646 was divergent to the Mozambican clade in the NSP1 and NSP4 trees, but clustered close to the monophyletic clade in the NSP2 tree. HGM0544 clustered with the main Mozambican clade in the NSP1, NSP2 and NSP4 trees. In the NSP5 tree, HJM1646 and HGM0544, along with 33 other Mozambican strains, formed a monophyletic clade that was interspersed with global strains (Figure 2i). MAN0033 clustered distinctly to the rest of the Mozambican strains and was closely related to contemporary Malawian strains isolated in 2012 across these trees (Figure 2e–i). The five G12P[6] [19] Mozambican strains fell within neighboring clusters to MAN0033 in the VP6 and NSP2 trees (Figure 2d,f).

#### *2.3. Evolutionary Analysis of VP7 and VP4 Genes*

A randomly subsampled dataset of 378 G1 genes that were representative of global strains temporally and genetically were analyzed (Supplementary Figure S1). The Mozambican strains detected between 2012 and 2017 shared a common ancestral strain circulating in 2009.9 (95% HPD 2008.1–2010.9). Of the Mozambican strains characterized in this study, 34 clustered within the same lineage and shared a most recent common ancestor in 2011.7 (95% HPD 2011.1–2012.0) and diverged from the closest related Indian strains around the same time. MAN0033 and HJM1646 clustered in the other major lineage present in the tree. MAN0033 and closely related Malawian strains shared a common ancestor in 2010.1 (95% HPD 2008.5–2010.9). These variants, circulating in Malawi, Zambia and Mozambique, diverged from a group of Indian G1P[8] strains in 2001.3 (95% HPD 1998.4–2003.5). HJM1646 was divergent to the other G1 strains from Mozambique in this lineage and shared its most recent common ancestor with Indian strains in 2012.9 (95% HPD 2012.4–2013.0) (Supplementary Figure S1).

212

**Figure 2.** Phylogenetic trees based on the ORF nucleotide sequences of the (**a**) VP1, (**b**) VP2, (**c**) VP3, (**d**) VP6, (**e**) NSP1, (**f**) NSP2, (**g**) NSP3, (**h**) NSP4 and (**i**) NSP5

genes of G1P[8] strains circulating in Mozambique and global strains obtained from GenBank. The trees were constructed based on the maximum likelihood method implemented in MEGA X [29], using the best-fit nucleotide substitution model General Time Reversible (GTR+G+I) for VP3, GTR+G for VP2, NSP2 and NSP3, Hasegawa Kishino Yano (HKY+G+I) for VP6 and NSP1, HKY+G for VP1, NSP4 and NSP5/6, determined by JModelTest [30]. Bootstrap values (1000 replicates) ≥70% are shown with DS-1 serving as an out-group (not shown in the final tree). Scale bar indicates genetic distance expressed as the number of nucleotide substitutions per site. Pre-vaccine Mozambican strains are indicated by blue squares, post-vaccine by red circles, the Rotarix® vaccine strain by a green triangle and Mozambican strains from a previous study [19] are indicated by black triangles.

A subsampled dataset of 235 P[8] genes, representative of global strains temporally and genetically, was also analyzed. Thirty-three Mozambican strains characterized in this study clustered within the same lineage and shared a common ancestor in 2013.0 (95% HPD 2012.1–2013.6) and diverged from the closest related Indian strains around 2011.6 (95% HPD 2011.2–2011.9) (Supplementary Figure S2). The most recent common ancestor of HGM0544 and HJM1646 (that was moderately divergent to the rest of the Mozambican strains in the major clade), was estimated to be 2014.1 (95% HPD 2013.0–2014.9). Clustering in a separate lineage to the other Mozambican G1P[8] strains, MAN0033 diverged from the closest Malawian strain in 2010.9 (95% HPD 2009.5–2011.7) (Supplementary Figure S2).

#### *2.4. Comparative Analysis of Neutralizing Antigenic Epitopes of the VP7 and VP4 Genes of Mozambican Strains and the Rotarix*® *Vaccine Strain*

The rotavirus VP7 protein consists of two antigenic epitopes, 7-1 and 7-2, with 7-1 subdivided into 7-1a and 7-1b [31]. The comparative analysis of the VP7 antigenic epitopes between Mozambican strains and the Rotarix® vaccine strain revealed amino acid substitutions in all three antigenic sites. However, most of the amino acid substitutions were observed in antigenic region 7-2. A total of 30 strains shared conserved amino acid differences at positions N147D and 25 strains at M217I. Sporadic mutations were observed in MAN0033 (unvaccinated) and HJM1646 (fully vaccinated) (S123N, K291R and M217T). The HJM1646 strain contained an additional amino acid substitution at N96S (Figure 3).


**Figure 3.** The alignment of amino acids corresponding to three VP7 antigenic epitopes (7-1a, 7-1b and 7-2). The amino acid sequence of Rotarix® is the reference strain and the conserved residues between the Rotarix® to Mozambican strains are indicated by dots (.) and residues that differ are in bold.

Activation of the protein VP4 requires proteolytic cleavage to produce the VP8\* and VP5\* subunits. These regions contain four (8-1 to 8-4) and five (5-1 to 5-5) antigenic epitopes, respectively [32,33]. The amino acid substitutions between the Rotarix® vaccine strain and Mozambican P[8] strains were concentrated in the 8-1 and 8-3 epitopes. There were five conserved amino acid substitutions, at positions E150D, N195D/G, S125N, S131R and N135D. Sporadic mutations were observed in MAN0033 (N195S and N113D), HJM0338 (S146N) and HGM1789 (P114T) (Figure 4).


**Figure 4.** The alignment of the amino acids corresponding to the VP4 antigenic epitopes (8-1, 8-2, 8-3, 8-4 for VP8\* and 5-1, 5-2, 5-3, 5-4, 5-5 for VP5). The amino acid sequence of Rotarix® is the reference strain and the conserved residues between the Rotarix® to Mozambican strains are indicated by dots (.) and residues that differ areinbold.

#### **3. Discussion**

In the present study, whole genome sequencing was performed for 36 G1P[8] RVA strains obtained from Mozambican children with gastroenteritis between 2012–2017 (12 from the pre-vaccine period and 24 from the post-vaccine period). This is the first study to perform whole genome analysis of G1P[8] strains in Mozambique, facilitating the description of genetic diversity and the origins of Mozambican strains.

Of the 36 strains characterized, 33 clustered within the same conserved Mozambican clade across all trees. Two strains, HGM0544 and HJM1646, showed varied patterns by clustering within and were distinct from the Mozambican clade across trees, suggesting these strains had undergone reassortment events. The strain, MAN0033, clustered distinctly from the rest of the Mozambican strains in all trees. This strain was closely related to a conserved group of Malawian G1P[8] strains suggesting that this strain may have been recently introduced from a neighboring country. No distinct clustering patterns were observed based on the year of isolation or vaccination status, which suggests that strains with limited sequence diversity may have circulated among children in the country over the five year period investigated (2012–2017). The homogeneous population of G1P[8] strains suggests that the introduction of Rotarix®has not resulted in a dramatic shift in the diversity of G1P[8] strains circulating in Mozambique. A similar finding was reported in South Africa where no distinct clustering was observed for strains from the pre- and post-vaccine introduction period [26]. Analysis of G1P[8] strains in Brazil over a 27 year period also did not detect any evidence of a selective pressure exerted by the mass introduction of Rotarix® [34]. In contrast, Australia and Belgium reported some unique clusters of G1P[8] strains following vaccine introduction, which may have been due to natural fluctuation or the first signs of vaccine-driven evolution [35]. In Rwanda, unique clusters of G1P[8] strains were identified following RotaTeq introduction [36]. Although neighboring countries reported the widespread (Malawi) and sporadic (South Africa) detection of G1P[8] strains that had undergone reassortment with DS-1 like strains [28,37], all Mozambican strains characterized in this study exhibited a typical Wa-like genetic backbone. Despite some reports of vaccine-derived G1P[8] strains detected in Australia and England, none of the strains identified in this study were derived from the Rotarix® vaccine [38,39].

Although there are seven recognized lineages described for global G1 sequences [40], the majority of Mozambican strains from this study clustered in lineage II, with only two strains clustering in lineage I. However, one strain that clustered in lineage I represented the oldest Mozambican strain sequenced in this study from 2012, which clustered with previously characterized G1P[8] Mozambican strains from 2011 [14]. This may suggest that lineage I strains were replaced in later years by G1 strains associated with lineage II [26,41]. The VP7 lineage I strains were detected in the south of Mozambique which may suggest geographical restriction in the circulation of strains. However, these results can be in part due to the short sampling period of this study. Of the four established lineages of the P[8] genotype [21], all Mozambican strains clustered in lineage III and shared a high level of genetic similarity, except strain MAN0033 which clustered in a distinct sub-lineage.

Maximum likelihood phylogenetic analysis showed that the majority of the Mozambican strains, with the exception of MAN0033, were most closely related to a conserved group of Indian strains across most genes. Even the two reassortant strains (HJM1646 and HGM0544) were most closely related to Indian strains. This suggests that there may have been multiple, contemporary introductions of diverse strains from a similar origin, perhaps India, into Mozambique. This was further supported by the results of the Bayesian analysis, where the time to the most recent common ancestor for the VP7 and VP4 genes of the main Mozambican clade were 2011.7 and 2013.0, respectively, and which had diverged from the closest Indian strain in 2011.7 and 2011.6, respectively. This suggests that the strains became endemic shortly after being introduced and became the dominant variant circulating in the population. These Indian strains were submitted directly to the GenBank database and no associated manuscripts were found, so it is unclear if these strains were associated with any particular outbreak or were detected as part of routine surveillance.

Overall, the VP7 and VP4 antigenic epitopes exhibited conserved substitutions among the Mozambican strains when compared to Rotarix®. The substitutions in the 7-2 VP7 epitope at position M217I, N147D and 8-1, 8-3 VP4 at positions E150D, N195D and S125N, S131R, N135D were observed in pre- and post-vaccine introduction strains, suggesting that these substitutions are not due to the vaccine introduction.

The main limitation of the study was the limited number of strains successfully sequenced, and that fewer strains were sequenced from the pre-vaccine period.

There is a need to expand the whole genome analysis to other strains detected in Mozambique such as G3, G9 in combination with P[4], P[6] and P[8] genotypes reported previously [18] in order to evaluate the possible influence of vaccine introduction on other rotavirus genotypes.

Mozambique has introduced the Rotarix® vaccine, however cases of rotavirus infection associated with G1P[8] strains resulting in hospitalization of children are still being reported. The present analysis showed that G1P[8] strains detected in the post-vaccine period did not undergo significant mutations in the epitope regions that could result in vaccine escape. However, the short post vaccine period analyzed (three years) may have influenced these results, as it may be too early to see major genetic changes associated with vaccine pressure. These results highlight the need for future studies to understand host factors such as the role of histo-blood group antigen status, nutritional status and enteric co-infections that can influence the vaccine effectiveness in Mozambique.

#### **4. Materials and Methods**

#### *4.1. Ethics Approval*

The ViNaDia Protocol was approved by the National Health Bioethics Committee of Mozambique (CNBS) under number (IRB00002657, reference Nr: 348/CNBS/13). Participants' anonymity and confidentiality were guaranteed.

#### *4.2. Sample Collection*

Forty-three fecal samples, collected between 2012 and 2017 that were positive for RVA by ELISA (Prospect EIA rotavirus, Basingstoke, UK) and identified as genotype G1P[8] by multiplex RT-PCR according to described protocols [42,43], were selected for sequencing according to year of isolation, location of collection (region of Mozambique) and the child's vaccination status. The samples were obtained from children <5 years of age, hospitalized with acute gastroenteritis, and collected at five sentinel sites of the National Diarrheal Surveillance (ViNaDia), which are Hospital Geral de Mavalane (HGM), Hospital Geral Jose Macamo (HJM), Hospital Central da Beira (HCB), Hospital Geral de Quelimane (HGQ) and Hospital Central de Nampula (HCN), and from a previous study of Centro de Investigação em Saúde da Manhiça (CISM) in Mozambique between 2012 and 2013 (Supplementary Figure S3) [15]. Clinical information was collected through a structured questionnaire from ViNaDia which included metadata such as age, gender, site and vaccination status.

#### *4.3. RNA Extraction and cDNA Synthesis*

Total RNA was extracted from stool samples with TRI-reagent (Sigma, Darmstadt, Germany) and single-stranded RNA was precipitated with lithium chloride. The self-priming PC3-T7 loop primer (Integrated DNA Technologies, Coralville, IA, USA) was ligated to dsRNA in order to obtain full-length sequences and cDNA was synthesized using the Maxima H Minus double-stranded cDNA kit (Thermo Fisher Scientific, Massachusetts, MA, USA) as previously described [20,44].

#### *4.4. Next Generation Sequencing*

The whole genome sequencing was performed using an Illumina MiSeq sequencing platform (Illumina, Inc. San Diego, CA, USA) at the Next Generation Sequencing Unit at the University of the Free State (NGS-UFS). Sequencing was completed using the Nextera XT DNA Library Preparation Kit (Illumina, Inc., San Diego, CA, USA) using protocols previously describe [19].

#### *4.5. Data Analyses*

A de novo assembly was performed for all samples using CLC Bio Genomics Workbench (12.0.3; Qiagen, Aarhus, Denmark); all contigs with an average coverage above 100 were identified on the Nucleotide Basic Local Alignment Search Tool (BLASTn at the National Center for Biotechnology information (NCBI). References were chosen based on the Blastn results for reference mapping and extraction of consensus sequences for each segment. The genotyping tools, Virus Pathogenic database and analysis resource (ViPR) [45] and RotaC v2.0 [46], were used to determine the genotype of each gene. The sequences were submitted to GenBank and accession numbers MT737379-MT737772 were assigned.

#### *4.6. Phylogenetic Analysis*

Multiple nucleotide sequence alignments with strains obtained from GenBank [24] were made with multiple sequence alignment program (MAFFT v7.450) [47] on Geneious prime v2020.0.3 and Multiple Sequence Comparison by Log Expectation (MUSCLE) [48] alignment available in Molecular Evolutionary Genetic Analysis X (MEGA X) [29]. The optimal nucleotide substitution model for phylogenetic analysis was selected based upon the Akaike information criterion (corrected) (AICc) ranking implemented in the model selection algorithm available on JModelTest [30] and the models selected for each segment were: Tamura-3 parameter (T92+G+I) [49] for VP7; General Time Reversible (GTR+G+I) [50] for VP3; GTR+G for VP2, VP4, NSP2 and NSP3; Hasegawa Kishino Yano (HKY+G+I) for VP6 and NSP1; and HKY+G for VP1, NSP4 and NSP5/6 [51]. The maximum-likelihood trees were generated using MEGA X [29] using 1000 bootstrap replicates to estimate branch support. Pairwise distance matrix nucleotides were obtained in MEGA X using the p-distance algorithm [29]. Amino acid sequences of the VP7 and VP4 Mozambican strains were aligned and epitopes were identified and compared to those of the vaccine strain Rotarix® (A41CB052A, with accession numbers JN849114 and JN849113 for VP7 and VP4) using MEGA X [23,31].

#### *4.7. Evolutionary Analysis*

Maximum likelihood trees were generated using the Randomized Accelerated Maximum Likelihood (RAxML) program (v2.0.0) [52], applying the nucleotide substitution model GTR+G. The trees were used as the input for TempEst 1.5.3 to plot root-to-tip genetic distances, and sequences not conforming to a linear evolutionary pattern were discarded [53]. Time-measured evolutionary histories were reconstructed using the Bayesian Evolutionary Analysis Sampling Trees (BEAST) Program package (v 1.7.5) [54]. The nucleotide substitution model Hasegawa-Kishino-Yano model (HKY+G) for VP7 and GTR+G for VP4 were selected based on AICc raking in jModelTest [30].

The parameters applied included a relaxed uncorrelated lognormal molecular clock to account for varied evolutionary rates among lineages and a coalescent Gaussian Markov random field (GMRF) Bayesian Skyride tree prior. Three independent Markov chain Monte Carlo (MCMC) chains were run for 200 million generations with sampling every 20,000 generations, with the first 10% discarded as burn-in. Convergence and mixing of the chains was assessed using Tracer (v1.7.1) and all parameters yielded effective sample sizes ≥ 200 [55]. The Maximum Clade Credibility (MCC) trees were summarized using TreeAnnotator (v1.10.4) [54]. The time-ordered MCC trees were visualized in FigTree (v1.4.4) (http://tree.bio.ed.ac.uk/software/figtree/).

#### **5. Conclusions**

This study provides important insights into the whole genome sequences of G1P[8] strains in Mozambique. Whilst similar strains were detected prior to and following vaccine introduction, multiple introductions of diverse strains from India highlight the importance of continuously monitoring the strains detected in Mozambique to determine if the strains are evolving by vaccine-induced selection or by natural evolutionary pressures.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2076-0817/9/12/1026/s1: Table S1: Genome assembly of Mozambican Wa-like G1P[8] strains from 2012 to 2017. The percentage identity was determined with BLASTn; Table S2: Nucleotide identities of the Mozambican and Rotarix® vaccine strains; Table S3: Mozambican G1P[8] strains. Figure S1: A simplified maximum clade credibility trees (MCC) for the G1 VP7 strains characterized between 1978 and 2017; Figure S2: A simplified maximum clade credibility trees (MCC) for the P[8] VP4 strains characterized between 2000 and 2017; Figure S3: Mozambique Map with the geographical location of study sentinel sites.

**Author Contributions:** Conceptualization: N.d.D., B.M. and E.D.J.; methodology: B.M., C.M.D., A.S. and E.D.J.; validation: N.d.D., C.M.D. and H.G.O.; formal analysis: B.M., E.D.J., C.M.D. and A.S.; investigation: B.M., E.D.J., J.J.C., J.L., A.C., A.F.L.B., M.C., I.C.-M. and S.S.B.; resources: N.d.D.; data curation: C.M.D. and A.S.; writing—original draft preparation: B.M. and E.D.J.; writing—review and editing: B.M., E.D.J., C.M.D., N.d.D., A.S., H.G.O., A.F.L.B. and A.C.; visualization: B.M., E.D.J., C.M.D., N.d.D., A.S. and H.G.O.; supervision: C.M.D. and N.d.D.; project administration: J.J.C. and funding acquisition: N.d.D. and H.G.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** The study was supported by European Foundation Initiative into African Research in Neglected Tropical Diseases that awarded a senior fellowship to ND (EFINTD; 89539); World Health Organization (WHO); Deutsche Forschungsgemeinschaft (DFG; JO369/5-1) to ND and HGO which supported scholarship of BM, AC, AFLB, AS and SSB; Fundo Nacional de Investigação (FNI) fellowships to BM, AC, and JC. EDJ Ph.D. was supported by Calouste Gulbenkian Foundation. CMD is supported through the Australian National Health and Medical Research Council with an Early Career Fellowship (1113269).

**Acknowledgments:** We want to thank the caretakers who consented for their children to be enrolled in the surveillance. For their efforts with recruitment, data collection and shipment of specimens to Maputo, we would like to thank Elda Anapakala, Esperança Guimarães, Júlia Sambo, Diocreciano Bero, Lena Manhique, Judite Salência, Félix Gundane, Aunésia Marurele, Délcio Muteto, Angelina Pereira, Mulaja Kabeya Étienne, Celso Gabriel, Titos Maulate, Julieta Ernesto, Francisca Ricardo, Siasa Mendes, Hércio Simbine, Susete de Carvalho, Marcos Joaquim, Elvira Sarguene, Fernando Vilanculos, Felicidade Martins, Dulce Graça, Edma Samuel, Vivaldo Pedro, Lúcia Matabel, Maria Safrina, Natércia Abreu, Vanessa da Silva, Nazareth Mabutana, Carlos Guilamba and Celina Nhamuave and Rui Cossa for providing the map illustration. For technical ICT support we thank Stephanus Riekert.

**Conflicts of Interest:** CMD has served on an advisory board for GSK (2019), all payments were paid directly to an administrative fund held by Murdoch Children's Research Institute. All other authors declare no conflicts of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Pathogens* Editorial Office E-mail: pathogens@mdpi.com www.mdpi.com/journal/pathogens

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-0365-2590-7